Begining of BC cleaning.

"The path is straight, but the slope is strong!"

src/main.cpp

@@ -13,7 +13,7 @@
 
 #include <Connectivity1D.hpp>
 #include <Mesh.hpp>
-#include <FacesBoundaryCondition.hpp>
+#include <BoundaryCondition.hpp>
 #include <AcousticSolver.hpp>
 
 #include <TinyVector.hpp>
@@ -127,23 +127,25 @@ int main(int argc, char *argv[])
     MeshType mesh(connectivity);
     MeshDataType mesh_data(mesh);
 
-    std::vector<FacesBoundaryCondition> faces_bc_list;
+    std::vector<BoundaryConditionHandler> bc_list;
     { // quite dirty!
-      FacesBoundaryCondition bc0(FacesBoundaryCondition::symmetry, 0,
-			       std::vector<unsigned int>({0u}));
-      faces_bc_list.push_back(bc0);
-      // FacesBoundaryCondition bc1(FacesBoundaryCondition::symmetry, 0,
-      // 				 std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
-      FacesBoundaryCondition bc1(FacesBoundaryCondition::pressure, 0.125,
+      SymmetryBoundaryCondition* sym_bc0
+	= new SymmetryBoundaryCondition(std::vector<unsigned int>({0u}));
+      std::shared_ptr<SymmetryBoundaryCondition> bc0(sym_bc0);
+      bc_list.push_back(BoundaryConditionHandler(bc0));
+
+      PressureBoundaryCondition* pres_bc1
+	= new PressureBoundaryCondition(0.1,
 					std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
-      faces_bc_list.push_back(bc1);
+      std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1);
+      bc_list.push_back(BoundaryConditionHandler(bc1));
     }
 
     UnknownsType unknowns(mesh_data);
 
     unknowns.initializeSod();
 
-    AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, faces_bc_list);
+    AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
 
     typedef TinyVector<MeshType::dimension> Rd;
 

src/scheme/AcousticSolver.hpp

@@ -11,7 +11,7 @@
 #include <Mesh.hpp>
 #include <MeshData.hpp>
 #include <FiniteVolumesEulerUnknowns.hpp>
-#include <FacesBoundaryCondition.hpp>
+#include <BoundaryCondition.hpp>
 
 template<typename MeshData>
 class AcousticSolver
@@ -22,7 +22,7 @@ class AcousticSolver
   MeshData& m_mesh_data;
   const MeshType& m_mesh;
   const typename MeshType::Connectivity& m_connectivity;
-  const std::vector<FacesBoundaryCondition>& m_faces_boundary_condition_list;
+  const std::vector<BoundaryConditionHandler>& m_boundary_condition_list;
 
   constexpr static size_t dimension = MeshType::dimension;
 
@@ -140,19 +140,21 @@ private:
 
   void applyBoundaryConditions()
   {
-    for (const auto& face_bc : m_faces_boundary_condition_list) {
-      switch (face_bc.type()) {
-      case FacesBoundaryCondition::normal_velocity: {
+    for (const auto& handler : m_boundary_condition_list) {
+      switch (handler.boundaryCondition().type()) {
+      case BoundaryCondition::normal_velocity: {
 	std::cerr << "normal_velocity NIY\n";
 	std::exit(0);
 	break;
       }
-      case FacesBoundaryCondition::velocity: {
+      case BoundaryCondition::velocity: {
 	std::cerr << "velocity NIY\n";
 	std::exit(0);
 	break;
       }
-      case FacesBoundaryCondition::pressure: {
+      case BoundaryCondition::pressure: {
+	const PressureBoundaryCondition& pressure_bc
+	  = dynamic_cast<const PressureBoundaryCondition&>(handler.boundaryCondition());
 	const Kokkos::View<const unsigned int**> face_cells
 	  = m_connectivity.faceCells();
 	const Kokkos::View<const unsigned short**> face_cell_local_face
@@ -160,18 +162,20 @@ private:
 	const Kokkos::View<const Rd**> Cjr
 	  = m_mesh_data.Cjr();
 
-	Kokkos::parallel_for(face_bc.numberOfFaces(), KOKKOS_LAMBDA(const int& l_number) {
+	Kokkos::parallel_for(pressure_bc.numberOfFaces(), KOKKOS_LAMBDA(const int& l_number) {
 	    // quite ugly: node/faces are melt in a bad way... at least works in 1d...
-	    const int l = face_bc.faceList()[l_number];
+	    const int l = pressure_bc.faceList()[l_number];
 	    assert(m_connectivity.faceNbCells() == 1);
 	    const unsigned int j = face_cells(l,0);
 	    const unsigned int L = face_cell_local_face(l,0);
 
-	    m_br(l) -= face_bc.value()*Cjr(j,L);
+	    m_br(l) -= pressure_bc.value()*Cjr(j,L);
 	  });
 	break;
       }
-      case FacesBoundaryCondition::symmetry: {
+      case BoundaryCondition::symmetry: {
+	const SymmetryBoundaryCondition& symmetry_bc
+	  = dynamic_cast<const SymmetryBoundaryCondition&>(handler.boundaryCondition());
 	const Rdd I = identity;
 	const Kokkos::View<const unsigned int**> face_cells
 	  = m_connectivity.faceCells();
@@ -180,9 +184,9 @@ private:
 	const Kokkos::View<const Rd**> njr
 	  = m_mesh_data.njr();
 
-	Kokkos::parallel_for(face_bc.numberOfFaces(), KOKKOS_LAMBDA(const int& l_number) {
+	Kokkos::parallel_for(symmetry_bc.numberOfFaces(), KOKKOS_LAMBDA(const int& l_number) {
 	    // quite ugly: node/faces are melt in a bad way... at least works in 1d...
-	    const int l = face_bc.faceList()[l_number];
+	    const int l = symmetry_bc.faceList()[l_number];
 	    assert(m_connectivity.faceNbCells() == 1);
 	    const unsigned int j = face_cells(l,0);
 	    const unsigned int L = face_cell_local_face(l,0);
@@ -282,11 +286,11 @@ private:
 public:
   AcousticSolver(MeshData& mesh_data,
 		 UnknownsType& unknowns,
-		 const std::vector<FacesBoundaryCondition>& faces_bc_list)
+		 const std::vector<BoundaryConditionHandler>& bc_list)
     : m_mesh_data(mesh_data),
       m_mesh(mesh_data.mesh()),
       m_connectivity(m_mesh.connectivity()),
-      m_faces_boundary_condition_list(faces_bc_list),
+      m_boundary_condition_list(bc_list),
       m_br("br", m_mesh.numberOfNodes()),
       m_Ajr("Ajr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
       m_Ar("Ar", m_mesh.numberOfNodes()),

src/scheme/BoundaryCondition.hpp

+#ifndef BOUNDARY_CONDITION_HANDLER_HPP
+#define BOUNDARY_CONDITION_HANDLER_HPP
+
+#include <vector>
+#include <memory>
+
+class BoundaryCondition
+{
+public:
+  enum Type {
+    velocity,
+    normal_velocity,
+    pressure,
+    symmetry
+  };
+
+  const Type m_type;
+protected:
+
+  BoundaryCondition(const Type& type)
+    : m_type(type)
+  {
+    ;
+  }
+
+public:
+  const Type& type() const
+  {
+    return m_type;
+  }
+
+  virtual ~BoundaryCondition() = default;
+};
+
+class PressureBoundaryCondition
+  : public BoundaryCondition
+{
+private:
+  const double m_value;
+
+  const size_t m_number_of_faces;
+  Kokkos::View<unsigned int*> m_face_list;
+
+public:
+  double value() const
+  {
+    return m_value;
+  }
+
+  const size_t& numberOfFaces() const
+  {
+    return m_number_of_faces;
+  }
+
+  const Kokkos::View<const unsigned int*> faceList() const
+  {
+    return m_face_list;
+  }
+
+  PressureBoundaryCondition(const double& value,
+			    const std::vector<unsigned int>& faces)
+    : BoundaryCondition(BoundaryCondition::pressure),
+      m_value(value),
+      m_number_of_faces(faces.size()),
+      m_face_list("faces_list", m_number_of_faces)
+  {
+    Kokkos::parallel_for(m_number_of_faces, KOKKOS_LAMBDA(const int& f){
+	m_face_list[f]=faces[f];
+      });
+  }
+
+  ~PressureBoundaryCondition() = default;
+};
+
+class SymmetryBoundaryCondition
+  : public BoundaryCondition
+{
+private:
+  const size_t m_number_of_faces;
+  Kokkos::View<unsigned int*> m_face_list;
+
+public:
+  const size_t& numberOfFaces() const
+  {
+    return m_number_of_faces;
+  }
+
+  const Kokkos::View<const unsigned int*> faceList() const
+  {
+    return m_face_list;
+  }
+
+  SymmetryBoundaryCondition(const std::vector<unsigned int>& faces)
+    : BoundaryCondition(BoundaryCondition::symmetry),
+      m_number_of_faces(faces.size()),
+      m_face_list("faces_list", m_number_of_faces)
+  {
+    Kokkos::parallel_for(m_number_of_faces, KOKKOS_LAMBDA(const int& f){
+	m_face_list[f]=faces[f];
+      });
+  }
+
+  ~SymmetryBoundaryCondition() = default;
+};
+
+class BoundaryConditionHandler
+{
+private:
+  std::shared_ptr<BoundaryCondition> m_boundary_condition;
+
+public:
+  const BoundaryCondition& boundaryCondition() const
+  {
+    return *m_boundary_condition;
+  }
+  
+  KOKKOS_INLINE_FUNCTION
+  BoundaryConditionHandler& operator=(BoundaryConditionHandler&) = default;
+
+  KOKKOS_INLINE_FUNCTION
+  BoundaryConditionHandler(const BoundaryConditionHandler&) = default;
+
+  KOKKOS_INLINE_FUNCTION
+  BoundaryConditionHandler(BoundaryConditionHandler&&) = default;
+
+  KOKKOS_INLINE_FUNCTION
+  BoundaryConditionHandler(std::shared_ptr<BoundaryCondition> boundary_condition)
+    : m_boundary_condition(boundary_condition)
+  {
+  }
+
+  ~BoundaryConditionHandler() = default;
+};
+
+#endif // BOUNDARY_CONDITION_HANDLER_HPP

src/scheme/FacesBoundaryCondition.hpp

-#ifndef FACES_BOUNDARY_CONDITION_HPP
-#define FACES_BOUNDARY_CONDITION_HPP
-
-#include <vector>
-
-class FacesBoundaryCondition
-{
-public:
-  enum Type
-  {
-    velocity,
-    normal_velocity,
-    pressure,
-    symmetry
-  };
-
-private:
-  const Type m_type;
-  const double m_value;
-  const size_t m_number_of_faces;
-  Kokkos::View<unsigned int*> m_face_list;
-
-public:
-  const size_t& numberOfFaces() const
-  {
-    return m_number_of_faces;
-  }
-
-  const Kokkos::View<const unsigned int*> faceList() const
-  {
-    return m_face_list;
-  }
-  
-  const Type& type() const
-  {
-    return m_type;
-  }
-
-  double value() const
-  {
-    return m_value;
-  }
-
-  KOKKOS_INLINE_FUNCTION
-  FacesBoundaryCondition(const Type& type,
-			 const double& value,
-			 const std::vector<unsigned int>& faces)
-    : m_type(type),
-      m_value(value),
-      m_number_of_faces(faces.size()),
-      m_face_list("face_list", m_number_of_faces)
-  {
-    Kokkos::parallel_for(m_number_of_faces, KOKKOS_LAMBDA(const int& f){
-	m_face_list[f]=faces[f];
-      });
-  }
-  ~FacesBoundaryCondition() = default;
-};
-
-#endif // FACES_BOUNDARY_CONDITION_HPP