tFETools.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tFETools.hh (25000B)
---
1 // Copyright (C) 2009--2011, 2013, 2014, 2015, 2016, 2017 Jed Brown nd Ed Bueler and Constantine Khroulev and David Maxwell
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 // The following are macros (instead of inline functions) so that error handling
20 // is less cluttered. They should be replaced with empty macros when in
21 // optimized mode.
22
23 #ifndef _FETOOLS_H_
24 #define _FETOOLS_H_
25
26 #include <cassert>
27
28 #include <vector>
29
30 #include <petscmat.h>
31
32 #include "pism/util/Vector2.hh"
33
34 namespace pism {
35 class IceModelVec;
36 class IceModelVec2S;
37 class IceModelVec2Int;
38 class IceModelVec2V;
39 class IceGrid;
40 namespace fem {
41
42 //! Hard-wired maximum number of points a quadrature can use. This is used as the size of arrays
43 //! storing quadrature point values. Some of the entries in such an array may not be used.
44 static const unsigned int MAX_QUADRATURE_SIZE = 9;
45
46 //! \brief Classes for implementing the Finite Element Method on an IceGrid.
47 /*! \file FETools.hh We assume that the reader has a basic understanding of the finite element method. The
48 following is a reminder of the method that also gives the background for the how to implement it
49 on an IceGrid with the tools in this module.
50
51 The IceGrid domain \f$\Omega\f$ is decomposed into a grid of rectangular physical elements indexed
52 by indices (i,j):
53
54 ~~~
55 (0,1) (1,1)
56 ---------
57 | | |
58 ---------
59 | | |
60 ---------
61 (0,0) (1,0)
62 ~~~
63
64 The index of an element corresponds with the index of its lower-left vertex in the grid.
65
66 The reference element is the square \f$[-1,1]\times[-1,1]\f$. For each physical element
67 \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
68 this implementation, the rectangles in the domain are all congruent, and the maps F_{ij} are all
69 the same up to a translation.
70
71 On the reference element, vertices are ordered as follows:
72
73 ~~~
74 3 o---------o 2
75 | |
76 | |
77 | |
78 0 o---------o 1
79 ~~~
80
81 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
82 vertex \f$k\f$, and equals 0 at the remaining vertices.
83
84 For each node \f$(i,j)\f$ in the physical domain there is a basis function \f$\psi_{ij}\f$ that equals 1 at
85 vertex \f$(i,j)\f$, and equals zero at all other vertices, and that on element \f$E_{i'j'}\f$
86 can be written as \f$\phi_k\circ F_{i',j'}^{-1}\f$ for some index \f$k\f$:
87 \f[
88 \psi_{ij}\Big|_{E_{i'j'}} = \phi_k\circ F_{i',j'}^{-1}.
89 \f]
90
91 The hat functions \f$\psi_{ij}\f$ form a basis of the function space of piecewise-linear functions.
92 A (scalar) piecewise-linear function \f$v=v(x,y)\f$ on the domain is a linear combination
93 \f[
94 v = \sum_{i,j} c_{ij}\psi_{ij}.
95 \f]
96
97 Let \f$G(w,v)\f$ denote the weak form of the PDE under consideration. For example, for the
98 scalar Poisson equation \f$-\Delta w = f\f$,
99 \f[
100 G(w,v) = \int_\Omega \nabla w \cdot \nabla v -f v \;dx.
101 \f]
102 The SSA weak form is more complicated, in particular because there is dimension 2 at each vertex,
103 corresponding to x- and y-velocity components, but it is in many ways like the Poisson weak form.
104
105 In the continuous problem we seek to find a trial function \f$w\f$ such that \f$G(w,v)=0\f$ for
106 all suitable test functions \f$v\f$.
107
108 In the discrete problem, we seek a finite element function
109 \f$w_h\f$ such that \f$G(w_h,\psi_{ij})=0\f$ for all suitable indices \f$(i,j)\f$. For realistic
110 problems, the integral defining \f$G\f$ cannot be evaluated exactly, but is approximated with some
111 \f$G_h\f$ that arises from numerical quadrature rule: integration on an element \f$E\f$ is
112 approximated with
113 \f[
114 \int_E f dx \approx \sum_{q} f(x_q) j_q
115 \f]
116 for certain points \f$x_q\f$ and weights \f$j_q\f$ (specific details are found in Quadrature).
117
118 The unknown \f$w_h\f$ is represented by an IceVec, \f$w_h=\sum_{ij} c_{ij} \psi_{ij}\f$ where
119 \f$c_{ij}\f$ are the coefficients of the vector. The solution of the finite element problem
120 requires the following computations:
121
122 -# Evaluation of the residuals \f$r_{ij} = G_h(w_h,\psi_{ij})\f$
123 -# Evaluation of the Jacobian matrix
124 \f[
125 J_{(ij)\;(kl)}=\frac{d}{dc_{kl}} G_h(w_h,\psi_{ij}).
126 \f]
127
128 Computations of these 'integrals' are done by adding up the contributions from each element (an
129 ElementIterator helps with iterating over the elements). For a fixed element, there is a locally
130 defined trial function \f$\hat w_h\f$ (with 4 degrees of freedom in the scalar case) and 4 local
131 test functions \f$\hat\psi_k\f$, one for each vertex.
132
133 The contribution to the integrals proceeds as follows (for concreteness
134 in the case of computing the residual):
135
136 - Extract from the global degrees of freedom \f$c\f$ defining \f$w_h\f$ the local degrees of
137 freedom \f$d\f$ defining \f$\hat w_h\f$. (ElementMap)
138
139 - Evaluate the local trial function \f$w_h\f$ (values and derivatives as needed) at the quadrature
140 points \f$x_q\f$ (Quadrature)
141
142 - Evaluate the local test functions (again values and derivatives) at the quadrature points.
143 (Quadrature)
144
145 - Obtain the quadrature weights \f$j_q\f$ for the element (Quadrature)
146
147 - Compute the values of the integrand \f$G(\hat w_h,\psi_k)\f$ at each quadrature point (call
148 these \f$g_{qk}\f$) and form the weighted sums \f$y_k=\sum_{q} j_q g_{qk}\f$.
149
150 - Each sum \f$y_k\f$ is the contribution of the current element to a residual entry \f$r_{ij}\f$,
151 where local degree of freedom \f$k\f$ corresponds with global degree of freedom \f$(i,j)\f$. The
152 local contibutions now need to be added to the global residual vector (ElementMap).
153
154 Computation of the Jacobian is similar, except that there are now multiple integrals per element
155 (one for each local degree of freedom of \f$\hat w_h\f$).
156
157 All of the above is a simplified description of what happens in practice. The complications below
158 treated by the following classes, and discussed in their documentation:
159
160 - Ghost elements (as well as periodic boundary conditions): ElementIterator
161 - Dirichlet data: ElementMap
162 - Vector valued functions: (ElementMap, Quadrature)
163
164 The classes in this module are not intended to be a fancy finite element package. Their purpose is
165 to clarify the steps that occur in the computation of residuals and Jacobians in SSAFEM, and to
166 isolate and consolidate the hard steps so that they are not scattered about the code.
167 */
168
169 //! Struct for gathering the value and derivative of a function at a point.
170 /*! Germ in meant in the mathematical sense, sort of. */
171 struct Germ
172 {
173 //! Function value.
174 double val;
175 //! Function deriviative with respect to x.
176 double dx;
177 //! Function derivative with respect to y.
178 double dy;
179 };
180
181 //! Coordinates of a quadrature point, in the (xi, eta) coordinate space (i.e. on the reference
182 //! element).
183 struct QuadPoint {
184 double xi;
185 double eta;
186 };
187
188 //! Q0 element information.
189 // FIXME: not sure if Q0 is the right notation here.
190 namespace q0 {
191 //! Number of shape functions on a Q0 element.
192 const int n_chi = 4;
193 //! Evaluate a piecewise-constant shape function and its derivatives.
194 Germ chi(unsigned int k, const QuadPoint &p);
195 } // end of namespace q0
196
197 //! Q1 element information.
198 namespace q1 {
199 //! Number of shape functions on a Q1 element.
200 const int n_chi = 4;
201 //! Number of sides per element.
202 const int n_sides = 4;
203 //! Evaluate a Q1 shape function and its derivatives with respect to xi and eta.
204 Germ chi(unsigned int k, const QuadPoint &p);
205
206 //! Nodes incident to a side. Used to extract nodal values and add contributions.
207 const unsigned int incident_nodes[n_sides][2] = {
208 {0, 1}, {1, 2}, {2, 3}, {3, 0}
209 };
210
211 //! Outward-pointing normal vectors to sides of a Q1 element aligned with the x and y axes.
212 const Vector2* outward_normals();
213
214 }
215
216 //! @brief P1 element embedded in a Q1 element.
217 //! function at node 2).
218 namespace p1 {
219 //! Evaluate a P1 shape function and its derivatives with respect to xi and eta.
220 Germ chi(unsigned int k, const QuadPoint &p);
221
222 //! Number of sides per element.
223 const int n_sides = 3;
224
225 //! Nodes incident to a side. Used to extract nodal values and add contributions.
226 const unsigned int incident_nodes[q1::n_chi][n_sides][2] = {
227 {{0, 1}, {1, 3}, {3, 0}},
228 {{0, 1}, {1, 2}, {2, 0}},
229 {{1, 2}, {2, 3}, {3, 1}},
230 {{2, 3}, {3, 0}, {0, 2}}
231 };
232 } // end of namespace p1
233
234
235 // define Germs, which is an array of q1::n_chi "Germ"s
236 typedef Germ Germs[q1::n_chi];
237
238 //! The mapping from global to local degrees of freedom.
239 /*! Computations of residual and Jacobian entries in the finite element method are done by iterating
240 of the elements and adding the various contributions from each element. To do this, degrees of
241 freedom from the global vector of unknowns must be remapped to element-local degrees of freedom
242 for the purposes of local computation, (and the results added back again to the global residual
243 and Jacobian arrays).
244
245 An ElementMap mediates the transfer between element-local and global degrees of freedom. In this very
246 concrete implementation, the global degrees of freedom are either scalars (double's) or vectors
247 (Vector2's), one per node in the IceGrid, and the local degrees of freedom on the element are
248 q1::n_chi (%i.e. four) scalars or vectors, one for each vertex of the element.
249
250 The ElementMap is also (perhaps awkwardly) overloaded to mediate transfering locally computed
251 contributions to residual and Jacobian matrices to their global counterparts.
252
253 See also: \link FETools.hh FiniteElement/IceGrid background material\endlink.
254 */
255 class ElementMap {
256 public:
257 ElementMap(const IceGrid &grid);
258 ~ElementMap();
259
260 /*! @brief Extract nodal values for the element (`i`,`j`) from global array `x_global`
261 into the element-local array `result`.
262 */
263 template<typename T>
264 void nodal_values(T const* const* x_global, T* result) const {
265 for (unsigned int k = 0; k < q1::n_chi; ++k) {
266 int i = 0, j = 0;
267 local_to_global(k, i, j);
268 result[k] = x_global[j][i]; // note the indexing order
269 }
270 }
271
272 /*! @brief Extract nodal values for the element (`i`,`j`) from global IceModelVec `x_global`
273 into the element-local array `result`.
274 */
275 template<class C, typename T>
276 void nodal_values(const C& x_global, T* result) const {
277 for (unsigned int k = 0; k < q1::n_chi; ++k) {
278 int i = 0, j = 0;
279 local_to_global(k, i, j);
280 result[k] = x_global(i, j); // note the indexing order
281 }
282 }
283
284 /*! @brief Get nodal values of an integer mask. */
285 void nodal_values(const IceModelVec2Int &x_global, int *result) const;
286
287 /*! @brief Add the values of element-local contributions `y` to the global vector `y_global`. */
288 /*! The element-local array `local` should be an array of Nk values.
289 *
290 * Use this to add residual contributions.
291 */
292 template<typename T>
293 void add_contribution(const T *local, T** y_global) const {
294 for (unsigned int k = 0; k < fem::q1::n_chi; k++) {
295 if (m_row[k].k == 1) {
296 // skip rows marked as "invalid"
297 continue;
298 }
299 const int
300 i = m_row[k].i,
301 j = m_row[k].j;
302 y_global[j][i] += local[k]; // note the indexing order
303 }
304 }
305
306 template<class C, typename T>
307 void add_contribution(const T *local, C& y_global) const {
308 for (unsigned int k = 0; k < fem::q1::n_chi; k++) {
309 if (m_row[k].k == 1) {
310 // skip rows marked as "invalid"
311 continue;
312 }
313 const int
314 i = m_row[k].i,
315 j = m_row[k].j;
316 y_global(i, j) += local[k]; // note the indexing order
317 }
318 }
319
320 void reset(int i, int j);
321
322 void mark_row_invalid(int k);
323 void mark_col_invalid(int k);
324
325 //! Convert a local degree of freedom index `k` to a global degree of freedom index (`i`,`j`).
326 void local_to_global(int k, int &i, int &j) const {
327 i = m_i + m_i_offset[k];
328 j = m_j + m_j_offset[k];
329 }
330
331 /*! @brief Add Jacobian contributions. */
332 void add_contribution(const double *K, Mat J) const;
333
334 private:
335 //! Constant for marking invalid row/columns.
336 //!
337 //! Has to be negative because MatSetValuesBlockedStencil supposedly ignores negative indexes.
338 //! Seems like it has to be negative and below the smallest allowed index, i.e. -2 and below with
339 //! the stencil of width 1 (-1 *is* an allowed index). We use -2^30 and *don't* use PETSC_MIN_INT,
340 //! because PETSC_MIN_INT depends on PETSc's configuration flags.
341 static const int m_invalid_dof = -1073741824;
342 static const int m_i_offset[q1::n_chi];
343 static const int m_j_offset[q1::n_chi];
344
345 //! Indices of the current element.
346 int m_i, m_j;
347
348 //! Stencils for updating entries of the Jacobian matrix.
349 MatStencil m_row[q1::n_chi], m_col[q1::n_chi];
350
351 // the grid (for marking ghost DOFs as "invalid")
352 const IceGrid &m_grid;
353 };
354
355 //! Manages iterating over element indices.
356 /*! When computing residuals and Jacobians, there is a loop over all the elements
357 in the IceGrid, and computations are done on each element. The IceGrid
358 has an underlying Petsc DA, and our processor does not own all of the nodes in the grid.
359 Therefore we should not perform computation on all of the elements. In general,
360 an element will have ghost (o) and real (*) vertices:
361 \verbatim
362 o---*---*---*---o
363 | | | | |
364 o---*---*---*---o
365 | | | | |
366 o---o---o---o---o
367 \endverbatim
368 The strategy is to do computations on this processor on every element that has
369 a vertex that is owned by this processor. But we only update entries in the
370 global residual and Jacobian matrices if the corresponding row corresponds to a
371 vertex that we own. In the worst case, where each vertex of an element is owned by
372 a different processor, the computations for that element will be repeated four times,
373 once for each processor.
374
375 This same strategy also correctly deals with periodic boundary conditions. The way Petsc deals
376 with periodic boundaries can be thought of as using a kind of ghost. So the rule still works:
377 compute on all elements containg a real vertex, but only update rows corresponding to that real
378 vertex.
379
380 The calculation of what elements to index over needs to account for ghosts and the
381 presence or absense of periodic boundary conditions in the IceGrid. The ElementIterator performs
382 that computation for you (see ElementIterator::xs and friends).
383
384 See also: \link FETools.hh FiniteElement/IceGrid background material\endlink.
385 */
386 class ElementIterator {
387 public:
388 ElementIterator(const IceGrid &g);
389
390 /*!\brief The total number of elements to be iterated over. Useful for creating per-element storage.*/
391 int element_count() {
392 return xm*ym;
393 }
394
395 /*!\brief Convert an element index (`i`,`j`) into a flattened (1-d) array index, with the first
396 element (`i`, `j`) to be iterated over corresponding to flattened index 0. */
397 int flatten(int i, int j) {
398 return (i-xs) + (j-ys)*xm;
399 }
400
401 //! x-coordinate of the first element to loop over.
402 int xs;
403 //! total number of elements to loop over in the x-direction.
404 int xm;
405
406 //! y-coordinate of the first element to loop over.
407 int ys;
408 //! total number of elements to loop over in the y-direction.
409 int ym;
410
411 //! x-index of the first local element.
412 int lxs;
413 //! total number local elements in x direction.
414 int lxm;
415
416 //! y-index of the first local element.
417 int lys;
418 //! total number local elements in y direction.
419 int lym;
420 };
421
422 class Quadrature {
423 public:
424 Quadrature(unsigned int N);
425 virtual ~Quadrature();
426
427 //! Quadrature size (the number of points).
428 unsigned int n() const {
429 return m_Nq;
430 }
431
432 //! Vaues of `q1::n_chi` trial functions and their derivatives with respect to `x` and `y`, for
433 //! each of `n()` quadrature points.
434 const Germs* test_function_values() const {
435 return m_germs;
436 }
437
438 //! Determinant of the Jacobian of the map from the reference element to the physical element,
439 //! evaluated at quadrature points and multiplied by corresponding quadrature weights.
440 const double* weights() const {
441 return m_W;
442 }
443
444 Germ test_function_values(unsigned int q, unsigned int k) const;
445 double weights(unsigned int q) const;
446
447 protected:
448 //! Number of quadrature points.
449 const unsigned int m_Nq;
450
451 double *m_W;
452
453 Germs* m_germs;
454
455 //! Jacobian of the map from the reference element to a physical element.
456 double m_J[2][2];
457
458 // pointer to a 2D shape function
459 typedef Germ (*ShapeFunction2)(unsigned int k, const QuadPoint &p);
460
461 void initialize(ShapeFunction2 f,
462 unsigned int n_chi,
463 const QuadPoint *points,
464 const double *weights);
465 };
466
467 //! Quadratures on a Q1 element with sides aligned along coordinate axes.
468 /**
469 This code isolates the computation of the Jacobian of the map from the reference element to the
470 physical element and its inverse, which are used to convert partial derivatives with respect to
471 xi, eta to partial derivatives with respect to x and y.
472 */
473 class UniformQxQuadrature : public Quadrature {
474 protected:
475 UniformQxQuadrature(unsigned int size, double dx, double dy, double L);
476 };
477
478 //! Numerical integration of finite element functions.
479 /*! The core of the finite element method is the computation of integrals over elements.
480 For nonlinear problems, or problems with non-constant coefficients (%i.e. any real problem)
481 the integration has to be done approximately:
482 \f[
483 \int_E f(x)\; dx \approx \sum_q f(x_q) w_q
484 \f]
485 for certain quadrature points \f$x_q\f$ and weights \f$w_q\f$. A quadrature is used
486 to evaluate finite element functions at quadrature points, and to compute weights \f$w_q\f$
487 for a given element.
488
489 In this concrete implementation, the reference element \f$R\f$ is the square
490 \f$[-1,1]\times[-1,1]\f$. On a given element, nodes (o) and quadrature points (*)
491 are ordered as follows:
492
493 ~~~
494 3 o------------------o 2
495 | 3 2 |
496 | * * |
497 | |
498 | |
499 | * * |
500 | 0 1 |
501 0 o------------------o 1
502 ~~~
503
504 There are four quad points per element, which occur at \f$x,y=\pm 1/\sqrt{3}\f$. This corresponds
505 to the tensor product of Gaussian integration on an interval that is exact for cubic functions on
506 the interval.
507
508 Integration on a physical element can be thought of as being done by change of variables. The
509 quadrature weights need to be modified, and the Quadrature takes care of this. Because all
510 elements in an IceGrid are congruent, the quadrature weights are the same for each element, and
511 are computed upon initialization.
512
513 See also: \link FETools.hh FiniteElement/IceGrid background material\endlink.
514 */
515 class Q1Quadrature4 : public UniformQxQuadrature {
516 public:
517 Q1Quadrature4(double dx, double dy, double L=1.0);
518 private:
519 static const unsigned int m_size = 4;
520 };
521
522 //! The 9-point 2D Gaussian quadrature on the square [-1,1]*[-1,1].
523 class Q1Quadrature9 : public UniformQxQuadrature {
524 public:
525 Q1Quadrature9(double dx, double dy, double L=1.0);
526 private:
527 static const unsigned int m_size = 9;
528 };
529
530 //! The 16-point 2D Gaussian quadrature on the square [-1,1]*[-1,1].
531 class Q1Quadrature16 : public UniformQxQuadrature {
532 public:
533 Q1Quadrature16(double dx, double dy, double L=1.0);
534 private:
535 static const unsigned int m_size = 16;
536 };
537
538 class Q0Quadrature1e4 : public UniformQxQuadrature {
539 public:
540 Q0Quadrature1e4(double dx, double dy, double L=1.0);
541 private:
542 static const unsigned int m_size_1d = 100;
543 static const unsigned int m_size = m_size_1d * m_size_1d;
544 };
545
546 class Q1Quadrature1e4 : public UniformQxQuadrature {
547 public:
548 Q1Quadrature1e4(double dx, double dy, double L=1.0);
549 private:
550 static const unsigned int m_size_1d = 100;
551 static const unsigned int m_size = m_size_1d * m_size_1d;
552 };
553
554 //! Quadratures on a P1 element.
555 class P1Quadrature : public Quadrature {
556 protected:
557 P1Quadrature(unsigned int size, unsigned int N,
558 double dx, double dy, double L);
559 };
560
561 class P1Quadrature3 : public P1Quadrature {
562 public:
563 P1Quadrature3(unsigned int N, double dx, double dy, double L);
564 private:
565 static const unsigned int m_size = 3;
566 };
567
568 /*! @brief Compute the values at the quadrature points of a finite-element function.*/
569 /*! There should be room for Q.N() values in the output array `result`. */
570 template <typename T>
571 void quadrature_point_values(Quadrature &Q, const T *x, T *result) {
572 const Germs *test = Q.test_function_values();
573 const unsigned int n = Q.n();
574 for (unsigned int q = 0; q < n; q++) {
575 result[q] = 0.0;
576 for (unsigned int k = 0; k < q1::n_chi; k++) {
577 result[q] += test[q][k].val * x[k];
578 }
579 }
580 }
581
582 /*! @brief Compute the values and partial derivatives at the quadrature points of a finite-element
583 function.*/
584 /*! There should be room for Q.N() values in the output vectors `vals`, `dx`, and `dy`.
585 Each element of `dx` is the derivative of the vector-valued finite-element function in the x
586 direction, and similarly for `dy`.
587 */
588 template <typename T>
589 void quadrature_point_values(Quadrature &Q, const T *x, T *vals, T *dx, T *dy) {
590 const Germs *test = Q.test_function_values();
591 const unsigned int n = Q.n();
592 for (unsigned int q = 0; q < n; q++) {
593 vals[q] = 0.0;
594 dx[q] = 0.0;
595 dy[q] = 0.0;
596 for (unsigned int k = 0; k < q1::n_chi; k++) {
597 const Germ &psi = test[q][k];
598 vals[q] += psi.val * x[k];
599 dx[q] += psi.dx * x[k];
600 dy[q] += psi.dy * x[k];
601 }
602 }
603 }
604
605 //* Parts shared by scalar and 2D vector Dirichlet data classes.
606 class DirichletData {
607 public:
608 void constrain(ElementMap &element);
609 operator bool() {
610 return m_indices != NULL;
611 }
612 protected:
613 DirichletData();
614 ~DirichletData();
615
616 void init(const IceModelVec2Int *indices, const IceModelVec *values, double weight = 1.0);
617 void finish(const IceModelVec *values);
618
619 const IceModelVec2Int *m_indices;
620 double m_indices_e[q1::n_chi];
621 double m_weight;
622 };
623
624 class DirichletData_Scalar : public DirichletData {
625 public:
626 DirichletData_Scalar(const IceModelVec2Int *indices, const IceModelVec2S *values,
627 double weight = 1.0);
628 ~DirichletData_Scalar();
629
630 void enforce(const ElementMap &element, double* x_e);
631 void enforce_homogeneous(const ElementMap &element, double* x_e);
632 void fix_residual(double const *const *const x_global, double **r_global);
633 void fix_residual_homogeneous(double **r_global);
634 void fix_jacobian(Mat J);
635 protected:
636 const IceModelVec2S *m_values;
637 };
638
639 class DirichletData_Vector : public DirichletData {
640 public:
641 DirichletData_Vector(const IceModelVec2Int *indices, const IceModelVec2V *values,
642 double weight);
643 ~DirichletData_Vector();
644
645 void enforce(const ElementMap &element, Vector2* x_e);
646 void enforce_homogeneous(const ElementMap &element, Vector2* x_e);
647 void fix_residual(Vector2 const *const *const x_global, Vector2 **r_global);
648 void fix_residual_homogeneous(Vector2 **r);
649 void fix_jacobian(Mat J);
650 protected:
651 const IceModelVec2V *m_values;
652 };
653
654 //! 2-point quadrature for sides of Q1 quadrilateral elements.
655 class BoundaryQuadrature2 {
656 public:
657 BoundaryQuadrature2(double dx, double dy, double L);
658
659 inline unsigned int n() const;
660
661 inline double weights(unsigned int side) const;
662
663 inline const Germ& germ(unsigned int side,
664 unsigned int func,
665 unsigned int pt) const;
666 private:
667 //! Number of quadrature points per side.
668 static const unsigned int m_Nq = 2;
669 Germ m_germs[q1::n_sides][m_Nq][q1::n_chi];
670 double m_W[q1::n_sides];
671 };
672
673 inline unsigned int BoundaryQuadrature2::n() const {
674 return m_Nq;
675 }
676
677 inline double BoundaryQuadrature2::weights(unsigned int side) const {
678 assert(side < q1::n_sides);
679 return m_W[side];
680 }
681
682 //! @brief Return the "germ" (value and partial derivatives) of a basis function @f$ \chi_k @f$
683 //! evaluated at the point `pt` on the side `side` of an element.
684 inline const Germ& BoundaryQuadrature2::germ(unsigned int side,
685 unsigned int q,
686 unsigned int k) const {
687 assert(side < q1::n_sides);
688 assert(k < q1::n_chi);
689 assert(q < 2);
690
691 return m_germs[side][q][k];
692 }
693
694 } // end of namespace fem
695 } // end of namespace pism
696
697 #endif/* _FETOOLS_H_*/