From 2569d3e72cf81ab46fcf5ba25d4a6df0e4d549d9 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Fri, 29 Jun 2018 19:33:02 +0200
Subject: [PATCH] Added 3d Cjr calculation
- added a bunch of 3d connectivity structures
- not yet validated
---
src/mesh/Connectivity3D.hpp | 119 +++++++++++++++++++++++++-----------
src/mesh/MeshData.hpp | 58 ++++++++++++++++--
2 files changed, 138 insertions(+), 39 deletions(-)
diff --git a/src/mesh/Connectivity3D.hpp b/src/mesh/Connectivity3D.hpp
index 94e65508b..8abc9b7b0 100644
--- a/src/mesh/Connectivity3D.hpp
+++ b/src/mesh/Connectivity3D.hpp
@@ -14,6 +14,7 @@
#include <RefNodeList.hpp>
#include <RefFaceList.hpp>
+#include <tuple>
class Connectivity3D
{
@@ -33,7 +34,8 @@ private:
Kokkos::View<double*> m_inv_cell_nb_nodes;
Kokkos::View<const unsigned short*> m_cell_nb_faces;
- Kokkos::View<unsigned int**> m_cell_faces;
+ Kokkos::View<const unsigned int**> m_cell_faces;
+ Kokkos::View<const bool**> m_cell_faces_is_reversed;
Kokkos::View<unsigned short*> m_node_nb_cells;
Kokkos::View<unsigned int**> m_node_cells;
@@ -49,11 +51,14 @@ private:
size_t m_max_nb_node_per_cell;
size_t m_max_nb_face_per_cell;
+ size_t m_max_nb_node_per_face;
struct Face
{
+ bool reversed = false;
+
KOKKOS_INLINE_FUNCTION
- std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list) const
+ std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
{
const auto min_id = std::min_element(node_list.begin(), node_list.end());
const size_t shift = std::distance(min_id, node_list.begin());
@@ -63,6 +68,7 @@ private:
for (size_t i=0; i<node_list.size(); ++i) {
const size_t reverse_shift = node_list.size()-shift;
rotated_node_list[(reverse_shift+node_list.size()-i)%node_list.size()] = node_list[i];
+ reversed = true;
}
} else {
for (size_t i=0; i<node_list.size(); ++i) {
@@ -113,8 +119,8 @@ private:
{
Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", m_number_of_cells);
- typedef std::pair<unsigned int, unsigned short> CellFaceId;
- std::map<Face, std::vector<CellFaceId>> face_cells_map;
+ typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo;
+ std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
for (unsigned int j=0; j<m_number_of_cells; ++j) {
const unsigned short cell_nb_nodes = m_cell_nb_nodes[j];
switch (cell_nb_nodes) {
@@ -126,35 +132,46 @@ private:
}
case 8: { // hexahedron
// face 0
- face_cells_map[Face({m_cell_nodes(j,0),
- m_cell_nodes(j,1),
- m_cell_nodes(j,2),
- m_cell_nodes(j,3)})].push_back(std::make_pair(j,0));
+ Face f0({m_cell_nodes(j,0),
+ m_cell_nodes(j,1),
+ m_cell_nodes(j,2),
+ m_cell_nodes(j,3)});
+ face_cells_map[f0].push_back(std::make_tuple(j, 0, f0.reversed));
+
// face 1
- face_cells_map[Face({m_cell_nodes(j,4),
- m_cell_nodes(j,5),
- m_cell_nodes(j,6),
- m_cell_nodes(j,7)})].push_back(std::make_pair(j,1));
+ Face f1({m_cell_nodes(j,4),
+ m_cell_nodes(j,7),
+ m_cell_nodes(j,6),
+ m_cell_nodes(j,5)});
+ face_cells_map[f1].push_back(std::make_tuple(j, 1, f1.reversed));
+
// face 2
- face_cells_map[Face({m_cell_nodes(j,0),
- m_cell_nodes(j,3),
- m_cell_nodes(j,7),
- m_cell_nodes(j,4)})].push_back(std::make_pair(j,2));
+ Face f2({m_cell_nodes(j,0),
+ m_cell_nodes(j,3),
+ m_cell_nodes(j,7),
+ m_cell_nodes(j,4)});
+ face_cells_map[f2].push_back(std::make_tuple(j, 2, f2.reversed));
+
// face 3
- face_cells_map[Face({m_cell_nodes(j,1),
- m_cell_nodes(j,2),
- m_cell_nodes(j,6),
- m_cell_nodes(j,5)})].push_back(std::make_pair(j,3));
+ Face f3({m_cell_nodes(j,1),
+ m_cell_nodes(j,5),
+ m_cell_nodes(j,6),
+ m_cell_nodes(j,2)});
+ face_cells_map[f3].push_back(std::make_tuple(j, 3, f3.reversed));
+
// face 4
- face_cells_map[Face({m_cell_nodes(j,0),
- m_cell_nodes(j,1),
- m_cell_nodes(j,5),
- m_cell_nodes(j,4)})].push_back(std::make_pair(j,4));
+ Face f4({m_cell_nodes(j,0),
+ m_cell_nodes(j,4),
+ m_cell_nodes(j,5),
+ m_cell_nodes(j,1)});
+ face_cells_map[f4].push_back(std::make_tuple(j, 4, f4.reversed));
+
// face 5
- face_cells_map[Face({m_cell_nodes(j,3),
- m_cell_nodes(j,2),
- m_cell_nodes(j,6),
- m_cell_nodes(j,7)})].push_back(std::make_pair(j,5));
+ Face f5({m_cell_nodes(j,3),
+ m_cell_nodes(j,2),
+ m_cell_nodes(j,6),
+ m_cell_nodes(j,7)});
+ face_cells_map[f5].push_back(std::make_tuple(j, 5, f5.reversed));
cell_nb_faces[j] = 6;
break;
@@ -178,9 +195,9 @@ private:
m_face_nb_nodes = face_nb_nodes;
std::cerr << __FILE__ << ':' << __LINE__ << ':'
- << rang::fg::red << " max_nb_node_per_face forced to 4" << rang::fg::reset << '\n';
- const size_t max_nb_node_per_face = 4;
- Kokkos::View<unsigned int**> face_nodes("face_nodes", m_number_of_faces, max_nb_node_per_face);
+ << rang::fg::red << " m_max_nb_node_per_face forced to 4" << rang::fg::reset << '\n';
+ m_max_nb_node_per_face = 4;
+ Kokkos::View<unsigned int**> face_nodes("face_nodes", m_number_of_faces, m_max_nb_node_per_face);
{
int l=0;
for (const auto& face_cells_vector : face_cells_map) {
@@ -212,7 +229,7 @@ private:
for (const auto& face_cells_vector : face_cells_map) {
const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- unsigned int cell_number = cells_vector[lj].first;
+ const unsigned int cell_number = std::get<0>(cells_vector[lj]);
face_cells(l,lj) = cell_number;
}
++l;
@@ -229,8 +246,9 @@ private:
for (const auto& face_cells_vector : face_cells_map) {
const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- unsigned int cell_number = cells_vector[lj].first;
- unsigned short cell_local_face = cells_vector[lj].second;
+ unsigned int cell_number;
+ unsigned short cell_local_face;
+ std::tie(cell_number, cell_local_face, std::ignore) = cells_vector[lj];
cell_faces(cell_number,cell_local_face) = l;
}
++l;
@@ -238,6 +256,22 @@ private:
}
m_cell_faces = cell_faces;
+ Kokkos::View<bool**> cell_faces_is_reversed("cell_faces_is_reversed", m_number_of_faces, m_max_nb_face_per_cell);
+ {
+ int l=0;
+ for (const auto& face_cells_vector : face_cells_map) {
+ const auto& cells_vector = face_cells_vector.second;
+ for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
+ unsigned int cell_number;
+ bool reversed;
+ std::tie(cell_number, std::ignore, reversed);
+ cell_faces_is_reversed(cell_number,lj) = reversed;
+ }
+ ++l;
+ }
+ }
+ m_cell_faces_is_reversed = cell_faces_is_reversed;
+
Kokkos::View<unsigned short**> face_cell_local_face("face_cell_local_face",
m_number_of_faces, m_max_nb_node_per_cell);
{
@@ -245,7 +279,7 @@ private:
for (const auto& face_cells_vector : face_cells_map) {
const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- unsigned short cell_local_face = cells_vector[lj].second;
+ unsigned short cell_local_face = std::get<1>(cells_vector[lj]);
face_cell_local_face(l,lj) = cell_local_face;
}
++l;
@@ -300,11 +334,21 @@ private:
return m_number_of_cells;
}
+ const size_t& maxNbFacePerCell() const
+ {
+ return m_max_nb_face_per_cell;
+ }
+
const size_t& maxNbNodePerCell() const
{
return m_max_nb_node_per_cell;
}
+ const size_t& maxNbNodePerFace() const
+ {
+ return m_max_nb_node_per_face;
+ }
+
const Kokkos::View<const unsigned int**> cellNodes() const
{
return m_cell_nodes;
@@ -315,6 +359,11 @@ private:
return m_cell_faces;
}
+ const Kokkos::View<const bool**> cellFacesIsReversed() const
+ {
+ return m_cell_faces_is_reversed;
+ }
+
const Kokkos::View<const unsigned short*> nodeNbCells() const
{
return m_node_nb_cells;
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index f65f3aa09..59411b7f9 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -4,6 +4,8 @@
#include <Kokkos_Core.hpp>
#include <TinyVector.hpp>
+#include <map>
+
template <typename M>
class MeshData
{
@@ -116,21 +118,69 @@ private:
const Kokkos::View<const Rd*> xr = m_mesh.xr();
- std::cerr << "Cjr are not computed in 3D!\n";
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::View<Rd**> Nlr("Nlr", m_mesh.connectivity().numberOfFaces(), m_mesh.connectivity().maxNbNodePerFace());
+
+ const Kokkos::View<const unsigned int**>& face_nodes
+ = m_mesh.connectivity().faceNodes();
+ const Kokkos::View<const unsigned short*>& face_nb_nodes
+ = m_mesh.connectivity().faceNbNodes();
+ Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
+ const size_t nb_nodes = face_nb_nodes[l];
+ std::vector<Rd> dxr(nb_nodes);
+ for (size_t r=0; r<nb_nodes; ++r) {
+ dxr[r] = xr[face_nodes(l,(r+1)%nb_nodes)] - xr[face_nodes(l,(r+nb_nodes-1)%nb_nodes)];
+ }
+ const double inv_12_nb_nodes = 1./(12.*nb_nodes);
+ for (size_t r=0; r<nb_nodes; ++r) {
+ Rd Nr = zero;
+ const Rd two_dxr = 2*dxr[r];
+ for (size_t s=0; s<nb_nodes; ++s) {
+ Nr += crossProduct((two_dxr - dxr[s]), xr[face_nodes(l,s)]);
+ }
+ Nr *= inv_12_nb_nodes;
+ Nr -= (1./6.)*crossProduct(dxr[r], xr[face_nodes(l,r)]);
+ Nlr(l,r) = Nr;
+ }
+ });
+
+ const Kokkos::View<const unsigned int**> cell_faces
+ = m_mesh.connectivity().cellFaces();
+ const Kokkos::View<const bool**> cell_faces_is_reversed
+ = m_mesh.connectivity().cellFacesIsReversed();
+
+ const Kokkos::View<const unsigned short*> cell_nb_faces
+ = m_mesh.connectivity().cellNbFaces();
+
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
for (int R=0; R<cell_nb_nodes[j]; ++R) {
m_Cjr(j,R) = zero;
}
+ std::map<unsigned int, unsigned short> node_id_to_local;
+ for (size_t R=0; R<cell_nb_nodes[j]; ++R) {
+ node_id_to_local[cell_nodes(j,R)] = R;
+ }
+ for (size_t L=0; L<cell_nb_faces[j]; ++L) {
+ const size_t l = cell_faces(j, L);
+ if (cell_faces_is_reversed(j, L)) {
+ for (size_t rl = 0; rl<face_nb_nodes[l]; ++rl) {
+ m_Cjr(j, node_id_to_local[face_nodes(l,rl)]) -= Nlr(l, rl);
+ }
+ } else {
+ for (size_t rl = 0; rl<face_nb_nodes[l]; ++rl) {
+ m_Cjr(j, node_id_to_local[face_nodes(l,rl)]) += Nlr(l, rl);
+ }
+ }
+ }
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
for (int R=0; R<cell_nb_nodes[j]; ++R) {
const Rd& Cjr = m_Cjr(j,R);
m_ljr(j,R) = l2Norm(Cjr);
}
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
for (int R=0; R<cell_nb_nodes[j]; ++R) {
const Rd& Cjr = m_Cjr(j,R);
const double inv_ljr = 1./m_ljr(j,R);
--
GitLab