Add Blatter::nodal_parameter_values()

This method isolates the code setting nodal values of 2D parameters. It can be re-implemented in a derived class to add verification tests that use periodic boundary conditions.
pism · Oct 14, 2020 · 0306ea6 · 0306ea6
1 parent fd9d6d5
commit 0306ea6
Show file tree

Hide file tree

Showing 4 changed files with 57 additions and 20 deletions.
diff --git a/src/stressbalance/blatter/Blatter.cc b/src/stressbalance/blatter/Blatter.cc
@@ -543,6 +543,40 @@ void Blatter::init_ice_hardness(const Inputs &inputs) {
   } // end of the loop over grid points
 }
 
+/*!
+ * Get values of 2D parameters at element nodes.
+ *
+ * This method is re-implemented by derived classes that use periodic boundary conditions.
+ */
+void Blatter::nodal_parameter_values(const fem::Q1Element3 &element,
+                                     Parameters **P,
+                                     int i,
+                                     int j,
+                                     int *node_type,
+                                     double *bottom_elevation,
+                                     double *ice_thickness,
+                                     double *surface_elevation,
+                                     double *sea_level) const {
+
+  for (int n = 0; n < fem::q13d::n_chi; ++n) {
+    auto I = element.local_to_global(i, j, 0, n);
+
+    auto p = P[I.j][I.i];
+
+    node_type[n]         = p.node_type;
+    bottom_elevation[n]  = p.bed;
+    ice_thickness[n]     = p.thickness;
+
+    if (surface_elevation) {
+      surface_elevation[n] = p.bed + p.thickness;
+    }
+
+    if (sea_level) {
+      sea_level[n]         = p.sea_level;
+    }
+  }
+}
+
 Blatter::~Blatter() {
   // empty
 }

diff --git a/src/stressbalance/blatter/Blatter.hh b/src/stressbalance/blatter/Blatter.hh
@@ -105,6 +105,16 @@ protected:
 
   bool partially_submerged_face(int face, const double *z, const double *sea_level);
 
+  virtual void nodal_parameter_values(const fem::Q1Element3 &element,
+                                      Parameters **P,
+                                      int i,
+                                      int j,
+                                      int *node_type,
+                                      double *bottom,
+                                      double *thickness,
+                                      double *surface,
+                                      double *sea_level) const;
+
   virtual bool neumann_bc_face(int face, const int *node_type);
 
   virtual bool dirichlet_node(const DMDALocalInfo &info, const fem::Element3::GlobalIndex& I);

diff --git a/src/stressbalance/blatter/jacobian.cc b/src/stressbalance/blatter/jacobian.cc
@@ -237,15 +237,13 @@ void Blatter::compute_jacobian(DMDALocalInfo *petsc_info,
   for (int j = info.gys; j < info.gys + info.gym - 1; j++) {
     for (int i = info.gxs; i < info.gxs + info.gxm - 1; i++) {
 
-      for (int n = 0; n < Nk; ++n) {
-        auto I = element.local_to_global(i, j, 0, n);
-
-        auto p = P[I.j][I.i];
-
-        bottom_elevation[n] = p.bed;
-        ice_thickness[n]    = p.thickness;
-        node_type[n]        = p.node_type;
-      }
+      // Initialize 2D geometric info at element nodes
+      nodal_parameter_values(element, P, i, j,
+                             node_type,
+                             bottom_elevation,
+                             ice_thickness,
+                             NULL,
+                             NULL);
 
       // skip ice-free (exterior) columns
       if (exterior_element(node_type)) {

diff --git a/src/stressbalance/blatter/residual.cc b/src/stressbalance/blatter/residual.cc
@@ -291,17 +291,12 @@ void Blatter::compute_residual(DMDALocalInfo *petsc_info,
     for (int i = info.gxs; i < info.gxs + info.gxm - 1; i++) {
 
       // Initialize 2D geometric info at element nodes
-      for (int n = 0; n < Nk; ++n) {
-        auto I = element.local_to_global(i, j, 0, n);
-
-        auto p = P[I.j][I.i];
-
-        bottom_elevation[n]  = p.bed;
-        ice_thickness[n]     = p.thickness;
-        node_type[n]         = p.node_type;
-        sea_level[n]         = p.sea_level;
-        surface_elevation[n] = p.bed + p.thickness;
-      }
+      nodal_parameter_values(element, P, i, j,
+                             node_type,
+                             bottom_elevation,
+                             ice_thickness,
+                             surface_elevation,
+                             sea_level);
 
       // skip ice-free (exterior) elements
       if (exterior_element(node_type)) {