PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
FEM.hh
Go to the documentation of this file.
1 /* Copyright (C) 2020, 2021 PISM Authors
2  *
3  * This file is part of PISM.
4  *
5  * PISM is free software; you can redistribute it and/or modify it under the
6  * terms of the GNU General Public License as published by the Free Software
7  * Foundation; either version 3 of the License, or (at your option) any later
8  * version.
9  *
10  * PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13  * details.
14  *
15  * You should have received a copy of the GNU General Public License
16  * along with PISM; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 #ifndef PISM_FEM_H
20 #define PISM_FEM_H
21 
22 
23 //! \brief Classes for implementing the Finite Element Method on an IceGrid.
24 /*! @file FEM.hh We assume that the reader has a basic understanding of the finite element method. The
25  following is a reminder of the method that also gives the background for the how to implement it
26  on an IceGrid with the tools in this module.
27 
28  The IceGrid domain \f$\Omega\f$ is decomposed into a grid of rectangular physical elements indexed
29  by indices (i,j):
30 
31  ~~~
32  (0,1) (1,1)
33  ---------
34  | | |
35  ---------
36  | | |
37  ---------
38  (0,0) (1,0)
39  ~~~
40 
41  The index of an element corresponds with the index of its lower-left vertex in the grid.
42 
43  The reference element is the square \f$[-1,1]\times[-1,1]\f$. For each physical element
44  \f$E_{ij}\f$, there is an map \f$F_{ij}\f$ from the reference element \f$R\f$ to \f$E_{ij}\f$. In
45  this implementation, the rectangles in the domain are all congruent, and the maps F_{ij} are all
46  the same up to a translation.
47 
48  On the reference element, vertices are ordered as follows:
49 
50  ~~~
51  3 o---------o 2
52  | |
53  | |
54  | |
55  0 o---------o 1
56  ~~~
57 
58  For each vertex \f$k\in\{0,1,2,3\}\f$ there is an element basis function \f$\phi_k\f$ that is bilinear, equals 1 at
59  vertex \f$k\f$, and equals 0 at the remaining vertices.
60 
61  For each node \f$(i,j)\f$ in the physical domain there is a basis function \f$\psi_{ij}\f$ that equals 1 at
62  vertex \f$(i,j)\f$, and equals zero at all other vertices, and that on element \f$E_{i'j'}\f$
63  can be written as \f$\phi_k\circ F_{i',j'}^{-1}\f$ for some index \f$k\f$:
64  \f[
65  \psi_{ij}\Big|_{E_{i'j'}} = \phi_k\circ F_{i',j'}^{-1}.
66  \f]
67 
68  The hat functions \f$\psi_{ij}\f$ form a basis of the function space of piecewise-linear functions.
69  A (scalar) piecewise-linear function \f$v=v(x,y)\f$ on the domain is a linear combination
70  \f[
71  v = \sum_{i,j} c_{ij}\psi_{ij}.
72  \f]
73 
74  Let \f$G(w,v)\f$ denote the weak form of the PDE under consideration. For example, for the
75  scalar Poisson equation \f$-\Delta w = f\f$,
76  \f[
77  G(w,v) = \int_\Omega \nabla w \cdot \nabla v -f v \;dx.
78  \f]
79  The SSA weak form is more complicated, in particular because there is dimension 2 at each vertex,
80  corresponding to x- and y-velocity components, but it is in many ways like the Poisson weak form.
81 
82  In the continuous problem we seek to find a trial function \f$w\f$ such that \f$G(w,v)=0\f$ for
83  all suitable test functions \f$v\f$.
84 
85  In the discrete problem, we seek a finite element function
86  \f$w_h\f$ such that \f$G(w_h,\psi_{ij})=0\f$ for all suitable indices \f$(i,j)\f$. For realistic
87  problems, the integral defining \f$G\f$ cannot be evaluated exactly, but is approximated with some
88  \f$G_h\f$ that arises from numerical quadrature rule: integration on an element \f$E\f$ is
89  approximated with
90  \f[
91  \int_E f dx \approx \sum_{q} f(x_q) j_q
92  \f]
93  for certain points \f$x_q\f$ and weights \f$j_q\f$ (specific details are found in Quadrature).
94 
95  The unknown \f$w_h\f$ is represented by an IceVec, \f$w_h=\sum_{ij} c_{ij} \psi_{ij}\f$ where
96  \f$c_{ij}\f$ are the coefficients of the vector. The solution of the finite element problem
97  requires the following computations:
98 
99  -# Evaluation of the residuals \f$r_{ij} = G_h(w_h,\psi_{ij})\f$
100  -# Evaluation of the Jacobian matrix
101  \f[
102  J_{(ij)\;(kl)}=\frac{d}{dc_{kl}} G_h(w_h,\psi_{ij}).
103  \f]
104 
105  Computations of these 'integrals' are done by adding up the contributions from each element (an
106  ElementIterator helps with iterating over the elements). For a fixed element, there is a locally
107  defined trial function \f$\hat w_h\f$ (with 4 degrees of freedom in the scalar case) and 4 local
108  test functions \f$\hat\psi_k\f$, one for each vertex.
109 
110  The contribution to the integrals proceeds as follows (for concreteness
111  in the case of computing the residual):
112 
113  - Extract from the global degrees of freedom \f$c\f$ defining \f$w_h\f$ the local degrees of
114  freedom \f$d\f$ defining \f$\hat w_h\f$. (Element)
115 
116  - Evaluate the local trial function \f$w_h\f$ (values and derivatives as needed) at the quadrature
117  points \f$x_q\f$ (Quadrature)
118 
119  - Evaluate the local test functions (again values and derivatives) at the quadrature points.
120  (Quadrature)
121 
122  - Obtain the quadrature weights \f$j_q\f$ for the element (Quadrature)
123 
124  - Compute the values of the integrand \f$G(\hat w_h,\psi_k)\f$ at each quadrature point (call
125  these \f$g_{qk}\f$) and form the weighted sums \f$y_k=\sum_{q} j_q g_{qk}\f$.
126 
127  - Each sum \f$y_k\f$ is the contribution of the current element to a residual entry \f$r_{ij}\f$,
128  where local degree of freedom \f$k\f$ corresponds with global degree of freedom \f$(i,j)\f$. The
129  local contibutions now need to be added to the global residual vector (Element).
130 
131  Computation of the Jacobian is similar, except that there are now multiple integrals per element
132  (one for each local degree of freedom of \f$\hat w_h\f$).
133 
134  All of the above is a simplified description of what happens in practice. The complications below
135  treated by the following classes, and discussed in their documentation:
136 
137  - Ghost elements (as well as periodic boundary conditions): ElementIterator
138  - Dirichlet data: Element
139  - Vector valued functions: (Element, Quadrature)
140 
141  The classes in this module are not intended to be a fancy finite element package. Their purpose is
142  to clarify the steps that occur in the computation of residuals and Jacobians in SSAFEM, and to
143  isolate and consolidate the hard steps so that they are not scattered about the code.
144 */
145 
146 namespace pism {
147 namespace fem {
148 
149 //! Struct for gathering the value and derivative of a function at a point.
150 /*! Germ in meant in the mathematical sense, sort of. */
151 struct Germ {
152  //! Function value.
153  double val;
154  //! Function deriviative with respect to x.
155  double dx;
156  //! Function derivative with respect to y.
157  double dy;
158  //! Function derivative with respect to z.
159  double dz;
160 };
161 
162 //! Coordinates of a quadrature point, in the (xi, eta) coordinate space (i.e. on the
163 //! reference element).
164 struct QuadPoint {
165  double xi;
166  double eta;
167  double zeta;
168 };
169 
170 //! Hard-wired maximum number of points a quadrature can use. This is used as the size of
171 //! arrays storing quadrature point values. Some of the entries in such an array may not
172 //! be used.
173 const unsigned int MAX_QUADRATURE_SIZE = 9;
174 
175 //! 1D (linear) elements
176 namespace linear {
177 const int n_chi = 2;
178 Germ chi(unsigned int k, const QuadPoint &pt);
179 } // end of namespace linear
180 
181 //! Q0 element information
182 // FIXME: not sure if Q0 is the right notation here.
183 namespace q0 {
184 const int n_chi = 4;
185 const int n_sides = 4;
186 Germ chi(unsigned int k, const QuadPoint &p);
187 } // end of namespace q0
188 
189 //! Q1 element information
190 namespace q1 {
191 const int n_chi = 4;
192 const int n_sides = 4;
193 Germ chi(unsigned int k, const QuadPoint &p);
194 } // end of namespace q1
195 
196 //! P1 element information
197 namespace p1 {
198 Germ chi(unsigned int k, const QuadPoint &p);
199 const int n_chi = 3;
200 const int n_sides = 3;
201 } // end of namespace p1
202 
206 
207 ElementType element_type(int node_type[q1::n_chi]);
208 
209 //! Q1 element information.
210 namespace q13d {
211 
212 //! Number of shape functions on a Q1 element.
213 const int n_chi = 8;
214 //! Number of sides per element.
215 const int n_faces = 6;
216 //! Evaluate a Q1 shape function and its derivatives with respect to xi and eta.
217 Germ chi(unsigned int k, const QuadPoint &p);
218 
219 /*! Nodes incident to a side. Used to extract nodal values and add contributions.
220  *
221  * The order of faces is used in Q1Element3Face::reset()
222  */
223 const unsigned int incident_nodes[n_faces][4] =
224  {{3, 0, 4, 7}, // 0 - left, xi = -1
225  {1, 2, 6, 5}, // 1 - right, xi = +1
226  {0, 1, 5, 4}, // 2 - front, eta = -1
227  {2, 3, 7, 6}, // 3 - back, eta = +1
228  {0, 1, 2, 3}, // 4 - bottom, zeta = -1
229  {4, 5, 6, 7} // 5 - top, zeta = +1
230 };
231 
237  FACE_TOP = 5};
238 
239 } // end of namespace q13d
240 
241 } // end of namespace fem
242 } // end of namespace pism
243 
244 #include "DirichletData.hh"
245 
246 #include "Element.hh"
247 
248 #include "ElementIterator.hh"
249 
250 #include "Quadrature.hh"
251 
252 #endif /* PISM_FEM_H */
Germ chi(unsigned int k, const QuadPoint &pt)
Linear basis functions on the interval [-1, -1].
Definition: FEM.cc:29
const int n_chi
Definition: FEM.hh:177
const int n_chi
Definition: FEM.hh:199
Germ chi(unsigned int k, const QuadPoint &pt)
P1 basis functions on the reference element with nodes (0,0), (1,0), (0,1).
Definition: FEM.cc:88
const int n_sides
Definition: FEM.hh:200
Germ chi(unsigned int k, const QuadPoint &pt)
Definition: FEM.cc:46
const int n_chi
Definition: FEM.hh:184
const int n_sides
Definition: FEM.hh:185
const int n_faces
Number of sides per element.
Definition: FEM.hh:215
const int n_chi
Number of shape functions on a Q1 element.
Definition: FEM.hh:213
Germ chi(unsigned int k, const QuadPoint &p)
Evaluate a Q1 shape function and its derivatives with respect to xi and eta.
Definition: FEM.cc:138
const unsigned int incident_nodes[n_faces][4]
Definition: FEM.hh:223
const int n_sides
Definition: FEM.hh:192
const int n_chi
Definition: FEM.hh:191
Germ chi(unsigned int k, const QuadPoint &pt)
Q1 basis functions on the reference element with nodes (-1,-1), (1,-1), (1,1), (-1,...
Definition: FEM.cc:70
const unsigned int MAX_QUADRATURE_SIZE
Definition: FEM.hh:173
ElementType
Definition: FEM.hh:203
@ ELEMENT_Q
Definition: FEM.hh:203
@ ELEMENT_P2
Definition: FEM.hh:204
@ ELEMENT_P1
Definition: FEM.hh:204
@ ELEMENT_EXTERIOR
Definition: FEM.hh:205
@ ELEMENT_P0
Definition: FEM.hh:204
@ ELEMENT_P3
Definition: FEM.hh:204
ElementType element_type(int node_type[q1::n_chi])
Definition: FEM.cc:104
static const double k
Definition: exactTestP.cc:45
double dy
Function derivative with respect to y.
Definition: FEM.hh:157
double val
Function value.
Definition: FEM.hh:153
double dz
Function derivative with respect to z.
Definition: FEM.hh:159
double dx
Function deriviative with respect to x.
Definition: FEM.hh:155
Struct for gathering the value and derivative of a function at a point.
Definition: FEM.hh:151