diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 3f453aed9f58cd88b1cdb4345b0cc95e5f4c9888..71f1d1c2d787ad1069fe57687125b3599c917567 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -102,6 +102,8 @@ Connectivity(const ConnectivityDescriptor& descriptor)
face_owner = convert_to_array(descriptor.face_owner_vector);
m_face_owner = face_owner;
}
+
+ m_ref_face_list = descriptor.refFaceList();
}
}
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index d38c0e1630ae54485974713d08f054922deb08fb..b2c7bfec713c1f1cbbbfe85f6ce2fbe35584917f 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -39,6 +39,7 @@
class ConnectivityDescriptor
{
private:
+ std::vector<RefFaceList> m_ref_face_list;
public:
std::vector<std::vector<unsigned int>> cell_by_node_vector;
@@ -60,6 +61,16 @@ class ConnectivityDescriptor
std::vector<int> face_owner_vector;
std::vector<int> node_owner_vector;
+ const std::vector<RefFaceList>& refFaceList() const
+ {
+ return m_ref_face_list;
+ }
+
+ void addRefFaceList(const RefFaceList& ref_face_list)
+ {
+ m_ref_face_list.push_back(ref_face_list);
+ }
+
ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
@@ -396,6 +407,7 @@ class Connectivity final
return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces);
}
+ [[deprecated("Use connectivity descriptor")]]
void addRefFaceList(const RefFaceList& ref_face_list)
{
m_ref_face_list.push_back(ref_face_list);
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 2484795d1d333b9fa82f37a39f47b0af24e57541..bc7c9cdb37c5d5c464c15a2e868721c8d369e3b9 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -950,6 +950,161 @@ void GmshReader::_dispatch()
}
}
+ // Getting references
+ Array<const size_t> number_of_ref_face_list_per_proc
+ = parallel::allGather(mesh.connectivity().numberOfRefFaceList());
+
+ const size_t number_of_face_list_sender
+ = [&] () {
+ size_t number_of_face_list_sender=0;
+ for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+ number_of_face_list_sender
+ += (number_of_ref_face_list_per_proc[i_rank] > 0);
+ }
+ return number_of_face_list_sender;
+ }();
+
+ if (number_of_face_list_sender > 0) {
+ if (number_of_face_list_sender > 1) {
+ perr() << __FILE__ << ':' << __LINE__ << ": "
+ << rang::fgB::red
+ <<"need to check that knowing procs know the same ref_face_lists!"
+ << rang::fg::reset << '\n';
+ }
+
+ if (number_of_face_list_sender < parallel::size()) {
+ const size_t sender_rank
+ = [&] () {
+ size_t i_rank = 0;
+ for (; i_rank < parallel::size(); ++i_rank) {
+ if (number_of_ref_face_list_per_proc[i_rank] > 0) {
+ break;
+ }
+ }
+ return i_rank;
+ }();
+
+ Assert(number_of_face_list_sender < parallel::size());
+
+ // sending references tags
+ Array<RefId::TagNumberType> ref_tag_list{number_of_ref_face_list_per_proc[sender_rank]};
+ if (parallel::rank() == sender_rank){
+ for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
+ ++i_ref_face_list) {
+ auto ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
+ ref_tag_list[i_ref_face_list] = ref_face_list.refId().tagNumber();
+ }
+ }
+ parallel::broadcast(ref_tag_list, sender_rank);
+
+ // sending references name size
+ Array<size_t> ref_name_size_list{number_of_ref_face_list_per_proc[sender_rank]};
+ if (parallel::rank() == sender_rank){
+ for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
+ ++i_ref_face_list) {
+ auto ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
+ ref_name_size_list[i_ref_face_list] = ref_face_list.refId().tagName().size();
+ }
+ }
+ parallel::broadcast(ref_name_size_list, sender_rank);
+
+ // sending references name size
+ Array<RefId::TagNameType::value_type> ref_name_cat{Sum{ref_name_size_list}};
+ if (parallel::rank() == sender_rank){
+ size_t i_char=0;
+ for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
+ ++i_ref_face_list) {
+ auto ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
+ for (auto c : ref_face_list.refId().tagName()) {
+ ref_name_cat[i_char++] = c;
+ }
+ }
+ }
+ parallel::broadcast(ref_name_cat, sender_rank);
+
+ std::vector<RefId> ref_id_list
+ = [&] () {
+ std::vector<RefId> ref_id_list;
+ ref_id_list.reserve(ref_name_size_list.size());
+ size_t begining=0;
+ for (size_t i_ref=0; i_ref < ref_name_size_list.size(); ++i_ref) {
+ const size_t size = ref_name_size_list[i_ref];
+ ref_id_list.emplace_back(ref_tag_list[i_ref],
+ std::string{&(ref_name_cat[begining]), size});
+ begining += size;
+ }
+ return ref_id_list;
+ } ();
+
+ using block_type = int32_t;
+ constexpr size_t block_size = sizeof(block_type);
+ const size_t nb_block = ref_id_list.size()/block_size + (ref_id_list.size()%block_size != 0);
+ for (size_t i_block=0; i_block<nb_block; ++i_block) {
+ FaceValue<block_type> face_references(mesh.connectivity());
+ face_references.fill(0);
+
+ if (mesh.connectivity().numberOfRefFaceList() > 0) {
+ const size_t max_i_ref = std::min(ref_id_list.size(), block_size*(i_block+1));
+ for (size_t i_ref=block_size*i_block, i=0; i_ref<max_i_ref; ++i_ref, ++i) {
+ block_type ref_bit{1<<i};
+ auto ref_face_list = mesh.connectivity().refFaceList(i_ref);
+
+ const auto& face_list = ref_face_list.faceList();
+ for (size_t i_face=0; i_face<face_list.size(); ++i_face) {
+ const FaceId& face_id = face_list[i_face];
+ face_references[face_id] |= ref_bit;
+ }
+ }
+ }
+
+ std::vector<Array<const block_type>> send_face_refs_by_proc(parallel::size());
+ for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+ Array<block_type> send_face_refs(nb_face_to_send_by_proc[i_rank]);
+ const Array<const FaceId> send_face_id = send_face_id_by_proc[i_rank];
+ parallel_for(send_face_id.size(), PASTIS_LAMBDA(const size_t& l) {
+ const FaceId& face_id = send_face_id[l];
+ send_face_refs[l] = face_references[face_id];
+ });
+ send_face_refs_by_proc[i_rank] = send_face_refs;
+ }
+
+ std::vector<Array<block_type>> recv_face_refs_by_proc(parallel::size());
+ for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+ recv_face_refs_by_proc[i_rank] = Array<block_type>(nb_face_to_recv_by_proc[i_rank]);
+ }
+ parallel::exchange(send_face_refs_by_proc, recv_face_refs_by_proc);
+
+ std::vector<block_type> face_refs(new_descriptor.face_number_vector.size());
+ for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+ for (size_t r=0; r<recv_face_owner_by_proc[i_rank].size(); ++r) {
+ const FaceId& face_id = recv_face_id_correspondance_by_proc[i_rank][r];
+ face_refs[face_id] = recv_face_refs_by_proc[i_rank][r];
+ }
+ }
+
+ const size_t max_i_ref = std::min(ref_id_list.size(), block_size*(i_block+1));
+ for (size_t i_ref=block_size*i_block, i=0; i_ref<max_i_ref; ++i_ref, ++i) {
+ block_type ref_bit{1<<i};
+
+ std::vector<FaceId> face_id_vector;
+
+ for (size_t i_face=0; i_face<face_refs.size(); ++i_face) {
+ const FaceId face_id{i_face};
+ if (face_refs[face_id] & ref_bit) {
+ face_id_vector.push_back(face_id);
+ }
+ }
+
+ Array<const FaceId> face_id_array = convert_to_array(face_id_vector);
+
+ new_descriptor.addRefFaceList(RefFaceList(ref_id_list[i_ref], face_id_array));
+ }
+
+ pout() << __FILE__ << ':' << __LINE__ << ": remains to build lists\n";
+ }
+ }
+ }
+
}
using ConnectivityType = Connectivity<Dimension>;
@@ -982,94 +1137,6 @@ void GmshReader::_dispatch()
<< rang::fg::reset << '\n';
}
- Array<const size_t> number_of_ref_face_list_per_proc
- = parallel::allGather(mesh.connectivity().numberOfRefFaceList());
-
- const size_t number_of_face_list_sender
- = [&] () {
- size_t number_of_face_list_sender=0;
- for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
- number_of_face_list_sender
- += (number_of_ref_face_list_per_proc[i_rank] > 0);
- }
- return number_of_face_list_sender;
- }();
-
- if (number_of_face_list_sender > 0) {
- if (number_of_face_list_sender > 1) {
- perr() << __FILE__ << ':' << __LINE__ << ": "
- << rang::fgB::red
- <<"need to check that knowing procs know the same ref_face_lists!"
- << rang::fg::reset << '\n';
- }
-
- if (number_of_face_list_sender < parallel::size()) {
- const size_t sender_rank
- = [&] () {
- size_t i_rank = 0;
- for (; i_rank < parallel::size(); ++i_rank) {
- if (number_of_ref_face_list_per_proc[i_rank] > 0) {
- break;
- }
- }
- return i_rank;
- }();
-
- Assert(number_of_face_list_sender < parallel::size());
-
- // sending references tags
- Array<RefId::TagNumberType> ref_tag_list{number_of_ref_face_list_per_proc[sender_rank]};
- if (parallel::rank() == sender_rank){
- for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
- ++i_ref_face_list) {
- auto ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
- ref_tag_list[i_ref_face_list] = ref_face_list.refId().tagNumber();
- }
- }
- parallel::broadcast(ref_tag_list, sender_rank);
-
- // sending references name size
- Array<size_t> ref_name_size_list{number_of_ref_face_list_per_proc[sender_rank]};
- if (parallel::rank() == sender_rank){
- for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
- ++i_ref_face_list) {
- auto ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
- ref_name_size_list[i_ref_face_list] = ref_face_list.refId().tagName().size();
- }
- }
- parallel::broadcast(ref_name_size_list, sender_rank);
-
- // sending references name size
- Array<RefId::TagNameType::value_type> ref_name_cat{Sum{ref_name_size_list}};
- if (parallel::rank() == sender_rank){
- size_t i_char=0;
- for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
- ++i_ref_face_list) {
- auto ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
- for (auto c : ref_face_list.refId().tagName()) {
- ref_name_cat[i_char++] = c;
- }
- }
- }
- parallel::broadcast(ref_name_cat, sender_rank);
-
- std::vector<RefId> ref_id_list
- = [&] () {
- std::vector<RefId> ref_id_list;
- ref_id_list.reserve(ref_name_size_list.size());
- size_t begining=0;
- for (size_t i_ref=0; i_ref < ref_name_size_list.size(); ++i_ref) {
- const size_t size = ref_name_size_list[i_ref];
- ref_id_list.emplace_back(ref_tag_list[i_ref],
- std::string{&(ref_name_cat[begining]), size});
- begining += size;
- }
- return ref_id_list;
- } ();
-
- }
- }
-
// Finally build the mesh
NodeValue<Rd> new_xr(*p_connectivity);
for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index a52a91c99d61c028fdc9f634c896d0ce06aa15dd..3d470179ffd541130593e1a308ee7a0907f2738b 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -14,6 +14,8 @@
#include <IConnectivity.hpp>
#include <iostream>
+#include <Messenger.hpp>
+
template <size_t dimension>
class MeshNodeBoundary
{
@@ -212,6 +214,21 @@ _getNormal(const MeshType& mesh)
}
}
+ Array<R2> xmin_array = parallel::allGather(xmin);
+ Array<R2> xmax_array = parallel::allGather(xmax);
+ for (size_t i=0; i<xmin_array.size(); ++i) {
+ const R2& x = xmin_array[i];
+ if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
+ xmin = x;
+ }
+ }
+ for (size_t i=0; i<xmax_array.size(); ++i) {
+ const R2& x = xmax_array[i];
+ if ((x[0]>xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
+ xmax = x;
+ }
+ }
+
if (xmin == xmax) {
perr() << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
std::exit(1);
@@ -360,23 +377,35 @@ _getOutgoingNormal(const MeshType& mesh)
const R2 normal = this->_getNormal(mesh);
- const NodeValue<const R2>& xr = mesh.xr();
- const auto& cell_to_node_matrix
- = mesh.connectivity().cellToNodeMatrix();
+ double max_height = 0;
- const auto& node_to_cell_matrix
- = mesh.connectivity().nodeToCellMatrix();
+ if (m_node_list.size()>0) {
+ const NodeValue<const R2>& xr = mesh.xr();
+ const auto& cell_to_node_matrix
+ = mesh.connectivity().cellToNodeMatrix();
+
+ const auto& node_to_cell_matrix
+ = mesh.connectivity().nodeToCellMatrix();
+
+ const NodeId r0 = m_node_list[0];
+ const CellId j0 = node_to_cell_matrix[r0][0];
+ const auto& j0_nodes = cell_to_node_matrix[j0];
+ for (size_t r=0; r<j0_nodes.size(); ++r) {
+ const double height = (xr[j0_nodes[r]]-xr[r0], normal);
+ if (std::abs(height) > std::abs(max_height)) {
+ max_height = height;
+ }
+ }
+ }
- const NodeId r0 = m_node_list[0];
- const CellId j0 = node_to_cell_matrix[r0][0];
- const auto& j0_nodes = cell_to_node_matrix[j0];
- double max_height = 0;
- for (size_t r=0; r<j0_nodes.size(); ++r) {
- const double height = (xr[j0_nodes[r]]-xr[r0], normal);
+ Array<double> max_height_array = parallel::allGather(max_height);
+ for (size_t i=0; i<max_height_array.size(); ++i) {
+ const double height = max_height_array[i];
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
}
+
if (max_height > 0) {
return -normal;
} else {
diff --git a/src/mesh/RefFaceList.hpp b/src/mesh/RefFaceList.hpp
index ab66943936c091735c404b287d587a15d81a1a9f..cbc7d3b997983c0bdc03d9513684cde8bd669873 100644
--- a/src/mesh/RefFaceList.hpp
+++ b/src/mesh/RefFaceList.hpp
@@ -7,8 +7,8 @@
class RefFaceList
{
private:
- const RefId m_ref_id;
- const Array<const FaceId> m_face_id_list;
+ RefId m_ref_id;
+ Array<const FaceId> m_face_id_list;
public:
const RefId& refId() const