Add Blatter::residual_surface()

This method is needed to implement the top-surface residual contributions in MMS-based
verification tests.

src/stressbalance/blatter/Blatter.hh

@@ -159,6 +159,10 @@ protected:
                               const Vector2 *u_nodal,
                               Vector2 *residual);
 
+  virtual void residual_surface(const fem::Q1Element3 &element,
+                                const fem::Q1Element3Face &face,
+                                Vector2 *residual);
+
   virtual void residual_lateral(const fem::Q1Element3 &element,
                                 const fem::Q1Element3Face &face,
                                 const double *z_nodal,

src/stressbalance/blatter/residual.cc

@@ -150,6 +150,21 @@ void Blatter::residual_basal(const fem::Q1Element3 &element,
   }
 }
 
+/*!
+ * Top surface contribution to the residual.
+ *
+ * Used by verification tests ONLY.
+ */
+void Blatter::residual_surface(const fem::Q1Element3 &element,
+                               const fem::Q1Element3Face &face,
+                               Vector2 *residual) {
+  (void) element;
+  (void) face;
+  (void) residual;
+  // In normal circumstances the top surface contribution is zero (natural BCs apply).
+}
+
+
 /*!
  * Computes the residual contribution of lateral boundary conditions.
  *
@@ -371,6 +386,14 @@ void Blatter::compute_residual(DMDALocalInfo *petsc_info,
           residual_lateral(element, *face, z, sea_level, R_nodal);
         } // end of the loop over element faces
 
+        // top boundary (verification tests only)
+        if (k == info.mz - 2) {
+          // face 5 is the top face in fem::q13d::incident_nodes
+          m_face4.reset(5, z);
+
+          residual_surface(element, m_face4, R_nodal);
+        }
+
         element.add_contribution(R_nodal, R);
       } // end of the loop over i
     } // end of the loop over j