PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
Classes | Namespaces | Enumerations | Functions | Variables
FEM.hh File Reference

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]
 

Detailed Description

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):

(0,1) (1,1)
---------
| | |
---------
| | |
---------
(0,0) (1,0)

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:

3 o---------o 2
| |
| |
| |
0 o---------o 1

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:

  1. Evaluation of the residuals \(r_{ij} = G_h(w_h,\psi_{ij})\)
  2. Evaluation of the Jacobian matrix

    \[ 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.