PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
Element.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_ELEMENT_H
20 #define PISM_ELEMENT_H
21 
22 #include <vector>
23 #include <cassert> // assert
24 
25 #include <petscdmda.h> // DMDALocalInfo
26 
27 #include "FEM.hh"
28 #include "pism/util/Vector2.hh"
29 #include "pism/util/petscwrappers/Mat.hh" // Mat, MatStencil
30 
31 
32 namespace pism {
33 
34 class IceGrid;
35 class IceModelVec2Int;
36 
37 struct Vector3 {
38  double x, y, z;
39 };
40 
41 namespace fem {
42 
43 class Quadrature;
44 
45 //! The mapping from global to local degrees of freedom.
46 /*! Computations of residual and Jacobian entries in the finite element method are done by
47  iterating of the elements and adding the various contributions from each element. To do
48  this, degrees of freedom from the global vector of unknowns must be remapped to
49  element-local degrees of freedom for the purposes of local computation, (and the results
50  added back again to the global residual and Jacobian arrays).
51 
52  An Element mediates the transfer between element-local and global degrees of freedom and
53  provides values of shape functions at quadrature points. In this implementation, the
54  global degrees of freedom are either scalars (double's) or vectors (Vector2's), one per
55  node in the IceGrid, and the local degrees of freedom on the element are q1::n_chi or
56  p1::n_chi (%i.e. four or three) scalars or vectors, one for each vertex of the element.
57 
58  In addition to this, the Element transfers locally computed contributions to residual
59  and Jacobian matrices to their global counterparts.
60 */
61 class Element {
62 public:
63  ~Element() = default;
64 
65  int n_chi() const {
66  return m_n_chi;
67  }
68 
69  /*!
70  * `chi(q, k)` returns values and partial derivatives of the `k`-th shape function at a
71  * quadrature point `q`.
72  */
73  const Germ& chi(unsigned int q, unsigned int k) const {
74  assert(q < m_Nq);
75  assert(k < m_n_chi);
76  return m_germs[q * m_n_chi + k];
77  }
78 
79  //! Number of quadrature points
80  int n_pts() const {
81  return m_Nq;
82  }
83 
84  //! Weight of the quadrature point `q`
85  double weight(unsigned int q) const {
86  assert(q < m_Nq);
87  return m_weights[q];
88  }
89 
90  /*! @brief Given nodal values, compute the values at quadrature points.*/
91  //! The output array `result` should have enough elements to hold values at all
92  //! quadrature points.
93  template <typename T>
94  void evaluate(const T *x, T *result) const {
95  for (unsigned int q = 0; q < m_Nq; q++) {
96  result[q] = 0.0;
97  for (unsigned int k = 0; k < m_n_chi; k++) {
98  result[q] += m_germs[q * m_n_chi + k].val * x[k];
99  }
100  }
101  }
102 
103 protected:
104  Element(const IceGrid &grid, int Nq, int n_chi, int block_size);
105  Element(const DMDALocalInfo &grid_info, int Nq, int n_chi, int block_size);
106 
107  /*! @brief Add Jacobian contributions. */
108  void add_contribution(const double *K, Mat J) const;
109 
110  void mark_row_invalid(unsigned int k);
111  void mark_col_invalid(unsigned int k);
112 
113  DMDALocalInfo m_grid;
114 
115  // pointer to a shape function
116  typedef Germ (*ShapeFunction)(unsigned int k, const QuadPoint &p);
117 
118  void initialize(const double J[3][3],
119  ShapeFunction f,
120  unsigned int n_chi,
121  const std::vector<QuadPoint>& points,
122  const std::vector<double>& W);
123 
124  //! grid offsets used to extract nodal values from the grid and add contributions from
125  //! an element to the residual and the Jacobian.
126  std::vector<int> m_i_offset;
127  std::vector<int> m_j_offset;
128 
129  //! Number of nodes (and therefore the number of shape functions) in this particular
130  //! type of elements.
131  const unsigned int m_n_chi;
132 
133  //! Indices of the current element.
134  int m_i, m_j;
135 
136  //! Number of quadrature points
137  const unsigned int m_Nq;
138 
139  const int m_block_size;
140 
141  std::vector<Germ> m_germs;
142 
143  //! Stencils for updating entries of the Jacobian matrix.
144  std::vector<MatStencil> m_row, m_col;
145 
146  //! Constant for marking invalid row/columns.
147  //!
148  //! Has to be negative because MatSetValuesBlockedStencil supposedly ignores negative indexes.
149  //! Seems like it has to be negative and below the smallest allowed index, i.e. -2 and below with
150  //! the stencil of width 1 (-1 *is* an allowed index). We use -2^30 and *don't* use PETSC_MIN_INT,
151  //! because PETSC_MIN_INT depends on PETSc's configuration flags.
152  static const int m_invalid_dof = -1073741824;
153 
154  //! Quadrature weights for a particular physical element
155  std::vector<double> m_weights;
156 private:
158  : m_n_chi(0),
159  m_Nq(0),
160  m_block_size(0) {
161  // empty
162  }
163 };
164 
165 class Element2 : public Element {
166 public:
167 
168  void reset(int i, int j);
169 
170  //! Convert a local degree of freedom index `k` to a global degree of freedom index (`i`,`j`).
171  void local_to_global(unsigned int k, int &i, int &j) const {
172  i = m_i + m_i_offset[k];
173  j = m_j + m_j_offset[k];
174  }
175 
176  Vector2 normal(unsigned int side) const {
177  assert(side < m_n_chi);
178  return m_normals[side];
179  }
180 
181  unsigned int n_sides() const {
182  return n_chi();
183  }
184 
185  double side_length(unsigned int side) const {
186  assert(side < m_n_chi);
187  return m_side_lengths[side];
188  }
189 
190  using Element::evaluate;
191  /*! @brief Given nodal values, compute the values and partial derivatives at the
192  * quadrature points.*/
193  //! Output arrays should have enough elements to hold values at all quadrature points.`
194  template <typename T>
195  void evaluate(const T *x, T *vals, T *dx, T *dy) {
196  for (unsigned int q = 0; q < m_Nq; q++) {
197  vals[q] = 0.0;
198  dx[q] = 0.0;
199  dy[q] = 0.0;
200  for (unsigned int k = 0; k < m_n_chi; k++) {
201  const Germ &psi = m_germs[q * m_n_chi + k];
202  vals[q] += psi.val * x[k];
203  dx[q] += psi.dx * x[k];
204  dy[q] += psi.dy * x[k];
205  }
206  }
207  }
208 
209  /*! @brief Get nodal values of an integer mask. */
210  void nodal_values(const IceModelVec2Int &x_global, int *result) const;
211 
212  /*! @brief Extract nodal values for the element (`i`,`j`) from global array `x_global`
213  into the element-local array `result`.
214  */
215  template<typename T>
216  void nodal_values(T const* const* x_global, T* result) const {
217  for (unsigned int k = 0; k < m_n_chi; ++k) {
218  int i = 0, j = 0;
219  local_to_global(k, i, j);
220  result[k] = x_global[j][i]; // note the indexing order
221  }
222  }
223 
225  /*! @brief Add the values of element-local contributions `y` to the global vector `y_global`. */
226  /*! The element-local array `local` should be an array of Nk values.
227  *
228  * Use this to add residual contributions.
229  */
230  template<typename T>
231  void add_contribution(const T *local, T** y_global) const {
232  for (unsigned int k = 0; k < m_n_chi; k++) {
233  if (m_row[k].i == m_invalid_dof) {
234  // skip rows marked as "invalid"
235  continue;
236  }
237  const int
238  i = m_row[k].i,
239  j = m_row[k].j;
240  y_global[j][i] += local[k]; // note the indexing order
241  }
242  }
243 
246 
247 protected:
248  Element2(const IceGrid &grid, int Nq, int n_chi, int block_size);
249  Element2(const DMDALocalInfo &grid_info, int Nq, int n_chi, int block_size);
250 
251  std::vector<Vector2> m_normals;
252 
253  std::vector<double> m_side_lengths;
254 };
255 
256 //! Q1 element with sides parallel to X and Y axes
257 class Q1Element2 : public Element2 {
258 public:
259  Q1Element2(const IceGrid &grid, const Quadrature &quadrature);
260  Q1Element2(const DMDALocalInfo &grid_info,
261  double dx, double dy,
262  const Quadrature &quadrature);
263 };
264 
265 //! P1 element embedded in Q1Element2
266 class P1Element2 : public Element2 {
267 public:
268  P1Element2(const IceGrid &grid, const Quadrature &quadrature, int type);
269 };
270 
271 class Element3 : public Element {
272 public:
273  using Element::evaluate;
274  /*! @brief Given nodal values, compute the values and partial derivatives at the
275  * quadrature points.*/
276  //! Output arrays should have enough elements to hold values at all quadrature points.`
277  template <typename T>
278  void evaluate(const T *x, T *vals, T *dx, T *dy, T *dz) const {
279  for (unsigned int q = 0; q < m_Nq; q++) {
280  vals[q] = 0.0;
281  dx[q] = 0.0;
282  dy[q] = 0.0;
283  dz[q] = 0.0;
284  for (unsigned int k = 0; k < m_n_chi; k++) {
285  const Germ &psi = m_germs[q * m_n_chi + k];
286  vals[q] += psi.val * x[k];
287  dx[q] += psi.dx * x[k];
288  dy[q] += psi.dy * x[k];
289  dz[q] += psi.dz * x[k];
290  }
291  }
292  }
293 
294  struct GlobalIndex {
295  int i, j, k;
296  };
297 
298  GlobalIndex local_to_global(int i, int j, int k, unsigned int n) const {
299  return {i + m_i_offset[n],
300  j + m_j_offset[n],
301  k + m_k_offset[n]};
302  }
303 
304  GlobalIndex local_to_global(unsigned int n) const {
305  return {m_i + m_i_offset[n],
306  m_j + m_j_offset[n],
307  m_k + m_k_offset[n]};
308  }
309 
310  /*! @brief Extract nodal values for the element (`i`,`j`,`k`) from global array `x_global`
311  into the element-local array `result`.
312  */
313  template<typename T>
314  void nodal_values(T const* const* const* x_global, T* result) const {
315  for (unsigned int n = 0; n < m_n_chi; ++n) {
316  auto I = local_to_global(n);
317  result[n] = x_global[I.j][I.i][I.k]; // note the STORAGE_ORDER
318  }
319  }
320 
322  /*! @brief Add the values of element-local contributions `y` to the global vector `y_global`. */
323  /*! The element-local array `local` should be an array of Nk values.
324  *
325  * Use this to add residual contributions.
326  */
327  template<typename T>
328  void add_contribution(const T *local, T*** y_global) const {
329  for (unsigned int n = 0; n < m_n_chi; n++) {
330  if (m_row[n].i == m_invalid_dof) {
331  // skip rows marked as "invalid"
332  continue;
333  }
334  auto I = local_to_global(n);
335  y_global[I.j][I.i][I.k] += local[n]; // note the STORAGE_ORDER
336  }
337  }
338 protected:
339  Element3(const DMDALocalInfo &grid_info, int Nq, int n_chi, int block_size);
340  Element3(const IceGrid &grid, int Nq, int n_chi, int block_size);
341 
342  std::vector<int> m_k_offset;
343 
344  int m_k;
345 };
346 
347 //! @brief 3D Q1 element
348 /* Regular grid in the x and y directions, scaled in the z direction.
349  *
350  */
351 class Q1Element3 : public Element3 {
352 public:
353  Q1Element3(const DMDALocalInfo &grid,
354  const Quadrature &quadrature,
355  double dx,
356  double dy,
357  double x_min,
358  double y_min);
359  Q1Element3(const IceGrid &grid, const Quadrature &quadrature);
360 
361  void reset(int i, int j, int k, const double *z);
362 
363  // return the x coordinate of node n
364  double x(int n) const {
365  return m_x_min + m_dx * (m_i + m_i_offset[n]);
366  }
367 
368  // return the y coordinate of node n
369  double y(int n) const {
370  return m_y_min + m_dy * (m_j + m_j_offset[n]);
371  }
372 
373  // return the z coordinate of node n
374  double z(int n) const {
375  return m_z_nodal[n];
376  }
377 
380 private:
381  double m_dx;
382  double m_dy;
383  double m_x_min;
384  double m_y_min;
386 
387  // values of shape functions and their derivatives with respect to xi,eta,zeta at all
388  // quadrature points
389  std::vector<Germ> m_chi;
390 
391  // quadrature points (on the reference element)
392  std::vector<QuadPoint> m_points;
393 
394  // quadrature weights (on the reference element)
395  std::vector<double> m_w;
396 };
397 
398 
400 public:
401  Q1Element3Face(double dx, double dy, const Quadrature &quadrature);
402 
403  //! Number of quadrature points
404  int n_pts() const {
405  return m_Nq;
406  }
407  double weight(unsigned int q) const {
408  assert(q < m_Nq);
409  return m_weights[q];
410  }
411 
412  const double& chi(unsigned int q, unsigned int k) const {
413  assert(q < m_Nq);
414  assert(k < m_n_chi);
415  return m_chi[q * m_n_chi + k];
416  }
417 
418  // outward pointing unit normal vector to an element face at the quadrature point q
419  const Vector3& normal(unsigned int q) const {
420  return m_normals[q];
421  }
422 
423  // NB: here z contains nodal z coordinates for *all* nodes of the element
424  void reset(int face, const double *z);
425 
426  template <typename T>
427  void evaluate(const T *x, T *result) const {
428  for (unsigned int q = 0; q < m_Nq; q++) {
429  result[q] = 0.0;
430  for (unsigned int k = 0; k < m_n_chi; k++) {
431  result[q] += m_chi[q * m_n_chi + k] * x[k];
432  }
433  }
434  }
435 protected:
436  // grid spacing
437  double m_dx;
438  double m_dy;
439 
440  // values of shape functions at all quadrature points and their derivatives with respect
441  // to xi, eta, zeta
442  std::vector<double> m_chi;
443 
444  // quadrature points (on the reference element corresponding to a face)
445  std::vector<QuadPoint> m_points;
446 
447  std::vector<Vector3> m_normals;
448 
449  // quadrature weights (on the reference element corresponding to a face)
450  std::vector<double> m_w;
451 
452  // quadrature weights (on the face of a physical element)
453  std::vector<double> m_weights;
454 
455  // number of quadrature points
456  const unsigned int m_n_chi;
457  const unsigned int m_Nq;
458 };
459 
460 } // end of namespace fem
461 } // end of namespace pism
462 
463 #endif /* PISM_ELEMENT_H */
Classes for implementing the Finite Element Method on an IceGrid.
Describes the PISM grid and the distribution of data across processors.
Definition: IceGrid.hh:228
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
Definition: iceModelVec.hh:389
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2.hh:29
void evaluate(const T *x, T *vals, T *dx, T *dy)
Given nodal values, compute the values and partial derivatives at the quadrature points.
Definition: Element.hh:195
void local_to_global(unsigned int k, int &i, int &j) const
Convert a local degree of freedom index k to a global degree of freedom index (i,j).
Definition: Element.hh:171
void reset(int i, int j)
Initialize the Element to element (i, j) for the purposes of inserting into global residual and Jacob...
Definition: Element.cc:187
std::vector< double > m_side_lengths
Definition: Element.hh:253
double side_length(unsigned int side) const
Definition: Element.hh:185
void add_contribution(const T *local, T **y_global) const
Add the values of element-local contributions y to the global vector y_global.
Definition: Element.hh:231
std::vector< Vector2 > m_normals
Definition: Element.hh:251
void nodal_values(T const *const *x_global, T *result) const
Extract nodal values for the element (i,j) from global array x_global into the element-local array re...
Definition: Element.hh:216
Element2(const IceGrid &grid, int Nq, int n_chi, int block_size)
Definition: Element.cc:166
Vector2 normal(unsigned int side) const
Definition: Element.hh:176
unsigned int n_sides() const
Definition: Element.hh:181
void nodal_values(const IceModelVec2Int &x_global, int *result) const
Get nodal values of an integer mask.
Definition: Element.cc:176
Element3(const DMDALocalInfo &grid_info, int Nq, int n_chi, int block_size)
Definition: Element.cc:376
GlobalIndex local_to_global(unsigned int n) const
Definition: Element.hh:304
void nodal_values(T const *const *const *x_global, T *result) const
Extract nodal values for the element (i,j,k) from global array x_global into the element-local array ...
Definition: Element.hh:314
std::vector< int > m_k_offset
Definition: Element.hh:342
void add_contribution(const T *local, T ***y_global) const
Add the values of element-local contributions y to the global vector y_global.
Definition: Element.hh:328
GlobalIndex local_to_global(int i, int j, int k, unsigned int n) const
Definition: Element.hh:298
void evaluate(const T *x, T *vals, T *dx, T *dy, T *dz) const
Given nodal values, compute the values and partial derivatives at the quadrature points.
Definition: Element.hh:278
const Germ & chi(unsigned int q, unsigned int k) const
Definition: Element.hh:73
void evaluate(const T *x, T *result) const
Given nodal values, compute the values at quadrature points.
Definition: Element.hh:94
int n_pts() const
Number of quadrature points.
Definition: Element.hh:80
int n_chi() const
Definition: Element.hh:65
std::vector< int > m_i_offset
Definition: Element.hh:126
void initialize(const double J[3][3], ShapeFunction f, unsigned int n_chi, const std::vector< QuadPoint > &points, const std::vector< double > &W)
Initialize shape function values and quadrature weights of a 2D physical element.
Definition: Element.cc:144
std::vector< double > m_weights
Quadrature weights for a particular physical element.
Definition: Element.hh:155
Germ(* ShapeFunction)(unsigned int k, const QuadPoint &p)
Definition: Element.hh:116
double weight(unsigned int q) const
Weight of the quadrature point q
Definition: Element.hh:85
std::vector< MatStencil > m_col
Definition: Element.hh:144
const unsigned int m_n_chi
Definition: Element.hh:131
void add_contribution(const double *K, Mat J) const
Add Jacobian contributions.
Definition: Element.cc:232
const int m_block_size
Definition: Element.hh:139
std::vector< Germ > m_germs
Definition: Element.hh:141
void mark_row_invalid(unsigned int k)
Mark that the row corresponding to local degree of freedom k should not be updated when inserting int...
Definition: Element.cc:212
const unsigned int m_Nq
Number of quadrature points.
Definition: Element.hh:137
static const int m_invalid_dof
Definition: Element.hh:152
DMDALocalInfo m_grid
Definition: Element.hh:113
std::vector< MatStencil > m_row
Stencils for updating entries of the Jacobian matrix.
Definition: Element.hh:144
void mark_col_invalid(unsigned int k)
Mark that the column corresponding to local degree of freedom k should not be updated when inserting ...
Definition: Element.cc:218
std::vector< int > m_j_offset
Definition: Element.hh:127
int m_i
Indices of the current element.
Definition: Element.hh:134
The mapping from global to local degrees of freedom.
Definition: Element.hh:61
P1Element2(const IceGrid &grid, const Quadrature &quadrature, int type)
Definition: Element.cc:296
P1 element embedded in Q1Element2.
Definition: Element.hh:266
Q1Element2(const IceGrid &grid, const Quadrature &quadrature)
Definition: Element.cc:240
Q1 element with sides parallel to X and Y axes.
Definition: Element.hh:257
const unsigned int m_Nq
Definition: Element.hh:457
Q1Element3Face(double dx, double dy, const Quadrature &quadrature)
Definition: Element.cc:517
std::vector< Vector3 > m_normals
Definition: Element.hh:447
double weight(unsigned int q) const
Definition: Element.hh:407
const double & chi(unsigned int q, unsigned int k) const
Definition: Element.hh:412
std::vector< double > m_w
Definition: Element.hh:450
std::vector< double > m_chi
Definition: Element.hh:442
void evaluate(const T *x, T *result) const
Definition: Element.hh:427
std::vector< QuadPoint > m_points
Definition: Element.hh:445
const Vector3 & normal(unsigned int q) const
Definition: Element.hh:419
std::vector< double > m_weights
Definition: Element.hh:453
const unsigned int m_n_chi
Definition: Element.hh:456
int n_pts() const
Number of quadrature points.
Definition: Element.hh:404
void reset(int face, const double *z)
Definition: Element.cc:534
void reset(int i, int j, int k, const double *z)
Definition: Element.cc:451
std::vector< QuadPoint > m_points
Definition: Element.hh:392
double y(int n) const
Definition: Element.hh:369
Q1Element3(const DMDALocalInfo &grid, const Quadrature &quadrature, double dx, double dy, double x_min, double y_min)
Definition: Element.cc:390
double m_z_nodal[q13d::n_chi]
Definition: Element.hh:385
std::vector< double > m_w
Definition: Element.hh:395
double z(int n) const
Definition: Element.hh:374
std::vector< Germ > m_chi
Definition: Element.hh:389
double x(int n) const
Definition: Element.hh:364
3D Q1 element
Definition: Element.hh:351
Numerical integration of finite element functions.
Definition: Quadrature.hh:62
Contains parameters of an input file grid.
Definition: IceGrid.hh:62
#define n
Definition: exactTestM.c:37
const int n_chi
Number of shape functions on a Q1 element.
Definition: FEM.hh:213
static double K(double psi_x, double psi_y, double speed, double epsilon)
static const double k
Definition: exactTestP.cc:45
const int J[]
Definition: ssafd_code.cc:34
const int I[]
Definition: ssafd_code.cc:24
double y
Definition: Element.hh:38
double z
Definition: Element.hh:38
double x
Definition: Element.hh:38
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