More work on FEM code

Now I'm ready to write the 3D Q1 element code.
pism · Oct 14, 2020 · bf3cc68 · bf3cc68
1 parent 859054c
commit bf3cc68
Show file tree

Hide file tree

Showing 11 changed files with 330 additions and 176 deletions.
diff --git a/src/inverse/functional/IPGroundedIceH1NormFunctional.cc b/src/inverse/functional/IPGroundedIceH1NormFunctional.cc
@@ -202,7 +202,7 @@ void IPGroundedIceH1NormFunctional2S::gradientAt(IceModelVec2S &x, IceModelVec2S
                                    m_cH1*(dxdx_qq*m_element.chi(q, k).dx + dxdy_qq*m_element.chi(q, k).dy));
         } // k
       } // q
-      m_element.add_contribution(gradient_e, gradient);
+      m_element.add_contribution(gradient_e, gradient.array());
     } // j
   } // i
 }

diff --git a/src/inverse/functional/IPTotalVariationFunctional.cc b/src/inverse/functional/IPTotalVariationFunctional.cc
@@ -126,7 +126,7 @@ void IPTotalVariationFunctional2S::gradientAt(IceModelVec2S &x, IceModelVec2S &g
             *(dxdx_qq*m_element.chi(q, k).dx + dxdy_qq*m_element.chi(q, k).dy);
         } // k
       } // q
-      m_element.add_contribution(gradient_e, gradient);
+      m_element.add_contribution(gradient_e, gradient.array());
     } // j
   } // i
 }

diff --git a/src/inverse/functional/IP_H1NormFunctional.cc b/src/inverse/functional/IP_H1NormFunctional.cc
@@ -174,7 +174,7 @@ void IP_H1NormFunctional2S::gradientAt(IceModelVec2S &x, IceModelVec2S &gradient
                                    m_cH1*(dxdx_qq*m_element.chi(q, k).dx + dxdy_qq*m_element.chi(q, k).dy));
         } // k
       } // q
-      m_element.add_contribution(gradient_e, gradient);
+      m_element.add_contribution(gradient_e, gradient.array());
     } // j
   } // i
 }

diff --git a/src/inverse/functional/IP_L2NormFunctional.cc b/src/inverse/functional/IP_L2NormFunctional.cc
@@ -147,7 +147,7 @@ void IP_L2NormFunctional2S::gradientAt(IceModelVec2S &x, IceModelVec2S &gradient
           gradient_e[k] += 2*W*x_qq*m_element.chi(q, k).val;
         } // k
       } // q
-      m_element.add_contribution(gradient_e, gradient);
+      m_element.add_contribution(gradient_e, gradient.array());
     } // j
   } // i
 }
@@ -279,7 +279,7 @@ void IP_L2NormFunctional2V::gradientAt(IceModelVec2V &x, IceModelVec2V &gradient
           gradient_e[k].v += gcommon*x_qq.v;
         } // k
       } // q
-      m_element.add_contribution(gradient_e, gradient);
+      m_element.add_contribution(gradient_e, gradient.array());
     } // j
   } // i
 }

diff --git a/src/stressbalance/ssa/SSAFEM.cc b/src/stressbalance/ssa/SSAFEM.cc
@@ -616,7 +616,7 @@ void SSAFEM::cache_residual_cfbc(const Inputs &inputs) {
         int node_type[Nk];
         m_q1_element.nodal_values(m_node_type, node_type);
 
-        fem::Element *E = nullptr;
+        fem::Element2 *E = nullptr;
         {
           auto type = fem::element_type(node_type);
 
@@ -708,7 +708,7 @@ void SSAFEM::cache_residual_cfbc(const Inputs &inputs) {
 
         } // loop over element sides
 
-        E->add_contribution(I.data(), m_boundary_integral);
+        E->add_contribution(I.data(), m_boundary_integral.array());
 
       } // i-loop
     } // j-loop
@@ -773,7 +773,7 @@ void SSAFEM::compute_local_function(Vector2 const *const *const velocity_global,
     for (int j = ys; j < ys + ym; j++) {
       for (int i = xs; i < xs + xm; i++) {
 
-        fem::Element *E = nullptr;
+        fem::Element2 *E = nullptr;
         {
           m_q1_element.reset(i, j);
 
@@ -987,7 +987,7 @@ void SSAFEM::compute_local_jacobian(Vector2 const *const *const velocity_global,
     for (int j = ys; j < ys + ym; j++) {
       for (int i = xs; i < xs + xm; i++) {
 
-        fem::Element *E = nullptr;
+        fem::Element2 *E = nullptr;
         {
           m_q1_element.reset(i, j);
 

diff --git a/src/util/fem/DirichletData.cc b/src/util/fem/DirichletData.cc
@@ -74,7 +74,7 @@ void DirichletData::finish(const IceModelVec *values) {
 }
 
 //! @brief Constrain `element`, i.e. ensure that quadratures do not contribute to Dirichlet nodes by marking corresponding rows and columns as "invalid".
-void DirichletData::constrain(Element &element) {
+void DirichletData::constrain(Element2 &element) {
   element.nodal_values(m_indices->array(), m_indices_e);
   auto n_chi = element.n_chi();
   for (int k = 0; k < n_chi; k++) {
@@ -95,7 +95,7 @@ DirichletData_Scalar::DirichletData_Scalar(const IceModelVec2Int *indices,
   init(indices, m_values, weight);
 }
 
-void DirichletData_Scalar::enforce(const Element &element, double* x_nodal) {
+void DirichletData_Scalar::enforce(const Element2 &element, double* x_nodal) {
   assert(m_values != NULL);
 
   element.nodal_values(m_indices->array(), m_indices_e);
@@ -108,11 +108,11 @@ void DirichletData_Scalar::enforce(const Element &element, double* x_nodal) {
   }
 }
 
-void DirichletData_Scalar::enforce_homogeneous(const Element &element, double* x_nodal) {
+void DirichletData_Scalar::enforce_homogeneous(const Element2 &element, double* x_nodal) {
   element.nodal_values(m_indices->array(), m_indices_e);
   for (int k = 0; k < element.n_chi(); k++) {
     if (m_indices_e[k] > 0.5) { // Dirichlet node
-      x_nodal[k] = 0.;
+      x_nodal[k] = 0.0;
     }
   }
 }
@@ -191,7 +191,7 @@ DirichletData_Vector::DirichletData_Vector(const IceModelVec2Int *indices,
   init(indices, m_values, weight);
 }
 
-void DirichletData_Vector::enforce(const Element &element, Vector2* x_nodal) {
+void DirichletData_Vector::enforce(const Element2 &element, Vector2* x_nodal) {
   assert(m_values != NULL);
 
   element.nodal_values(m_indices->array(), m_indices_e);
@@ -204,12 +204,11 @@ void DirichletData_Vector::enforce(const Element &element, Vector2* x_nodal) {
   }
 }
 
-void DirichletData_Vector::enforce_homogeneous(const Element &element, Vector2* x_nodal) {
+void DirichletData_Vector::enforce_homogeneous(const Element2 &element, Vector2* x_nodal) {
   element.nodal_values(m_indices->array(), m_indices_e);
   for (int k = 0; k < element.n_chi(); k++) {
     if (m_indices_e[k] > 0.5) { // Dirichlet node
-      x_nodal[k].u = 0.0;
-      x_nodal[k].v = 0.0;
+      x_nodal[k] = 0.0;
     }
   }
 }

diff --git a/src/util/fem/DirichletData.hh b/src/util/fem/DirichletData.hh
@@ -32,12 +32,12 @@ class IceModelVec2V;
 
 namespace fem {
 
-class Element;
+class Element2;
 
 //* Parts shared by scalar and 2D vector Dirichlet data classes.
 class DirichletData {
 public:
-  void constrain(Element &element);
+  void constrain(Element2 &element);
   operator bool() {
     return m_indices != NULL;
   }
@@ -59,8 +59,8 @@ public:
                        double weight = 1.0);
   ~DirichletData_Scalar();
 
-  void enforce(const Element &element, double* x_e);
-  void enforce_homogeneous(const Element &element, double* x_e);
+  void enforce(const Element2 &element, double* x_e);
+  void enforce_homogeneous(const Element2 &element, double* x_e);
   void fix_residual(double const *const *const x_global, double **r_global);
   void fix_residual_homogeneous(double **r_global);
   void fix_jacobian(Mat J);
@@ -74,8 +74,8 @@ public:
                        double weight);
   ~DirichletData_Vector();
 
-  void enforce(const Element &element, Vector2* x_e);
-  void enforce_homogeneous(const Element &element, Vector2* x_e);
+  void enforce(const Element2 &element, Vector2* x_e);
+  void enforce_homogeneous(const Element2 &element, Vector2* x_e);
   void fix_residual(Vector2 const *const *const x_global, Vector2 **r_global);
   void fix_residual_homogeneous(Vector2 **r);
   void fix_jacobian(Mat J);