tFETools.cc - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tFETools.cc (25284B)
---
1 // Copyright (C) 2009--2011, 2013, 2014, 2015, 2016, 2017, 2018 Jed Brown and 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 // Utility functions used by the SSAFEM code.
20
21 #include <cassert> // assert
22 #include <cstring> // memset
23 #include <cstdlib> // malloc, free
24
25 #include "FETools.hh"
26 #include "pism/util/IceGrid.hh"
27 #include "pism/util/iceModelVec.hh"
28 #include "pism/util/pism_utilities.hh"
29 #include "pism/util/Logger.hh"
30
31 #include "pism/util/error_handling.hh"
32
33 namespace pism {
34
35 //! FEM (Finite Element Method) utilities
36 namespace fem {
37
38 namespace q0 {
39 /*!
40 * Piecewise-constant shape functions.
41 */
42 Germ chi(unsigned int k, const QuadPoint &pt) {
43 assert(k < q0::n_chi);
44
45 Germ result;
46
47 if ((k == 0 and pt.xi <= 0.0 and pt.eta <= 0.0) or
48 (k == 1 and pt.xi > 0.0 and pt.eta <= 0.0) or
49 (k == 2 and pt.xi > 0.0 and pt.eta > 0.0) or
50 (k == 3 and pt.xi <= 0.0 and pt.eta > 0.0)) {
51 result.val = 1.0;
52 } else {
53 result.val = 0.0;
54 }
55
56 result.dx = 0.0;
57 result.dy = 0.0;
58
59 return result;
60 }
61 } // end of namespace q0
62
63 namespace q1 {
64
65 //! Q1 basis functions on the reference element with nodes (-1,-1), (1,-1), (1,1), (-1,1).
66 Germ chi(unsigned int k, const QuadPoint &pt) {
67 assert(k < q1::n_chi);
68
69 const double xi[] = {-1.0, 1.0, 1.0, -1.0};
70 const double eta[] = {-1.0, -1.0, 1.0, 1.0};
71
72 Germ result;
73
74 result.val = 0.25 * (1.0 + xi[k] * pt.xi) * (1.0 + eta[k] * pt.eta);
75 result.dx = 0.25 * xi[k] * (1.0 + eta[k] * pt.eta);
76 result.dy = 0.25 * eta[k] * (1.0 + xi[k] * pt.xi);
77
78 return result;
79 }
80
81
82 const Vector2* outward_normals() {
83 static const Vector2 n[] = {Vector2( 0.0, -1.0), // south
84 Vector2( 1.0, 0.0), // east
85 Vector2( 0.0, 1.0), // north
86 Vector2(-1.0, 0.0)}; // west
87 return n;
88 }
89
90 } // end of namespace q1
91
92 namespace p1 {
93
94 //! P1 basis functions on the reference element with nodes (0,0), (1,0), (0,1).
95 Germ chi(unsigned int k, const QuadPoint &pt) {
96 assert(k < q1::n_chi);
97 Germ result;
98
99 switch (k) {
100 case 0:
101 result.val = 1.0 - pt.xi - pt.eta;
102 result.dx = -1.0;
103 result.dy = -1.0;
104 break;
105 case 1:
106 result.val = pt.xi;
107 result.dx = 1.0;
108 result.dy = 0.0;
109 break;
110 case 2:
111 result.val = pt.eta;
112 result.dx = 0.0;
113 result.dy = 1.0;
114 break;
115 default: // the fourth (dummy) basis function
116 result.val = 0.0;
117 result.dx = 0.0;
118 result.dy = 0.0;
119 }
120 return result;
121 }
122 } // end of namespace p1
123
124
125 ElementIterator::ElementIterator(const IceGrid &g) {
126 // Start by assuming ghost elements exist in all directions.
127 // Elements are indexed by their lower left vertex. If there is a ghost
128 // element on the right, its i-index will be the same as the maximum
129 // i-index of a non-ghost vertex in the local grid.
130 xs = g.xs() - 1; // Start at ghost to the left.
131 int xf = g.xs() + g.xm() - 1; // End at ghost to the right.
132 ys = g.ys() - 1; // Start at ghost at the bottom.
133 int yf = g.ys() + g.ym() - 1; // End at ghost at the top.
134
135 lxs = g.xs();
136 int lxf = lxs + g.xm() - 1;
137 lys = g.ys();
138 int lyf = lys + g.ym() - 1;
139
140 // Now correct if needed. The only way there will not be ghosts is if the
141 // grid is not periodic and we are up against the grid boundary.
142
143 if (!(g.periodicity() & X_PERIODIC)) {
144 // Leftmost element has x-index 0.
145 if (xs < 0) {
146 xs = 0;
147 }
148 // Rightmost vertex has index g.Mx-1, so the rightmost element has index g.Mx-2
149 if (xf > (int)g.Mx() - 2) {
150 xf = g.Mx() - 2;
151 lxf = g.Mx() - 2;
152 }
153 }
154
155 if (!(g.periodicity() & Y_PERIODIC)) {
156 // Bottom element has y-index 0.
157 if (ys < 0) {
158 ys = 0;
159 }
160 // Topmost vertex has index g.My - 1, so the topmost element has index g.My - 2
161 if (yf > (int)g.My() - 2) {
162 yf = g.My() - 2;
163 lyf = g.My() - 2;
164 }
165 }
166
167 // Tally up the number of elements in each direction
168 xm = xf - xs + 1;
169 ym = yf - ys + 1;
170 lxm = lxf - lxs + 1;
171 lym = lyf - lys + 1;
172 }
173
174 ElementMap::ElementMap(const IceGrid &grid)
175 : m_grid(grid) {
176 reset(0, 0);
177 }
178
179 ElementMap::~ElementMap() {
180 // empty
181 }
182
183 void ElementMap::nodal_values(const IceModelVec2Int &x_global, int *result) const {
184 for (unsigned int k = 0; k < q1::n_chi; ++k) {
185 const int
186 ii = m_i + m_i_offset[k],
187 jj = m_j + m_j_offset[k];
188 result[k] = x_global.as_int(ii, jj);
189 }
190 }
191
192 /*!@brief Initialize the ElementMap to element (`i`, `j`) for the purposes of inserting into
193 global residual and Jacobian arrays. */
194 void ElementMap::reset(int i, int j) {
195 m_i = i;
196 m_j = j;
197
198 for (unsigned int k = 0; k < fem::q1::n_chi; ++k) {
199 m_col[k].i = i + m_i_offset[k];
200 m_col[k].j = j + m_j_offset[k];
201 m_col[k].k = 0;
202
203 m_row[k].i = m_col[k].i;
204 m_row[k].j = m_col[k].j;
205 m_row[k].k = m_col[k].k;
206 }
207
208 // We do not ever sum into rows that are not owned by the local rank.
209 for (unsigned int k = 0; k < fem::q1::n_chi; k++) {
210 int pism_i = m_row[k].i, pism_j = m_row[k].j;
211 if (pism_i < m_grid.xs() || m_grid.xs() + m_grid.xm() - 1 < pism_i ||
212 pism_j < m_grid.ys() || m_grid.ys() + m_grid.ym() - 1 < pism_j) {
213 mark_row_invalid(k);
214 }
215 }
216 }
217
218 /*!@brief Mark that the row corresponding to local degree of freedom `k` should not be updated
219 when inserting into the global residual or Jacobian arrays. */
220 void ElementMap::mark_row_invalid(int k) {
221 m_row[k].i = m_row[k].j = m_invalid_dof;
222 // We are solving a 2D system, so MatStencil::k is not used. Here we
223 // use it to mark invalid rows.
224 m_row[k].k = 1;
225 }
226
227 /*!@brief Mark that the column corresponding to local degree of freedom `k` should not be updated
228 when inserting into the global Jacobian arrays. */
229 void ElementMap::mark_col_invalid(int k) {
230 m_col[k].i = m_col[k].j = m_invalid_dof;
231 // We are solving a 2D system, so MatStencil::k is not used. Here we
232 // use it to mark invalid columns.
233 m_col[k].k = 1;
234 }
235
236 //! Add the contributions of an element-local Jacobian to the global Jacobian vector.
237 /*! The element-local Jacobian should be given as a row-major array of
238 * Nk*Nk values in the scalar case or (2Nk)*(2Nk) values in the
239 * vector valued case.
240 *
241 * Note that MatSetValuesBlockedStencil ignores negative indexes, so
242 * values in K corresponding to locations marked using
243 * mark_row_invalid() and mark_col_invalid() are ignored. (Just as they
244 * should be.)
245 */
246 void ElementMap::add_contribution(const double *K, Mat J) const {
247 PetscErrorCode ierr = MatSetValuesBlockedStencil(J,
248 fem::q1::n_chi, m_row,
249 fem::q1::n_chi, m_col,
250 K, ADD_VALUES);
251 PISM_CHK(ierr, "MatSetValuesBlockedStencil");
252 }
253
254 const int ElementMap::m_i_offset[4] = {0, 1, 1, 0};
255 const int ElementMap::m_j_offset[4] = {0, 0, 1, 1};
256
257 Quadrature::Quadrature(unsigned int N)
258 : m_Nq(N) {
259
260 m_W = (double*) malloc(m_Nq * sizeof(double));
261 if (m_W == NULL) {
262 throw std::runtime_error("Failed to allocate a Quadrature instance");
263 }
264
265 m_germs = (Germs*) malloc(m_Nq * q1::n_chi * sizeof(Germ));
266 if (m_germs == NULL) {
267 free(m_W);
268 throw std::runtime_error("Failed to allocate a Quadrature instance");
269 }
270 }
271
272 Germ Quadrature::test_function_values(unsigned int q, unsigned int k) const {
273 return m_germs[q][k];
274 }
275
276 double Quadrature::weights(unsigned int q) const {
277 return m_W[q];
278 }
279
280 UniformQxQuadrature::UniformQxQuadrature(unsigned int size, double dx, double dy, double scaling)
281 : Quadrature(size) {
282 // We use uniform Cartesian coordinates, so the Jacobian is constant and diagonal on every
283 // element.
284 //
285 // Note that the reference element is [-1,1]^2, hence the extra factor of 1/2.
286 m_J[0][0] = 0.5 * dx / scaling;
287 m_J[0][1] = 0.0;
288 m_J[1][0] = 0.0;
289 m_J[1][1] = 0.5 * dy / scaling;
290 }
291
292 Quadrature::~Quadrature() {
293 free(m_W);
294 m_W = NULL;
295
296 free(m_germs);
297 m_germs = NULL;
298 }
299
300 //! Build quadrature points and weights for a tensor product quadrature based on a 1D quadrature
301 //! rule. Uses the same 1D quadrature in both directions.
302 /**
303 @param[in] n 1D quadrature size (the resulting quadrature has size n*n)
304 @param[in] points1 1D quadrature points
305 @param[in] weights1 1D quadrature weights
306 @param[out] points resulting 2D quadrature points
307 @param[out] weights resulting 2D quadrature weights
308 */
309 static void tensor_product_quadrature(unsigned int n,
310 const double *points1,
311 const double *weights1,
312 QuadPoint *points,
313 double *weights) {
314 unsigned int q = 0;
315 for (unsigned int j = 0; j < n; ++j) {
316 for (unsigned int i = 0; i < n; ++i) {
317 points[q].xi = points1[i];
318 points[q].eta = points1[j];
319
320 weights[q] = weights1[i] * weights1[j];
321
322 ++q;
323 }
324 }
325 }
326
327 //! Determinant of a square matrix of size 2.
328 static double determinant(const double J[2][2]) {
329 return J[0][0] * J[1][1] - J[1][0] * J[0][1];
330 }
331
332 //! Compute the inverse of a two by two matrix.
333 static void invert(const double A[2][2], double A_inv[2][2]) {
334 const double det_A = determinant(A);
335
336 assert(det_A != 0.0);
337
338 A_inv[0][0] = A[1][1] / det_A;
339 A_inv[0][1] = -A[0][1] / det_A;
340 A_inv[1][0] = -A[1][0] / det_A;
341 A_inv[1][1] = A[0][0] / det_A;
342 }
343
344 //! Compute derivatives with respect to x,y using J^{-1} and derivatives with respect to xi, eta.
345 static Germ apply_jacobian_inverse(const double J_inv[2][2], const Germ &f) {
346 Germ result;
347 result.val = f.val;
348 result.dx = f.dx * J_inv[0][0] + f.dy * J_inv[0][1];
349 result.dy = f.dx * J_inv[1][0] + f.dy * J_inv[1][1];
350 return result;
351 }
352
353 //! Two-by-two Gaussian quadrature on a rectangle.
354 Q1Quadrature4::Q1Quadrature4(double dx, double dy, double L)
355 : UniformQxQuadrature(m_size, dx, dy, L) {
356
357 // coordinates and weights of the 2-point 1D Gaussian quadrature
358 const double
359 A = 1.0 / sqrt(3.0),
360 points2[2] = {-A, A},
361 weights2[2] = {1.0, 1.0};
362
363 QuadPoint points[m_size];
364 double W[m_size];
365
366 tensor_product_quadrature(2, points2, weights2, points, W);
367
368 initialize(q1::chi, q1::n_chi, points, W);
369 }
370
371 Q1Quadrature9::Q1Quadrature9(double dx, double dy, double L)
372 : UniformQxQuadrature(m_size, dx, dy, L) {
373 // The quadrature points on the reference square.
374
375 const double
376 A = 0.0,
377 B = sqrt(0.6),
378 points3[3] = {-B, A, B};
379
380 const double
381 w1 = 5.0 / 9.0,
382 w2 = 8.0 / 9.0,
383 weights3[3] = {w1, w2, w1};
384
385 QuadPoint points[m_size];
386 double W[m_size];
387
388 tensor_product_quadrature(3, points3, weights3, points, W);
389
390 initialize(q1::chi, q1::n_chi, points, W);
391 }
392
393 Q1Quadrature16::Q1Quadrature16(double dx, double dy, double L)
394 : UniformQxQuadrature(m_size, dx, dy, L) {
395 // The quadrature points on the reference square.
396 const double
397 A = sqrt(3.0 / 7.0 - (2.0 / 7.0) * sqrt(6.0 / 5.0)), // smaller magnitude
398 B = sqrt(3.0 / 7.0 + (2.0 / 7.0) * sqrt(6.0 / 5.0)), // larger magnitude
399 points4[4] = {-B, -A, A, B};
400
401 // The weights w_i for Gaussian quadrature on the reference element with these quadrature points
402 const double
403 w1 = (18.0 + sqrt(30.0)) / 36.0, // larger
404 w2 = (18.0 - sqrt(30.0)) / 36.0, // smaller
405 weights4[4] = {w2, w1, w1, w2};
406
407 QuadPoint points[m_size];
408 double W[m_size];
409
410 tensor_product_quadrature(4, points4, weights4, points, W);
411
412 initialize(q1::chi, q1::n_chi, points, W);
413 }
414
415
416 //! @brief 1e4-point (100x100) uniform (*not* Gaussian) quadrature for integrating discontinuous
417 //! functions.
418 Q0Quadrature1e4::Q0Quadrature1e4(double dx, double dy, double L)
419 : UniformQxQuadrature(m_size, dx, dy, L) {
420 assert(m_size_1d * m_size_1d == m_size);
421
422 double xi[m_size_1d], w[m_size_1d];
423 const double dxi = 2.0 / m_size_1d;
424 for (unsigned int k = 0; k < m_size_1d; ++k) {
425 xi[k] = -1.0 + dxi*(k + 0.5);
426 w[k] = 2.0 / m_size_1d;
427 }
428
429 QuadPoint points[m_size];
430 double W[m_size];
431
432 tensor_product_quadrature(m_size_1d, xi, w, points, W);
433
434 initialize(q0::chi, q0::n_chi, points, W);
435 }
436
437 //! @brief 1e4-point (100x100) uniform (*not* Gaussian) quadrature for integrating discontinuous
438 //! functions.
439 Q1Quadrature1e4::Q1Quadrature1e4(double dx, double dy, double L)
440 : UniformQxQuadrature(m_size, dx, dy, L) {
441 assert(m_size_1d * m_size_1d == m_size);
442
443 double xi[m_size_1d], w[m_size_1d];
444 const double dxi = 2.0 / m_size_1d;
445 for (unsigned int k = 0; k < m_size_1d; ++k) {
446 xi[k] = -1.0 + dxi*(k + 0.5);
447 w[k] = 2.0 / m_size_1d;
448 }
449
450 QuadPoint points[m_size];
451 double W[m_size];
452
453 tensor_product_quadrature(m_size_1d, xi, w, points, W);
454
455 initialize(q1::chi, q1::n_chi, points, W);
456 }
457
458 //! Initialize shape function values and weights of a 2D quadrature.
459 /** Assumes that the Jacobian does not depend on coordinates of the current quadrature point.
460 */
461 void Quadrature::initialize(ShapeFunction2 f,
462 unsigned int n_chi,
463 const QuadPoint *points,
464 const double *W) {
465
466 double J_inv[2][2];
467 invert(m_J, J_inv);
468
469 for (unsigned int q = 0; q < m_Nq; q++) {
470 for (unsigned int k = 0; k < n_chi; k++) {
471 m_germs[q][k] = apply_jacobian_inverse(J_inv, f(k, points[q]));
472 }
473 }
474
475 const double J_det = determinant(m_J);
476 for (unsigned int q = 0; q < m_Nq; q++) {
477 m_W[q] = J_det * W[q];
478 }
479 }
480
481 /** Create a quadrature on a P1 element aligned with coordinate axes and embedded in a Q1 element.
482
483 There are four possible P1 elements in a Q1 element. The argument `N` specifies which one,
484 numbering them by the node at the right angle in the "reference" element (0,0) -- (1,0) --
485 (0,1).
486 */
487 P1Quadrature::P1Quadrature(unsigned int size, unsigned int N,
488 double dx, double dy, double L)
489 : Quadrature(size) {
490
491 // Compute the Jacobian. The nodes of the selected triangle are
492 // numbered, the unused node is marked with an "X". In all triangles
493 // nodes are numbered in the counter-clockwise direction.
494 switch (N) {
495 case 0:
496 /*
497 2------X
498 | |
499 | |
500 0------1
501 */
502 m_J[0][0] = dx / L;
503 m_J[0][1] = 0.0;
504 m_J[1][0] = 0.0;
505 m_J[1][1] = dy / L;
506 break;
507 case 1:
508 /*
509 X------1
510 | |
511 | |
512 2------0
513 */
514 m_J[0][0] = 0.0;
515 m_J[0][1] = dy / L;
516 m_J[1][0] = -dx / L;
517 m_J[1][1] = 0.0;
518 break;
519 case 2:
520 /*
521 1------0
522 | |
523 | |
524 X------2
525 */
526 m_J[0][0] = -dx / L;
527 m_J[0][1] = 0.0;
528 m_J[1][0] = 0.0;
529 m_J[1][1] = -dy / L;
530 break;
531 case 3:
532 /*
533 0------2
534 | |
535 | |
536 1------X
537 */
538 m_J[0][0] = 0.0;
539 m_J[0][1] = -dy / L;
540 m_J[1][0] = dx / L;
541 m_J[1][1] = 0.0;
542 break;
543 }
544 }
545
546 //! Permute shape functions stored in `f` *in place*.
547 static void permute(const unsigned int p[q1::n_chi], Germ f[q1::n_chi]) {
548 Germ tmp[q1::n_chi];
549
550 // store permuted values to avoid overwriting f too soon
551 for (unsigned int k = 0; k < q1::n_chi; ++k) {
552 tmp[k] = f[p[k]];
553 }
554
555 // copy back into f
556 for (unsigned int k = 0; k < q1::n_chi; ++k) {
557 f[k] = tmp[k];
558 }
559 }
560
561 P1Quadrature3::P1Quadrature3(unsigned int N,
562 double dx, double dy, double L)
563 : P1Quadrature(m_size, N, dx, dy, L) {
564
565 const double
566 one_over_six = 1.0 / 6.0,
567 two_over_three = 2.0 / 3.0;
568
569 const QuadPoint points[3] = {{two_over_three, one_over_six},
570 {one_over_six, two_over_three},
571 {one_over_six, one_over_six}};
572
573 const double W[3] = {one_over_six, one_over_six, one_over_six};
574
575 // Note that we use q1::n_chi here.
576 initialize(p1::chi, q1::n_chi, points, W);
577
578 // Permute shape function values according to N, the index of this triangle in the Q1 element.
579 const unsigned int
580 X = 3, // index of the dummy shape function
581 p[4][4] = {{0, 1, X, 2},
582 {2, 0, 1, X},
583 {X, 2, 0, 1},
584 {1, X, 2, 0}};
585 for (unsigned int q = 0; q < m_size; ++q) {
586 permute(p[N], m_germs[q]);
587 }
588 }
589
590 DirichletData::DirichletData()
591 : m_indices(NULL), m_weight(1.0) {
592 for (unsigned int k = 0; k < q1::n_chi; ++k) {
593 m_indices_e[k] = 0;
594 }
595 }
596
597 DirichletData::~DirichletData() {
598 finish(NULL);
599 }
600
601 void DirichletData::init(const IceModelVec2Int *indices,
602 const IceModelVec *values,
603 double weight) {
604 m_weight = weight;
605
606 if (indices != NULL) {
607 indices->begin_access();
608 m_indices = indices;
609 }
610
611 if (values != NULL) {
612 values->begin_access();
613 }
614 }
615
616 void DirichletData::finish(const IceModelVec *values) {
617 if (m_indices != NULL) {
618 MPI_Comm com = m_indices->grid()->ctx()->com();
619 try {
620 m_indices->end_access();
621 m_indices = NULL;
622 } catch (...) {
623 handle_fatal_errors(com);
624 }
625 }
626
627 if (values != NULL) {
628 MPI_Comm com = values->grid()->ctx()->com();
629 try {
630 values->end_access();
631 } catch (...) {
632 handle_fatal_errors(com);
633 }
634 }
635 }
636
637 //! @brief Constrain `element`, i.e. ensure that quadratures do not contribute to Dirichlet nodes by marking corresponding rows and columns as "invalid".
638 void DirichletData::constrain(ElementMap &element) {
639 element.nodal_values(*m_indices, m_indices_e);
640 for (unsigned int k = 0; k < q1::n_chi; k++) {
641 if (m_indices_e[k] > 0.5) { // Dirichlet node
642 // Mark any kind of Dirichlet node as not to be touched
643 element.mark_row_invalid(k);
644 element.mark_col_invalid(k);
645 }
646 }
647 }
648
649 // Scalar version
650
651 DirichletData_Scalar::DirichletData_Scalar(const IceModelVec2Int *indices,
652 const IceModelVec2S *values,
653 double weight)
654 : m_values(values) {
655 init(indices, m_values, weight);
656 }
657
658 void DirichletData_Scalar::enforce(const ElementMap &element, double* x_nodal) {
659 assert(m_values != NULL);
660
661 element.nodal_values(*m_indices, m_indices_e);
662 for (unsigned int k = 0; k < q1::n_chi; k++) {
663 if (m_indices_e[k] > 0.5) { // Dirichlet node
664 int i = 0, j = 0;
665 element.local_to_global(k, i, j);
666 x_nodal[k] = (*m_values)(i, j);
667 }
668 }
669 }
670
671 void DirichletData_Scalar::enforce_homogeneous(const ElementMap &element, double* x_nodal) {
672 element.nodal_values(*m_indices, m_indices_e);
673 for (unsigned int k = 0; k < q1::n_chi; k++) {
674 if (m_indices_e[k] > 0.5) { // Dirichlet node
675 x_nodal[k] = 0.;
676 }
677 }
678 }
679
680 void DirichletData_Scalar::fix_residual(double const *const *const x_global, double **r_global) {
681 assert(m_values != NULL);
682
683 const IceGrid &grid = *m_indices->grid();
684
685 // For each node that we own:
686 for (Points p(grid); p; p.next()) {
687 const int i = p.i(), j = p.j();
688
689 if ((*m_indices)(i, j) > 0.5) {
690 // Enforce explicit dirichlet data.
691 r_global[j][i] = m_weight * (x_global[j][i] - (*m_values)(i,j));
692 }
693 }
694 }
695
696 void DirichletData_Scalar::fix_residual_homogeneous(double **r_global) {
697 const IceGrid &grid = *m_indices->grid();
698
699 // For each node that we own:
700 for (Points p(grid); p; p.next()) {
701 const int i = p.i(), j = p.j();
702
703 if ((*m_indices)(i, j) > 0.5) {
704 // Enforce explicit dirichlet data.
705 r_global[j][i] = 0.0;
706 }
707 }
708 }
709
710 void DirichletData_Scalar::fix_jacobian(Mat J) {
711 const IceGrid &grid = *m_indices->grid();
712
713 // Until now, the rows and columns correspoinding to Dirichlet data
714 // have not been set. We now put an identity block in for these
715 // unknowns. Note that because we have takes steps to not touching
716 // these columns previously, the symmetry of the Jacobian matrix is
717 // preserved.
718
719 const double identity = m_weight;
720 ParallelSection loop(grid.com);
721 try {
722 for (Points p(grid); p; p.next()) {
723 const int i = p.i(), j = p.j();
724
725 if ((*m_indices)(i, j) > 0.5) {
726 MatStencil row;
727 row.j = j;
728 row.i = i;
729 PetscErrorCode ierr = MatSetValuesBlockedStencil(J, 1, &row, 1, &row, &identity,
730 ADD_VALUES);
731 PISM_CHK(ierr, "MatSetValuesBlockedStencil"); // this may throw
732 }
733 }
734 } catch (...) {
735 loop.failed();
736 }
737 loop.check();
738 }
739
740 DirichletData_Scalar::~DirichletData_Scalar() {
741 finish(m_values);
742 m_values = NULL;
743 }
744
745 // Vector version
746
747 DirichletData_Vector::DirichletData_Vector(const IceModelVec2Int *indices,
748 const IceModelVec2V *values,
749 double weight)
750 : m_values(values) {
751 init(indices, m_values, weight);
752 }
753
754 void DirichletData_Vector::enforce(const ElementMap &element, Vector2* x_nodal) {
755 assert(m_values != NULL);
756
757 element.nodal_values(*m_indices, m_indices_e);
758 for (unsigned int k = 0; k < q1::n_chi; k++) {
759 if (m_indices_e[k] > 0.5) { // Dirichlet node
760 int i = 0, j = 0;
761 element.local_to_global(k, i, j);
762 x_nodal[k] = (*m_values)(i, j);
763 }
764 }
765 }
766
767 void DirichletData_Vector::enforce_homogeneous(const ElementMap &element, Vector2* x_nodal) {
768 element.nodal_values(*m_indices, m_indices_e);
769 for (unsigned int k = 0; k < q1::n_chi; k++) {
770 if (m_indices_e[k] > 0.5) { // Dirichlet node
771 x_nodal[k].u = 0.0;
772 x_nodal[k].v = 0.0;
773 }
774 }
775 }
776
777 void DirichletData_Vector::fix_residual(Vector2 const *const *const x_global, Vector2 **r_global) {
778 assert(m_values != NULL);
779
780 const IceGrid &grid = *m_indices->grid();
781
782 // For each node that we own:
783 for (Points p(grid); p; p.next()) {
784 const int i = p.i(), j = p.j();
785
786 if ((*m_indices)(i, j) > 0.5) {
787 // Enforce explicit dirichlet data.
788 r_global[j][i] = m_weight * (x_global[j][i] - (*m_values)(i, j));
789 }
790 }
791 }
792
793 void DirichletData_Vector::fix_residual_homogeneous(Vector2 **r_global) {
794 const IceGrid &grid = *m_indices->grid();
795
796 // For each node that we own:
797 for (Points p(grid); p; p.next()) {
798 const int i = p.i(), j = p.j();
799
800 if ((*m_indices)(i, j) > 0.5) {
801 // Enforce explicit dirichlet data.
802 r_global[j][i].u = 0.0;
803 r_global[j][i].v = 0.0;
804 }
805 }
806 }
807
808 void DirichletData_Vector::fix_jacobian(Mat J) {
809 const IceGrid &grid = *m_indices->grid();
810
811 // Until now, the rows and columns correspoinding to Dirichlet data
812 // have not been set. We now put an identity block in for these
813 // unknowns. Note that because we have takes steps to not touching
814 // these columns previously, the symmetry of the Jacobian matrix is
815 // preserved.
816
817 const double identity[4] = {m_weight, 0,
818 0, m_weight};
819 ParallelSection loop(grid.com);
820 try {
821 for (Points p(grid); p; p.next()) {
822 const int i = p.i(), j = p.j();
823
824 if ((*m_indices)(i, j) > 0.5) {
825 MatStencil row;
826 row.j = j;
827 row.i = i;
828 PetscErrorCode ierr = MatSetValuesBlockedStencil(J, 1, &row, 1, &row, identity,
829 ADD_VALUES);
830 PISM_CHK(ierr, "MatSetValuesBlockedStencil"); // this may throw
831 }
832 }
833 } catch (...) {
834 loop.failed();
835 }
836 loop.check();
837 }
838
839 DirichletData_Vector::~DirichletData_Vector() {
840 finish(m_values);
841 m_values = NULL;
842 }
843
844 BoundaryQuadrature2::BoundaryQuadrature2(double dx, double dy, double L) {
845
846 const double J[2][2] = {{0.5 * dx / L, 0.0},
847 {0.0, 0.5 * dy / L}};
848
849 // The inverse of the Jacobian.
850 double J_inv[2][2];
851 invert(J, J_inv);
852
853 // Note that all quadrature weights are 1.0 (and so they are implicitly included below).
854 //
855 // bottom
856 m_W[0] = J[0][0];
857 // right
858 m_W[1] = J[1][1];
859 // top
860 m_W[2] = J[0][0];
861 // left
862 m_W[3] = J[1][1];
863
864 const double C = 1.0 / sqrt(3);
865 const QuadPoint points[q1::n_sides][m_Nq] = {
866 {{ -C, -1.0}, { C, -1.0}}, // South
867 {{ 1.0, -C}, { 1.0, C}}, // East
868 {{ -C, 1.0}, { C, 1.0}}, // North
869 {{-1.0, -C}, {-1.0, C}} // West
870 };
871
872 memset(m_germs, 0, q1::n_sides*m_Nq*q1::n_chi*sizeof(Germ));
873
874 for (unsigned int side = 0; side < q1::n_sides; ++side) {
875 for (unsigned int q = 0; q < m_Nq; ++q) {
876 for (unsigned int k = 0; k < q1::n_chi; ++k) {
877 Germ phi = q1::chi(k, points[side][q]);
878
879 m_germs[side][q][k] = apply_jacobian_inverse(J_inv, phi);
880 }
881 }
882 }
883 }
884
885 } // end of namespace fem
886 } // end of namespace pism