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