Loading...
Searching...
No Matches
voloctree_mapper_example_00003.cpp
1/*---------------------------------------------------------------------------*\
2 *
3 * bitpit
4 *
5 * Copyright (C) 2015-2021 OPTIMAD engineering Srl
6 *
7 * -------------------------------------------------------------------------
8 * License
9 * This file is part of bitpit.
10 *
11 * bitpit is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU Lesser General Public License v3 (LGPL)
13 * as published by the Free Software Foundation.
14 *
15 * bitpit is distributed in the hope that it will be useful, but WITHOUT
16 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 * License for more details.
19 *
20 * You should have received a copy of the GNU Lesser General Public License
21 * along with bitpit. If not, see <http://www.gnu.org/licenses/>.
22 *
23\*---------------------------------------------------------------------------*/
24
32#include <array>
33#if BITPIT_ENABLE_MPI==1
34#include <mpi.h>
35#endif
36
37#include "bitpit_POD.hpp"
38#include "bitpit_patchkernel.hpp"
39#include "bitpit_voloctree.hpp"
40
41using namespace bitpit;
42using namespace pod;
43
47void run()
48{
49
50 log::cout() << " >> 2D octree patch" << "\n";
51
55 double x_0 = 10.;
56 double y_0 = 20.;
57 double z_0 = 30.;
58 double l = 1.5;
59
60 std::unique_ptr<PabloUniform> treePointer = std::unique_ptr<PabloUniform>(new PabloUniform(x_0, y_0, z_0, l, 2));
61 PabloUniform &octree = *treePointer;
62
63 std::cout << " Origin : ( " << octree.getX0() << ", " << octree.getY0() << ", " << octree.getZ0() << " )" << std::endl;
64 std::cout << " Length : " << octree.getL() << std::endl;
65
67 octree.adaptGlobalRefine();
68 octree.adaptGlobalRefine();
69 octree.adaptGlobalRefine();
70 octree.adaptGlobalRefine();
71
76 VolOctree *patch_2D_original = new VolOctree(std::move(treePointer), &treePointer);
77
78 patch_2D_original->initializeAdjacencies();
79 patch_2D_original->initializeInterfaces();
80
81 patch_2D_original->update();
82
83#if BITPIT_ENABLE_MPI==1
85 patch_2D_original->partition(true);
86#endif
87
88 std::vector<uint64_t> refineList;
89 refineList.push_back( 7);
90 refineList.push_back( 13);
91 refineList.push_back( 15);
92 refineList.push_back( 26);
93 refineList.push_back( 27);
94 refineList.push_back( 31);
95 refineList.push_back( 37);
96 refineList.push_back( 39);
97 refineList.push_back( 49);
98 refineList.push_back( 50);
99 refineList.push_back( 51);
100 refineList.push_back( 53);
101 refineList.push_back( 55);
102 refineList.push_back( 63);
103 refineList.push_back( 78);
104 refineList.push_back(100);
105 refineList.push_back(102);
106 refineList.push_back(105);
107 refineList.push_back(108);
108 refineList.push_back(109);
109 refineList.push_back(110);
110 refineList.push_back(135);
111 refineList.push_back(141);
112 refineList.push_back(143);
113 refineList.push_back(146);
114 refineList.push_back(147);
115 refineList.push_back(151);
116 refineList.push_back(153);
117 refineList.push_back(154);
118 refineList.push_back(155);
119 refineList.push_back(157);
120 refineList.push_back(159);
121 refineList.push_back(165);
122 refineList.push_back(167);
123 refineList.push_back(183);
124 refineList.push_back(198);
125 refineList.push_back(204);
126 refineList.push_back(206);
127 refineList.push_back(225);
128 refineList.push_back(228);
129 refineList.push_back(229);
130 refineList.push_back(230);
131
132
133 for (uint64_t ind : refineList) {
134#if BITPIT_ENABLE_MPI==1
135 int owner = patch_2D_original->getTree().getOwnerRank(ind);
136 if (patch_2D_original->getRank() == owner){
137 uint32_t lind = patch_2D_original->getTree().getLocalIdx(ind, owner);
138#else
139 uint32_t lind = patch_2D_original->getTree().getLocalIdx(ind);
140#endif
141 VolOctree::OctantInfo octinfo(lind, true);
142 long id = patch_2D_original->getOctantId(octinfo);
143 patch_2D_original->markCellForRefinement(id);
144#if BITPIT_ENABLE_MPI==1
145 }
146#endif
147 }
148
149 patch_2D_original->update(false);
150
151 patch_2D_original->getVTK().setName("mesh_original.0");
152#if BITPIT_ENABLE_MPI==1
153 patch_2D_original->setVTKWriteTarget(PatchKernel::WriteTarget::WRITE_TARGET_CELLS_INTERNAL);
154#endif
155 patch_2D_original->write();
156
160 std::unique_ptr<PabloUniform> treePointer2 = std::unique_ptr<PabloUniform>(new PabloUniform(x_0, y_0, z_0, l, 2));
161 PabloUniform &octree2 = *treePointer2;
162
164 octree2.adaptGlobalRefine();
165 octree2.adaptGlobalRefine();
166 octree2.adaptGlobalRefine();
167 octree2.adaptGlobalRefine();
168
170 VolOctree *patch_2D = new VolOctree(std::move(treePointer2), &treePointer2);
171 patch_2D->update(false);
172#if BITPIT_ENABLE_MPI==1
173 patch_2D->setVTKWriteTarget(PatchKernel::WriteTarget::WRITE_TARGET_CELLS_INTERNAL);
174#endif
175
177 for (int k = 0; k < 10; ++k) {
178 long nCells = patch_2D->getCellCount();
179 log::cout() << std::endl;
180 log::cout() << ">> Marking the cells to adapt... " << std::endl;
181
182 for (int i = 0; i < 100; ++i) {
183 long cellId;
184#if BITPIT_ENABLE_MPI==1
185 if (patch_2D->getRank() == 0)
186#endif
187 cellId = rand() % nCells * 2;
188#if BITPIT_ENABLE_MPI==1
189 MPI_Bcast(&cellId, 1, MPI_LONG, 0, MPI_COMM_WORLD);
190#endif
191
192 if (!patch_2D->getCells().exists(cellId)) {
193 continue;
194 }
195
196 patch_2D->markCellForRefinement(cellId);
197 }
198
199 for (int i = 0; i < 55; ++i) {
200 long cellId;
201#if BITPIT_ENABLE_MPI==1
202 if (patch_2D->getRank() == 0)
203#endif
204 cellId = rand() % nCells * 2;
205#if BITPIT_ENABLE_MPI==1
206 MPI_Bcast(&cellId, 1, MPI_LONG, 0, MPI_COMM_WORLD);
207#endif
208
209 if (!patch_2D->getCells().exists(cellId)) {
210 continue;
211 }
212
213 patch_2D->markCellForCoarsening(cellId);
214 for (auto neighId : patch_2D->findCellNeighs(cellId)) {
215 patch_2D->markCellForCoarsening(neighId);
216 }
217 }
218
219 log::cout() << std::endl;
220 log::cout() << ">> Initial number of cells... " << nCells << std::endl;
221
222 patch_2D->initializeAdjacencies();
223 patch_2D->initializeInterfaces();
224
225 patch_2D->update(false);
226
227 nCells = patch_2D->getCellCount();
228 log::cout() << ">> Final number of cells... " << nCells << std::endl;
229 }
230
231#if BITPIT_ENABLE_MPI==1
233 patch_2D->partition(true);
234#endif
235
237 log::cout() << "Cell count: " << patch_2D->getCellCount() << std::endl;
238 log::cout() << "Vertex count: " << patch_2D->getVertexCount() << std::endl;
239
240 patch_2D->getVTK().setName("mesh_random.0");
241 patch_2D->write();
242
244#if BITPIT_ENABLE_MPI==1
245 VolOctreeMapper mapobject(patch_2D, patch_2D_original, patch_2D->getCommunicator());
246#else
247 VolOctreeMapper mapobject(patch_2D, patch_2D_original);
248#endif
249
251 mapobject.initialize(true);
252
254 PiercedStorage<double> data(1, &patch_2D_original->getCells());
255 std::vector<double> vdata(patch_2D_original->getInternalCellCount());
256 int count = 0;
257 for (Cell & cell : patch_2D_original->getCells()){
258 if (cell.isInterior()){
259 long id = cell.getId();
260 data[id] = double(patch_2D_original->getCellLevel(id));
261 vdata[count] = data[id];
262 count++;
263 }
264 }
265
266 patch_2D_original->getVTK().setName("mesh_original.0");
267 patch_2D_original->getVTK().addData("data", VTKFieldType::SCALAR, VTKLocation::CELL, vdata);
268#if BITPIT_ENABLE_MPI==1
269 patch_2D_original->setVTKWriteTarget(PatchKernel::WriteTarget::WRITE_TARGET_CELLS_INTERNAL);
270#endif
271 patch_2D_original->write();
272
274 PiercedStorage<double> data2(1, &patch_2D->getCells());
275 data2.fill(0.0);
276 {
277
278 const PiercedStorage<mapping::Info> & mapper = mapobject.getMapping();
279
280#if BITPIT_ENABLE_MPI==1
281 //Communicate needed data for mapping
282 std::map<int, std::map<long, double> > datarec;
283 {
284 std::unordered_map<int, std::vector<long> > rankIDrec = mapobject.getReceivedMappedIds();
285 std::unordered_map<int, std::vector<long> > rankIDsend = mapobject.getSentMappedIds();
286
287 //build send buffers
288 MPI_Comm comm = MPI_COMM_WORLD;
289 DataCommunicator dataCommunicator(comm);
290 MPI_Barrier(comm);
291 std::size_t bytes = uint8_t(sizeof(double));
292 for (const auto &val : rankIDsend){
293 int rank = val.first;
294 //set size
295 std::size_t buffSize = val.second.size() * bytes;
296 dataCommunicator.setSend(rank,buffSize);
297 //fill buffer with octants
298 SendBuffer &sendBuffer = dataCommunicator.getSendBuffer(rank);
299 for (long ID : val.second){
300 sendBuffer << data[ID];
301 }
302 }
303
304 dataCommunicator.discoverRecvs();
305 dataCommunicator.startAllRecvs();
306 dataCommunicator.startAllSends();
307
308 for (const auto &val : rankIDrec){
309 int rank = val.first;
310 dataCommunicator.waitRecv(rank);
311 RecvBuffer & recvBuffer = dataCommunicator.getRecvBuffer(rank);
312 for (long ID : val.second){
313 recvBuffer >> datarec[rank][ID];
314 }
315 }
316 dataCommunicator.waitAllSends();
317
318 }
319#endif
320
321 std::vector<double> vdata2(patch_2D->getInternalCellCount());
322 count = 0;
323 for (Cell & cell : patch_2D->getCells()){
324 if (cell.isInterior()){
325 long id = cell.getId();
326 if (mapper[id].type == mapping::Type::TYPE_RENUMBERING){
327#if BITPIT_ENABLE_MPI==1
328 if (mapper[id].ranks[0] == patch_2D->getRank()){
329#endif
330 data2[id] = data[mapper[id].ids[0]];
331 vdata2[count] = data2[id];
332#if BITPIT_ENABLE_MPI==1
333 }
334 else{
335 data2[id] = datarec[mapper[id].ranks[0]][mapper[id].ids[0]];
336 vdata2[count] = data2[id];
337 }
338#endif
339 }
340 else if (mapper[id].type == mapping::Type::TYPE_COARSENING){
341 data2[id] = 0.0;
342 int n = mapper[id].ids.size();
343 int i = 0;
344 for (long idd : mapper[id].ids){
345#if BITPIT_ENABLE_MPI==1
346 if (mapper[id].ranks[i] == patch_2D->getRank()){
347#endif
348 data2[id] += data[idd] / n;
349#if BITPIT_ENABLE_MPI==1
350 }
351 else{
352 data2[id] += datarec[mapper[id].ranks[i]][idd] / n;
353 }
354#endif
355 i++;
356 }
357 vdata2[count] = data2[id];
358 }
359 else if (mapper[id].type == mapping::Type::TYPE_REFINEMENT){
360#if BITPIT_ENABLE_MPI==1
361 if (mapper[id].ranks[0] == patch_2D->getRank()){
362#endif
363 data2[id] = data[mapper[id].ids[0]];
364 vdata2[count] = data2[id];
365#if BITPIT_ENABLE_MPI==1
366 }
367 else{
368 data2[id] = datarec[mapper[id].ranks[0]][mapper[id].ids[0]];
369 vdata2[count] = data2[id];
370 }
371#endif
372 }
373 count++;
374 }
375 }
376
377 patch_2D->getVTK().addData("data", VTKFieldType::SCALAR, VTKLocation::CELL, vdata2);
378 patch_2D->getVTK().setName("mesh_random.1");
379 patch_2D->write();
380
381 }
382
384 PiercedStorage<double> data2inv(1, &patch_2D_original->getCells());
385 data2inv.fill(0.0);
386 {
387
388 const PiercedStorage<mapping::Info> & mapper = mapobject.getInverseMapping();
389
390#if BITPIT_ENABLE_MPI==1
391 //Communicate needed data for mapping
392 std::map<int, std::map<long, double> > datarec;
393 {
394 // recover id to be rec/sent
395 std::map<int, std::set<long> > rankIDsend;
396
397 const PiercedStorage<mapping::Info> & _mapper = mapobject.getMapping();
398
399 for (Cell & cell : patch_2D->getCells()){
400 long ID = cell.getId();
401 auto info = _mapper[ID];
402 int i = 0;
403 for (int rank : info.ranks){
404 if (rank != patch_2D->getRank()){
405 rankIDsend[rank].insert(ID);
406 }
407 i++;
408 }
409 }
410
411 //build send buffers
412 MPI_Comm comm = MPI_COMM_WORLD;
413 DataCommunicator dataCommunicator(comm);
414 MPI_Barrier(comm);
415 std::size_t bytes = uint8_t(sizeof(double) + sizeof(long));
416 for (const auto &val : rankIDsend){
417 int rank = val.first;
418 //set size
419 std::size_t buffSize = val.second.size() * bytes;
420 dataCommunicator.setSend(rank,buffSize);
421 //fill buffer with octants
422 SendBuffer &sendBuffer = dataCommunicator.getSendBuffer(rank);
423 for (long ID : val.second){
424 sendBuffer << ID;
425 sendBuffer << data2[ID];
426 }
427 }
428
429 dataCommunicator.discoverRecvs();
430 dataCommunicator.startAllRecvs();
431 dataCommunicator.startAllSends();
432
433 std::vector<int> recvRanks = dataCommunicator.getRecvRanks();
434 std::sort(recvRanks.begin(),recvRanks.end());
435
436 for (int rank : recvRanks){
437 dataCommunicator.waitRecv(rank);
438 RecvBuffer & recvBuffer = dataCommunicator.getRecvBuffer(rank);
439 long nof = recvBuffer.getSize() / bytes;
440 for (long ii = 0; ii < nof; ii++){
441 long ID;
442 recvBuffer >> ID;
443 recvBuffer >> datarec[rank][ID];
444 }
445 }
446 dataCommunicator.waitAllSends();
447
448 }
449#endif
450
451 std::vector<double> vdata2inv(patch_2D_original->getInternalCellCount());
452 count = 0;
453 for (Cell & cell : patch_2D_original->getCells()){
454 if (cell.isInterior()){
455 long id = cell.getId();
456 if (mapper[id].type == mapping::Type::TYPE_RENUMBERING){
457#if BITPIT_ENABLE_MPI==1
458 if (mapper[id].ranks[0] == patch_2D_original->getRank()){
459#endif
460 data2inv[id] = data2[mapper[id].ids[0]];
461 vdata2inv[count] = data2inv[id];
462#if BITPIT_ENABLE_MPI==1
463 }
464 else{
465 data2inv[id] = datarec[mapper[id].ranks[0]][mapper[id].ids[0]];
466 vdata2inv[count] = data2inv[id];
467 }
468#endif
469 }
470 else if (mapper[id].type == mapping::Type::TYPE_COARSENING){
471 data2inv[id] = 0.0;
472 int n = mapper[id].ids.size();
473 int i = 0;
474 for (long idd : mapper[id].ids){
475#if BITPIT_ENABLE_MPI==1
476 if (mapper[id].ranks[i] == patch_2D_original->getRank()){
477#endif
478 data2inv[id] += data2[idd] / n;
479#if BITPIT_ENABLE_MPI==1
480 }
481 else{
482 data2inv[id] += datarec[mapper[id].ranks[i]][idd] / n;
483 }
484#endif
485 i++;
486 }
487 vdata2inv[count] = data2inv[id];
488 }
489 else if (mapper[id].type == mapping::Type::TYPE_REFINEMENT){
490#if BITPIT_ENABLE_MPI==1
491 if (mapper[id].ranks[0] == patch_2D_original->getRank()){
492#endif
493 data2inv[id] = data2[mapper[id].ids[0]];
494 vdata2inv[count] = data2inv[id];
495#if BITPIT_ENABLE_MPI==1
496 }
497 else{
498 data2inv[id] = datarec[mapper[id].ranks[0]][mapper[id].ids[0]];
499 vdata2inv[count] = data2inv[id];
500 }
501#endif
502 }
503 count++;
504 }
505 }
506
507 patch_2D_original->getVTK().addData("data", VTKFieldType::SCALAR, VTKLocation::CELL, vdata2inv);
508 patch_2D_original->getVTK().setName("mesh_original.1");
509 patch_2D_original->write();
510
511 }
512}
513
518int main(int argc, char *argv[])
519{
520#if BITPIT_ENABLE_MPI==1
521 MPI_Init(&argc,&argv);
522#endif
523
525 try {
526 run();
527 } catch (const std::exception &exception) {
528 log::cout() << exception.what();
529 exit(1);
530 }
531
532#if BITPIT_ENABLE_MPI==1
533 MPI_Finalize();
534#endif
535
536}
bitpit::Cell
The Cell class defines the cells.
Definition cell.hpp:42
bitpit::DataCommunicator
The DataCommunicator class provides the infrastructure needed to exchange data among processes.
Definition communications.hpp:42
bitpit::PabloUniform
PABLO Uniform is an example of user class derived from ParaTree to map ParaTree in a uniform (square/...
Definition PabloUniform.hpp:62
bitpit::ParaTree::getLocalIdx
uint32_t getLocalIdx(uint64_t gidx) const
Definition ParaTree.cpp:1616
bitpit::ParaTree::adaptGlobalRefine
bool adaptGlobalRefine(bool mapper_flag=false)
Definition ParaTree.cpp:3573
bitpit::ParaTree::getOwnerRank
int getOwnerRank(uint64_t globalIdx) const
Definition ParaTree.cpp:3797
bitpit::PatchKernel::initializeAdjacencies
void initializeAdjacencies(AdjacenciesBuildStrategy strategy=ADJACENCIES_AUTOMATIC)
Definition patch_kernel.cpp:5853
bitpit::PatchKernel::markCellForRefinement
void markCellForRefinement(long id)
Definition patch_kernel.cpp:890
bitpit::PatchKernel::update
std::vector< adaption::Info > update(bool trackAdaption=true, bool squeezeStorage=false)
Definition patch_kernel.cpp:625
bitpit::PatchKernel::getCells
PiercedVector< Cell > & getCells()
Definition patch_kernel.cpp:2551
bitpit::PatchKernel::getVertexCount
virtual long getVertexCount() const
Definition patch_kernel.cpp:1644
bitpit::PatchKernel::getVTK
VTKUnstructuredGrid & getVTK()
Definition patch_kernel.cpp:7784
bitpit::PatchKernel::partition
std::vector< adaption::Info > partition(MPI_Comm communicator, const std::unordered_map< long, int > &cellRanks, bool trackPartitioning, bool squeezeStorage=false, std::size_t haloSize=1)
Definition patch_kernel_parallel.cpp:1317
bitpit::PatchKernel::markCellForCoarsening
void markCellForCoarsening(long id)
Definition patch_kernel.cpp:904
bitpit::PatchKernel::getCommunicator
const MPI_Comm & getCommunicator() const
Definition patch_kernel_parallel.cpp:128
bitpit::PatchKernel::initializeInterfaces
void initializeInterfaces(InterfacesBuildStrategy strategy=INTERFACES_AUTOMATIC)
Definition patch_kernel.cpp:6349
bitpit::PatchKernel::setVTKWriteTarget
void setVTKWriteTarget(WriteTarget targetCells)
Definition patch_kernel.cpp:7814
bitpit::PatchKernel::getCellCount
virtual long getCellCount() const
Definition patch_kernel.cpp:2521
bitpit::PatchKernel::write
void write(VTKWriteMode mode=VTKWriteMode::DEFAULT)
Definition patch_kernel.cpp:1179
bitpit::PatchKernel::getInternalCellCount
long getInternalCellCount() const
Definition patch_kernel.cpp:2531
bitpit::PiercedStorage
Metafunction for generating a pierced storage.
Definition piercedStorage.hpp:146
bitpit::RecvBuffer
Buffer to be used for receive communications.
Definition communications_buffers.hpp:107
bitpit::SendBuffer
Buffer to be used for send communications.
Definition communications_buffers.hpp:91
bitpit::VTK::addData
VTKField & addData(VTKField &&field)
Definition VTK.cpp:281
bitpit::VTK::setName
void setName(const std::string &)
Definition VTK.cpp:119
bitpit::VolOctreeMapper
The VolOctreeMapper is the class to map two meshes of class VolOctree.
Definition voloctree_mapper.hpp:41
bitpit::VolOctree
The VolOctree defines a Octree patch.
Definition voloctree.hpp:37
bitpit::VolOctree::getCellLevel
int getCellLevel(long id) const
Definition voloctree.cpp:906
bitpit::VolOctree::getTree
PabloUniform & getTree()
Gets a reference to the octree associated with the patch.
Definition voloctree.cpp:795
bitpit::VolOctree::getOctantId
long getOctantId(const OctantInfo &octantInfo) const
Definition voloctree.cpp:827
bitpit::log::cout
Logger & cout(log::Level defaultSeverity, log::Visibility defaultVisibility)
Definition logger.cpp:1705
bitpit::log::info
Logger & info(log::Visibility defaultVisibility)
Definition logger.cpp:1847
--- layout: doxygen_footer ---
Generated on Wed Jun 19 2024 10:21:04 for bitpit by   doxygen 1.10.0