|
PISM, A Parallel Ice Sheet Model
stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
|
Classes for implementing the Finite Element Method on an Grid. More...
#include "pism/util/fem/DirichletData.hh"#include "pism/util/fem/Element.hh"#include "pism/util/fem/ElementIterator.hh"#include "pism/util/fem/Quadrature.hh"Go to the source code of this file.
Classes | |
| struct | pism::fem::Germ |
| Struct for gathering the value and derivative of a function at a point. More... | |
| struct | pism::fem::QuadPoint |
Namespaces | |
| pism | |
| pism::fem | |
| pism::fem::linear | |
| 1D (linear) elements | |
| pism::fem::q0 | |
| Q0 element information. | |
| pism::fem::q1 | |
| Q1 element information. | |
| pism::fem::p1 | |
| P1 element information. | |
| pism::fem::q13d | |
| Q1 element information. | |
Enumerations | |
| enum | pism::fem::ElementType { pism::fem::ELEMENT_Q = -1 , pism::fem::ELEMENT_P0 = 0 , pism::fem::ELEMENT_P1 = 1 , pism::fem::ELEMENT_P2 = 2 , pism::fem::ELEMENT_P3 = 3 , pism::fem::ELEMENT_EXTERIOR } |
| enum | pism::fem::q13d::ElementFace { pism::fem::q13d::FACE_LEFT = 0 , pism::fem::q13d::FACE_RIGHT = 1 , pism::fem::q13d::FACE_FRONT = 2 , pism::fem::q13d::FACE_BACK = 3 , pism::fem::q13d::FACE_BOTTOM = 4 , pism::fem::q13d::FACE_TOP = 5 } |
Functions | |
| Germ | pism::fem::linear::chi (unsigned int k, const QuadPoint &pt) |
| Linear basis functions on the interval [-1, -1]. More... | |
| Germ | pism::fem::q0::chi (unsigned int k, const QuadPoint &pt) |
| Germ | pism::fem::q1::chi (unsigned int k, const QuadPoint &pt) |
| Q1 basis functions on the reference element with nodes (-1,-1), (1,-1), (1,1), (-1,1). More... | |
| Germ | pism::fem::p1::chi (unsigned int k, const QuadPoint &pt) |
| P1 basis functions on the reference element with nodes (0,0), (1,0), (0,1). More... | |
| ElementType | pism::fem::element_type (int node_type[q1::n_chi]) |
| Germ | pism::fem::q13d::chi (unsigned int k, const QuadPoint &p) |
| Evaluate a Q1 shape function and its derivatives with respect to xi and eta. More... | |
Variables | |
| const unsigned int | pism::fem::MAX_QUADRATURE_SIZE = 9 |
| const int | pism::fem::linear::n_chi = 2 |
| const int | pism::fem::q0::n_chi = 4 |
| const int | pism::fem::q0::n_sides = 4 |
| const int | pism::fem::q1::n_chi = 4 |
| const int | pism::fem::q1::n_sides = 4 |
| const int | pism::fem::p1::n_chi = 3 |
| const int | pism::fem::p1::n_sides = 3 |
| const int | pism::fem::q13d::n_chi = 8 |
| Number of shape functions on a Q1 element. More... | |
| const int | pism::fem::q13d::n_faces = 6 |
| Number of sides per element. More... | |
| const unsigned int | pism::fem::q13d::incident_nodes [n_faces][4] |
Classes for implementing the Finite Element Method on an Grid.
We assume that the reader has a basic understanding of the finite element method. The following is a reminder of the method that also gives the background for the how to implement it on an Grid with the tools in this module.
The Grid domain \(\Omega\) is decomposed into a grid of rectangular physical elements indexed by indices (i,j):
The index of an element corresponds with the index of its lower-left vertex in the grid.
The reference element is the square \([-1,1]\times[-1,1]\). For each physical element \(E_{ij}\), there is an map \(F_{ij}\) from the reference element \(R\) to \(E_{ij}\). In this implementation, the rectangles in the domain are all congruent, and the maps F_{ij} are all the same up to a translation.
On the reference element, vertices are ordered as follows:
For each vertex \(k\in\{0,1,2,3\}\) there is an element basis function \(\phi_k\) that is bilinear, equals 1 at vertex \(k\), and equals 0 at the remaining vertices.
For each node \((i,j)\) in the physical domain there is a basis function \(\psi_{ij}\) that equals 1 at vertex \((i,j)\), and equals zero at all other vertices, and that on element \(E_{i'j'}\) can be written as \(\phi_k\circ F_{i',j'}^{-1}\) for some index \(k\):
\[ \psi_{ij}\Big|_{E_{i'j'}} = \phi_k\circ F_{i',j'}^{-1}. \]
The hat functions \(\psi_{ij}\) form a basis of the function space of piecewise-linear functions. A (scalar) piecewise-linear function \(v=v(x,y)\) on the domain is a linear combination
\[ v = \sum_{i,j} c_{ij}\psi_{ij}. \]
Let \(G(w,v)\) denote the weak form of the PDE under consideration. For example, for the scalar Poisson equation \(-\Delta w = f\),
\[ G(w,v) = \int_\Omega \nabla w \cdot \nabla v -f v \;dx. \]
The SSA weak form is more complicated, in particular because there is dimension 2 at each vertex, corresponding to x- and y-velocity components, but it is in many ways like the Poisson weak form.
In the continuous problem we seek to find a trial function \(w\) such that \(G(w,v)=0\) for all suitable test functions \(v\).
In the discrete problem, we seek a finite element function \(w_h\) such that \(G(w_h,\psi_{ij})=0\) for all suitable indices \((i,j)\). For realistic problems, the integral defining \(G\) cannot be evaluated exactly, but is approximated with some \(G_h\) that arises from numerical quadrature rule: integration on an element \(E\) is approximated with
\[ \int_E f dx \approx \sum_{q} f(x_q) j_q \]
for certain points \(x_q\) and weights \(j_q\) (specific details are found in Quadrature).
The unknown \(w_h\) is represented by an IceVec, \(w_h=\sum_{ij} c_{ij} \psi_{ij}\) where \(c_{ij}\) are the coefficients of the vector. The solution of the finite element problem requires the following computations:
\[ J_{(ij)\;(kl)}=\frac{d}{dc_{kl}} G_h(w_h,\psi_{ij}). \]
Computations of these 'integrals' are done by adding up the contributions from each element (an ElementIterator helps with iterating over the elements). For a fixed element, there is a locally defined trial function \(\hat w_h\) (with 4 degrees of freedom in the scalar case) and 4 local test functions \(\hat\psi_k\), one for each vertex.
The contribution to the integrals proceeds as follows (for concreteness in the case of computing the residual):
Computation of the Jacobian is similar, except that there are now multiple integrals per element (one for each local degree of freedom of \(\hat w_h\)).
All of the above is a simplified description of what happens in practice. The complications below treated by the following classes, and discussed in their documentation:
The classes in this module are not intended to be a fancy finite element package. Their purpose is to clarify the steps that occur in the computation of residuals and Jacobians in SSAFEM, and to isolate and consolidate the hard steps so that they are not scattered about the code.
Definition in file FEM.hh.
1.9.1