Classes for implementing the Finite Element Method on an IceGrid. More...

#include <petscmat.h>
#include "iceModelVec.hh"
#include "flowlaws.hh"

Classes
struct	FEFunctionGerm
	Struct for gathering the value and derivative of a function at a point. More...
struct	FEVector2Germ
	Struct for gathering the value and derivative of a vector valued function at a point. More...
class	FEShapeQ1
	Computation of Q1 shape function values. More...
class	FEDOFMap
	The mapping from global to local degrees of freedom. More...
class	FEElementMap
	Manages iterating over element indices. More...
class	FEQuadrature
	Numerical integration of finite element functions. More...
Functions
int	PismIntMask (PetscScalar maskvalue)
	Legacy code that needs to vanish.

Detailed Description

Classes for implementing the Finite Element Method on an IceGrid.

We assume that the reader has a basic understanding of the finite element method at the level of ?????. The following is a reminder of the method that also gives the background for the how to implement it on an IceGrid with the tools in this module.

The IceGrid 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$ 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 that equals 1 at vertex , equals zero at all other vertices, and that on element $(i,j)$ can be written as $\phi_k\circ F_{i,j}^{-1}$ for some index $k$ .

A (scalar) finite element function $f$ on the domain is then a linear combination

$f_h = \sum_{i,j} c_{ij}\psi_{ij}.$

Let $G(w,\psi)$ denote the weak form of the PDE under consideration. For example, for the scalar Poisson equation $-\Delta w = f$ ,

$G(w,\psi) = \int_\Omega \nabla w \cdot \nabla \psi -f\psi\;dx.$

In the continuous problem we seek to find a trial function $w$ such that $G(w,\psi)=0$ for all suitable test functions $\psi$ . 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 realistice problems, the integral given by $G$ cannot be evaluated exactly, but is approximated with some $G_h$ that arises from numerical quadrature quadrature rule: integration on an element $E$ is approximated with

$\int_E f dx \approx \sum_{q} f(x_q) w_q$

for certain points $x_q$ and weights $j_q$ (specific details are found in FEQuadrature).

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:

Evaluation of the residuals $r_{ij} = G_h(w_h,\psi_{ij})$
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 FEElementMap 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):

Extract from the global degrees of freedom defining the local degrees of freedom defining $\hat w_h$ . (FEDOFMap)
Evaluate the local trial function (values and derivatives as needed) at the quadrature points (FEQuadrature)
Evaluate the local test functions (again values and derivatives) at the quadrature points. (FEQuadrature)
Obtain the quadrature weights $j_q$ for the element (FEQuadrature)
Compute the values of the integrand $g(\hat w_h,\psi_k)$ at each quadrature point (call these $g_{qk}$ ) and form the weighted sums $y_k=\sum_{q} j_q g_{qk}$ .
Each sum is the contribution of the current element to a residual entry $r_{ij}$ , where local degree of freedom corresponds with global degree of freedom . The local contibutions now need to be added to the global residual vector (FEDOFMap).

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:

Ghost elements (as well as periodic boundary conditions): FEElementMap
Dirichlet data: FEDOFMap
Vector valued functions: (FEDOFMap, FEQuadrature)

It should be mentioned that 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 consolodate the hard steps so that they are not scattered about the code.

Definition in file FETools.hh.

Function Documentation

int PismIntMask ( PetscScalar maskvalue )

Legacy code that needs to vanish.

Todo:: Make it go away.

Definition at line 415 of file FETools.cc.

Referenced by SSAFEM_Forward::assemble_T_rhs(), SSAFEM_Forward::assemble_TStarA_rhs(), and SSAFEM::FixDirichletValues().

src/base/stressbalance/FETools.hh File Reference

Classes

Functions

Detailed Description

Function Documentation