Add the third dimension to QuadPoint and Germ

src/util/fem/FEM.cc

@@ -34,7 +34,7 @@ Germ chi(unsigned int k, const QuadPoint &pt) {
 
   return {0.5 * (1.0 + xis[k] * pt.xi),
           0.5 * xis[k],
-          0.0};                         // unused
+          0.0, 0.0};            // unused
 }
 
 } // end of namespace linear
@@ -76,7 +76,8 @@ Germ chi(unsigned int k, const QuadPoint &pt) {
 
   return {0.25 * (1.0 + xi[k] * pt.xi) * (1.0 + eta[k] * pt.eta),
           0.25 * xi[k] * (1.0 + eta[k] * pt.eta),
-          0.25 * eta[k] * (1.0 + xi[k] * pt.xi)};
+          0.25 * eta[k] * (1.0 + xi[k] * pt.xi),
+          0.0};                 // unused
 }
 
 } // end of namespace q1
@@ -90,11 +91,11 @@ Germ chi(unsigned int k, const QuadPoint &pt) {
   switch (k) {
   default:
   case 0:
-    return {1.0 - pt.xi - pt.eta, -1.0, -1.0};
+    return {1.0 - pt.xi - pt.eta, -1.0, -1.0, 0.0};
   case 1:
-    return {pt.xi, 1.0, 0.0};
+    return {pt.xi, 1.0, 0.0, 0.0};
   case 2:
-    return {pt.eta, 0.0, 1.0};
+    return {pt.eta, 0.0, 1.0, 0.0};
  }
 }
 

src/util/fem/FEM.hh

@@ -155,13 +155,16 @@ struct Germ {
   double dx;
   //! Function derivative with respect to y.
   double dy;
+  //! Function derivative with respect to z.
+  double dz;
 };
 
 //! Coordinates of a quadrature point, in the (xi, eta) coordinate space (i.e. on the
 //! reference element).
 struct QuadPoint {
   double xi;
   double eta;
+  double zeta;
 };
 
 //! Hard-wired maximum number of points a quadrature can use. This is used as the size of

src/util/fem/Quadrature.cc

@@ -61,7 +61,7 @@ Gaussian2::Gaussian2(double D) {
   // coordinates and weights of the 2-point 1D Gaussian quadrature
   double A = 1.0 / std::sqrt(3.0);
 
-  m_points  = {{-A, 0.0}, {A, 0.0}};
+  m_points  = {{-A, 0.0, 0.0}, {A, 0.0, 0.0}};
   m_weights = {0.5 * D, 0.5 * D};
 }
 
@@ -139,9 +139,9 @@ P1Quadrature3::P1Quadrature3() {
     one_over_six   = 1.0 / 6.0,
     two_over_three = 2.0 / 3.0;
 
-  m_points = {{two_over_three, one_over_six},
-              {one_over_six,   two_over_three},
-              {one_over_six,   one_over_six}};
+  m_points = {{two_over_three, one_over_six, 0.0},
+              {one_over_six,   two_over_three, 0.0},
+              {one_over_six,   one_over_six, 0.0}};
 
   m_weights = {one_over_six, one_over_six, one_over_six};
 }