Added quite dirty BC code. Hopefully on the road to 2d/3d ...

- pressure and symmetry by now
- needs a cleaner implementation before velocity implementation

src/main.cpp

@@ -13,7 +13,7 @@
 
 #include <Connectivity1D.hpp>
 #include <Mesh.hpp>
-#include <BoundaryCondition.hpp>
+#include <FacesBoundaryCondition.hpp>
 #include <AcousticSolver.hpp>
 
 #include <TinyVector.hpp>
@@ -132,7 +132,9 @@ int main(int argc, char *argv[])
       FacesBoundaryCondition bc0(FacesBoundaryCondition::symmetry, 0,
 			       std::vector<unsigned int>({0u}));
       faces_bc_list.push_back(bc0);
-      FacesBoundaryCondition bc1(FacesBoundaryCondition::symmetry, 0,
+      // FacesBoundaryCondition bc1(FacesBoundaryCondition::symmetry, 0,
+      // 				 std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
+      FacesBoundaryCondition bc1(FacesBoundaryCondition::pressure, 0.125,
 				 std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
       faces_bc_list.push_back(bc1);
     }
@@ -141,7 +143,7 @@ int main(int argc, char *argv[])
 
     unknowns.initializeSod();
 
-    AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns);
+    AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, faces_bc_list);
 
     typedef TinyVector<MeshType::dimension> Rd;
 

src/scheme/AcousticSolver.hpp

@@ -11,6 +11,7 @@
 #include <Mesh.hpp>
 #include <MeshData.hpp>
 #include <FiniteVolumesEulerUnknowns.hpp>
+#include <FacesBoundaryCondition.hpp>
 
 template<typename MeshData>
 class AcousticSolver
@@ -21,6 +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;
 
   constexpr static size_t dimension = MeshType::dimension;
 
@@ -136,6 +138,68 @@ private:
     return m_br;
   }
 
+  void applyBoundaryConditions()
+  {
+    for (const auto& face_bc : m_faces_boundary_condition_list) {
+      switch (face_bc.type()) {
+      case FacesBoundaryCondition::normal_velocity: {
+	std::cerr << "normal_velocity NIY\n";
+	std::exit(0);
+	break;
+      }
+      case FacesBoundaryCondition::velocity: {
+	std::cerr << "velocity NIY\n";
+	std::exit(0);
+	break;
+      }
+      case FacesBoundaryCondition::pressure: {
+	const Kokkos::View<const unsigned int**> face_cells
+	  = m_connectivity.faceCells();
+	const Kokkos::View<const unsigned short**> face_cell_local_face
+	  = m_connectivity.faceCellLocalFace();
+	const Kokkos::View<const Rd**> Cjr
+	  = m_mesh_data.Cjr();
+
+	Kokkos::parallel_for(face_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];
+	    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);
+	  });
+	break;
+      }
+      case FacesBoundaryCondition::symmetry: {
+	const Rdd I = identity;
+	const Kokkos::View<const unsigned int**> face_cells
+	  = m_connectivity.faceCells();
+	const Kokkos::View<const unsigned short**> face_cell_local_face
+	  = m_connectivity.faceCellLocalFace();
+	const Kokkos::View<const Rd**> njr
+	  = m_mesh_data.njr();
+
+	Kokkos::parallel_for(face_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];
+	    assert(m_connectivity.faceNbCells() == 1);
+	    const unsigned int j = face_cells(l,0);
+	    const unsigned int L = face_cell_local_face(l,0);
+
+	    const Rd& n = njr(j,L);
+	    const Rdd nxn = tensorProduct(n,n);
+	    const Rdd P = I-nxn;
+
+	    m_Ar(l) = P*m_Ar(l)*P + nxn;
+	    m_br(l) = P*m_br(l);
+	  });
+	break;
+      }
+      }
+    }
+  }
+  
   Kokkos::View<Rd*>
   computeUr(const Kokkos::View<const Rdd*>& Ar,
 	    const Kokkos::View<const Rd*>& br) {
@@ -144,8 +208,6 @@ private:
     Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
 	m_ur[r]=invAr(r)*br(r);
       });
-    m_ur[0]=zero;
-    m_ur[m_mesh.numberOfNodes()-1]=zero;
 
     return m_ur;
   }
@@ -200,6 +262,8 @@ private:
     const Kokkos::View<const Rdd*> Ar = computeAr(Ajr);
     const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj);
 
+    this->applyBoundaryConditions();
+
     Kokkos::View<Rd*> ur = m_ur;
     Kokkos::View<Rd**> Fjr = m_Fjr;
     ur  = computeUr(Ar, br);
@@ -217,10 +281,12 @@ private:
 
 public:
   AcousticSolver(MeshData& mesh_data,
-		 UnknownsType& unknowns)
+		 UnknownsType& unknowns,
+		 const std::vector<FacesBoundaryCondition>& faces_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_br("br", m_mesh.numberOfNodes()),
       m_Ajr("Ajr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
       m_Ar("Ar", m_mesh.numberOfNodes()),