Feature/glace boundary conditions

• Add some of the missing boundary conditions for the acoustic solver
  • pressure
  • velocity
• Add Eucclhyd variant
• Create the language/algorithm data structure. This might be the place to define hard-coded (that is using C++) algorithms

src/mesh/MeshData.hpp

@@ -32,16 +32,44 @@ class MeshData : public IMeshData
 
  private:
   const MeshType& m_mesh;
-  std::shared_ptr<NodeValuePerCell<const Rd>> m_Cjr;
-  std::shared_ptr<NodeValuePerCell<const double>> m_ljr;
-  std::shared_ptr<NodeValuePerCell<const Rd>> m_njr;
+  NodeValuePerCell<const Rd> m_Cjr;
+  NodeValuePerCell<const double> m_ljr;
+  NodeValuePerCell<const Rd> m_njr;
   FaceValue<const Rd> m_xl;
-  std::shared_ptr<CellValue<const Rd>> m_xj;
-  std::shared_ptr<CellValue<const double>> m_Vj;
-  std::shared_ptr<NodeValuePerFace<const Rd>> m_Nlr;
+  CellValue<const Rd> m_xj;
+  CellValue<const double> m_Vj;
+  FaceValue<const double> m_ll;
+  NodeValuePerFace<const Rd> m_Nlr;
 
   PUGS_INLINE
   void
+  _compute_ll()
+  {
+    if constexpr (Dimension == 1) {
+      static_assert(Dimension != 1, "ll does not make sense in 1d");
+    } else {
+      const auto& nlr = this->Nlr();
+
+      FaceValue<double> ll{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;
+          for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+            lenght += l2Norm(nlr(face_id, i_node));
+          }
+
+          ll[face_id] = lenght;
+        });
+
+      m_ll = ll;
+    }
+  }
+
+  PUGS_INLINE void
   _computeFaceIsobarycenter()
   {   // Computes vertices isobarycenter
     if constexpr (Dimension == 1) {
@@ -78,7 +106,7 @@ class MeshData : public IMeshData
           const auto& cell_nodes = cell_to_node_matrix[j];
           xj[j]                  = 0.5 * (xr[cell_nodes[0]] + xr[cell_nodes[1]]);
         });
-      m_xj = std::make_shared<CellValue<const Rd>>(xj);
+      m_xj = xj;
     } else {
       const NodeValue<const Rd>& xr = m_mesh.xr();
 
@@ -95,7 +123,7 @@ class MeshData : public IMeshData
           }
           xj[j] = inv_cell_nb_nodes[j] * X;
         });
-      m_xj = std::make_shared<CellValue<const Rd>>(xj);
+      m_xj = xj;
     }
   }
 
@@ -119,7 +147,7 @@ class MeshData : public IMeshData
         }
         Vj[j] = inv_Dimension * sum_cjr_xr;
       });
-    m_Vj = std::make_shared<CellValue<const double>>(Vj);
+    m_Vj = Vj;
   }
 
   PUGS_INLINE
@@ -147,7 +175,7 @@ class MeshData : public IMeshData
           Nlr(l, 0) = Nr;
           Nlr(l, 1) = Nr;
         });
-      m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr);
+      m_Nlr = Nlr;
     } else {
       const NodeValue<const Rd>& xr = m_mesh.xr();
 
@@ -174,7 +202,7 @@ class MeshData : public IMeshData
             Nlr(l, r) = Nr;
           }
         });
-      m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr);
+      m_Nlr = Nlr;
     }
   }
 
@@ -189,7 +217,7 @@ class MeshData : public IMeshData
           Cjr(j, 0) = -1;
           Cjr(j, 1) = 1;
         });
-      m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+      m_Cjr = Cjr;
     } else if constexpr (Dimension == 2) {
       const NodeValue<const Rd>& xr   = m_mesh.xr();
       const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
@@ -205,7 +233,7 @@ class MeshData : public IMeshData
             Cjr(j, R)       = Rd{half_xrp_xrm[1], -half_xrp_xrm[0]};
           }
         });
-      m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+      m_Cjr = Cjr;
     } else if (Dimension == 3) {
       auto Nlr = this->Nlr();
 
@@ -252,7 +280,7 @@ class MeshData : public IMeshData
           }
         });
 
-      m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+      m_Cjr = Cjr;
     }
   }
 
@@ -265,13 +293,13 @@ class MeshData : public IMeshData
       NodeValuePerCell<double> ljr(m_mesh.connectivity());
       parallel_for(
         ljr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = 1; });
-      m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
+      m_ljr = ljr;
 
     } else {
       NodeValuePerCell<double> ljr(m_mesh.connectivity());
       parallel_for(
         Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm(Cjr[jr]); });
-      m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
+      m_ljr = ljr;
     }
   }
 
@@ -289,19 +317,18 @@ class MeshData : public IMeshData
       NodeValuePerCell<Rd> njr(m_mesh.connectivity());
       parallel_for(
         Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / ljr[jr]) * Cjr[jr]; });
-      m_njr = std::make_shared<NodeValuePerCell<const Rd>>(njr);
+      m_njr = njr;
     }
   }
 
   void
   _checkCellVolume()
   {
-    Assert(m_Vj);
-    auto& Vj = *m_Vj;
+    Assert(m_Vj.isBuilt());
 
     bool is_valid = [&] {
       for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
-        if (Vj[j] <= 0) {
+        if (m_Vj[j] <= 0) {
           return false;
         }
       }
@@ -321,44 +348,54 @@ class MeshData : public IMeshData
     return m_mesh;
   }
 
+  PUGS_INLINE
+  FaceValue<const double>
+  ll()
+  {
+    if (not m_ll.isBuilt()) {
+      this->_compute_ll();
+    }
+    return m_ll;
+  }
+
   PUGS_INLINE
   NodeValuePerFace<const Rd>
   Nlr()
   {
-    if (not m_Nlr) {
+    if (not m_Nlr.isBuilt()) {
       this->_computeNlr();
     }
-    return *m_Nlr;
+    return m_Nlr;
   }
 
   PUGS_INLINE
   NodeValuePerCell<const Rd>
   Cjr()
   {
-    if (not m_Cjr) {
+    if (not m_Cjr.isBuilt()) {
       this->_computeCjr();
     }
-    return *m_Cjr;
+    return m_Cjr;
   }
 
   PUGS_INLINE
   NodeValuePerCell<const double>
   ljr()
   {
-    if (not m_ljr) {
+    if (not m_ljr.isBuilt()) {
       this->_compute_ljr();
     }
-    return *m_ljr;
+    return m_ljr;
   }
 
   PUGS_INLINE
   NodeValuePerCell<const Rd>
   njr()
   {
-    if (not m_njr) {
+    if (not m_njr.isBuilt()) {
       this->_compute_njr();
     }
-    return *m_njr;
+    return m_njr;
   }
 
   PUGS_INLINE
@@ -375,21 +412,21 @@ class MeshData : public IMeshData
   CellValue<const Rd>
   xj()
   {
-    if (not m_xj) {
+    if (not m_xj.isBuilt()) {
       this->_computeCellIsobarycenter();
     }
-    return *m_xj;
+    return m_xj;
   }
 
   PUGS_INLINE
   CellValue<const double>
   Vj()
   {
-    if (not m_Vj) {
+    if (not m_Vj.isBuilt()) {
       this->_computeCellVolume();
       this->_checkCellVolume();
     }
-    return *m_Vj;
+    return m_Vj;
   }
 
  private: