PISM, A Parallel Ice Sheet Model  stable v2.0.4 committed by Constantine Khrulev on 2022-05-25 12:02:27 -0800
Element.cc
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 #include <cassert> // assert
20 #include <cmath> // std::sqrt()
21 
22 #include "Element.hh"
23 #include "pism/util/IceGrid.hh"
24 #include "pism/util/iceModelVec.hh"
25 #include "pism/util/error_handling.hh"
26 #include "pism/util/petscwrappers/DM.hh"
27 
28 namespace pism {
29 namespace fem {
30 
31 //! Determinant of a 3x3 matrix
32 static double det(const double a[3][3]) {
33  return (a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1]) -
34  a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0]) +
35  a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]));
36 }
37 
38 //! Cross product of two 3D vectors
39 static Vector3 cross(const Vector3 &a, const Vector3 &b) {
40  return {a.y * b.z - a.z * b.y,
41  a.z * b.x - a.x * b.z,
42  a.x * b.y - a.y * b.x};
43 }
44 
45 // extract a row of a 3x3 matrix
46 static Vector3 row(const double A[3][3], size_t k) {
47  return {A[k][0], A[k][1], A[k][2]};
48 }
49 
50 // extract a column of a 3x3 matrix
51 static Vector3 column(const double A[3][3], size_t k) {
52  return {A[0][k], A[1][k], A[2][k]};
53 }
54 
55 // dot product of a vector and a [dx, dy, dz] vector in Germ
56 static double dot(const Vector3 &v, const Germ &a) {
57  return a.dx * v.x + a.dy * v.y + a.dz + v.z;
58 }
59 
60 //! Invert a 3x3 matrix
61 static void invert(const double A[3][3], double result[3][3]) {
62  const double det_A = det(A);
63 
64  assert(det_A != 0.0);
65 
66  Vector3
67  x0 = column(A, 0),
68  x1 = column(A, 1),
69  x2 = column(A, 2),
70  a = cross(x1, x2),
71  b = cross(x2, x0),
72  c = cross(x0, x1);
73 
74  double A_cofactor[3][3] = {{a.x, a.y, a.z},
75  {b.x, b.y, b.z},
76  {c.x, c.y, c.z}};
77 
78  // inverse(A) = 1/det(A) * transpose(A_cofactor)
79  for (int i = 0; i < 3; ++i) {
80  for (int j = 0; j < 3; ++j) {
81  // note the transpose on the RHS
82  result[i][j] = A_cofactor[j][i] / det_A;
83  }
84  }
85 }
86 
87 //! Compute derivatives with respect to x,y using J^{-1} and derivatives with respect to xi, eta.
88 static Germ multiply(const double A[3][3], const Germ &v) {
89  // FIXME: something is not right here
90  return {v.val,
91  dot(row(A, 0), v),
92  dot(row(A, 1), v),
93  dot(row(A, 2), v)};
94 }
95 
96 static void set_to_identity(double A[3][3]) {
97  A[0][0] = 1.0;
98  A[0][1] = 0.0;
99  A[0][2] = 0.0;
100  A[1][0] = 0.0;
101  A[1][1] = 1.0;
102  A[1][2] = 0.0;
103  A[2][0] = 0.0;
104  A[2][1] = 0.0;
105  A[2][2] = 1.0;
106 }
107 
108 Element::Element(const IceGrid &grid, int Nq, int n_chi, int block_size)
109  : m_n_chi(n_chi),
110  m_Nq(Nq),
111  m_block_size(block_size) {
112 
113  // get sub-domain information from the grid:
114  auto da = grid.get_dm(1, 0); // dof = 1, stencil_width = 0
115  PetscErrorCode ierr = DMDAGetLocalInfo(*da, &m_grid);
116  if (ierr != 0) {
117  throw std::runtime_error("Failed to allocate an Element instance");
118  }
119  // reset da: we don't want to end up depending on it
120  m_grid.da = NULL;
121 
122  m_germs.resize(Nq * n_chi);
123  m_row.resize(block_size);
124  m_col.resize(block_size);
125 }
126 
127 Element::Element(const DMDALocalInfo &grid_info, int Nq, int n_chi, int block_size)
128  : m_n_chi(n_chi),
129  m_Nq(Nq),
130  m_block_size(block_size) {
131 
132  m_grid = grid_info;
133  // reset da: we don't want to end up depending on it
134  m_grid.da = NULL;
135 
136  m_germs.resize(Nq * n_chi);
137  m_row.resize(block_size);
138  m_col.resize(block_size);
139 }
140 
141 //! Initialize shape function values and quadrature weights of a 2D physical element.
142 /** Assumes that the Jacobian does not depend on coordinates of the current quadrature point.
143  */
144 void Element::initialize(const double J[3][3],
145  ShapeFunction f,
146  unsigned int n_chi,
147  const std::vector<QuadPoint>& points,
148  const std::vector<double>& W) {
149 
150  double J_inv[3][3];
151  invert(J, J_inv);
152 
153  for (unsigned int q = 0; q < m_Nq; q++) {
154  for (unsigned int k = 0; k < n_chi; k++) {
155  m_germs[q * m_n_chi + k] = multiply(J_inv, f(k, points[q]));
156  }
157  }
158 
159  m_weights.resize(m_Nq);
160  const double J_det = det(J);
161  for (unsigned int q = 0; q < m_Nq; q++) {
162  m_weights[q] = J_det * W[q];
163  }
164 }
165 
166 Element2::Element2(const IceGrid &grid, int Nq, int n_chi, int block_size)
167  : Element(grid, Nq, n_chi, block_size) {
168  // empty
169 }
170 
171 Element2::Element2(const DMDALocalInfo &grid_info, int Nq, int n_chi, int block_size)
172  : Element(grid_info, Nq, n_chi, block_size) {
173  // empty
174 }
175 
176 void Element2::nodal_values(const IceModelVec2Int &x_global, int *result) const {
177  for (unsigned int k = 0; k < m_n_chi; ++k) {
178  const int
179  ii = m_i + m_i_offset[k],
180  jj = m_j + m_j_offset[k];
181  result[k] = x_global.as_int(ii, jj);
182  }
183 }
184 
185 /*!@brief Initialize the Element to element (`i`, `j`) for the purposes of inserting into
186  global residual and Jacobian arrays. */
187 void Element2::reset(int i, int j) {
188  m_i = i;
189  m_j = j;
190 
191  for (unsigned int k = 0; k < m_n_chi; ++k) {
192  m_col[k].i = i + m_i_offset[k];
193  m_col[k].j = j + m_j_offset[k];
194  m_col[k].k = 0;
195  m_col[k].c = 0;
196  }
197  m_row = m_col;
198 
199  // We never sum into rows that are not owned by the local rank.
200  for (unsigned int k = 0; k < m_n_chi; k++) {
201  int pism_i = 0, pism_j = 0;
202  local_to_global(k, pism_i, pism_j);
203  if (pism_i < m_grid.xs or m_grid.xs + m_grid.xm - 1 < pism_i or
204  pism_j < m_grid.ys or m_grid.ys + m_grid.ym - 1 < pism_j) {
206  }
207  }
208 }
209 
210 /*!@brief Mark that the row corresponding to local degree of freedom `k` should not be updated
211  when inserting into the global residual or Jacobian arrays. */
212 void Element::mark_row_invalid(unsigned int k) {
213  m_row[k].i = m_row[k].j = m_row[k].k = m_invalid_dof;
214 }
215 
216 /*!@brief Mark that the column corresponding to local degree of freedom `k` should not be updated
217  when inserting into the global Jacobian arrays. */
218 void Element::mark_col_invalid(unsigned int k) {
219  m_col[k].i = m_col[k].j = m_col[k].k = m_invalid_dof;
220 }
221 
222 //! Add the contributions of an element-local Jacobian to the global Jacobian matrix.
223 /*! The element-local Jacobian should be given as a row-major array of
224  * Nk*Nk values in the scalar case or (2Nk)*(2Nk) values in the
225  * vector valued case.
226  *
227  * Note that MatSetValuesBlockedStencil ignores negative indexes, so
228  * values in K corresponding to locations marked using
229  * mark_row_invalid() and mark_col_invalid() are ignored. (Just as they
230  * should be.)
231  */
232 void Element::add_contribution(const double *K, Mat J) const {
233  PetscErrorCode ierr = MatSetValuesBlockedStencil(J,
234  m_block_size, m_row.data(),
235  m_block_size, m_col.data(),
236  K, ADD_VALUES);
237  PISM_CHK(ierr, "MatSetValuesBlockedStencil");
238 }
239 
240 Q1Element2::Q1Element2(const IceGrid &grid, const Quadrature &quadrature)
241  : Element2(grid, quadrature.weights().size(), q1::n_chi, q1::n_chi) {
242 
243  double dx = grid.dx();
244  double dy = grid.dy();
245 
246  m_i_offset = {0, 1, 1, 0};
247  m_j_offset = {0, 0, 1, 1};
248 
249  // south, east, north, west
250  m_normals = {{0.0, -1.0}, {1.0, 0.0}, {0.0, 1.0}, {-1.0, 0.0}};
251 
252  m_side_lengths = {dx, dy, dx, dy};
253 
254  // compute the Jacobian
255 
256  double J[3][3];
258  J[0][0] = 0.5 * dx;
259  J[0][1] = 0.0;
260  J[1][0] = 0.0;
261  J[1][1] = 0.5 * dy;
262 
263  // initialize germs and quadrature weights for the quadrature on this physical element
264  initialize(J, q1::chi, q1::n_chi, quadrature.points(), quadrature.weights());
265  reset(0, 0);
266 }
267 
268 Q1Element2::Q1Element2(const DMDALocalInfo &grid_info,
269  double dx,
270  double dy,
271  const Quadrature &quadrature)
272  : Element2(grid_info, quadrature.weights().size(), q1::n_chi, q1::n_chi) {
273 
274  m_i_offset = {0, 1, 1, 0};
275  m_j_offset = {0, 0, 1, 1};
276 
277  // south, east, north, west
278  m_normals = {{0.0, -1.0}, {1.0, 0.0}, {0.0, 1.0}, {-1.0, 0.0}};
279 
280  m_side_lengths = {dx, dy, dx, dy};
281 
282  // compute the Jacobian
283 
284  double J[3][3];
286  J[0][0] = 0.5 * dx;
287  J[0][1] = 0.0;
288  J[1][0] = 0.0;
289  J[1][1] = 0.5 * dy;
290 
291  // initialize germs and quadrature weights for the quadrature on this physical element
292  initialize(J, q1::chi, q1::n_chi, quadrature.points(), quadrature.weights());
293  reset(0, 0);
294 }
295 
296 P1Element2::P1Element2(const IceGrid &grid, const Quadrature &quadrature, int type)
297  : Element2(grid, quadrature.weights().size(), p1::n_chi, q1::n_chi) {
298 
299  double dx = grid.dx();
300  double dy = grid.dy();
301 
302  // outward pointing normals for all sides of a Q1 element with sides aligned with X and
303  // Y axes
304  const Vector2
305  n01( 0.0, -1.0), // south
306  n12( 1.0, 0.0), // east
307  n23( 0.0, 1.0), // north
308  n30(-1.0, 0.0); // west
309 
310  // "diagonal" sides
311  Vector2
312  n13( 1.0, dx / dy), // 1-3 diagonal, outward for element 0
313  n20(-1.0, dx / dy); // 2-0 diagonal, outward for element 1
314 
315  // normalize
316  n13 /= n13.magnitude();
317  n20 /= n20.magnitude();
318 
319  // coordinates of nodes of a Q1 element this P1 element is embedded in (up to
320  // translation)
321  Vector2 p[] = {{0, 0}, {dx, 0}, {dx, dy}, {0, dy}};
322  std::vector<Vector2> pts;
323 
324  switch (type) {
325  case 0:
326  m_i_offset = {0, 1, 0};
327  m_j_offset = {0, 0, 1};
328  m_normals = {n01, n13, n30};
329  pts = {p[0], p[1], p[3]};
330  break;
331  case 1:
332  m_i_offset = {1, 1, 0};
333  m_j_offset = {0, 1, 0};
334  m_normals = {n12, n20, n01};
335  pts = {p[1], p[2], p[0]};
336  break;
337  case 2:
338  m_i_offset = {1, 0, 1};
339  m_j_offset = {1, 1, 0};
340  m_normals = {n23, -1.0 * n13, n12};
341  pts = {p[2], p[3], p[1]};
342  break;
343  case 3:
344  default:
345  m_i_offset = {0, 0, 1};
346  m_j_offset = {1, 0, 1};
347  m_normals = {n30, -1.0 * n20, n23};
348  pts = {p[3], p[0], p[2]};
349  break;
350  }
351 
352  reset(0, 0);
353  // make sure add_contribution() ignores entries corresponding to the 4-th node
354  mark_row_invalid(3);
355  mark_col_invalid(3);
356 
357  m_side_lengths.resize(n_sides());
358  // compute side lengths
359  for (unsigned int k = 0; k < n_sides(); ++k) {
360  m_side_lengths[k] = (pts[k] - pts[(k + 1) % 3]).magnitude();
361  }
362 
363  // compute the Jacobian
364  double J[3][3];
366  J[0][0] = pts[1].u - pts[0].u;
367  J[0][1] = pts[1].v - pts[0].v;
368  J[1][0] = pts[2].u - pts[0].u;
369  J[1][1] = pts[2].v - pts[0].v;
370 
371  // initialize germs and quadrature weights for the quadrature on this physical element
372  initialize(J, p1::chi, p1::n_chi, quadrature.points(), quadrature.weights());
373  reset(0, 0);
374 }
375 
376 Element3::Element3(const DMDALocalInfo &grid_info, int Nq, int n_chi, int block_size)
377  : Element(grid_info, Nq, n_chi, block_size) {
378  m_i = 0;
379  m_j = 0;
380  m_k = 0;
381 }
382 
383 Element3::Element3(const IceGrid &grid, int Nq, int n_chi, int block_size)
384  : Element(grid, Nq, n_chi, block_size) {
385  m_i = 0;
386  m_j = 0;
387  m_k = 0;
388 }
389 
390 Q1Element3::Q1Element3(const DMDALocalInfo &grid_info,
391  const Quadrature &quadrature,
392  double dx,
393  double dy,
394  double x_min,
395  double y_min)
396  : Element3(grid_info, quadrature.weights().size(), q13d::n_chi, q13d::n_chi),
397  m_dx(dx),
398  m_dy(dy),
399  m_x_min(x_min),
400  m_y_min(y_min),
401  m_points(quadrature.points()),
402  m_w(quadrature.weights()) {
403 
404  m_weights.resize(m_Nq);
405 
406  m_i_offset = {0, 1, 1, 0, 0, 1, 1, 0};
407  m_j_offset = {0, 0, 1, 1, 0, 0, 1, 1};
408  m_k_offset = {0, 0, 0, 0, 1, 1, 1, 1};
409 
410  // store values of shape functions on the reference element
411  m_chi.resize(m_Nq * m_n_chi);
412  for (unsigned int q = 0; q < m_Nq; q++) {
413  for (unsigned int n = 0; n < m_n_chi; n++) {
414  m_chi[q * m_n_chi + n] = q13d::chi(n, m_points[q]);
415  }
416  }
417  m_germs = m_chi;
418 }
419 
420 Q1Element3::Q1Element3(const IceGrid &grid, const Quadrature &quadrature)
421  : Element3(grid, quadrature.weights().size(), q13d::n_chi, q13d::n_chi),
422  m_dx(grid.dx()),
423  m_dy(grid.dy()),
424  m_points(quadrature.points()),
425  m_w(quadrature.weights()) {
426 
427  m_weights.resize(m_Nq);
428 
429  m_i_offset = {0, 1, 1, 0, 0, 1, 1, 0};
430  m_j_offset = {0, 0, 1, 1, 0, 0, 1, 1};
431  m_k_offset = {0, 0, 0, 0, 1, 1, 1, 1};
432 
433  // store values of shape functions on the reference element
434  m_chi.resize(m_Nq * m_n_chi);
435  for (unsigned int q = 0; q < m_Nq; q++) {
436  for (unsigned int n = 0; n < m_n_chi; n++) {
437  m_chi[q * m_n_chi + n] = q13d::chi(n, m_points[q]);
438  }
439  }
440 }
441 
442 
443 /*! Initialize the element `i,j,k`.
444  *
445  *
446  * @param[in] i i-index of the lower left node
447  * @param[in] j j-index of the lower left node
448  * @param[in] k k-index of the lower left node
449  * @param[in] z z-coordinates of the nodes of this element
450  */
451 void Q1Element3::reset(int i, int j, int k, const double *z) {
452  // Record i,j,k corresponding to the current element:
453  m_i = i;
454  m_j = j;
455  m_k = k;
456 
457  // store z coordinates
458  for (unsigned int n = 0; n < m_n_chi; ++n) {
459  m_z_nodal[n] = z[n];
460  }
461 
462  // Set row and column info used to add contributions:
463  for (unsigned int n = 0; n < m_n_chi; ++n) {
464  m_col[n].i = k + m_k_offset[n]; // x -- z (STORAGE_ORDER)
465  m_col[n].j = i + m_i_offset[n]; // y -- x (STORAGE_ORDER)
466  m_col[n].k = j + m_j_offset[n]; // z -- y (STORAGE_ORDER)
467  m_col[n].c = 0;
468  }
469  m_row = m_col;
470 
471  // Mark rows that we don't own as invalid:
472  for (unsigned int n = 0; n < m_n_chi; n++) {
473  auto I = local_to_global(n);
474  if (I.i < m_grid.xs or m_grid.xs + m_grid.xm - 1 < I.i or
475  I.j < m_grid.ys or m_grid.ys + m_grid.ym - 1 < I.j) {
477  }
478  }
479 
480  // Compute J^{-1} and use it to compute m_germs and m_weights:
481  for (unsigned int q = 0; q < m_Nq; q++) {
482 
483  Vector3 dz{0.0, 0.0, 0.0};
484  for (unsigned int n = 0; n < m_n_chi; ++n) {
485  auto &chi = m_chi[q * m_n_chi + n];
486  dz.x += chi.dx * z[n];
487  dz.y += chi.dy * z[n];
488  dz.z += chi.dz * z[n];
489  }
490 
491  double J[3][3] = {{m_dx / 2.0, 0.0, dz.x},
492  { 0.0, m_dy / 2.0, dz.y},
493  { 0.0, 0.0, dz.z}};
494 
495  double J_det = J[0][0] * J[1][1] * J[2][2];
496 
497  assert(J_det != 0.0);
498 
499  double J_inv[3][3] = {{1.0 / J[0][0], 0.0, -J[0][2] / (J[0][0] * J[2][2])},
500  {0.0, 1.0 / J[1][1], -J[1][2] / (J[1][1] * J[2][2])},
501  {0.0, 0.0, 1.0 / J[2][2]}};
502 
503  m_weights[q] = J_det * m_w[q];
504 
505  for (unsigned int n = 0; n < m_n_chi; n++) {
506  auto &chi = m_chi[q * m_n_chi + n];
507  // FIXME: I should be able to use multiply() defined above, but there must be a bug
508  // there...
509  m_germs[q * m_n_chi + n] = {chi.val,
510  J_inv[0][0] * chi.dx + J_inv[0][1] * chi.dy + J_inv[0][2] * chi.dz,
511  J_inv[1][0] * chi.dx + J_inv[1][1] * chi.dy + J_inv[1][2] * chi.dz,
512  J_inv[2][0] * chi.dx + J_inv[2][1] * chi.dy + J_inv[2][2] * chi.dz};
513  }
514  }
515 }
516 
517 Q1Element3Face::Q1Element3Face(double dx, double dy, const Quadrature &quadrature)
518  : m_dx(dx),
519  m_dy(dy),
520  m_points(quadrature.points()),
521  m_w(quadrature.weights()),
522  m_n_chi(q13d::n_chi),
523  m_Nq(m_w.size()) {
524 
525  // Note: here I set m_n_chi to q13d::n_chi (8) while each face has only four basis
526  // functions that are not zero on it. One could make evaluate() cheaper by omitting
527  // these zeros. (This would make the code somewhat more complex.)
528 
529  m_chi.resize(m_Nq * m_n_chi);
530  m_weights.resize(m_Nq);
531  m_normals.resize(m_Nq);
532 }
533 
534 void Q1Element3Face::reset(int face, const double *z) {
535 
536  // Turn coordinates of a 2D quadrature point on [-1,1]*[-1,1] into coordinates on a face
537  // of the cube [-1,1]*[-1,1]*[-1,1].
538  for (unsigned int q = 0; q < m_Nq; q++) {
539  auto pt = m_points[q];
540  QuadPoint P;
541 
542  // This expresses parameterizations of faces of the reference element: for example,
543  // face 0 is parameterized by (s, t) -> [-1, s, t].
544  switch (face) {
546  P = {-1.0, pt.xi, pt.eta};
547  break;
549  P = { 1.0, pt.xi, pt.eta};
550  break;
552  P = {pt.xi, -1.0, pt.eta};
553  break;
555  P = {pt.xi, 1.0, pt.eta};
556  break;
558  P = {pt.xi, pt.eta, -1.0};
559  break;
560  case fem::q13d::FACE_TOP:
561  P = {pt.xi, pt.eta, 1.0};
562  break;
563  }
564 
565  // Compute dz/dxi, dz/deta and dz/dzeta.
566  Vector3 dz{0.0, 0.0, 0.0};
567  for (unsigned int n = 0; n < m_n_chi; ++n) {
568  // Note: chi(n, point) for a particular face does not depend on the element geometry
569  // and could be computed in advance (in the constructor). On the other hand, I
570  // expect that the number of faces in a Neumann boundary is relatively small and
571  // this optimization is not likely to pay off.
572  auto chi = q13d::chi(n, P);
573 
574  m_chi[q * m_n_chi + n] = chi.val;
575 
576  dz.x += chi.dx * z[n];
577  dz.y += chi.dy * z[n];
578  dz.z += chi.dz * z[n];
579  }
580 
581  // Use the magnitude of the normal to a face to turn quadrature weights on
582  // [-1,1]*[-1,1] into quadrature weights on a face of this physical element.
583 
584  m_weights[q] = m_w[q];
585 
586  // Sign (-1 or +1) used to set normal orientation This relies of the order of faces in
587  // fem::q13d::incident_nodes and above (see the code setting QuadPoint P).
588  double sign = 2 * (face % 2) - 1;
589 
590  switch (face) {
593  m_weights[q] *= 0.5 * m_dy * dz.z;
594  m_normals[q] = {sign, 0.0, 0.0};
595  break;
598  m_weights[q] *= 0.5 * m_dx * dz.z;
599  m_normals[q] = {0.0, sign, 0.0};
600  break;
602  case fem::q13d::FACE_TOP:
603  {
604  double
605  a = 0.5 * m_dy * dz.x,
606  b = 0.5 * m_dx * dz.y,
607  c = 0.25 * m_dx * m_dy,
608  M = std::sqrt(a * a + b * b + c * c);
609  m_weights[q] *= M;
610 
611  M *= sign;
612  m_normals[q] = {a / M, b / M, c / M};
613  }
614  break;
615  }
616  } // end of the loop over quadrature points
617 }
618 
619 
620 } // end of namespace fem
621 } // end of namespace pism
double dx() const
Horizontal grid spacing.
Definition: IceGrid.cc:983
double dy() const
Horizontal grid spacing.
Definition: IceGrid.cc:988
std::shared_ptr< petsc::DM > get_dm(int dm_dof, int stencil_width) const
Get a PETSc DM ("distributed array manager") object for given dof (number of degrees of freedom per g...
Definition: IceGrid.cc:841
Describes the PISM grid and the distribution of data across processors.
Definition: IceGrid.hh:228
int as_int(int i, int j) const
A simple class "hiding" the fact that the mask is stored as floating-point scalars (instead of intege...
Definition: iceModelVec.hh:389
double magnitude() const
Magnitude.
Definition: Vector2.hh:40
This class represents a 2D vector field (such as ice velocity) at a certain grid point.
Definition: Vector2.hh:29
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
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
std::vector< Vector2 > m_normals
Definition: Element.hh:251
Element2(const IceGrid &grid, int Nq, int n_chi, int block_size)
Definition: Element.cc:166
unsigned int n_sides() const
Definition: Element.hh:181
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
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
std::vector< int > m_k_offset
Definition: Element.hh:342
GlobalIndex local_to_global(int i, int j, int k, unsigned int n) const
Definition: Element.hh:298
const Germ & chi(unsigned int q, unsigned int k) const
Definition: Element.hh:73
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
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
Q1Element2(const IceGrid &grid, const Quadrature &quadrature)
Definition: Element.cc:240
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
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
std::vector< QuadPoint > m_points
Definition: Element.hh:445
std::vector< double > m_weights
Definition: Element.hh:453
const unsigned int m_n_chi
Definition: Element.hh:456
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
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
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
std::vector< Germ > m_chi
Definition: Element.hh:389
const std::vector< double > & weights() const
Definition: Quadrature.cc:29
const std::vector< QuadPoint > & points() const
Definition: Quadrature.cc:25
Numerical integration of finite element functions.
Definition: Quadrature.hh:62
Contains parameters of an input file grid.
Definition: IceGrid.hh:62
#define PISM_CHK(errcode, name)
#define n
Definition: exactTestM.c:37
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
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 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
static double det(const double a[3][3])
Determinant of a 3x3 matrix.
Definition: Element.cc:32
static double dot(const Vector3 &v, const Germ &a)
Definition: Element.cc:56
static void invert(const double A[3][3], double result[3][3])
Invert a 3x3 matrix.
Definition: Element.cc:61
static Vector3 row(const double A[3][3], size_t k)
Definition: Element.cc:46
static void set_to_identity(double A[3][3])
Definition: Element.cc:96
static Germ multiply(const double A[3][3], const Germ &v)
Compute derivatives with respect to x,y using J^{-1} and derivatives with respect to xi,...
Definition: Element.cc:88
static Vector3 column(const double A[3][3], size_t k)
Definition: Element.cc:51
static Vector3 cross(const Vector3 &a, const Vector3 &b)
Cross product of two 3D vectors.
Definition: Element.cc:39
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