Add Blatter::m_face{4,100}

"Face" elements don't use the information about domain distribution (only nodal values for
the current element), so they can be created once and kept as data members.

(This was not true in the code that supported coarsening in horizontal directions: dx and
dy were different on different multigrid levels.)

src/stressbalance/blatter/Blatter.cc

@@ -245,7 +245,13 @@ bool Blatter::neumann_bc_face(int face, const int *node_type) {
  * @param[in] coarsening_factor grid coarsening factor
  */
 Blatter::Blatter(IceGrid::ConstPtr grid, int Mz, int n_levels, int coarsening_factor)
-  : ShallowStressBalance(grid) {
+  : ShallowStressBalance(grid),
+    m_face4(grid->dx(), grid->dy(), fem::Q1Quadrature4()),    // 4-point Gaussian quadrature
+    m_face100(grid->dx(), grid->dy(), fem::Q1QuadratureN(10)) // 100-point quadrature for grounding lines
+{
+
+  assert(m_face4.n_pts() <= m_Nq);
+  assert(m_face100.n_pts() <= m_Nq);
 
   auto pism_da = grid->get_dm(1, 0);
 

src/stressbalance/blatter/Blatter.hh

@@ -93,6 +93,9 @@ protected:
 
   Vector2 m_work2[m_n_work][m_Nq];
 
+  fem::Q1Element3Face m_face4;
+  fem::Q1Element3Face m_face100;
+
   void init_impl();
 
   void define_model_state_impl(const File &output) const;

src/stressbalance/blatter/jacobian.cc

@@ -194,18 +194,12 @@ void Blatter::compute_jacobian(DMDALocalInfo *petsc_info,
 
   fem::Q1Element3 element(info, dx, dy, fem::Q13DQuadrature8());
 
-  fem::Q1Element3Face
-    face4(dx, dy, fem::Q1Quadrature4()),     // 4-point Gaussian quadrature
-    face100(dx, dy, fem::Q1QuadratureN(10)); // 100-point quadrature for grounding lines
-
   // Maximum number of nodes per element
   const int Nk = fem::q13d::n_chi;
   assert(element.n_chi() <= Nk);
 
   // Maximum number of quadrature points per element or face
   assert(element.n_pts() <= m_Nq);
-  assert(face4.n_pts() <= m_Nq);
-  assert(face100.n_pts() <= m_Nq);
 
   // scalar quantities
   double x[Nk], y[Nk], z[Nk];
@@ -292,7 +286,7 @@ void Blatter::compute_jacobian(DMDALocalInfo *petsc_info,
             floatation[n]         = P[I.j][I.i].floatation;
           }
 
-          fem::Q1Element3Face *face = grounding_line(floatation) ? &face100 : &face4;
+          fem::Q1Element3Face *face = grounding_line(floatation) ? &m_face100 : &m_face4;
 
           // face 4 is the bottom face in fem::q13d::incident_nodes
           face->reset(4, z);

src/stressbalance/blatter/residual.cc

@@ -257,13 +257,8 @@ void Blatter::compute_residual(DMDALocalInfo *petsc_info,
     dy = m_grid->dy();
 
   fem::Q1Element3 element(info, dx, dy, fem::Q13DQuadrature8());
-  fem::Q1Element3Face
-    face4(dx, dy, fem::Q1Quadrature4()),     // 4-point Gaussian quadrature
-    face100(dx, dy, fem::Q1QuadratureN(10)); // 100-point quadrature for grounding lines
-                                             // and partially-submerged faces
+
   assert(element.n_pts() <= m_Nq);
-  assert(face4.n_pts() <= m_Nq);
-  assert(face100.n_pts() <= m_Nq);
 
   // Number of nodes per element.
   const int Nk = fem::q13d::n_chi;
@@ -364,7 +359,7 @@ void Blatter::compute_residual(DMDALocalInfo *petsc_info,
           }
 
           // use an N*N-point equally-spaced quadrature at grounding lines
-          fem::Q1Element3Face *face = grounding_line(floatation) ? &face100 : &face4;
+          fem::Q1Element3Face *face = grounding_line(floatation) ? &m_face100 : &m_face4;
           // face 4 is the bottom face in fem::q13d::incident_nodes
           face->reset(4, z);
 
@@ -376,7 +371,7 @@ void Blatter::compute_residual(DMDALocalInfo *petsc_info,
         for (int f = 0; f < 4 and neumann_bc_face(f, node_type); ++f) {
           // use an N*N-point equally-spaced quadrature at for partially-submerged faces
           fem::Q1Element3Face *face = (partially_submerged_face(f, z, sea_level) ?
-                                       &face100 : &face4);
+                                       &m_face100 : &m_face4);
           face->reset(f, z);
 
           residual_lateral(element, *face, z, sea_level, R_nodal);