From 32ba54e6253fa319b8692b0c255147e00c7fd34e Mon Sep 17 00:00:00 2001 From: Stephane Del Pino <stephane.delpino44@gmail.com> Date: Thu, 21 Jun 2018 12:34:40 +0200 Subject: [PATCH] Normal to flat boundaries are now computed --- src/main.cpp | 4 +- src/mesh/MeshBoundary.hpp | 121 +++++++++++++++++++++++++++++--------- 2 files changed, 94 insertions(+), 31 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index 187eed116..fea066cf5 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -242,7 +242,7 @@ int main(int argc, char *argv[]) std::vector<BoundaryConditionHandler> bc_list; { #warning Should extract boundaries from mesh itself to avoid inconsistencies - std::map<RefId, MeshBoundary<MeshType::dimension>> ref_boundary; + std::map<RefId, MeshNodeBoundary<MeshType::dimension>> ref_boundary; for (const auto& bc_descriptor : bc_descriptor_list) { switch (bc_descriptor->type()) { case BoundaryConditionDescriptor::Type::symmetry: { @@ -254,7 +254,7 @@ 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] = MeshFlatBoundary<MeshType::dimension>(mesh, ref_face_list.faceList()); + ref_boundary[ref] = MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list.faceList()); } } break; diff --git a/src/mesh/MeshBoundary.hpp b/src/mesh/MeshBoundary.hpp index 46a5f55d0..be48b13a6 100644 --- a/src/mesh/MeshBoundary.hpp +++ b/src/mesh/MeshBoundary.hpp @@ -1,24 +1,23 @@ -#ifndef MESH_BOUNDARY_HPP -#define MESH_BOUNDARY_HPP +#ifndef MESH_NODE_BOUNDARY_HPP +#define MESH_NODE_BOUNDARY_HPP #include <Kokkos_Core.hpp> #include <TinyVector.hpp> #include <iostream> -#warning Must change rewrite that to compute associated node lists template <size_t dimension> -class MeshBoundary +class MeshNodeBoundary { protected: Kokkos::View<const unsigned int*> m_node_list; public: - MeshBoundary& operator=(const MeshBoundary&) = default; - MeshBoundary& operator=(MeshBoundary&&) = default; + MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default; + MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default; template <typename MeshType> - MeshBoundary(const MeshType& mesh, - const Kokkos::View<const unsigned int*>& face_list) + MeshNodeBoundary(const MeshType& mesh, + const Kokkos::View<const unsigned int*>& face_list) { static_assert(dimension == MeshType::dimension); const Kokkos::View<const unsigned short*> face_nb_cells = mesh.connectivity().faceNbCells(); @@ -53,48 +52,80 @@ class MeshBoundary m_node_list = node_list; } - MeshBoundary() = default; - MeshBoundary(const MeshBoundary&) = default; - MeshBoundary(MeshBoundary&&) = default; - virtual ~MeshBoundary() = default; + MeshNodeBoundary() = default; + MeshNodeBoundary(const MeshNodeBoundary&) = default; + MeshNodeBoundary(MeshNodeBoundary&&) = default; + virtual ~MeshNodeBoundary() = default; }; template <size_t dimension> -class MeshFlatBoundary - : public MeshBoundary<dimension> +class MeshFlatNodeBoundary + : public MeshNodeBoundary<dimension> { typedef TinyVector<dimension, double> Rd; private: const Rd m_outgoing_normal; + template <typename MeshType> + inline Rd _getNormal(const MeshType& mesh); + + template <typename MeshType> + inline void _checkBoundaryIsFlat(const TinyVector<2,double>& normal, + const TinyVector<2,double>& xmin, + const TinyVector<2,double>& xmax, + const MeshType& mesh) const; + template <typename MeshType> inline Rd _getOutgoingNormal(const MeshType& mesh); public: - MeshFlatBoundary& operator=(const MeshFlatBoundary&) = default; - MeshFlatBoundary& operator=(MeshFlatBoundary&&) = default; + MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default; + MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default; template <typename MeshType> - MeshFlatBoundary(const MeshType& mesh, - const Kokkos::View<const unsigned int*>& face_list) - : MeshBoundary<dimension>(mesh, face_list), + MeshFlatNodeBoundary(const MeshType& mesh, + const Kokkos::View<const unsigned int*>& face_list) + : MeshNodeBoundary<dimension>(mesh, face_list), m_outgoing_normal(_getOutgoingNormal(mesh)) { ; } - MeshFlatBoundary() = default; - MeshFlatBoundary(const MeshFlatBoundary&) = default; - MeshFlatBoundary(MeshFlatBoundary&&) = default; - virtual ~MeshFlatBoundary() = default; - + MeshFlatNodeBoundary() = default; + MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default; + MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default; + virtual ~MeshFlatNodeBoundary() = default; }; +template<> +template <typename MeshType> +void MeshFlatNodeBoundary<2>:: +_checkBoundaryIsFlat(const TinyVector<2,double>& normal, + const TinyVector<2,double>& xmin, + const TinyVector<2,double>& xmax, + const MeshType& mesh) const +{ + static_assert(MeshType::dimension == 2); + typedef TinyVector<2,double> R2; + const R2 origin = 0.5*(xmin+xmax); + const double length = l2Norm(xmax-xmin); + + const Kokkos::View<const R2*> xr = mesh.xr(); + + Kokkos::parallel_for(m_node_list.extent(0), KOKKOS_LAMBDA(const size_t& r) { + const R2& x = xr[m_node_list[r]]; + if ((x-origin,normal)>1E-13*length) { + std::cerr << "this FlatBoundary is not flat!\n"; + std::exit(1); + } + }); +} + template <> template <typename MeshType> inline TinyVector<2,double> -MeshFlatBoundary<2>:: -_getOutgoingNormal(const MeshType& mesh) +MeshFlatNodeBoundary<2>:: +_getNormal(const MeshType& mesh) { static_assert(MeshType::dimension == 2); typedef TinyVector<2,double> R2; @@ -123,14 +154,46 @@ _getOutgoingNormal(const MeshType& mesh) } R2 dx = xmax-xmin; - dx *= 1./std::sqrt((dx,dx)); + dx *= 1./l2Norm(dx); R2 normal(-dx[1], dx[0]); - std::cout << "xmin=" << xmin << " xmax=" << xmax << " normal=" << normal << '\n'; + this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh); return normal; } +template <> +template <typename MeshType> +inline TinyVector<2,double> +MeshFlatNodeBoundary<2>:: +_getOutgoingNormal(const MeshType& mesh) +{ + static_assert(MeshType::dimension == 2); + typedef TinyVector<2,double> R2; + + const R2 normal = this->_getNormal(mesh); + + const Kokkos::View<const R2*>& 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_BOUNDARY_HPP +#endif // MESH_NODE_BOUNDARY_HPP -- GitLab