Skip to content
Snippets Groups Projects
Commit 2569d3e7 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Added 3d Cjr calculation

- added a bunch of 3d connectivity structures
- not yet validated
parent d16d62f7
No related branches found
No related tags found
No related merge requests found
...@@ -14,6 +14,7 @@ ...@@ -14,6 +14,7 @@
#include <RefNodeList.hpp> #include <RefNodeList.hpp>
#include <RefFaceList.hpp> #include <RefFaceList.hpp>
#include <tuple>
class Connectivity3D class Connectivity3D
{ {
...@@ -33,7 +34,8 @@ private: ...@@ -33,7 +34,8 @@ private:
Kokkos::View<double*> m_inv_cell_nb_nodes; Kokkos::View<double*> m_inv_cell_nb_nodes;
Kokkos::View<const unsigned short*> m_cell_nb_faces; 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 short*> m_node_nb_cells;
Kokkos::View<unsigned int**> m_node_cells; Kokkos::View<unsigned int**> m_node_cells;
...@@ -49,11 +51,14 @@ private: ...@@ -49,11 +51,14 @@ private:
size_t m_max_nb_node_per_cell; size_t m_max_nb_node_per_cell;
size_t m_max_nb_face_per_cell; size_t m_max_nb_face_per_cell;
size_t m_max_nb_node_per_face;
struct Face struct Face
{ {
bool reversed = false;
KOKKOS_INLINE_FUNCTION 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 auto min_id = std::min_element(node_list.begin(), node_list.end());
const size_t shift = std::distance(min_id, node_list.begin()); const size_t shift = std::distance(min_id, node_list.begin());
...@@ -63,6 +68,7 @@ private: ...@@ -63,6 +68,7 @@ private:
for (size_t i=0; i<node_list.size(); ++i) { for (size_t i=0; i<node_list.size(); ++i) {
const size_t reverse_shift = node_list.size()-shift; const size_t reverse_shift = node_list.size()-shift;
rotated_node_list[(reverse_shift+node_list.size()-i)%node_list.size()] = node_list[i]; rotated_node_list[(reverse_shift+node_list.size()-i)%node_list.size()] = node_list[i];
reversed = true;
} }
} else { } else {
for (size_t i=0; i<node_list.size(); ++i) { for (size_t i=0; i<node_list.size(); ++i) {
...@@ -113,8 +119,8 @@ private: ...@@ -113,8 +119,8 @@ private:
{ {
Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", m_number_of_cells); Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", m_number_of_cells);
typedef std::pair<unsigned int, unsigned short> CellFaceId; typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo;
std::map<Face, std::vector<CellFaceId>> face_cells_map; std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
for (unsigned int j=0; j<m_number_of_cells; ++j) { for (unsigned int j=0; j<m_number_of_cells; ++j) {
const unsigned short cell_nb_nodes = m_cell_nb_nodes[j]; const unsigned short cell_nb_nodes = m_cell_nb_nodes[j];
switch (cell_nb_nodes) { switch (cell_nb_nodes) {
...@@ -126,35 +132,46 @@ private: ...@@ -126,35 +132,46 @@ private:
} }
case 8: { // hexahedron case 8: { // hexahedron
// face 0 // face 0
face_cells_map[Face({m_cell_nodes(j,0), Face f0({m_cell_nodes(j,0),
m_cell_nodes(j,1), m_cell_nodes(j,1),
m_cell_nodes(j,2), m_cell_nodes(j,2),
m_cell_nodes(j,3)})].push_back(std::make_pair(j,0)); m_cell_nodes(j,3)});
face_cells_map[f0].push_back(std::make_tuple(j, 0, f0.reversed));
// face 1 // face 1
face_cells_map[Face({m_cell_nodes(j,4), Face f1({m_cell_nodes(j,4),
m_cell_nodes(j,5), m_cell_nodes(j,7),
m_cell_nodes(j,6), m_cell_nodes(j,6),
m_cell_nodes(j,7)})].push_back(std::make_pair(j,1)); m_cell_nodes(j,5)});
face_cells_map[f1].push_back(std::make_tuple(j, 1, f1.reversed));
// face 2 // face 2
face_cells_map[Face({m_cell_nodes(j,0), Face f2({m_cell_nodes(j,0),
m_cell_nodes(j,3), m_cell_nodes(j,3),
m_cell_nodes(j,7), m_cell_nodes(j,7),
m_cell_nodes(j,4)})].push_back(std::make_pair(j,2)); m_cell_nodes(j,4)});
face_cells_map[f2].push_back(std::make_tuple(j, 2, f2.reversed));
// face 3 // face 3
face_cells_map[Face({m_cell_nodes(j,1), Face f3({m_cell_nodes(j,1),
m_cell_nodes(j,2), m_cell_nodes(j,5),
m_cell_nodes(j,6), m_cell_nodes(j,6),
m_cell_nodes(j,5)})].push_back(std::make_pair(j,3)); m_cell_nodes(j,2)});
face_cells_map[f3].push_back(std::make_tuple(j, 3, f3.reversed));
// face 4 // face 4
face_cells_map[Face({m_cell_nodes(j,0), Face f4({m_cell_nodes(j,0),
m_cell_nodes(j,1), m_cell_nodes(j,4),
m_cell_nodes(j,5), m_cell_nodes(j,5),
m_cell_nodes(j,4)})].push_back(std::make_pair(j,4)); m_cell_nodes(j,1)});
face_cells_map[f4].push_back(std::make_tuple(j, 4, f4.reversed));
// face 5 // face 5
face_cells_map[Face({m_cell_nodes(j,3), Face f5({m_cell_nodes(j,3),
m_cell_nodes(j,2), m_cell_nodes(j,2),
m_cell_nodes(j,6), m_cell_nodes(j,6),
m_cell_nodes(j,7)})].push_back(std::make_pair(j,5)); m_cell_nodes(j,7)});
face_cells_map[f5].push_back(std::make_tuple(j, 5, f5.reversed));
cell_nb_faces[j] = 6; cell_nb_faces[j] = 6;
break; break;
...@@ -178,9 +195,9 @@ private: ...@@ -178,9 +195,9 @@ private:
m_face_nb_nodes = face_nb_nodes; m_face_nb_nodes = face_nb_nodes;
std::cerr << __FILE__ << ':' << __LINE__ << ':' std::cerr << __FILE__ << ':' << __LINE__ << ':'
<< rang::fg::red << " max_nb_node_per_face forced to 4" << rang::fg::reset << '\n'; << rang::fg::red << " m_max_nb_node_per_face forced to 4" << rang::fg::reset << '\n';
const size_t max_nb_node_per_face = 4; m_max_nb_node_per_face = 4;
Kokkos::View<unsigned int**> face_nodes("face_nodes", m_number_of_faces, max_nb_node_per_face); Kokkos::View<unsigned int**> face_nodes("face_nodes", m_number_of_faces, m_max_nb_node_per_face);
{ {
int l=0; int l=0;
for (const auto& face_cells_vector : face_cells_map) { for (const auto& face_cells_vector : face_cells_map) {
...@@ -212,7 +229,7 @@ private: ...@@ -212,7 +229,7 @@ private:
for (const auto& face_cells_vector : face_cells_map) { for (const auto& face_cells_vector : face_cells_map) {
const auto& cells_vector = face_cells_vector.second; const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) { 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; face_cells(l,lj) = cell_number;
} }
++l; ++l;
...@@ -229,8 +246,9 @@ private: ...@@ -229,8 +246,9 @@ private:
for (const auto& face_cells_vector : face_cells_map) { for (const auto& face_cells_vector : face_cells_map) {
const auto& cells_vector = face_cells_vector.second; const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) { for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
unsigned int cell_number = cells_vector[lj].first; unsigned int cell_number;
unsigned short cell_local_face = cells_vector[lj].second; 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; cell_faces(cell_number,cell_local_face) = l;
} }
++l; ++l;
...@@ -238,6 +256,22 @@ private: ...@@ -238,6 +256,22 @@ private:
} }
m_cell_faces = cell_faces; 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", Kokkos::View<unsigned short**> face_cell_local_face("face_cell_local_face",
m_number_of_faces, m_max_nb_node_per_cell); m_number_of_faces, m_max_nb_node_per_cell);
{ {
...@@ -245,7 +279,7 @@ private: ...@@ -245,7 +279,7 @@ private:
for (const auto& face_cells_vector : face_cells_map) { for (const auto& face_cells_vector : face_cells_map) {
const auto& cells_vector = face_cells_vector.second; const auto& cells_vector = face_cells_vector.second;
for (unsigned short lj=0; lj<cells_vector.size(); ++lj) { 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; face_cell_local_face(l,lj) = cell_local_face;
} }
++l; ++l;
...@@ -300,11 +334,21 @@ private: ...@@ -300,11 +334,21 @@ private:
return m_number_of_cells; return m_number_of_cells;
} }
const size_t& maxNbFacePerCell() const
{
return m_max_nb_face_per_cell;
}
const size_t& maxNbNodePerCell() const const size_t& maxNbNodePerCell() const
{ {
return m_max_nb_node_per_cell; 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 const Kokkos::View<const unsigned int**> cellNodes() const
{ {
return m_cell_nodes; return m_cell_nodes;
...@@ -315,6 +359,11 @@ private: ...@@ -315,6 +359,11 @@ private:
return m_cell_faces; 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 const Kokkos::View<const unsigned short*> nodeNbCells() const
{ {
return m_node_nb_cells; return m_node_nb_cells;
......
...@@ -4,6 +4,8 @@ ...@@ -4,6 +4,8 @@
#include <Kokkos_Core.hpp> #include <Kokkos_Core.hpp>
#include <TinyVector.hpp> #include <TinyVector.hpp>
#include <map>
template <typename M> template <typename M>
class MeshData class MeshData
{ {
...@@ -116,11 +118,59 @@ private: ...@@ -116,11 +118,59 @@ private:
const Kokkos::View<const Rd*> xr = m_mesh.xr(); const Kokkos::View<const Rd*> xr = m_mesh.xr();
std::cerr << "Cjr are not computed in 3D!\n"; 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) { Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
for (int R=0; R<cell_nb_nodes[j]; ++R) { for (int R=0; R<cell_nb_nodes[j]; ++R) {
m_Cjr(j,R) = zero; 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) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment