diff --git a/src/main.cpp b/src/main.cpp index 56f458abf21399fb097aa0dc1948936b225d174a..e1f9260baa1daab57d7615faf26980b13ea221be 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -238,7 +238,6 @@ int main(int argc, char *argv[]) timer.reset(); MeshDataType mesh_data(mesh); - std::vector<BoundaryConditionHandler> bc_list; { #warning Should extract boundaries from mesh itself to avoid inconsistencies @@ -254,7 +253,10 @@ int main(int argc, char *argv[]) const RefId& ref = ref_face_list.refId(); if (ref == sym_bc_descriptor.boundaryDescriptor()) { std::cout << "Found mesh boundary to impose " << sym_bc_descriptor << '\n'; - ref_boundary[ref] = MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list.faceList()); + SymmetryBoundaryCondition<MeshType::dimension>* sym_bc + = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list)); + std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc); + bc_list.push_back(BoundaryConditionHandler(bc)); } } break; @@ -265,48 +267,6 @@ int main(int argc, char *argv[]) } } } - for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList(); - ++i_ref_face_list) { - const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list); - const RefId& ref = ref_face_list.refId(); - - TinyVector<2> normal(0,0); - if ((ref.tagName()== std::string("XMIN")) or (ref.tagName()=="XMAX")) { - normal = TinyVector<2>(1,0); - } else if ((ref.tagName()=="YMIN") or (ref.tagName()=="YMAX")) { - normal = TinyVector<2>(0,1); - } else { - std::cout << rang::fg::red << "ignoring " << ref << rang::style::reset << '\n'; - break; - } - - std::set<unsigned int> node_id_set; - const Kokkos::View<const unsigned int**> face_nodes = mesh.connectivity().faceNodes(); - - const Kokkos::View<const unsigned int*>& face_id_list = ref_face_list.faceList(); - - for (unsigned int l=0; l<face_id_list.extent(0); ++l) { - for (unsigned short r=0; r<2; ++r) { - node_id_set.insert(face_nodes(face_id_list[l],r)); - } - } - - Kokkos::View<unsigned int*> node_id_list("node_id_list", node_id_set.size()); - { - int r=0; - for (auto node_id : node_id_set) { - node_id_list[r] = node_id; - ++r; - } - } - - RefNodeList ref_node_list(ref, node_id_list); - - SymmetryBoundaryCondition<MeshType::dimension>* sym_bc - = new SymmetryBoundaryCondition<MeshType::dimension>(ref_node_list, normal); - std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc); - bc_list.push_back(BoundaryConditionHandler(bc)); - } } UnknownsType unknowns(mesh_data); diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp index be48b13a694609eca1707d02975221b6c5d593d7..5a2e3c227d71b514478f485f44909054d7628de0 100644 --- a/src/mesh/MeshNodeBoundary.hpp +++ b/src/mesh/MeshNodeBoundary.hpp @@ -4,6 +4,9 @@ #include <Kokkos_Core.hpp> #include <TinyVector.hpp> +#include <RefNodeList.hpp> +#include <RefFaceList.hpp> + #include <iostream> template <size_t dimension> @@ -15,13 +18,19 @@ class MeshNodeBoundary MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default; MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default; + const Kokkos::View<const unsigned int*>& nodeList() const + { + return m_node_list; + } + template <typename MeshType> MeshNodeBoundary(const MeshType& mesh, - const Kokkos::View<const unsigned int*>& face_list) + const RefFaceList& ref_face_list) { static_assert(dimension == MeshType::dimension); const Kokkos::View<const unsigned short*> face_nb_cells = mesh.connectivity().faceNbCells(); + const Kokkos::View<const unsigned int*>& face_list = ref_face_list.faceList(); Kokkos::parallel_for(face_list.extent(0), KOKKOS_LAMBDA(const int& l){ if (face_nb_cells[face_list[l]]>1) { std::cerr << "internal faces cannot be used to define mesh boundaries\n"; @@ -79,13 +88,18 @@ class MeshFlatNodeBoundary template <typename MeshType> inline Rd _getOutgoingNormal(const MeshType& mesh); public: + const Rd& outgoingNormal() const + { + return m_outgoing_normal; + } + MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default; MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default; template <typename MeshType> MeshFlatNodeBoundary(const MeshType& mesh, - const Kokkos::View<const unsigned int*>& face_list) - : MeshNodeBoundary<dimension>(mesh, face_list), + const RefFaceList& ref_face_list) + : MeshNodeBoundary<dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh)) { ; diff --git a/src/scheme/BoundaryCondition.hpp b/src/scheme/BoundaryCondition.hpp index 5486c62375f6b222d124a94701fbb6f47aaf46e7..c59eae3d347d6354a6772c3c2001bfbf412595d5 100644 --- a/src/scheme/BoundaryCondition.hpp +++ b/src/scheme/BoundaryCondition.hpp @@ -7,6 +7,7 @@ #include <Kokkos_Core.hpp> #include <RefNodeList.hpp> +#include <MeshNodeBoundary.hpp> class BoundaryCondition { @@ -84,29 +85,27 @@ public: typedef TinyVector<dimension, double> Rd; private: - const RefNodeList m_ref_node_list; - // Kokkos::View<unsigned int*> m_node_list; - const Rd m_outgoing_normal; + + const MeshFlatNodeBoundary<dimension> m_mesh_flat_node_boundary; public: const Rd& outgoingNormal() const { - return m_outgoing_normal; + return m_mesh_flat_node_boundary.outgoingNormal(); } + size_t numberOfNodes() const { - return m_ref_node_list.nodeList().extent(0); + return m_mesh_flat_node_boundary.nodeList().extent(0); } const Kokkos::View<const unsigned int*> nodeList() const { - return m_ref_node_list.nodeList(); + return m_mesh_flat_node_boundary.nodeList(); } - SymmetryBoundaryCondition(const RefNodeList& ref_node_list, - const Rd& outgoing_normal) + SymmetryBoundaryCondition(const MeshFlatNodeBoundary<dimension>& mesh_flat_node_boundary) : BoundaryCondition(BoundaryCondition::symmetry), - m_ref_node_list(ref_node_list), - m_outgoing_normal(outgoing_normal) + m_mesh_flat_node_boundary(mesh_flat_node_boundary) { ; }