From 49759b2895b05974c91f534d1e857aee44cdd861 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Thu, 3 May 2018 08:43:09 +0200
Subject: [PATCH] 2D symmetry boundary conditions almost working
- code is awful, will require a lot of clean-up...
---
src/main.cpp | 48 +++++++++++++++++++---------
src/mesh/Connectivity2D.hpp | 59 +++++++++++++++++++++++++++++++++++
src/mesh/GmshReader.cpp | 46 +++++++++++++++++++++++++--
src/scheme/AcousticSolver.hpp | 6 ++--
4 files changed, 140 insertions(+), 19 deletions(-)
diff --git a/src/main.cpp b/src/main.cpp
index 3cf9e8f8e..190ffae27 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -131,22 +131,41 @@ int main(int argc, char *argv[])
Kokkos::Timer timer;
timer.reset();
- MeshDataType mesh_data(mesh);
+ MeshDataType mesh_data(mesh);
std::vector<BoundaryConditionHandler> bc_list;
- // { // quite dirty!
- // SymmetryBoundaryCondition<MeshType::dimension>* sym_bc0
- // = new SymmetryBoundaryCondition<MeshType::dimension>(std::vector<unsigned int>({0u}),
- // TinyVector<1>(-1));
- // std::shared_ptr<SymmetryBoundaryCondition<1>> bc0(sym_bc0);
- // bc_list.push_back(BoundaryConditionHandler(bc0));
-
- // PressureBoundaryCondition* pres_bc1
- // = new PressureBoundaryCondition(0.1,
- // std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
- // std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1);
- // bc_list.push_back(BoundaryConditionHandler(bc1));
- // }
+ { // quite dirty!
+ for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) {
+ Connectivity2D::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary);
+ std::cout << i_boundary
+ << " -> "
+ << mesh.connectivity().nodesBoundary(i_boundary).second.extent(0)
+ << '\n';
+ unsigned int ref = nodes_boundary.first;
+ TinyVector<2> normal(0,0);
+ if ((ref == 3) or (ref == 4)) {
+ normal = TinyVector<2>(0,1);
+ } else {
+ normal = TinyVector<2>(1,0);
+ }
+ std::cout << "boundary=" << i_boundary << " ref=" << ref << " n=" << normal << '\n';
+ 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));
+ }
+
+ // PressureBoundaryCondition* pres_bc1
+ // = new PressureBoundaryCondition(0.1,
+ // std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
+ // std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1);
+ // bc_list.push_back(BoundaryConditionHandler(bc1));
+ }
UnknownsType unknowns(mesh_data);
@@ -178,7 +197,6 @@ int main(int argc, char *argv[])
if (t+dt>tmax) {
dt=tmax-t;
}
-
acoustic_solver.computeNextStep(t,dt, unknowns);
block_eos.updatePandCFromRhoE();
diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp
index fe7eb44d0..52ba08034 100644
--- a/src/mesh/Connectivity2D.hpp
+++ b/src/mesh/Connectivity2D.hpp
@@ -13,6 +13,10 @@ class Connectivity2D
public:
static constexpr size_t dimension = 2;
+ typedef std::pair<unsigned int, Kokkos::View<const unsigned int*>> FacesBoundary;
+ typedef std::pair<unsigned int, Kokkos::View<const unsigned int*>> NodesBoundary;
+ typedef std::vector<FacesBoundary> FacesBoundaryList;
+ typedef std::vector<NodesBoundary> NodesBoundaryList;
private:
const size_t m_number_of_cells;
size_t m_number_of_faces;
@@ -39,6 +43,9 @@ private:
size_t m_max_nb_node_per_cell;
+ FacesBoundaryList m_faces_boundary_list;
+ NodesBoundaryList m_nodes_boundary_list;
+
void _computeNodeCellConnectivities()
{
// Computes inefficiently node->cells connectivity [Version 0]
@@ -315,6 +322,11 @@ public:
return m_face_cells;
}
+ const Kokkos::View<const unsigned int**> faceNodes() const
+ {
+ return m_face_nodes;
+ }
+
const Kokkos::View<const unsigned short**> nodeCellLocalNode() const
{
return m_node_cell_local_node;
@@ -328,9 +340,56 @@ public:
unsigned int getFaceNumber(const unsigned int node0_id,
const unsigned int node1_id) const
{
+#warning rewrite
+ const unsigned int n0_id = std::min(node0_id, node1_id);
+ const unsigned int n1_id = std::max(node0_id, node1_id);
+ // worst code ever
+ for (unsigned int l=0; l<m_number_of_faces; ++l) {
+ if ((m_face_nodes(l,0) == n0_id) and (m_face_nodes(l,1) == n1_id)) {
+ return l;
+ }
+ }
+ std::cerr << "Face not found!\n";
+ std::exit(0);
return 0;
}
+ size_t numberOfFaceBoundaries() const
+ {
+ return m_faces_boundary_list.size();
+ }
+
+ const FacesBoundary& facesBoundary(const size_t& i) const
+ {
+ return m_faces_boundary_list[i];
+ }
+
+ void addFaceBoundary(const unsigned int& ref,
+ const std::vector<unsigned int>& face_id_vector)
+ {
+ Kokkos::View<unsigned int*> faces_boundary("faces_boundary", face_id_vector.size());
+ Kokkos::parallel_for(face_id_vector.size(), KOKKOS_LAMBDA(const int& l) {
+ faces_boundary[l] = face_id_vector[l];
+ });
+ m_faces_boundary_list.push_back(std::make_pair(ref, faces_boundary));
+ }
+
+ size_t numberOfNodeBoundaries() const
+ {
+ return m_nodes_boundary_list.size();
+ }
+
+ const NodesBoundary& nodesBoundary(const size_t& i) const
+ {
+ return m_nodes_boundary_list[i];
+ }
+
+ void addNodeBoundary(const unsigned int& ref,
+ const Kokkos::View<const unsigned int*> nodes_boundary)
+ {
+ m_nodes_boundary_list.push_back(std::make_pair(ref, nodes_boundary));
+ }
+
Connectivity2D(const Connectivity2D&) = delete;
Connectivity2D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 6577a60bd..8fe4ba632 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -7,6 +7,7 @@
#include <Connectivity2D.hpp>
#include <Mesh.hpp>
+#include <map>
template<typename T>
inline std::string stringify(const T & t)
@@ -800,10 +801,51 @@ GmshReader::__proceedData()
m_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes);
Connectivity2D& connectivity = *m_connectivity;
+ std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
for (unsigned int e=0; e<__edges.extent(0); ++e) {
- unsigned int edge_number = connectivity.getFaceNumber(__edges[e][0], __edges[e][1]);
+ const unsigned int edge_number = connectivity.getFaceNumber(__edges[e][0], __edges[e][1]);
+ const unsigned int ref = __edges_ref[e];
+ ref_faces_map[ref].push_back(edge_number);
}
- // __edges[edgeNumber]
+ for (const auto& ref_face_list : ref_faces_map) {
+ std::cout << ref_face_list.first << ":\n";
+ for (auto face_id : ref_face_list.second) {
+ std::cout << face_id << ' ';
+ }
+ std::cout << '\n';
+ }
+ for (const auto& ref_face_list : ref_faces_map) {
+ connectivity.addFaceBoundary(ref_face_list.first,
+ ref_face_list.second);
+ }
+
+ for (size_t i=0; i<connectivity.numberOfFaceBoundaries(); ++i) {
+ const Connectivity2D::FacesBoundary& faces_boundary
+ = connectivity.facesBoundary(i);
+
+ const unsigned int ref = faces_boundary.first;
+
+ std::set<unsigned int> node_id_set;
+ const Kokkos::View<const unsigned int*> face_ids
+ = faces_boundary.second;
+
+ for (unsigned int l=0; l<face_ids.extent(0); ++l) {
+ for (unsigned short r=0; r<2; ++r) {
+ node_id_set.insert(face_nodes(face_ids[l],r));
+ }
+ }
+
+ Kokkos::View<unsigned int*> node_id_list("node_ids", node_id_set.size());
+ {
+ int r=0;
+ for (auto node_id : node_id_set) {
+ node_id_list[r] = node_id;
+ ++r;
+ }
+ }
+ connectivity.addNodeBoundary(ref, node_id_list);
+ }
+// __edges[edgeNumber]
// = Edge(a, b);
// __edges_ref[edgeNumber] = __references[i];
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 03488f21c..edfe605c4 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -176,23 +176,25 @@ private:
case BoundaryCondition::symmetry: {
const SymmetryBoundaryCondition<dimension>& symmetry_bc
= dynamic_cast<const SymmetryBoundaryCondition<dimension>&>(handler.boundaryCondition());
-
const Rd& n = symmetry_bc.outgoingNormal();
+
const Rdd I = identity;
const Rdd nxn = tensorProduct(n,n);
const Rdd P = I-nxn;
Kokkos::parallel_for(symmetry_bc.numberOfNodes(), KOKKOS_LAMBDA(const int& r_number) {
const int r = symmetry_bc.nodeList()[r_number];
- assert(m_connectivity.nodeNbCells()[r] == 1);
+ // assert(m_connectivity.nodeNbCells()[r] == 1);
m_Ar(r) = P*m_Ar(r)*P + nxn;
m_br(r) = P*m_br(r);
+ std::cout << "A(" << r << ")=" << m_Ar(r) << '\n';
});
break;
}
}
}
+ std::exit(0);
}
Kokkos::View<Rd*>
--
GitLab