Fix face area calculation

Also add face (area) normal and face unit normal computation

src/mesh/MeshData.hpp

@@ -42,36 +42,75 @@ class MeshData : public IMeshData
   CellValue<const Rd> m_cell_iso_barycenter;
   CellValue<const double> m_Vj;
   CellValue<const double> m_sum_r_ljr;
+  FaceValue<const Rd> m_Nl;
+  FaceValue<const Rd> m_nl;
   FaceValue<const double> m_ll;
 
   PUGS_INLINE
   void
-  _compute_ll()
+  _computeNl()
   {
     if constexpr (Dimension == 1) {
-      static_assert(Dimension != 1, "ll does not make sense in 1d");
+      static_assert(Dimension != 1, "Nl does not make sense in 1d");
     } else {
       const auto& Nlr = this->Nlr();
 
-      FaceValue<double> ll{m_mesh.connectivity()};
+      FaceValue<Rd> Nl{m_mesh.connectivity()};
 
       const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
       parallel_for(
         m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
           const auto& face_nodes = face_to_node_matrix[face_id];
 
-          double lenght = 0;
+          Rd N = zero;
           for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-            lenght += l2Norm(Nlr(face_id, i_node));
+            N += Nlr(face_id, i_node);
           }
-
-          ll[face_id] = lenght;
+          Nl[face_id] = N;
         });
 
+      m_Nl = Nl;
+    }
+  }
+
+  PUGS_INLINE
+  void
+  _compute_ll()
+  {
+    if constexpr (Dimension == 1) {
+      static_assert(Dimension != 1, "ll does not make sense in 1d");
+    } else {
+      const auto& Nl = this->Nl();
+
+      FaceValue<double> ll{m_mesh.connectivity()};
+
+      parallel_for(
+        m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { ll[face_id] = L2Norm(Nl[face_id]); });
+
       m_ll = ll;
     }
   }
 
+  PUGS_INLINE
+  void
+  _compute_nl()
+  {
+    if constexpr (Dimension == 1) {
+      static_assert(Dimension != 1, "Nl does not make sense in 1d");
+    } else {
+      const auto& Nl = this->Nl();
+      const auto& ll = this->ll();
+
+      FaceValue<Rd> nl{m_mesh.connectivity()};
+
+      const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+      parallel_for(
+        m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) { nl[face_id] = (1. / ll[face_id]) * Nl[face_id]; });
+
+      m_Nl = Nl;
+    }
+  }
+
   PUGS_INLINE
   void
   _compute_nlr()
@@ -310,7 +349,7 @@ class MeshData : public IMeshData
             Rd Nr            = zero;
             const Rd two_dxr = 2 * dxr[r];
             for (size_t s = 0; s < nb_nodes; ++s) {
-              Nr += crossProduct((two_dxr - dxr[s]), xr[face_nodes[s]]);
+              Nr += crossProduct(two_dxr - dxr[s], xr[face_nodes[s]]);
             }
             Nr *= inv_12_nb_nodes;
             Nr -= (1. / 6.) * crossProduct(dxr[r], xr[face_nodes[r]]);
@@ -358,8 +397,7 @@ class MeshData : public IMeshData
       const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
 
       NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
-      parallel_for(
-        Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Cjr[jr] = zero; });
+      Cjr.fill(zero);
 
       parallel_for(
         m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
@@ -491,6 +529,16 @@ class MeshData : public IMeshData
     return m_mesh;
   }
 
+  PUGS_INLINE
+  FaceValue<const Rd>
+  Nl()
+  {
+    if (not m_Nl.isBuilt()) {
+      this->_computeNl();
+    }
+    return m_Nl;
+  }
+
   PUGS_INLINE
   FaceValue<const double>
   ll()
@@ -501,6 +549,16 @@ class MeshData : public IMeshData
     return m_ll;
   }
 
+  PUGS_INLINE
+  FaceValue<const Rd>
+  nl()
+  {
+    if (not m_nl.isBuilt()) {
+      this->_compute_nl();
+    }
+    return m_nl;
+  }
+
   PUGS_INLINE
   NodeValuePerFace<const Rd>
   Nlr()