From 9c1d99179f6e21b7f0e92f50eae1d333d3ed755f Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 19 Jun 2018 19:29:50 +0200
Subject: [PATCH] Added first version of mesh boundaries calculation for BC
---
src/main.cpp | 14 +++-
src/mesh/Connectivity2D.hpp | 5 ++
src/mesh/MeshBoundary.hpp | 124 ++++++++++++++++++++++++++++++++++++
3 files changed, 142 insertions(+), 1 deletion(-)
create mode 100644 src/mesh/MeshBoundary.hpp
diff --git a/src/main.cpp b/src/main.cpp
index 58166eda7..187eed116 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -23,6 +23,8 @@
#include <BoundaryConditionDescriptor.hpp>
+#include <MeshBoundary.hpp>
+
#include <GmshReader.hpp>
#include <CLI/CLI.hpp>
@@ -239,12 +241,22 @@ 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;
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);
- std::cout << sym_bc_descriptor << '\n';
+ 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();
+ 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());
+ }
+ }
break;
}
default: {
diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp
index 290e826aa..19ea67e0c 100644
--- a/src/mesh/Connectivity2D.hpp
+++ b/src/mesh/Connectivity2D.hpp
@@ -266,6 +266,11 @@ private:
return m_face_nb_cells;
}
+ const Kokkos::View<const unsigned short*> faceNbNodes() const
+ {
+ return m_face_nb_nodes;
+ }
+
const Kokkos::View<const unsigned int**> nodeCells() const
{
return m_node_cells;
diff --git a/src/mesh/MeshBoundary.hpp b/src/mesh/MeshBoundary.hpp
new file mode 100644
index 000000000..8e3961e52
--- /dev/null
+++ b/src/mesh/MeshBoundary.hpp
@@ -0,0 +1,124 @@
+#ifndef MESH_BOUNDARY_HPP
+#define MESH_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
+{
+ public:
+ MeshBoundary& operator=(const MeshBoundary&) = default;
+ MeshBoundary& operator=(MeshBoundary&&) = default;
+
+ template <typename MeshType>
+ MeshBoundary(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();
+
+ 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";
+ std::exit(1);
+ }
+ });
+ }
+
+ MeshBoundary() = default;
+ MeshBoundary(const MeshBoundary&) = default;
+ MeshBoundary(MeshBoundary&&) = default;
+ virtual ~MeshBoundary() = default;
+};
+
+
+template <size_t dimension>
+class MeshFlatBoundary
+ : public MeshBoundary<dimension>
+{
+ typedef TinyVector<dimension, double> Rd;
+ private:
+ const Rd m_outgoing_normal;
+
+ template <typename MeshType>
+ inline Rd _getOutgoingNormal(const MeshType& mesh,
+ const Kokkos::View<const unsigned int*>& face_list);
+ public:
+ MeshFlatBoundary& operator=(const MeshFlatBoundary&) = default;
+ MeshFlatBoundary& operator=(MeshFlatBoundary&&) = default;
+
+ template <typename MeshType>
+ MeshFlatBoundary(const MeshType& mesh,
+ const Kokkos::View<const unsigned int*>& face_list)
+ : MeshBoundary<dimension>(mesh, face_list),
+ m_outgoing_normal(_getOutgoingNormal(mesh, face_list))
+ {
+ ;
+ }
+
+ MeshFlatBoundary() = default;
+ MeshFlatBoundary(const MeshFlatBoundary&) = default;
+ MeshFlatBoundary(MeshFlatBoundary&&) = default;
+ virtual ~MeshFlatBoundary() = default;
+
+};
+
+template <>
+template <typename MeshType>
+inline TinyVector<2,double>
+MeshFlatBoundary<2>::
+_getOutgoingNormal(const MeshType& mesh,
+ const Kokkos::View<const unsigned int*>& face_list)
+{
+ static_assert(MeshType::dimension == 2);
+ typedef TinyVector<2,double> R2;
+
+ const Kokkos::View<const unsigned short*> face_nb_nodes = mesh.connectivity().faceNbNodes();
+ const Kokkos::View<const unsigned int**> face_nodes = mesh.connectivity().faceNodes();
+ const Kokkos::View<const R2*> xr = mesh.xr();
+
+ R2 xmin(std::numeric_limits<double>::max(),
+ std::numeric_limits<double>::max());
+
+ R2 xmax(-std::numeric_limits<double>::max(),
+ -std::numeric_limits<double>::max());
+
+ if (face_list.extent(0)==0) {
+ std::cerr << "face_list is empty! Cannot compute normal!\n";
+ std::exit(1);
+ }
+
+ for (size_t l=0; l<face_list.extent(0); ++l) {
+ const unsigned int face_number = face_list[l];
+ for (int r=0; r<face_nb_nodes(face_number); ++r) {
+ const R2& x = xr(face_nodes(face_number,r));
+ if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
+ xmin = x;
+ }
+ if ((x[0]>xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
+ xmax = x;
+ }
+ }
+ }
+
+ if (xmin == xmax) {
+ std::cerr << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
+ std::exit(1);
+ }
+
+ R2 dx = xmax-xmin;
+ dx *= 1./std::sqrt((dx,dx));
+
+ R2 normal(-dx[1], dx[0]);
+
+ std::cout << "xmin=" << xmin << " xmax=" << xmax << " normal=" << normal << '\n';
+
+ return normal;
+}
+
+
+#endif // MESH_BOUNDARY_HPP
--
GitLab