preCICE v3.1.1
Loading...
Searching...
No Matches
ReceivedPartition.cpp
Go to the documentation of this file.
2#include <algorithm>
3#include <map>
4#include <memory>
5#include <ostream>
6#include <utility>
7#include <vector>
8
10#include "com/Extra.hpp"
11#include "com/SharedPointer.hpp"
12#include "logging/LogMacros.hpp"
13#include "m2n/M2N.hpp"
14#include "mapping/Mapping.hpp"
16#include "mesh/BoundingBox.hpp"
17#include "mesh/Filter.hpp"
18#include "mesh/Mesh.hpp"
19#include "mesh/Vertex.hpp"
22#include "profiling/Event.hpp"
23#include "utils/IntraComm.hpp"
24#include "utils/algorithm.hpp"
25#include "utils/assertion.hpp"
26#include "utils/fmt.hpp"
27
29
30namespace precice::partition {
31
33 const mesh::PtrMesh &mesh, GeometricFilter geometricFilter, double safetyFactor, bool allowDirectAccess)
34 : Partition(mesh),
35 _geometricFilter(geometricFilter),
36 _bb(mesh->getDimensions()),
37 _dimensions(mesh->getDimensions()),
38 _safetyFactor(safetyFactor),
39 _allowDirectAccess(allowDirectAccess)
40{
41}
42
44{
46 PRECICE_ASSERT(_mesh->empty());
47
48 // for two-level initialization, receive mesh partitions
49 if (m2n().usesTwoLevelInitialization()) {
50 PRECICE_INFO("Receive mesh partitions for mesh {}", _mesh->getName());
51 Event e("partition.receiveMeshPartitions." + _mesh->getName(), profiling::Synchronize);
52
54 // Primary rank receives remote mesh's global number of vertices
55 int globalNumberOfVertices = -1;
56 m2n().getPrimaryRankCommunication()->receive(globalNumberOfVertices, 0);
57 _mesh->setGlobalNumberOfVertices(globalNumberOfVertices);
58 }
59
60 // each rank receives max/min global vertex indices from connected remote ranks
63 // each rank receives mesh partition from connected remote ranks
65
66 } else {
67 // for one-level initialization receive complete mesh on primary rank
68 PRECICE_INFO("Receive global mesh {}", _mesh->getName());
69 Event e("partition.receiveGlobalMesh." + _mesh->getName(), profiling::Synchronize);
70
72 // a ReceivedPartition can only have one communication, @todo nicer design
73 com::receiveMesh(*(m2n().getPrimaryRankCommunication()), 0, *_mesh);
74 _mesh->setGlobalNumberOfVertices(_mesh->nVertices());
75 }
76 }
77
78 // for both initialization concepts broadcast and set the global number of vertices
80 utils::IntraComm::getCommunication()->broadcast(_mesh->getGlobalNumberOfVertices());
81 }
83 int globalNumberOfVertices = -1;
84 utils::IntraComm::getCommunication()->broadcast(globalNumberOfVertices, 0);
85 PRECICE_ASSERT(globalNumberOfVertices >= 0);
86 _mesh->setGlobalNumberOfVertices(globalNumberOfVertices);
87 }
88}
89
91{
93
94 // handle coupling mode first (i.e. serial participant)
96 PRECICE_DEBUG("Handle partition data structures for serial participant");
97
99 // Prepare the bounding boxes
101 // Filter out vertices not laying in the bounding box
102 mesh::Mesh filteredMesh("FilteredMesh", _dimensions, mesh::Mesh::MESH_ID_UNDEFINED);
103 // To discuss: maybe check this somewhere in the ParticipantImpl, as we have now a similar check for the parallel case
104 PRECICE_CHECK(!_bb.empty(), "You are running this participant in serial mode and the bounding box on mesh \"{}\", is empty. Did you call setMeshAccessRegion with valid data?", _mesh->getName());
105 unsigned int nFilteredVertices = 0;
106 mesh::filterMesh(filteredMesh, *_mesh, [&](const mesh::Vertex &v) { if(!_bb.contains(v))
107 ++nFilteredVertices;
108 return _bb.contains(v); });
109
110 PRECICE_WARN_IF(nFilteredVertices > 0,
111 "{} vertices on mesh \"{}\" have been filtered out due to the defined bounding box in \"setMeshAccessRegion\" "
112 "in serial mode. Associated data values of the filtered vertices will be filled with zero values in order to provide valid data for other participants when reading data.",
113 nFilteredVertices, _mesh->getName());
114
115 _mesh->clear();
116 _mesh->addMesh(filteredMesh);
117 }
118
119 mesh::Mesh::VertexDistribution vertexDistribution;
120 int vertexCounter = 0;
121 for (mesh::Vertex &v : _mesh->vertices()) {
122 v.setOwner(true);
123 vertexDistribution[0].push_back(vertexCounter);
124 vertexCounter++;
125 }
126 PRECICE_ASSERT(_mesh->getVertexDistribution().empty());
127 _mesh->setVertexDistribution(std::move(vertexDistribution));
128 PRECICE_ASSERT(_mesh->getVertexOffsets().empty());
129 _mesh->setVertexOffsets({vertexCounter});
130 return;
131 }
132
133 // check to prevent false configuration
136 "The received mesh {} needs a mapping, either from it, to it, or both. Maybe you don't want to receive this mesh at all?",
137 _mesh->getName());
138 }
139
140 // To better understand steps (2) to (5), it is recommended to look at BU's thesis, especially Figure 69 on page 89
141 // for RBF-based filtering. https://mediatum.ub.tum.de/doc/1320661/document.pdf
142
143 // (1) Bounding-Box-Filter
145
146 // (2) Tag vertices 1st round (i.e. who could be owned by this rank)
147 PRECICE_DEBUG("Tag vertices for filtering: 1st round.");
148 // go to both meshes, vertex is tagged if already one mesh tags him
150
151 // (3) Define which vertices are owned by this rank
152 PRECICE_DEBUG("Create owner information.");
154
155 // (4) Tag vertices 2nd round (what should be filtered out)
156 PRECICE_DEBUG("Tag vertices for filtering: 2nd round.");
158
159 // (5) Filter mesh according to tag
160 PRECICE_INFO("Filter mesh {} by mappings", _mesh->getName());
161 Event e5("partition.filterMeshMappings" + _mesh->getName(), profiling::Synchronize);
162 mesh::Mesh filteredMesh("FilteredMesh", _dimensions, mesh::Mesh::MESH_ID_UNDEFINED);
163 mesh::filterMesh(filteredMesh, *_mesh, [&](const mesh::Vertex &v) { return v.isTagged(); });
164 PRECICE_DEBUG("Mapping filter, filtered from {} to {} vertices, {} to {} edges, and {} to {} triangles.",
165 _mesh->nVertices(), filteredMesh.nVertices(),
166 _mesh->edges().size(), filteredMesh.edges().size(),
167 _mesh->triangles().size(), filteredMesh.triangles().size());
168
169 _mesh->clear();
170 _mesh->addMesh(filteredMesh);
171 e5.stop();
172
173 // (6) Compute vertex distribution or local communication map
174 if (m2n().usesTwoLevelInitialization()) {
175
176 PRECICE_INFO("Compute communication map for mesh {}", _mesh->getName());
177 Event e6("partition.computeCommunicationMap." + _mesh->getName(), profiling::Synchronize);
178
179 // Fill two data structures: remoteCommunicationMap and this rank's communication map (_mesh->getCommunicationMap()).
180 // remoteCommunicationMap: connectedRank -> {remote local vertex index}
181 // _mesh->getCommunicationMap(): connectedRank -> {this rank's local vertex index}
182 // A vertex belongs to a specific connected rank if its global vertex ID lies within the ranks min and max.
183 mesh::Mesh::CommunicationMap remoteCommunicationMap;
184
185 for (size_t vertexIndex = 0; vertexIndex < _mesh->nVertices(); ++vertexIndex) {
186 for (size_t rankIndex = 0; rankIndex < _mesh->getConnectedRanks().size(); ++rankIndex) {
187 int globalVertexIndex = _mesh->vertex(vertexIndex).getGlobalIndex();
188 if (globalVertexIndex <= _remoteMaxGlobalVertexIDs[rankIndex] && globalVertexIndex >= _remoteMinGlobalVertexIDs[rankIndex]) {
189 int remoteRank = _mesh->getConnectedRanks()[rankIndex];
190 remoteCommunicationMap[remoteRank].push_back(globalVertexIndex - _remoteMinGlobalVertexIDs[rankIndex]); // remote local vertex index
191 _mesh->getCommunicationMap()[remoteRank].push_back(vertexIndex); // this rank's local vertex index
192 }
193 }
194 }
195
196 // communicate remote communication map to all remote connected ranks
197 m2n().scatterAllCommunicationMap(remoteCommunicationMap, *_mesh);
198
199 } else {
200
201 PRECICE_INFO("Feedback distribution for mesh {}", _mesh->getName());
202 Event e6("partition.feedbackMesh." + _mesh->getName(), profiling::Synchronize);
204 int numberOfVertices = _mesh->nVertices();
205 std::vector<VertexID> vertexIDs(numberOfVertices, -1);
206 for (int i = 0; i < numberOfVertices; i++) {
207 vertexIDs[i] = _mesh->vertex(i).getGlobalIndex();
208 }
209 PRECICE_DEBUG("Send partition feedback to primary rank");
210 utils::IntraComm::getCommunication()->sendRange(vertexIDs, 0);
211 } else { // Primary
212
213 mesh::Mesh::VertexDistribution vertexDistribution;
214 int numberOfVertices = _mesh->nVertices();
215 std::vector<VertexID> vertexIDs(numberOfVertices, -1);
216 for (int i = 0; i < numberOfVertices; i++) {
217 vertexIDs[i] = _mesh->vertex(i).getGlobalIndex();
218 }
219 vertexDistribution[0] = std::move(vertexIDs);
220
221 for (int secondaryRank : utils::IntraComm::allSecondaryRanks()) {
222 PRECICE_DEBUG("Receive partition feedback from the secondary rank {}", secondaryRank);
223 vertexDistribution[secondaryRank] = utils::IntraComm::getCommunication()->receiveRange(secondaryRank, com::AsVectorTag<VertexID>{});
224 }
225 PRECICE_ASSERT(_mesh->getVertexDistribution().empty());
226 _mesh->setVertexDistribution(std::move(vertexDistribution));
227 }
228 }
229
230 // (7) Compute vertex offsets
231 PRECICE_DEBUG("Compute vertex offsets");
233
234 // send number of vertices
235 PRECICE_DEBUG("Send number of vertices: {}", _mesh->nVertices());
236 int numberOfVertices = _mesh->nVertices();
237 utils::IntraComm::getCommunication()->send(numberOfVertices, 0);
238
239 // receive vertex offsets
240 mesh::Mesh::VertexOffsets vertexOffsets;
241 utils::IntraComm::getCommunication()->broadcast(vertexOffsets, 0);
242 PRECICE_DEBUG("My vertex offsets: {}", vertexOffsets);
243 PRECICE_ASSERT(_mesh->getVertexOffsets().empty());
244 _mesh->setVertexOffsets(std::move(vertexOffsets));
245
246 } else if (utils::IntraComm::isPrimary()) {
247
249 vertexOffsets[0] = _mesh->nVertices();
250
251 // receive number of secondary vertices and fill vertex offsets
252 for (int secondaryRank : utils::IntraComm::allSecondaryRanks()) {
253 int numberOfSecondaryRankVertices = -1;
254 utils::IntraComm::getCommunication()->receive(numberOfSecondaryRankVertices, secondaryRank);
255 PRECICE_ASSERT(numberOfSecondaryRankVertices >= 0);
256 vertexOffsets[secondaryRank] = numberOfSecondaryRankVertices + vertexOffsets[secondaryRank - 1];
257 }
258
259 // broadcast vertex offsets
260 PRECICE_DEBUG("My vertex offsets: {}", vertexOffsets);
261 utils::IntraComm::getCommunication()->broadcast(vertexOffsets);
262 PRECICE_ASSERT(_mesh->getVertexOffsets().empty());
263 _mesh->setVertexOffsets(std::move(vertexOffsets));
264 }
265}
266
267namespace {
268auto errorMeshFilteredOut(const std::string &meshName, const int rank)
269{
270 return fmt::format("The re-partitioning completely filtered out the mesh \"{0}\" received on rank {1} "
271 "at the coupling interface, although the provided mesh partition on this rank is "
272 "non-empty. Most probably, the coupling interfaces of your coupled participants do "
273 "not match geometry-wise. Please check your geometry setup again. Small overlaps or "
274 "gaps are no problem. If your geometry setup is correct and if you have very different "
275 "mesh resolutions on both sides, you may want to increase the safety-factor: "
276 "\"<receive-mesh mesh=\"{0} \" ... safety-factor=\"N\"/> (default value is 0.5) of the "
277 "decomposition strategy or disable the filtering completely: "
278 "\"<receive-mesh mesh=\"{0}\" ... geometric-filter=\"no-filter\" />",
279 meshName, rank);
280}
281} // namespace
282
284{
285 PRECICE_TRACE(static_cast<int>(_geometricFilter));
286
287 if (m2n().usesTwoLevelInitialization()) {
288 std::string msg = "The received mesh " + _mesh->getName() +
289 " cannot solely be filtered on the primary rank "
290 "(option \"filter-on-primary-rank\") if it is communicated by an m2n communication that uses "
291 "two-level initialization. Use \"filter-on-secondary-rank\" or \"no-filter\" instead.";
293 }
294
296
297 if (_geometricFilter == ON_PRIMARY_RANK) { // filter on primary rank and communicate reduced mesh then
298
299 PRECICE_ASSERT(not m2n().usesTwoLevelInitialization());
300 PRECICE_INFO("Pre-filter mesh {} by bounding box on primary rank", _mesh->getName());
301 Event e("partition.preFilterMesh." + _mesh->getName(), profiling::Synchronize);
302
304 PRECICE_DEBUG("Send bounding box to primary rank");
306 PRECICE_DEBUG("Receive filtered mesh");
308
310 PRECICE_CHECK(not _mesh->empty(), errorMeshFilteredOut(_mesh->getName(), utils::IntraComm::getRank()));
311 }
312
313 } else { // Primary
316
317 for (int secondaryRank : utils::IntraComm::allSecondaryRanks()) {
318 mesh::BoundingBox secondaryBB(_bb.getDimension());
319 com::receiveBoundingBox(*utils::IntraComm::getCommunication(), secondaryRank, secondaryBB);
320
321 PRECICE_DEBUG("From secondary rank {}, bounding mesh: {}", secondaryRank, secondaryBB);
322 mesh::Mesh secondaryMesh("SecondaryMesh", _dimensions, mesh::Mesh::MESH_ID_UNDEFINED);
323 mesh::filterMesh(secondaryMesh, *_mesh, [&secondaryBB](const mesh::Vertex &v) { return secondaryBB.contains(v); });
324 PRECICE_DEBUG("Send filtered mesh to secondary rank: {}", secondaryRank);
325 com::sendMesh(*utils::IntraComm::getCommunication(), secondaryRank, secondaryMesh);
326 }
327
328 // Now also filter the remaining primary mesh
329 mesh::Mesh filteredMesh("FilteredMesh", _dimensions, mesh::Mesh::MESH_ID_UNDEFINED);
330 mesh::filterMesh(filteredMesh, *_mesh, [&](const mesh::Vertex &v) { return _bb.contains(v); });
331 PRECICE_DEBUG("Primary rank mesh, filtered from {} to {} vertices, {} to {} edges, and {} to {} triangles.",
332 _mesh->nVertices(), filteredMesh.nVertices(),
333 _mesh->edges().size(), filteredMesh.edges().size(),
334 _mesh->triangles().size(), filteredMesh.triangles().size());
335 _mesh->clear();
336 _mesh->addMesh(filteredMesh);
337
339 PRECICE_CHECK(not _mesh->empty(), errorMeshFilteredOut(_mesh->getName(), utils::IntraComm::getRank()));
340 }
341 }
342 } else {
343 if (not m2n().usesTwoLevelInitialization()) {
344 PRECICE_INFO("Broadcast mesh {}", _mesh->getName());
345 Event e("partition.broadcastMesh." + _mesh->getName(), profiling::Synchronize);
346
349 } else { // Primary
352 }
353 }
355
356 PRECICE_INFO("Filter mesh {} by bounding box on secondary ranks", _mesh->getName());
357 Event e("partition.filterMeshBB." + _mesh->getName(), profiling::Synchronize);
358
359 mesh::Mesh filteredMesh("FilteredMesh", _dimensions, mesh::Mesh::MESH_ID_UNDEFINED);
360 mesh::filterMesh(filteredMesh, *_mesh, [&](const mesh::Vertex &v) { return _bb.contains(v); });
361
362 PRECICE_DEBUG("Bounding box filter, filtered from {} to {} vertices, {} to {} edges, and {} to {} triangles.",
363 _mesh->nVertices(), filteredMesh.nVertices(),
364 _mesh->edges().size(), filteredMesh.edges().size(),
365 _mesh->triangles().size(), filteredMesh.triangles().size());
366
367 _mesh->clear();
368 _mesh->addMesh(filteredMesh);
370 PRECICE_CHECK(not _mesh->empty(), errorMeshFilteredOut(_mesh->getName(), utils::IntraComm::getRank()));
371 }
372 } else {
374 }
375 }
376}
377
379{
381
382 _mesh->clear();
383 _mesh->clearPartitioning();
384 _boundingBoxPrepared = false;
387
388 // @todo handle coupling mode (i.e. serial participant)
389 // @todo treatment of multiple m2ns
390 PRECICE_ASSERT(_m2ns.size() == 1);
391 if (not m2n().usesTwoLevelInitialization())
392 return;
393
394 // receive and broadcast number of remote ranks
395 int numberOfRemoteRanks = -1;
397 m2n().getPrimaryRankCommunication()->receive(numberOfRemoteRanks, 0);
398 utils::IntraComm::getCommunication()->broadcast(numberOfRemoteRanks);
399 } else {
401 utils::IntraComm::getCommunication()->broadcast(numberOfRemoteRanks, 0);
402 }
403
404 // define and initialize remote bounding box map
405 mesh::Mesh::BoundingBoxMap remoteBBMap;
406 mesh::BoundingBox initialBB(_mesh->getDimensions());
407
408 for (int remoteRank = 0; remoteRank < numberOfRemoteRanks; remoteRank++) {
409 remoteBBMap.emplace(remoteRank, initialBB);
410 }
411
412 // receive and broadcast remote bounding box map
414 com::receiveBoundingBoxMap(*m2n().getPrimaryRankCommunication(), 0, remoteBBMap);
416 } else {
419 }
420
421 // prepare local bounding box
423
424 if (utils::IntraComm::isPrimary()) { // Primary
425 mesh::Mesh::CommunicationMap connectionMap; // local ranks -> {remote ranks}
426 std::vector<Rank> connectedRanksList; // local ranks with any connection
427
428 // connected ranks for primary rank
429 std::vector<Rank> connectedRanks;
430 for (auto &remoteBB : remoteBBMap) {
431 if (_bb.overlapping(remoteBB.second)) {
432 connectedRanks.push_back(remoteBB.first); // connected remote ranks for this rank
433 }
434 }
435 PRECICE_ASSERT(_mesh->getConnectedRanks().empty());
436 _mesh->setConnectedRanks(connectedRanks);
437 if (not connectedRanks.empty()) {
438 connectionMap.emplace(0, std::move(connectedRanks));
439 connectedRanksList.push_back(0);
440 }
441
442 // receive connected ranks from secondary ranks and add them to the connection map
443 for (int rank : utils::IntraComm::allSecondaryRanks()) {
444 std::vector<Rank> secondaryConnectedRanks = utils::IntraComm::getCommunication()->receiveRange(rank, com::AsVectorTag<Rank>{});
445 if (!secondaryConnectedRanks.empty()) {
446 connectedRanksList.push_back(rank);
447 connectionMap.emplace(rank, std::move(secondaryConnectedRanks));
448 }
449 }
450
451 // send connectionMap to other primary rank
452 m2n().getPrimaryRankCommunication()->sendRange(connectedRanksList, 0);
453 PRECICE_CHECK(not connectionMap.empty(),
454 "The mesh \"{}\" of this participant seems to have no partitions at the coupling interface. "
455 "Check that both mapped meshes are describing the same geometry. "
456 "If you deal with very different mesh resolutions, consider increasing the safety-factor in the <receive-mesh /> tag.",
457 _mesh->getName());
458 com::sendConnectionMap(*m2n().getPrimaryRankCommunication(), 0, connectionMap);
459 } else {
461
462 std::vector<Rank> connectedRanks;
463 for (const auto &remoteBB : remoteBBMap) {
464 if (_bb.overlapping(remoteBB.second)) {
465 connectedRanks.push_back(remoteBB.first);
466 }
467 }
468 PRECICE_ASSERT(_mesh->getConnectedRanks().empty());
469 _mesh->setConnectedRanks(connectedRanks);
470
471 // send connected ranks to primary rank
472 utils::IntraComm::getCommunication()->sendRange(connectedRanks, 0);
473 }
474}
475
477{
479
481 return;
482
483 PRECICE_DEBUG("Merge bounding boxes and increase by safety factor");
484
485 // Reset the BoundingBox
487
488 // Create BB around all "other" meshes
489 for (mapping::PtrMapping &fromMapping : _fromMappings) {
490 auto other_bb = fromMapping->getOutputMesh()->getBoundingBox();
491 _bb.expandBy(other_bb);
494 }
495 for (mapping::PtrMapping &toMapping : _toMappings) {
496 auto other_bb = toMapping->getInputMesh()->getBoundingBox();
497 _bb.expandBy(other_bb);
500 }
501
502 // Expand by user-defined bounding box in case a direct access is desired
503 if (_allowDirectAccess) {
504 auto &other_bb = _mesh->getBoundingBox();
505 _bb.expandBy(other_bb);
506
507 // The safety factor is for mapping based partitionings applied, as usual.
508 // For the direct access, however, we don't apply any safety factor scaling.
509 // If the user defines a safety factor and the partitioning is entirely based
510 // on the defined access region (setMeshAccessRegion), we raise a warning
511 // to inform the user
512 const float defaultSafetyFactor = 0.5;
514 utils::IntraComm::isPrimary() && !hasAnyMapping() && (_safetyFactor != defaultSafetyFactor),
515 "The received mesh \"{}\" was entirely partitioned based on the defined access region "
516 "(setMeshAccessRegion) and a safety-factor was defined. However, the safety factor "
517 "will be ignored in this case. You may want to modify the access region by modifying "
518 "the specified region in the function itself.",
519 _mesh->getName());
521 }
522}
523
525{
527 Event e("partition.createOwnerInformation." + _mesh->getName(), profiling::Synchronize);
528
529 /*
530 We follow different approaches for two-level and one-level methods. For 1LI, a centeralized
531 approach is followed, while the 2LI employs a local/parallel scheme for assigning vertices
532 to corresponding ranks.
533 */
534
535 if (m2n().usesTwoLevelInitialization()) {
536 /*
537 This function ensures that each vertex is owned by only a single rank and
538 is not shared among ranks. Initially, the vertices are checked against the
539 bounding box of each rank. If a vertex fits into only a single bounding box,
540 the vertex is assigned to that rank. If it fits to various bbs, the rank with the
541 lowest number of vertices gets ownership to keep the load as balanced as
542 possible.
543
544 Following steps are taken:
545
546 1- receive local bb map from primary rank
547 2- filter bb map to keep the connected ranks
548 3- own the vertices that only fit into this rank's bb
549 4- send number of owned vertices and the list of shared vertices to neighbors
550 5- for the remaining vertices: check if we have less vertices -> own it!
551 */
552
553 // #1: receive local bb map from primary rank
554 // Define and initialize localBBMap to save local bbs
555
557 for (Rank rank = 0; rank < utils::IntraComm::getSize(); rank++) {
558 localBBMap.emplace(rank, mesh::BoundingBox(_dimensions));
559 }
560
561 // Define a bb map to save the local connected ranks and respective boundingboxes
562 mesh::Mesh::BoundingBoxMap localConnectedBBMap;
563
564 // store global IDs and list of possible shared vertices
565
566 // Global IDs to be saved in:
567 std::vector<VertexID> sharedVerticesGlobalIDs;
568 // Local IDs to be saved in:
569 std::vector<VertexID> sharedVerticesLocalIDs;
570
571 // store possible shared vertices in a map to communicate with neighbors, map: rank -> vertex_global_id
572 mesh::Mesh::CommunicationMap sharedVerticesSendMap;
573
574 // receive list of possible shared vertices from neighboring ranks
575 mesh::Mesh::CommunicationMap sharedVerticesReceiveMap;
576
578
579 // Insert bounding box of primary ranks
580 localBBMap.at(0) = _bb;
581
582 // primary rank receives local bb from each secondary rank
583 for (int secondaryRank = 1; secondaryRank < utils::IntraComm::getSize(); secondaryRank++) {
584 com::receiveBoundingBox(*utils::IntraComm::getCommunication(), secondaryRank, localBBMap.at(secondaryRank));
585 }
586
587 // primary rank broadcast localBBMap to all secondary ranks
589 } else if (utils::IntraComm::isSecondary()) {
590 // secondary ranks send local bb to primary rank
592 // secondary ranks receive localBBMap from primary rank
594 }
595
596 // #2: filter bb map to keep the connected ranks
597 // remove the own bb from the map since we compare the own bb only with other ranks bb.
598 localBBMap.erase(utils::IntraComm::getRank());
599 // find and store local connected ranks
600 for (const auto &localBB : localBBMap) {
601 if (_bb.overlapping(localBB.second)) {
602 localConnectedBBMap.emplace(localBB.first, localBB.second);
603 }
604 }
605
606 // First, insert all comm partnerts, otherwise we end up in a deadlock further down
607 // when exchanging the vector sizes as vertices might be shared from the other
608 // connected rank(s), although the actual size we request here is zero
609 // i.e., we never tell the other ranks that we don't want any vertices
610 // See also test Integration/Parallel/TestBoundingBoxInitializationEmpty
611 for (auto &neighborRank : localConnectedBBMap)
612 sharedVerticesSendMap[neighborRank.first] = std::vector<VertexID>();
613
614 // #3: check vertices and keep only those that fit into the current rank's bb
615 const int numberOfVertices = _mesh->nVertices();
616 PRECICE_DEBUG("Tag vertices, number of vertices {}", numberOfVertices);
617 std::vector<int> tags(numberOfVertices, 1);
618 std::vector<VertexID> globalIDs(numberOfVertices, -1);
619 int ownedVerticesCount = 0; // number of vertices owned by this rank
620 for (int i = 0; i < numberOfVertices; i++) {
621 globalIDs[i] = _mesh->vertex(i).getGlobalIndex();
622 if (_mesh->vertex(i).isTagged()) {
623 bool vertexIsShared = false;
624 for (const auto &neighborRank : localConnectedBBMap) {
625 if (neighborRank.second.contains(_mesh->vertex(i))) {
626 vertexIsShared = true;
627 sharedVerticesSendMap[neighborRank.first].push_back(globalIDs[i]);
628 }
629 }
630
631 if (not vertexIsShared) {
632 tags[i] = 1;
633 ownedVerticesCount++;
634 } else {
635 sharedVerticesGlobalIDs.push_back(globalIDs[i]);
636 sharedVerticesLocalIDs.push_back(i);
637 }
638 } else {
639 tags[i] = 0;
640 }
641 }
642
643 // #4: Exchange number of already owned vertices with the neighbors
644
645 // to store receive requests.
646 std::vector<com::PtrRequest> vertexNumberRequests;
647
648 // Define and initialize for load balancing
649 std::map<int, int> neighborRanksVertexCount;
650 for (auto &neighborRank : localConnectedBBMap) {
651 neighborRanksVertexCount.emplace(neighborRank.first, 0);
652 }
653
654 // Asynchronous receive number of owned vertices from neighbor ranks
655 for (auto &neighborRank : localConnectedBBMap) {
656 auto request = utils::IntraComm::getCommunication()->aReceive(neighborRanksVertexCount.at(neighborRank.first), neighborRank.first);
657 vertexNumberRequests.push_back(request);
658 }
659
660 // Synchronous send number of owned vertices to neighbor ranks
661 for (auto &neighborRank : localConnectedBBMap) {
662 utils::IntraComm::getCommunication()->send(ownedVerticesCount, neighborRank.first);
663 }
664
665 // wait until all aReceives are complete.
666 for (auto &rqst : vertexNumberRequests) {
667 rqst->wait();
668 }
669
670 // Exchange list of shared vertices with the neighbor
671 // to store send requests.
672 std::vector<com::PtrRequest> vertexListRequests;
673
674 for (auto &receivingRank : sharedVerticesSendMap) {
675 int sendSize = receivingRank.second.size();
676 auto request = utils::IntraComm::getCommunication()->aSend(sendSize, receivingRank.first);
677 vertexListRequests.push_back(request);
678 if (sendSize != 0) {
679 auto request = utils::IntraComm::getCommunication()->aSend(span<const int>{receivingRank.second}, receivingRank.first);
680 vertexListRequests.push_back(request);
681 }
682 }
683
684 for (auto &neighborRank : sharedVerticesSendMap) {
685 int receiveSize = 0;
686 utils::IntraComm::getCommunication()->receive(receiveSize, neighborRank.first);
687 if (receiveSize != 0) {
688 std::vector<int> receivedSharedVertices(receiveSize, -1);
689 utils::IntraComm::getCommunication()->receive(span<int>{receivedSharedVertices}, neighborRank.first);
690 sharedVerticesReceiveMap.insert(std::make_pair(neighborRank.first, receivedSharedVertices));
691 }
692 }
693
694 // wait until all aSends are complete.
695 for (auto &rqst : vertexListRequests) {
696 rqst->wait();
697 }
698
699 // #5: Second round assignment according to the number of owned vertices
700 // In case that a vertex can be shared between two ranks, the rank with lower
701 // vertex count will own the vertex.
702 // If both ranks have same vertex count, the lower rank will own the vertex.
703
704 // To do so, we look at all vertices shared with all neighbors
705 for (auto &sharingRank : sharedVerticesReceiveMap) {
706 // First, check if we would change the ownership at all (by default initialization
707 // above, we would be considered as owner of the shared vertices)
708 // If this is fulfilled, we need to set the ownership to false
709 if ((ownedVerticesCount > neighborRanksVertexCount[sharingRank.first]) ||
710 (ownedVerticesCount == neighborRanksVertexCount[sharingRank.first] && utils::IntraComm::getRank() > sharingRank.first)) {
711 // In such a case, we need to find out the tags we need to switch to 'false', as we don't want
712 // to own them any longer. We compute the intersection of all globalIDs this rank shares with
713 // others and the globalIDs of the neighbor rank we are just considering.
714 // The set_intersection_indices gives us the indices in the sharedVerticesGlobalIDs, which we
715 // can use in the sharedVerticesLocalIDs to get the actual index in the 'tags' vector.
717 precice::utils::set_intersection_indices(sharedVerticesGlobalIDs.begin(), sharedVerticesGlobalIDs.begin(), sharedVerticesGlobalIDs.end(),
718 sharingRank.second.begin(), sharingRank.second.end(), std::back_inserter(res));
719 for (auto r : res)
720 tags[sharedVerticesLocalIDs[r]] = static_cast<int>(false);
721 }
722 }
723
725 PRECICE_DEBUG("{} of {} vertices of mesh {} have been filtered out on rank {} since they have no influence on the mapping.",
726 std::count(tags.begin(), tags.end(), 0), tags.size(), _mesh->getName(), utils::IntraComm::getRank());
727 // end of two-level initialization section
728 } else {
730 int numberOfVertices = _mesh->nVertices();
731 utils::IntraComm::getCommunication()->send(numberOfVertices, 0);
732
733 if (numberOfVertices != 0) {
734 PRECICE_DEBUG("Tag vertices, number of vertices {}", numberOfVertices);
735 std::vector<int> tags(numberOfVertices, -1);
736 std::vector<VertexID> globalIDs(numberOfVertices, -1);
737 bool atInterface = false;
738 for (int i = 0; i < numberOfVertices; i++) {
739 globalIDs[i] = _mesh->vertex(i).getGlobalIndex();
740 if (_mesh->vertex(i).isTagged()) {
741 tags[i] = 1;
742 atInterface = true;
743 } else {
744 tags[i] = 0;
745 }
746 }
747 PRECICE_DEBUG("My tags: {}", tags);
748 PRECICE_DEBUG("My global IDs: {}", globalIDs);
749 PRECICE_DEBUG("Send tags and global IDs");
750 utils::IntraComm::getCommunication()->sendRange(tags, 0);
751 utils::IntraComm::getCommunication()->sendRange(globalIDs, 0);
752 utils::IntraComm::getCommunication()->send(atInterface, 0);
753
754 PRECICE_DEBUG("Receive owner information");
756 PRECICE_DEBUG("My owner information: {}", ownerVec);
757 PRECICE_ASSERT(ownerVec.size() == static_cast<std::size_t>(numberOfVertices));
758 setOwnerInformation(ownerVec);
759 }
760 }
761
762 else if (utils::IntraComm::isPrimary()) {
763 // To temporary store which vertices already have an owner
764 std::vector<VertexID> globalOwnerVec(_mesh->getGlobalNumberOfVertices(), 0);
765 // The same per rank
767 // Global IDs per rank
769 // Tag information per rank
771
772 // Fill primary data
773 PRECICE_DEBUG("Tag vertices of primary rank");
774 bool primaryRankAtInterface = false;
775 secondaryOwnerVecs[0].resize(_mesh->nVertices());
776 secondaryGlobalIDs[0].resize(_mesh->nVertices());
777 secondaryTags[0].resize(_mesh->nVertices());
778 for (size_t i = 0; i < _mesh->nVertices(); i++) {
779 secondaryGlobalIDs[0][i] = _mesh->vertex(i).getGlobalIndex();
780 if (_mesh->vertex(i).isTagged()) {
781 primaryRankAtInterface = true;
782 secondaryTags[0][i] = 1;
783 } else {
784 secondaryTags[0][i] = 0;
785 }
786 }
787 PRECICE_DEBUG("My tags: {}", secondaryTags[0]);
788
789 // receive secondary data
790 Rank ranksAtInterface = 0;
791 if (primaryRankAtInterface)
792 ranksAtInterface++;
793
795 int localNumberOfVertices = -1;
796 utils::IntraComm::getCommunication()->receive(localNumberOfVertices, rank);
797 PRECICE_DEBUG("Rank {} has {} vertices.", rank, localNumberOfVertices);
798 secondaryOwnerVecs[rank].resize(localNumberOfVertices, 0);
799
800 if (localNumberOfVertices != 0) {
801 PRECICE_DEBUG("Receive tags from secondary rank {}", rank);
802 secondaryTags[rank] = utils::IntraComm::getCommunication()->receiveRange(rank, com::AsVectorTag<int>{});
803 secondaryGlobalIDs[rank] = utils::IntraComm::getCommunication()->receiveRange(rank, com::AsVectorTag<VertexID>{});
804 PRECICE_DEBUG("Rank {} has tags {}", rank, secondaryTags[rank]);
805 PRECICE_DEBUG("Rank {} has global IDs {}", rank, secondaryGlobalIDs[rank]);
806 bool atInterface = false;
807 utils::IntraComm::getCommunication()->receive(atInterface, rank);
808 if (atInterface)
809 ranksAtInterface++;
810 }
811 }
812
813 // Decide upon owners,
814 PRECICE_DEBUG("Decide owners, first round by rough load balancing");
815 // Provide a more descriptive error message if direct access was enabled
816 PRECICE_CHECK(!(ranksAtInterface == 0 && _allowDirectAccess),
817 "After repartitioning of mesh \"{}\" all ranks are empty. "
818 "Please check the dimensions of the provided bounding box "
819 "(in \"setMeshAccessRegion\") and verify that it covers vertices "
820 "in the mesh or check the definition of the provided meshes.",
821 _mesh->getName());
822 PRECICE_ASSERT(ranksAtInterface != 0);
823 int localGuess = _mesh->getGlobalNumberOfVertices() / ranksAtInterface; // Guess for a decent load balancing
824 // First round: every secondary rank gets localGuess vertices
825 for (Rank rank : utils::IntraComm::allRanks()) {
826 int counter = 0;
827 for (size_t i = 0; i < secondaryOwnerVecs[rank].size(); i++) {
828 // Vertex has no owner yet and rank could be owner
829 if (globalOwnerVec[secondaryGlobalIDs[rank][i]] == 0 && secondaryTags[rank][i] == 1) {
830 secondaryOwnerVecs[rank][i] = 1; // Now rank is owner
831 globalOwnerVec[secondaryGlobalIDs[rank][i]] = 1; // Vertex now has owner
832 counter++;
833 if (counter == localGuess)
834 break;
835 }
836 }
837 }
838
839 // Second round: distribute all other vertices in a greedy way
840 PRECICE_DEBUG("Decide owners, second round in greedy way");
841 for (Rank rank : utils::IntraComm::allRanks()) {
842 for (size_t i = 0; i < secondaryOwnerVecs[rank].size(); i++) {
843 if (globalOwnerVec[secondaryGlobalIDs[rank][i]] == 0 && secondaryTags[rank][i] == 1) {
844 secondaryOwnerVecs[rank][i] = 1;
845 globalOwnerVec[secondaryGlobalIDs[rank][i]] = rank + 1;
846 }
847 }
848 }
849
850 // Send information back to secondary ranks
852 if (not secondaryTags[rank].empty()) {
853 PRECICE_DEBUG("Send owner information to secondary rank {}", rank);
854 utils::IntraComm::getCommunication()->sendRange(secondaryOwnerVecs[rank], rank);
855 }
856 }
857 // Primary rank data
858 PRECICE_DEBUG("My owner information: {}", secondaryOwnerVecs[0]);
859 setOwnerInformation(secondaryOwnerVecs[0]);
860
861#ifndef NDEBUG
862 for (size_t i = 0; i < globalOwnerVec.size(); i++) {
863 PRECICE_DEBUG_IF(globalOwnerVec[i] == 0,
864 "The Vertex with global index {} of mesh: {} was completely filtered out, since it has no influence on any mapping.",
865 i, _mesh->getName());
866 }
867#endif
868 auto filteredVertices = std::count(globalOwnerVec.begin(), globalOwnerVec.end(), 0);
869 PRECICE_WARN_IF(filteredVertices,
870 "{} of {} vertices of mesh {} have been filtered out since they have no influence on the mapping.{}",
871 filteredVertices, _mesh->getGlobalNumberOfVertices(), _mesh->getName(),
872 _allowDirectAccess ? " Associated data values of the filtered vertices will be filled with zero values in order to "
873 "provide valid data for other participants when reading data."
874 : "");
875 }
876 }
877}
878
880{
881 for (const auto &fromMapping : _fromMappings) {
882 if (not fromMapping->getOutputMesh()->empty()) {
883 return true;
884 }
885 }
886 for (const auto &toMapping : _toMappings) {
887 if (not toMapping->getInputMesh()->empty()) {
888 return true;
889 }
890 }
891 return false;
892}
893
895{
896 return not(_fromMappings.empty() && _toMappings.empty());
897}
898
900{
901 // We want to have every vertex within the box if we access the mesh directly
902 if (_allowDirectAccess) {
903 _mesh->tagAll();
904 return;
905 }
906
907 for (const mapping::PtrMapping &fromMapping : _fromMappings) {
908 fromMapping->tagMeshFirstRound();
909 }
910 for (const mapping::PtrMapping &toMapping : _toMappings) {
911 toMapping->tagMeshFirstRound();
912 }
913}
914
916{
917 // We have already tagged every node in this case in the first round
918 if (_allowDirectAccess) {
919 return;
920 }
921
922 for (const mapping::PtrMapping &fromMapping : _fromMappings) {
923 fromMapping->tagMeshSecondRound();
924 }
925 for (const mapping::PtrMapping &toMapping : _toMappings) {
926 toMapping->tagMeshSecondRound();
927 }
928}
929
931{
932 size_t i = 0;
933 for (mesh::Vertex &vertex : _mesh->vertices()) {
934 PRECICE_ASSERT(i < ownerVec.size());
935 PRECICE_ASSERT(ownerVec[i] != -1);
936 vertex.setOwner(ownerVec[i] == 1);
937 i++;
938 }
939}
940
942{
943 PRECICE_ASSERT(_m2ns.size() == 1);
944 return *_m2ns[0];
945}
946
947} // namespace precice::partition
#define PRECICE_WARN_IF(condition,...)
Definition LogMacros.hpp:21
#define PRECICE_DEBUG(...)
Definition LogMacros.hpp:64
#define PRECICE_DEBUG_IF(condition,...)
Definition LogMacros.hpp:66
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_INFO(...)
Definition LogMacros.hpp:13
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
T at(T... args)
T back_inserter(T... args)
T begin(T... args)
M2N communication class. This layer is necessary since communication between two participants can be ...
Definition M2N.hpp:31
void scatterAllCommunicationMap(std::map< int, std::vector< int > > &localCommunicationMap, mesh::Mesh &mesh)
Scatters a communication map over connected ranks on remote participant (concerning the given mesh)
Definition M2N.cpp:277
void broadcastReceiveAll(std::vector< int > &itemToReceive, mesh::Mesh &mesh)
Receives an int per connected rank on remote participant (concerning the given mesh) @para[out] itemT...
Definition M2N.cpp:370
com::PtrCommunication getPrimaryRankCommunication()
Get the basic communication between the 2 primary ranks.
Definition M2N.cpp:195
void broadcastReceiveAllMesh(mesh::Mesh &mesh)
Receive mesh partitions per connected rank on remote participant (concerning the given mesh)
Definition M2N.cpp:379
An axis-aligned bounding box around a (partition of a) mesh.
bool contains(const Vertex &vertex) const
Checks if vertex in contained in _bb.
bool empty() const
Check if every dimension's length is equal to zero.
void expandBy(const BoundingBox &otherBB)
Expand bounding box using another bounding box.
bool overlapping(const BoundingBox &otherBB) const
Checks whether two bounding boxes are overlapping.
int getDimension() const
Getter dimension of the bounding box.
void scaleBy(double safetyFactor)
Increase the size of bounding box by safety margin.
Container and creator for meshes.
Definition Mesh.hpp:39
static constexpr MeshID MESH_ID_UNDEFINED
Use if the id of the mesh is not necessary.
Definition Mesh.hpp:58
std::size_t nVertices() const
Returns the number of vertices.
Definition Mesh.cpp:63
TriangleContainer & triangles()
Returns modifiable container holding all triangles.
Definition Mesh.cpp:78
EdgeContainer & edges()
Returns modifiable container holding all edges.
Definition Mesh.cpp:68
Vertex of a mesh.
Definition Vertex.hpp:16
bool isTagged() const
Definition Vertex.cpp:32
Abstract base class for partitions.
Definition Partition.hpp:29
std::vector< mapping::PtrMapping > _toMappings
Definition Partition.hpp:67
std::vector< mapping::PtrMapping > _fromMappings
Definition Partition.hpp:65
std::vector< m2n::PtrM2N > _m2ns
m2n connection to each connected participant
Definition Partition.hpp:70
void prepareBoundingBox()
Sets _bb to the union with the mesh from fromMapping resp. toMapping, also enlage by _safetyFactor.
ReceivedPartition(const mesh::PtrMesh &mesh, GeometricFilter geometricFilter, double safetyFactor, bool allowDirectAccess=false)
Constructor.
void compute() override
The partition is computed, i.e. the mesh re-partitioned if required and all data structures are set u...
std::vector< int > _remoteMaxGlobalVertexIDs
Max global vertex IDs of remote connected ranks.
void tagMeshFirstRound()
Tag mesh in first round according to all mappings.
std::vector< int > _remoteMinGlobalVertexIDs
Min global vertex IDs of remote connected ranks.
bool hasAnyMapping() const
Returns whether any mapping is defined.
void tagMeshSecondRound()
Tag mesh in second round according to all mappings.
m2n::M2N & m2n()
return the one m2n, a ReceivedPartition can only have one m2n
void compareBoundingBoxes() override
Intersections between bounding boxes around each rank are computed.
bool _boundingBoxPrepared
Is the local other (i.e. provided) bounding box already prepared (i.e. has prepareBoundingBox() been ...
GeometricFilter
Defines the type of geometric filter used.
@ ON_PRIMARY_RANK
Filter at primary rank and communicate only filtered mesh.
@ ON_SECONDARY_RANKS
Filter after communication on all secondary ranks.
@ NO_FILTER
No geometric filter used (e.g. for RBF mappings)
void setOwnerInformation(const std::vector< int > &ownerVec)
Helper function for 'createOwnerFunction' to set local owner information.
void communicate() override
The mesh is communicated between both primary ranks (if required)
void stop()
Stops a running event.
Definition Event.cpp:37
A C++ 11 implementation of the non-owning C++20 std::span type.
Definition span.hpp:284
PRECICE_SPAN_CONSTEXPR11 span< element_type, Count > first() const
Definition span.hpp:416
static int getSize()
Number of ranks. This includes ranks from both participants, e.g. minimal size is 2.
Definition IntraComm.cpp:47
static Rank getRank()
Current rank.
Definition IntraComm.cpp:42
static bool isPrimary()
True if this process is running the primary rank.
Definition IntraComm.cpp:52
static auto allSecondaryRanks()
Returns an iterable range over salve ranks [1, _size)
Definition IntraComm.hpp:37
static auto allRanks()
Returns an iterable range over all ranks [0, _size)
Definition IntraComm.hpp:43
static bool isParallel()
True if this process is running in parallel.
Definition IntraComm.cpp:62
static bool isSecondary()
True if this process is running a secondary rank.
Definition IntraComm.cpp:57
static com::PtrCommunication & getCommunication()
Intra-participant communication.
Definition IntraComm.hpp:31
T clear(T... args)
T count(T... args)
T emplace(T... args)
T empty(T... args)
T end(T... args)
T erase(T... args)
T insert(T... args)
T make_pair(T... args)
void receiveBoundingBoxMap(Communication &communication, int rankSender, mesh::Mesh::BoundingBoxMap &bbm)
Definition Extra.cpp:63
void broadcastReceiveBoundingBoxMap(Communication &communication, mesh::Mesh::BoundingBoxMap &bbm)
Definition Extra.cpp:73
void sendMesh(Communication &communication, int rankReceiver, const mesh::Mesh &mesh)
Definition Extra.cpp:8
void broadcastSendMesh(Communication &communication, const mesh::Mesh &mesh)
Definition Extra.cpp:18
void sendBoundingBox(Communication &communication, int rankReceiver, const mesh::BoundingBox &bb)
Definition Extra.cpp:48
void sendConnectionMap(Communication &communication, int rankReceiver, const mesh::Mesh::ConnectionMap &cm)
Definition Extra.cpp:28
void broadcastReceiveMesh(Communication &communication, mesh::Mesh &mesh)
Definition Extra.cpp:23
void receiveMesh(Communication &communication, int rankSender, mesh::Mesh &mesh)
Definition Extra.cpp:13
void broadcastSendBoundingBoxMap(Communication &communication, const mesh::Mesh::BoundingBoxMap &bbm)
Definition Extra.cpp:68
void receiveBoundingBox(Communication &communication, int rankSender, mesh::BoundingBox &bb)
Definition Extra.cpp:53
void filterMesh(Mesh &destination, const Mesh &source, UnaryPredicate p)
Definition Filter.hpp:17
contains the partitioning of distributed meshes.
Definition Partition.cpp:5
static constexpr SynchronizeTag Synchronize
Convenience instance of the SynchronizeTag.
Definition Event.hpp:21
void set_intersection_indices(InputIt1 ref1, InputIt1 first1, InputIt1 last1, InputIt2 first2, InputIt2 last2, OutputIt d_first)
This function is by and large the same as std::set_intersection(). The only difference is that we don...
Definition algorithm.hpp:34
int Rank
Definition Types.hpp:37
T push_back(T... args)
T resize(T... args)
T size(T... args)