diff --git a/src/main.cpp b/src/main.cpp index e1f9260baa1daab57d7615faf26980b13ea221be..6d1016fb58ac752c94cff745958c2ae8b0108c0a 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -131,6 +131,17 @@ int main(int argc, char *argv[]) switch (p_mesh->meshDimension()) { case 1: { + std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"}; + std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list; + for (const auto& sym_boundary_name : sym_boundary_name_list){ + std::shared_ptr<BoundaryDescriptor> boudary_descriptor + = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name)); + SymmetryBoundaryConditionDescriptor* sym_bc_descriptor + = new SymmetryBoundaryConditionDescriptor(boudary_descriptor); + + bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor)); + } + typedef Connectivity1D ConnectivityType; typedef Mesh<ConnectivityType> MeshType; typedef MeshData<MeshType> MeshDataType; @@ -143,27 +154,32 @@ int main(int argc, char *argv[]) MeshDataType mesh_data(mesh); std::vector<BoundaryConditionHandler> bc_list; - // { // quite dirty! - // for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) { - // ConnectivityType::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary); - // unsigned int ref = nodes_boundary.first; - // TinyVector<2> normal(0,0); - // if ((ref == 5) or (ref == 6)) { - // normal = TinyVector<2>(0,1); - // } else { - // normal = TinyVector<2>(1,0); - // } - // const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second; - // std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0)); - // for (size_t r=0; r<nodes_ids.extent(0); ++r) { - // node_boundary_vector[r] = nodes_ids[r]; - // } - // SymmetryBoundaryCondition<MeshType::dimension>* sym_bc - // = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal); - // std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc); - // bc_list.push_back(BoundaryConditionHandler(bc)); - // } - // } + { + for (const auto& bc_descriptor : bc_descriptor_list) { + switch (bc_descriptor->type()) { + case BoundaryConditionDescriptor::Type::symmetry: { + const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor + = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor); + for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefNodeList(); + ++i_ref_node_list) { + const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(i_ref_node_list); + const RefId& ref = ref_node_list.refId(); + if (ref == sym_bc_descriptor.boundaryDescriptor()) { + SymmetryBoundaryCondition<MeshType::dimension>* sym_bc + = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_node_list)); + std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc); + bc_list.push_back(BoundaryConditionHandler(bc)); + } + } + break; + } + default: { + std::cerr << "Unknown BCDescription\n"; + std::exit(1); + } + } + } + } UnknownsType unknowns(mesh_data); @@ -212,6 +228,15 @@ int main(int argc, char *argv[]) method_cost_map["AcousticSolverWithMesh"] = timer.seconds(); + { // gnuplot output for density + const Kokkos::View<const Rd*> xj = mesh_data.xj(); + const Kokkos::View<const double*> rhoj = unknowns.rhoj(); + std::ofstream fout("rho"); + for (size_t j=0; j<mesh.numberOfCells(); ++j) { + fout << xj[j][0] << ' ' << rhoj[j] << '\n'; + } + } + break; } case 2: { @@ -240,8 +265,6 @@ int main(int argc, char *argv[]) std::vector<BoundaryConditionHandler> bc_list; { -#warning Should extract boundaries from mesh itself to avoid inconsistencies - std::map<RefId, MeshNodeBoundary<MeshType::dimension>> ref_boundary; for (const auto& bc_descriptor : bc_descriptor_list) { switch (bc_descriptor->type()) { case BoundaryConditionDescriptor::Type::symmetry: { @@ -252,7 +275,6 @@ int main(int argc, char *argv[]) const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list); const RefId& ref = ref_face_list.refId(); if (ref == sym_bc_descriptor.boundaryDescriptor()) { - std::cout << "Found mesh boundary to impose " << sym_bc_descriptor << '\n'; 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); diff --git a/src/mesh/Connectivity1D.hpp b/src/mesh/Connectivity1D.hpp index 4e26a077827b4f9edf2a798cd5e3da8b039e6837..3bb13bab1b03191344ceca3570eea785f11e31b8 100644 --- a/src/mesh/Connectivity1D.hpp +++ b/src/mesh/Connectivity1D.hpp @@ -7,12 +7,17 @@ #include <TinyVector.hpp> #include <ConnectivityUtils.hpp> +#include <RefId.hpp> +#include <RefNodeList.hpp> + class Connectivity1D { public: static constexpr size_t dimension = 1; private: + std::vector<RefNodeList> m_ref_node_list; + size_t m_number_of_nodes; const size_t& m_number_of_faces = m_number_of_nodes; const size_t m_number_of_cells; @@ -41,8 +46,8 @@ private: Kokkos::View<unsigned int*[2]> cell_nodes("cell_nodes", number_of_cells); Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) { - cell_nodes(j,0) = j; - cell_nodes(j,1) = j+1; + cell_nodes(j,0) = j; + cell_nodes(j,1) = j+1; }); return cell_nodes; @@ -54,13 +59,28 @@ private: { Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", number_of_cells); Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) { - cell_nb_nodes[j] = 2; + cell_nb_nodes[j] = 2; }); return cell_nb_nodes; } public: + void addRefNodeList(const RefNodeList& ref_node_list) + { + m_ref_node_list.push_back(ref_node_list); + } + + size_t numberOfRefNodeList() const + { + return m_ref_node_list.size(); + } + + const RefNodeList& refNodeList(const size_t& i) const + { + return m_ref_node_list[i]; + } + const size_t& numberOfNodes() const { return m_number_of_nodes; @@ -146,40 +166,40 @@ public: { Assert(number_of_cells>0); Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) { - m_inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j]; + m_inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j]; }); ConnectivityUtils utils; utils.computeNodeCellConnectivity(m_cell_nodes, - m_cell_nb_nodes, - m_number_of_cells, - m_max_nb_node_per_cell, - m_number_of_nodes, - m_node_nb_cells, - m_node_cells, - m_node_cell_local_node); + m_cell_nb_nodes, + m_number_of_cells, + m_max_nb_node_per_cell, + m_number_of_nodes, + m_node_nb_cells, + m_node_cells, + m_node_cell_local_node); } Connectivity1D(const Kokkos::View<const unsigned short*> cell_nb_nodes, - const Kokkos::View<const unsigned int**> cell_nodes) + const Kokkos::View<const unsigned int**> cell_nodes) : m_number_of_cells (cell_nb_nodes.extent(0)), m_cell_nodes (cell_nodes), m_cell_nb_nodes (cell_nb_nodes), m_inv_cell_nb_nodes ("inv_cell_nb_nodes", m_number_of_cells) { Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) { - m_inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j]; + m_inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j]; }); ConnectivityUtils utils; utils.computeNodeCellConnectivity(m_cell_nodes, - m_cell_nb_nodes, - m_number_of_cells, - m_max_nb_node_per_cell, - m_number_of_nodes, - m_node_nb_cells, - m_node_cells, - m_node_cell_local_node); + m_cell_nb_nodes, + m_number_of_cells, + m_max_nb_node_per_cell, + m_number_of_nodes, + m_node_nb_cells, + m_node_cells, + m_node_cell_local_node); } ~Connectivity1D() diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp index 19ea67e0ce360669a75e922bc0807216ec66e7b5..1e518651375ba642758787c46610b283bf4698b6 100644 --- a/src/mesh/Connectivity2D.hpp +++ b/src/mesh/Connectivity2D.hpp @@ -19,10 +19,10 @@ class Connectivity2D public: static constexpr size_t dimension = 2; +private: std::vector<RefFaceList> m_ref_face_list; std::vector<RefNodeList> m_ref_node_list; -private: const size_t m_number_of_cells; size_t m_number_of_faces; size_t m_number_of_nodes; diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp index 4eb7dde8099d0a17efb79b655e751bc13b78edc0..5c35495d4cb1cd7cf58d30f52c1f0f227d2a42f0 100644 --- a/src/mesh/GmshReader.cpp +++ b/src/mesh/GmshReader.cpp @@ -841,8 +841,6 @@ GmshReader::__proceedData() connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list)); } - const Kokkos::View<const unsigned int**> face_nodes = connectivity.faceNodes(); - std::map<unsigned int, std::vector<unsigned int>> ref_points_map; for (unsigned int r=0; r<__points.extent(0); ++r) { const unsigned int point_number = __points[r]; @@ -885,7 +883,25 @@ GmshReader::__proceedData() cell_nb_nodes[j] = 2; } - std::shared_ptr<Connectivity1D> connectivity(new Connectivity1D(cell_nb_nodes, cell_nodes)); + std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_nb_nodes, cell_nodes)); + Connectivity1D& connectivity = *p_connectivity; + + std::map<unsigned int, std::vector<unsigned int>> ref_points_map; + for (unsigned int r=0; r<__points.extent(0); ++r) { + const unsigned int point_number = __points[r]; + const unsigned int& ref = __points_ref[r]; + ref_points_map[ref].push_back(point_number); + } + + for (const auto& ref_point_list : ref_points_map) { + Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size()); + for (size_t j=0; j<ref_point_list.second.size(); ++j) { + point_list[j]=ref_point_list.second[j]; + } + const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first); + connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list)); + } + typedef Mesh<Connectivity1D> MeshType; typedef TinyVector<1, double> Rd; @@ -894,7 +910,7 @@ GmshReader::__proceedData() xr[i][0] = __vertices[i][0]; } - m_mesh = std::shared_ptr<IMesh>(new MeshType(connectivity, xr)); + m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr)); } else { std::cerr << "*** using a dimension 0 mesh is forbidden!\n"; diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp index 5a2e3c227d71b514478f485f44909054d7628de0..c720228a455c28f3b1d6b08c55d7e8533e0320a5 100644 --- a/src/mesh/MeshNodeBoundary.hpp +++ b/src/mesh/MeshNodeBoundary.hpp @@ -61,6 +61,14 @@ class MeshNodeBoundary m_node_list = node_list; } + template <typename MeshType> + MeshNodeBoundary(const MeshType& mesh, + const RefNodeList& ref_node_list) + : m_node_list(ref_node_list.nodeList()) + { + static_assert(dimension == MeshType::dimension); + } + MeshNodeBoundary() = default; MeshNodeBoundary(const MeshNodeBoundary&) = default; MeshNodeBoundary(MeshNodeBoundary&&) = default; @@ -105,6 +113,15 @@ class MeshFlatNodeBoundary ; } + template <typename MeshType> + MeshFlatNodeBoundary(const MeshType& mesh, + const RefNodeList& ref_node_list) + : MeshNodeBoundary<dimension>(mesh, ref_node_list), + m_outgoing_normal(_getOutgoingNormal(mesh)) + { + ; + } + MeshFlatNodeBoundary() = default; MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default; MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default; @@ -135,6 +152,23 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal, }); } +template <> +template <typename MeshType> +inline TinyVector<1,double> +MeshFlatNodeBoundary<1>:: +_getNormal(const MeshType& mesh) +{ + static_assert(MeshType::dimension == 1); + typedef TinyVector<1,double> R; + + if (m_node_list.extent(0) != 1) { + std::cerr << "Node boundaries in 1D require to have exactly 1 node\n"; + std::exit(1); + } + + return R(1); +} + template <> template <typename MeshType> inline TinyVector<2,double> @@ -209,5 +243,37 @@ _getOutgoingNormal(const MeshType& mesh) } } +template <> +template <typename MeshType> +inline TinyVector<1,double> +MeshFlatNodeBoundary<1>:: +_getOutgoingNormal(const MeshType& mesh) +{ + static_assert(MeshType::dimension == 1); + typedef TinyVector<1,double> R; + + const R normal = this->_getNormal(mesh); + + const Kokkos::View<const R*>& xr = mesh.xr(); + const Kokkos::View<const unsigned int**>& node_cells = mesh.connectivity().nodeCells(); + const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes(); + const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes(); + + const size_t r0 = m_node_list[0]; + const size_t j0 = node_cells(r0,0); + double max_height = 0; + for (int r=0; r<cell_nb_nodes(j0); ++r) { + const double height = (xr(cell_nodes(j0, r))-xr(r0), normal); + if (std::abs(height) > std::abs(max_height)) { + max_height = height; + } + } + if (max_height > 0) { + return -normal; + } else { + return normal; + } +} + #endif // MESH_NODE_BOUNDARY_HPP