Loading...
Searching...
No Matches
pod_voloctree.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
25#include "pod_voloctree.hpp"
26
27namespace bitpit {
28
41# if BITPIT_ENABLE_MPI
45PODVolOctree::PODVolOctree(MPI_Comm comm):PODKernel(comm)
46# else
47PODVolOctree::PODVolOctree():PODKernel()
48# endif
49{
50}
51
58
64std::unique_ptr<VolumeKernel> PODVolOctree::createMesh()
65{
66#if BITPIT_ENABLE_MPI
67 return std::unique_ptr<VolumeKernel>(new VolOctree(getCommunicator()));
68#else
69 return std::unique_ptr<VolumeKernel>(new VolOctree());
70#endif
71}
72
79std::unique_ptr<VolumeMapper> PODVolOctree::_computeMapper(const VolumeKernel* mesh, bool fillInv)
80{
81 const VolOctree* meshPOD = static_cast<const VolOctree*>(getMesh());
82 const VolOctree* _mesh = static_cast<const VolOctree*>(mesh);
83
84 std::unique_ptr<VolumeMapper> mapper;
85#if BITPIT_ENABLE_MPI
86 mapper = std::unique_ptr<VolumeMapper>(new VolOctreeMapper(meshPOD, _mesh, getCommunicator()));
87#else
88 mapper = std::unique_ptr<VolumeMapper>(new VolOctreeMapper(meshPOD, _mesh));
89#endif
90
91 mapper->initialize(fillInv);
92
93 return mapper;
94}
95
103pod::PODField PODVolOctree::mapPODFieldToPOD(const pod::PODField & field, const std::unordered_set<long> * targetCells)
104{
105
106 //Check target cells
107 std::unordered_set<long> targetCellsStorage;
108 if (!targetCells) {
109 for (const Cell &cell : getMesh()->getCells())
110 targetCellsStorage.insert(cell.getId());
111
112 targetCells = &targetCellsStorage;
113 }
114
115 // Map data fields on pod mesh
116 std::size_t nsf = field.scalar->getFieldCount();
117 std::size_t nvf = field.vector->getFieldCount();
118
119 const PiercedStorage<mapping::Info> &mappingInfo = getMapper()->getMapping();
120
121 pod::PODField mappedField(nsf, nvf, getMesh(), &getMesh()->getCells());
122
123
124#if BITPIT_ENABLE_MPI
125 if ( getMapper()->checkPartition()){
126#endif
127
128 for (long id : *targetCells){
129 std::size_t rawIndex = m_meshPOD->getCells().getRawIndex(id);
130
131 double *datamappedS = mappedField.scalar->data(id);
132 std::array<double,3> *datamappedV = mappedField.vector->data(id);
133
134 if (mappingInfo[id].type == mapping::Type::TYPE_RENUMBERING){
135 bool dataB = field.mask->at(mappingInfo[id].ids[0]);
136 mappedField.mask->set(id, dataB);
137
138 double *dataS = field.scalar->data(mappingInfo[id].ids[0]);
139 for (std::size_t i = 0; i < nsf; i++){
140 double *dataSi = dataS + i;
141 double *datamappedSi = datamappedS + i;
142 (*datamappedSi) = (*dataSi);
143 }
144
145 std::array<double,3> *dataV = field.vector->data(mappingInfo[id].ids[0]);
146 for (std::size_t i = 0; i < nvf; i++) {
147 std::array<double,3> *dataVi = dataV + i;
148 std::array<double,3> *datamappedVi = datamappedV + i;
149 (*datamappedVi) = (*dataVi);
150 }
151 }
152 else if (mappingInfo[id].type == mapping::Type::TYPE_COARSENING){
153 mappedField.mask->set(id, true);
154
155 for (std::size_t i = 0; i < nsf; i++){
156 double *datamappedSi = datamappedS + i;
157 (*datamappedSi) = 0.0;
158 }
159 for (std::size_t i = 0; i < nvf; i++){
160 std::array<double,3> *datamappedVi = datamappedV + i;
161 (*datamappedVi) = std::array<double,3>{0.0, 0.0, 0.0};
162 }
163
164 bool dataB, dataMappedB = true;
165 double *dataS;
166 std::array<double,3> *dataV;
167 double volmapped = getRawCellVolume(rawIndex);
168 for (long idd : mappingInfo[id].ids){
169 dataB = field.mask->at(idd);
170 dataMappedB &= dataB;
171
172 dataS = field.scalar->data(idd);
173 double vol = field.mesh->evalCellVolume(idd);
174 for (std::size_t i = 0; i < nsf; i++){
175 double *dataSi = dataS + i;
176 double *datamappedSi = datamappedS + i;
177 (*datamappedSi) += (*dataSi) * vol / volmapped;
178 }
179 dataV = field.vector->data(idd);
180 for (std::size_t i = 0; i < nvf; i++) {
181 std::array<double,3> *dataVi = dataV + i;
182 std::array<double,3> *datamappedVi = datamappedV + i;
183 (*datamappedVi) += (*dataVi) * vol / volmapped;
184 }
185 }
186 mappedField.mask->set(id, dataMappedB);
187
188 }
189 else if (mappingInfo[id].type == mapping::Type::TYPE_REFINEMENT){
190 bool dataB = field.mask->at(mappingInfo[id].ids[0]);
191 double *dataS = field.scalar->data(mappingInfo[id].ids[0]);
192 std::array<double,3> *dataV = field.vector->data(mappingInfo[id].ids[0]);
193
194 mappedField.mask->set(id, dataB);
195
196 for (std::size_t i = 0; i < nsf; i++){
197 double *dataSi = dataS + i;
198 double *datamappedSi = datamappedS + i;
199 (*datamappedSi) = (*dataSi);
200 }
201 for (std::size_t i = 0; i < nvf; i++) {
202 std::array<double,3> *dataVi = dataV + i;
203 std::array<double,3> *datamappedVi = datamappedV + i;
204 (*datamappedVi) = (*dataVi);
205 }
206 }
207 }
208
209#if BITPIT_ENABLE_MPI
210 }
211 else{
212
213 //Mapping between partitioned meshes
214
215 //Communicate data and volumes
216 std::map<int, std::map<long, bool> > dataBrec;
217 std::map<int, std::map<long, std::vector<double> > > dataSrec;
218 std::map<int, std::map<long, std::vector<std::array<double, 3> > > > dataVrec;
219 std::map<int, std::map<long, double> > volrec;
220
221 communicatePODField(field, dataBrec, dataSrec, dataVrec, volrec);
222
223 for (long id : *targetCells){
224
225 double *datamappedS = mappedField.scalar->data(id);
226 std::array<double,3> *datamappedV = mappedField.vector->data(id);
227
228 if (mappingInfo[id].type == mapping::Type::TYPE_RENUMBERING){
229
230 int rank = mappingInfo[id].ranks[0];
231 if (rank != m_rank){
232 bool dataB = dataBrec[rank][mappingInfo[id].ids[0]];
233 mappedField.mask->set(id, dataB);
234
235 double *dataS = dataSrec[rank][mappingInfo[id].ids[0]].data();
236 for (std::size_t i = 0; i < nsf; i++){
237 double *dataSi = dataS + i;
238 double *datamappedSi = datamappedS + i;
239 (*datamappedSi) = (*dataSi);
240 }
241
242 std::array<double,3> *dataV = dataVrec[rank][mappingInfo[id].ids[0]].data();
243 for (std::size_t i = 0; i < nvf; i++) {
244 std::array<double,3> *dataVi = dataV + i;
245 std::array<double,3> *datamappedVi = datamappedV + i;
246 (*datamappedVi) = (*dataVi);
247 }
248
249 }
250 else{
251 bool dataB = field.mask->at(mappingInfo[id].ids[0]);
252 mappedField.mask->set(id, dataB);
253
254 double *dataS = field.scalar->data(mappingInfo[id].ids[0]);
255 for (std::size_t i = 0; i < nsf; i++){
256 double *dataSi = dataS + i;
257 double *datamappedSi = datamappedS + i;
258 (*datamappedSi) = (*dataSi);
259 }
260
261 std::array<double,3> *dataV = field.vector->data(mappingInfo[id].ids[0]);
262 for (std::size_t i = 0; i < nvf; i++) {
263 std::array<double,3> *dataVi = dataV + i;
264 std::array<double,3> *datamappedVi = datamappedV + i;
265 (*datamappedVi) = (*dataVi);
266 }
267 }//end if other rank
268
269 }
270 else if (mappingInfo[id].type == mapping::Type::TYPE_COARSENING){
271 mappedField.mask->set(id, true);
272
273 for (std::size_t i = 0; i < nsf; i++){
274 double *datamappedSi = datamappedS + i;
275 (*datamappedSi) = 0.0;
276 }
277 for (std::size_t i = 0; i < nvf; i++){
278 std::array<double,3> *datamappedVi = datamappedV + i;
279 (*datamappedVi) = std::array<double,3>{0.0, 0.0, 0.0};
280 }
281
282 bool dataB, dataMappedB = true;
283 double *dataS;
284 std::array<double,3> *dataV;
285 double vol;
286 double volmapped = getCellVolume(id);
287 int ii = 0;
288 for (long idd : mappingInfo[id].ids){
289 int rank = mappingInfo[id].ranks[ii];
290 if (rank != m_rank){
291 dataB = dataBrec[rank][idd];
292 dataS = dataSrec[rank][idd].data();
293 dataV = dataVrec[rank][idd].data();
294 vol = volrec[rank][idd];
295 }
296 else{
297 dataB = field.mask->at(idd);
298 dataS = field.scalar->data(idd);
299 dataV = field.vector->data(idd);
300 vol = field.mesh->evalCellVolume(idd);
301 }
302 dataMappedB &= dataB;
303
304 for (std::size_t i = 0; i < nsf; i++){
305 double *dataSi = dataS + i;
306 double *datamappedSi = datamappedS + i;
307 (*datamappedSi) += (*dataSi) * vol / volmapped;
308 }
309 for (std::size_t i = 0; i < nvf; i++) {
310 std::array<double,3> *dataVi = dataV + i;
311 std::array<double,3> *datamappedVi = datamappedV + i;
312 (*datamappedVi) += (*dataVi) * vol / volmapped;
313 }
314 ii++;
315 }
316 mappedField.mask->set(id, dataMappedB);
317
318 }
319 else if (mappingInfo[id].type == mapping::Type::TYPE_REFINEMENT){
320 bool dataB;
321 double *dataS;
322 std::array<double,3> *dataV;
323 int rank = mappingInfo[id].ranks[0];
324 if (rank != m_rank){
325 dataB = dataBrec[rank][mappingInfo[id].ids[0]];
326 dataS = dataSrec[rank][mappingInfo[id].ids[0]].data();
327 dataV = dataVrec[rank][mappingInfo[id].ids[0]].data();
328 }
329 else{
330 dataB = field.mask->at(mappingInfo[id].ids[0]);
331 dataS = field.scalar->data(mappingInfo[id].ids[0]);
332 dataV = field.vector->data(mappingInfo[id].ids[0]);
333 }
334
335 mappedField.mask->set(id, dataB);
336
337 for (std::size_t i = 0; i < nsf; i++){
338 double *dataSi = dataS + i;
339 double *datamappedSi = datamappedS + i;
340 (*datamappedSi) = (*dataSi);
341 }
342 for (std::size_t i = 0; i < nvf; i++) {
343 std::array<double,3> *dataVi = dataV + i;
344 std::array<double,3> *datamappedVi = datamappedV + i;
345 (*datamappedVi) = (*dataVi);
346 }
347 }
348 } // end target cells
349
350 }//end partitioned meshes
351#endif
352 return mappedField;
353
354}
355
363void PODVolOctree::mapPODFieldFromPOD(pod::PODField & field, const std::unordered_set<long> * targetCells,
364 const pod::PODField & mappedField)
365{
366
367 //Check target cells
368 std::unordered_set<long> targetCellsStorage;
369 if (!targetCells) {
370 for (const Cell &cell : field.mesh->getCells())
371 targetCellsStorage.insert(cell.getId());
372
373 targetCells = &targetCellsStorage;
374 }
375
376 // Map data fields from pod mesh
377 std::size_t nsf = field.scalar->getFieldCount();
378 std::size_t nvf = field.vector->getFieldCount();
379
380 const PiercedStorage<mapping::Info> &inverseMapping = getMapper()->getInverseMapping();
381
382#if BITPIT_ENABLE_MPI
383 if ( getMapper()->checkPartition()){
384#endif
385
386
387 for (long id : *targetCells){
388 double *dataS = field.scalar->data(id);
389 std::array<double,3> *dataV = field.vector->data(id);
390
391 if (inverseMapping[id].type == mapping::Type::TYPE_RENUMBERING){
392 bool datamappedB = mappedField.mask->at(inverseMapping[id].ids[0]);
393 field.mask->set(id, datamappedB);
394 double *datamappedS = mappedField.scalar->data(inverseMapping[id].ids[0]);
395 for (std::size_t i = 0; i < nsf; i++){
396 double *dataSi = dataS + i;
397 double *datamappedSi = datamappedS + i;
398 (*dataSi) = (*datamappedSi);
399 }
400 std::array<double,3> *datamappedV = mappedField.vector->data(inverseMapping[id].ids[0]);
401 for (std::size_t i = 0; i < nvf; i++) {
402 std::array<double,3> *dataVi = dataV + i;
403 std::array<double,3> *datamappedVi = datamappedV + i;
404 (*dataVi) = (*datamappedVi);
405 }
406 }
407 else if (inverseMapping[id].type == mapping::Type::TYPE_COARSENING){
408 bool dataB = true;
409 for (std::size_t i = 0; i < nsf; i++){
410 double *dataSi = dataS + i;
411 (*dataSi) = 0.0;
412 }
413 for (std::size_t i = 0; i < nvf; i++){
414 std::array<double,3> *dataVi = dataV + i;
415 (*dataVi) = std::array<double,3>{0.0, 0.0, 0.0};
416 }
417 bool datamappedB;
418 double *datamappedS;
419 std::array<double,3> *datamappedV;
420 double vol = field.mesh->evalCellVolume(id);
421 for (long idd : inverseMapping[id].ids){
422 std::size_t rawIndexIdd = m_meshPOD->getCells().getRawIndex(idd);
423 datamappedB = mappedField.mask->at(idd);
424 dataB &= datamappedB;
425 datamappedS = mappedField.scalar->data(idd);
426 double volmapped = getRawCellVolume(rawIndexIdd);
427 for (std::size_t i = 0; i < nsf; i++){
428 double *dataSi = dataS + i;
429 double *datamappedSi = datamappedS + i;
430 (*dataSi) += (*datamappedSi) * volmapped / vol;
431 }
432 datamappedV = mappedField.vector->data(idd);
433 for (std::size_t i = 0; i < nvf; i++) {
434 std::array<double,3> *dataVi = dataV + i;
435 std::array<double,3> *datamappedVi = datamappedV + i;
436 (*dataVi) += (*datamappedVi) * volmapped / vol;
437 }
438 }
439 field.mask->set(id, dataB);
440 }
441 else if (inverseMapping[id].type == mapping::Type::TYPE_REFINEMENT){
442 bool datamappedB = mappedField.mask->at(inverseMapping[id].ids[0]);
443 field.mask->set(id, datamappedB);
444 double *datamappedS = mappedField.scalar->data(inverseMapping[id].ids[0]);
445 for (std::size_t i = 0; i < nsf; i++){
446 double *dataSi = dataS + i;
447 double *datamappedSi = datamappedS + i;
448 (*dataSi) = (*datamappedSi);
449 }
450 std::array<double,3> *datamappedV = mappedField.vector->data(inverseMapping[id].ids[0]);
451 for (std::size_t i = 0; i < nvf; i++) {
452 std::array<double,3> *dataVi = dataV + i;
453 std::array<double,3> *datamappedVi = datamappedV + i;
454 (*dataVi) = (*datamappedVi);
455 }
456 }
457 }//end target cells
458
459#if BITPIT_ENABLE_MPI
460 }
461 else{
462
463 //Mapping between partitioned meshes
464
465 //Communicate data and volumes
466 std::map<int, std::map<long, bool> > dataBrec;
467 std::map<int, std::map<long, std::vector<double> > > dataSrec;
468 std::map<int, std::map<long, std::vector<std::array<double, 3> > > > dataVrec;
469 std::map<int, std::map<long, double> > volrec;
470
471 communicatePODFieldFromPOD(mappedField, dataBrec, dataSrec, dataVrec, volrec);
472
473 for (long id : *targetCells){
474 double *dataS = field.scalar->data(id);
475 std::array<double,3> *dataV = field.vector->data(id);
476
477 if (inverseMapping[id].type == mapping::Type::TYPE_RENUMBERING){
478 bool datamappedB;
479 double *datamappedS;
480 std::array<double,3> *datamappedV;
481 int rank = inverseMapping[id].ranks[0];
482 if (rank != m_rank){
483 datamappedB = dataBrec[rank][inverseMapping[id].ids[0]];
484 datamappedS = dataSrec[rank][inverseMapping[id].ids[0]].data();
485 datamappedV = dataVrec[rank][inverseMapping[id].ids[0]].data();
486 }
487 else{
488 datamappedB = mappedField.mask->at(inverseMapping[id].ids[0]);
489 datamappedS = mappedField.scalar->data(inverseMapping[id].ids[0]);
490 datamappedV = mappedField.vector->data(inverseMapping[id].ids[0]);
491 }//endif different rank
492
493 field.mask->set(id, datamappedB);
494 for (std::size_t i = 0; i < nsf; i++){
495 double *dataSi = dataS + i;
496 double *datamappedSi = datamappedS + i;
497 (*dataSi) = (*datamappedSi);
498 }
499 for (std::size_t i = 0; i < nvf; i++) {
500 std::array<double,3> *dataVi = dataV + i;
501 std::array<double,3> *datamappedVi = datamappedV + i;
502 (*dataVi) = (*datamappedVi);
503 }
504 }
505 else if (inverseMapping[id].type == mapping::Type::TYPE_COARSENING){
506 bool dataB = true;
507 for (std::size_t i = 0; i < nsf; i++){
508 double *dataSi = dataS + i;
509 (*dataSi) = 0.0;
510 }
511 for (std::size_t i = 0; i < nvf; i++){
512 std::array<double,3> *dataVi = dataV + i;
513 (*dataVi) = std::array<double,3>{0.0, 0.0, 0.0};
514 }
515 bool datamappedB;
516 double *datamappedS;
517 std::array<double,3> *datamappedV;
518 double volmapped;
519 double vol = field.mesh->evalCellVolume(id);
520 int ii=0;
521 for (long idd : inverseMapping[id].ids){
522 int rank = inverseMapping[id].ranks[ii];
523 if (rank != m_rank){
524 datamappedB = dataBrec[rank][idd];
525 datamappedS = dataSrec[rank][idd].data();
526 datamappedV = dataVrec[rank][idd].data();
527 volmapped = volrec[rank][idd];
528 }
529 else{
530 datamappedB = mappedField.mask->at(idd);
531 datamappedS = mappedField.scalar->data(idd);
532 datamappedV = mappedField.vector->data(idd);
533 volmapped = getCellVolume(idd);
534 }//endif different rank
535
536 dataB &= datamappedB;
537 for (std::size_t i = 0; i < nsf; i++){
538 double *dataSi = dataS + i;
539 double *datamappedSi = datamappedS + i;
540 (*dataSi) += (*datamappedSi) * volmapped / vol;
541 }
542 for (std::size_t i = 0; i < nvf; i++) {
543 std::array<double,3> *dataVi = dataV + i;
544 std::array<double,3> *datamappedVi = datamappedV + i;
545 (*dataVi) += (*datamappedVi) * volmapped / vol;
546 }
547 ii++;
548 }
549 field.mask->set(id, dataB);
550 }
551 else if (inverseMapping[id].type == mapping::Type::TYPE_REFINEMENT){
552 bool datamappedB;
553 double *datamappedS;
554 std::array<double,3> *datamappedV;
555 int rank = inverseMapping[id].ranks[0];
556 if (rank != m_rank){
557 datamappedB = dataBrec[rank][inverseMapping[id].ids[0]];
558 datamappedS = dataSrec[rank][inverseMapping[id].ids[0]].data();
559 datamappedV = dataVrec[rank][inverseMapping[id].ids[0]].data();
560 }
561 else{
562 datamappedB = mappedField.mask->at(inverseMapping[id].ids[0]);
563 datamappedS = mappedField.scalar->data(inverseMapping[id].ids[0]);
564 datamappedV = mappedField.vector->data(inverseMapping[id].ids[0]);
565 }//endif different rank
566
567 field.mask->set(id, datamappedB);
568 for (std::size_t i = 0; i < nsf; i++){
569 double *dataSi = dataS + i;
570 double *datamappedSi = datamappedS + i;
571 (*dataSi) = (*datamappedSi);
572 }
573 for (std::size_t i = 0; i < nvf; i++) {
574 std::array<double,3> *dataVi = dataV + i;
575 std::array<double,3> *datamappedVi = datamappedV + i;
576 (*dataVi) = (*datamappedVi);
577 }
578 }
579 }//end target cells
580
581 }//endif partitioned
582#endif
583
584}
585
596PiercedStorage<double> PODVolOctree::mapFieldsToPOD(const PiercedStorage<double> & fields, const VolumeKernel * mesh,
597 const std::unordered_set<long> * targetCells,
598 const std::vector<std::size_t> &scalarIds,
599 const std::vector<std::array<std::size_t, 3>> &vectorIds)
600{
601
602 //Check target cells
603 std::unordered_set<long> targetCellsStorage;
604 if (!targetCells) {
605 for (const Cell &cell : getMesh()->getCells())
606 targetCellsStorage.insert(cell.getId());
607
608 targetCells = &targetCellsStorage;
609 }
610
611 // Map data fields on pod mesh
612 std::size_t nsf = scalarIds.size();
613 std::size_t nvf = vectorIds.size();
614
615 const PiercedStorage<mapping::Info> &mappingInfo = getMapper()->getMapping();
616
617 PiercedStorage<double> mappedFields(fields.getFieldCount(), &getMesh()->getCells());
618
619#if BITPIT_ENABLE_MPI
620 if ( getMapper()->checkPartition()){
621#endif
622
623 for (long id : *targetCells){
624 std::size_t rawIndex = m_meshPOD->getCells().getRawIndex(id);
625 double *datamapped = mappedFields.data(id);
626 if (mappingInfo[id].type == mapping::Type::TYPE_RENUMBERING){
627 const double *data = fields.data(mappingInfo[id].ids[0]);
628 for (std::size_t i = 0; i < nsf; i++){
629 const double *datai = data + scalarIds[i];
630 double *datamappedi = datamapped + scalarIds[i];
631 (*datamappedi) = (*datai);
632 }
633 for (std::size_t i = 0; i < nvf; i++) {
634 for (std::size_t j = 0; j < 3; j++) {
635 const double *datai = data + vectorIds[i][j];
636 double *datamappedi = datamapped + vectorIds[i][j];
637 (*datamappedi) = (*datai);
638 }
639 }
640 }
641 else if (mappingInfo[id].type == mapping::Type::TYPE_COARSENING){
642 for (std::size_t i = 0; i < nsf; i++){
643 double *datamappedi = datamapped + scalarIds[i];
644 (*datamappedi) = 0.0;
645 }
646 for (std::size_t i = 0; i < nvf; i++){
647 for (std::size_t j = 0; j < 3; j++) {
648 double *datamappedi = datamapped + vectorIds[i][j];
649 (*datamappedi) = 0.0;
650 }
651 }
652 double volmapped = getRawCellVolume(rawIndex);
653 for (long idd : mappingInfo[id].ids){
654 const double *data = fields.data(idd);
655 double vol = mesh->evalCellVolume(idd);
656
657 for (std::size_t i = 0; i < nsf; i++){
658 const double *datai = data + scalarIds[i];
659 double *datamappedi = datamapped + scalarIds[i];
660 (*datamappedi) += (*datai) * vol / volmapped;
661 }
662 for (std::size_t i = 0; i < nvf; i++) {
663 for (std::size_t j = 0; j < 3; j++) {
664 const double *datai = data + vectorIds[i][j];
665 double *datamappedi = datamapped + vectorIds[i][j];
666 (*datamappedi) += (*datai) * vol / volmapped;
667 }
668 }
669 }
670 }
671 else if (mappingInfo[id].type == mapping::Type::TYPE_REFINEMENT){
672 const double *data = fields.data(mappingInfo[id].ids[0]);
673 for (std::size_t i = 0; i < nsf; i++){
674 const double *datai = data + scalarIds[i];
675 double *datamappedi = datamapped + scalarIds[i];
676 (*datamappedi) = (*datai);
677 }
678 for (std::size_t i = 0; i < nvf; i++) {
679 for (std::size_t j = 0; j < 3; j++) {
680 const double *datai = data + vectorIds[i][j];
681 double *datamappedi = datamapped + vectorIds[i][j];
682 (*datamappedi) = (*datai);
683 }
684 }
685 }
686 }//end targetcells
687
688#if BITPIT_ENABLE_MPI
689 }
690 else{
691
692 //Mapping between partitioned meshes
693
694 //Communicate data and volumes
695 std::map<int, std::map<long, std::vector<double> > > datarec;
696 std::map<int, std::map<long, double> > volrec;
697
698 communicateField(fields, mesh, datarec, volrec);
699
700 for (long id : *targetCells){
701 double *datamapped = mappedFields.data(id);
702 if (mappingInfo[id].type == mapping::Type::TYPE_RENUMBERING){
703 const double *data;
704 int rank = mappingInfo[id].ranks[0];
705 if (rank != m_rank){
706 data = datarec[rank][mappingInfo[id].ids[0]].data();
707 }
708 else{
709 data = fields.data(mappingInfo[id].ids[0]);
710 }//endif different rank
711
712 for (std::size_t i = 0; i < nsf; i++){
713 const double *datai = data + scalarIds[i];
714 double *datamappedi = datamapped + scalarIds[i];
715 (*datamappedi) = (*datai);
716 }
717 for (std::size_t i = 0; i < nvf; i++) {
718 for (std::size_t j = 0; j < 3; j++) {
719 const double *datai = data + vectorIds[i][j];
720 double *datamappedi = datamapped + vectorIds[i][j];
721 (*datamappedi) = (*datai);
722 }
723 }
724 }
725 else if (mappingInfo[id].type == mapping::Type::TYPE_COARSENING){
726 for (std::size_t i = 0; i < nsf; i++){
727 double *datamappedi = datamapped + scalarIds[i];
728 (*datamappedi) = 0.0;
729 }
730 for (std::size_t i = 0; i < nvf; i++){
731 for (std::size_t j = 0; j < 3; j++) {
732 double *datamappedi = datamapped + vectorIds[i][j];
733 (*datamappedi) = 0.0;
734 }
735 }
736 double volmapped = getCellVolume(id);
737 int ii=0;
738 for (long idd : mappingInfo[id].ids){
739 int rank = mappingInfo[id].ranks[ii];
740 const double *data;
741 double vol;
742 if (rank != m_rank){
743 data = datarec[rank][idd].data();
744 vol = volrec[rank][idd];
745 }
746 else{
747 data = fields.data(idd);
748 vol = mesh->evalCellVolume(idd);
749 }//endif different rank
750
751 for (std::size_t i = 0; i < nsf; i++){
752 const double *datai = data + scalarIds[i];
753 double *datamappedi = datamapped + scalarIds[i];
754 (*datamappedi) += (*datai) * vol / volmapped;
755 }
756 for (std::size_t i = 0; i < nvf; i++) {
757 for (std::size_t j = 0; j < 3; j++) {
758 const double *datai = data + vectorIds[i][j];
759 double *datamappedi = datamapped + vectorIds[i][j];
760 (*datamappedi) += (*datai) * vol / volmapped;
761 }
762 }
763 ii++;
764 }
765 }
766 else if (mappingInfo[id].type == mapping::Type::TYPE_REFINEMENT){
767 const double* data;
768 int rank = mappingInfo[id].ranks[0];
769 if (rank != m_rank){
770 data = datarec[rank][mappingInfo[id].ids[0]].data();
771 }
772 else{
773 data = fields.data(mappingInfo[id].ids[0]);
774 }//endif different rank
775 for (std::size_t i = 0; i < nsf; i++){
776 const double *datai = data + scalarIds[i];
777 double *datamappedi = datamapped + scalarIds[i];
778 (*datamappedi) = (*datai);
779 }
780 for (std::size_t i = 0; i < nvf; i++) {
781 for (std::size_t j = 0; j < 3; j++) {
782 const double *datai = data + vectorIds[i][j];
783 double *datamappedi = datamapped + vectorIds[i][j];
784 (*datamappedi) = (*datai);
785 }
786 }
787 }
788 }//end targetcells
789
790 }//endif partitioned mesh
791#endif
792
793 return mappedFields;
794
795}
796
807void PODVolOctree::mapFieldsFromPOD(PiercedStorage<double> & fields, const VolumeKernel * mesh,
808 const std::unordered_set<long> * targetCells, const PiercedStorage<double> & mappedFields,
809 const std::vector<std::size_t> &scalarIds, const std::vector<std::array<std::size_t, 3>> &vectorIds)
810{
811
812 //Check target cells
813 std::unordered_set<long> targetCellsStorage;
814 if (!targetCells) {
815 for (const Cell &cell : mesh->getCells())
816 targetCellsStorage.insert(cell.getId());
817
818 targetCells = &targetCellsStorage;
819 }
820
821 // Map data fields from pod mesh
822 std::size_t nsf = scalarIds.size();
823 std::size_t nvf = vectorIds.size();
824
825 const PiercedStorage<mapping::Info> &inverseMapping = getMapper()->getInverseMapping();
826
827#if BITPIT_ENABLE_MPI
828 if ( getMapper()->checkPartition()){
829#endif
830
831 for (long id : *targetCells){
832 double *data = fields.data(id);
833 if (inverseMapping[id].type == mapping::Type::TYPE_RENUMBERING){
834 const double *datamapped = mappedFields.data(inverseMapping[id].ids[0]);
835 for (std::size_t i = 0; i < nsf; i++){
836 double *datai = data + scalarIds[i];
837 const double *datamappedi = datamapped + scalarIds[i];
838 (*datai) = (*datamappedi);
839 }
840 for (std::size_t i = 0; i < nvf; i++) {
841 for (std::size_t j = 0; j < 3; j++) {
842 double *datai = data + vectorIds[i][j];
843 const double *datamappedi = datamapped + vectorIds[i][j];
844 (*datai) = (*datamappedi);
845 }
846 }
847 }
848 else if (inverseMapping[id].type == mapping::Type::TYPE_COARSENING){
849 for (std::size_t i = 0; i < nsf; i++){
850 double *datai = data + scalarIds[i];
851 (*datai) = 0.0;
852 }
853 for (std::size_t i = 0; i < nvf; i++){
854 for (std::size_t j = 0; j < 3; j++) {
855 double *datai = data + vectorIds[i][j];
856 (*datai) = 0.0;
857 }
858 }
859 double vol = mesh->evalCellVolume(id);
860 for (long idd : inverseMapping[id].ids){
861 std::size_t rawIndexIdd = m_meshPOD->getCells().getRawIndex(idd);
862 const double * datamapped = mappedFields.data(idd);
863 double volmapped = getRawCellVolume(rawIndexIdd);
864
865 for (std::size_t i = 0; i < nsf; i++){
866 double *datai = data + scalarIds[i];
867 const double *datamappedi = datamapped + scalarIds[i];
868 (*datai) += (*datamappedi) * volmapped / vol;
869 }
870 for (std::size_t i = 0; i < nvf; i++) {
871 for (std::size_t j = 0; j < 3; j++) {
872 double *datai = data + vectorIds[i][j];
873 const double *datamappedi = datamapped + vectorIds[i][j];
874 (*datai) += (*datamappedi) * volmapped / vol;
875 }
876 }
877 }
878 }
879 else if (inverseMapping[id].type == mapping::Type::TYPE_REFINEMENT){
880 const double *datamapped = mappedFields.data(inverseMapping[id].ids[0]);
881 for (std::size_t i = 0; i < nsf; i++){
882 double *datai = data + scalarIds[i];
883 const double *datamappedi = datamapped + scalarIds[i];
884 (*datai) = (*datamappedi);
885 }
886 for (std::size_t i = 0; i < nvf; i++) {
887 for (std::size_t j = 0; j < 3; j++) {
888 double *datai = data + vectorIds[i][j];
889 const double *datamappedi = datamapped + vectorIds[i][j];
890 (*datai) = (*datamappedi);
891 }
892 }
893 }
894 }//end target cells
895
896#if BITPIT_ENABLE_MPI
897 }
898 else{
899
900 //Mapping between partitioned meshes
901
902 //Communicate data and volumes
903 std::map<int, std::map<long, std::vector<double> > > datarec;
904 std::map<int, std::map<long, double> > volrec;
905
906 communicateFieldFromPOD(mappedFields, mesh, datarec, volrec);
907
908 for (long id : *targetCells){
909 double *data = fields.data(id);
910 if (inverseMapping[id].type == mapping::Type::TYPE_RENUMBERING){
911 const double *datamapped;
912 int rank = inverseMapping[id].ranks[0];
913 if (rank != m_rank){
914 datamapped = datarec[rank][inverseMapping[id].ids[0]].data();
915 }
916 else{
917 datamapped = mappedFields.data(inverseMapping[id].ids[0]);
918 }
919
920 for (std::size_t i = 0; i < nsf; i++){
921 double *datai = data + scalarIds[i];
922 const double *datamappedi = datamapped + scalarIds[i];
923 (*datai) = (*datamappedi);
924 }
925 for (std::size_t i = 0; i < nvf; i++) {
926 for (std::size_t j = 0; j < 3; j++) {
927 double *datai = data + vectorIds[i][j];
928 const double *datamappedi = datamapped + vectorIds[i][j];
929 (*datai) = (*datamappedi);
930 }
931 }
932 }
933 else if (inverseMapping[id].type == mapping::Type::TYPE_COARSENING){
934 for (std::size_t i = 0; i < nsf; i++){
935 double *datai = data + scalarIds[i];
936 (*datai) = 0.0;
937 }
938 for (std::size_t i = 0; i < nvf; i++){
939 for (std::size_t j = 0; j < 3; j++) {
940 double *datai = data + vectorIds[i][j];
941 (*datai) = 0.0;
942 }
943 }
944 double vol = mesh->evalCellVolume(id);
945 int ii = 0;
946 for (long idd : inverseMapping[id].ids){
947 const double *datamapped;
948 double volmapped;
949 int rank = inverseMapping[id].ranks[ii];
950 if (rank != m_rank){
951 datamapped = datarec[rank][idd].data();
952 volmapped = volrec[rank][idd];
953 }
954 else{
955 datamapped = mappedFields.data(idd);
956 volmapped = getCellVolume(idd);
957 }
958
959 for (std::size_t i = 0; i < nsf; i++){
960 double *datai = data + scalarIds[i];
961 const double *datamappedi = datamapped + scalarIds[i];
962 (*datai) += (*datamappedi) * volmapped / vol;
963 }
964 for (std::size_t i = 0; i < nvf; i++) {
965 for (std::size_t j = 0; j < 3; j++) {
966 double *datai = data + vectorIds[i][j];
967 const double *datamappedi = datamapped + vectorIds[i][j];
968 (*datai) += (*datamappedi) * volmapped / vol;
969 }
970 }
971 ii++;
972 }
973 }
974 else if (inverseMapping[id].type == mapping::Type::TYPE_REFINEMENT){
975 const double *datamapped;
976 int rank = inverseMapping[id].ranks[0];
977 if (rank != m_rank){
978 datamapped = datarec[rank][inverseMapping[id].ids[0]].data();
979 }
980 else{
981 datamapped = mappedFields.data(inverseMapping[id].ids[0]);
982 }
983 for (std::size_t i = 0; i < nsf; i++){
984 double *datai = data + scalarIds[i];
985 const double *datamappedi = datamapped + scalarIds[i];
986 (*datai) = (*datamappedi);
987 }
988 for (std::size_t i = 0; i < nvf; i++) {
989 for (std::size_t j = 0; j < 3; j++) {
990 double *datai = data + vectorIds[i][j];
991 const double *datamappedi = datamapped + vectorIds[i][j];
992 (*datai) = (*datamappedi);
993 }
994 }
995 }
996 }//end target cells
997
998
999
1000 }//endif partitioned
1001#endif
1002
1003
1004}
1005
1014PiercedStorage<bool> PODVolOctree::mapBoolFieldToPOD(const PiercedStorage<bool> & field, const VolumeKernel * mesh,
1015 const std::unordered_set<long> * targetCells)
1016{
1017
1018 // Map bool field (number of fields = 1) on pod mesh
1019 PiercedStorage<bool> mappedField(1, &getMesh()->getCells());
1020 mapBoolFieldToPOD(field, mesh, targetCells, mappedField);
1021 return mappedField;
1022
1023}
1024
1033void PODVolOctree::mapBoolFieldToPOD(const PiercedStorage<bool> & field, const VolumeKernel * mesh,
1034 const std::unordered_set<long> * targetCells, PiercedStorage<bool> & mappedField)
1035{
1036 BITPIT_UNUSED(mesh);
1037
1038 //Check target cells
1039 std::unordered_set<long> targetCellsStorage;
1040 if (!targetCells) {
1041 for (const Cell &cell : getMesh()->getCells())
1042 targetCellsStorage.insert(cell.getId());
1043
1044 targetCells = &targetCellsStorage;
1045 }
1046
1047 // Map bool field (number of fields = 1) on pod mesh
1048 const PiercedStorage<mapping::Info> &mappingInfo = getMapper()->getMapping();
1049
1050 mappedField.unsetKernel();
1051 mappedField.setStaticKernel(&getMesh()->getCells());
1052 mappedField.fill(false);
1053
1054#if BITPIT_ENABLE_MPI
1055 if ( getMapper()->checkPartition()){
1056#endif
1057
1058 for (long id : *targetCells){
1059 if (mappingInfo[id].type == mapping::Type::TYPE_RENUMBERING){
1060 bool dataB = field.at(mappingInfo[id].ids[0]);
1061 mappedField.set(id, dataB);
1062 }
1063 else if (mappingInfo[id].type == mapping::Type::TYPE_COARSENING){
1064 mappedField.set(id, false);
1065 bool dataB, dataMappedB = true;
1066 for (long idd : mappingInfo[id].ids){
1067 dataB = field.at(idd);
1068 dataMappedB &= dataB;
1069 }
1070 mappedField.set(id, dataMappedB);
1071 }
1072 else if (mappingInfo[id].type == mapping::Type::TYPE_REFINEMENT){
1073 bool dataB = field.at(mappingInfo[id].ids[0]);
1074 mappedField.set(id, dataB);
1075 }
1076 }
1077
1078#if BITPIT_ENABLE_MPI
1079 }
1080 else{
1081 //Mapping between partitioned meshes
1082
1083 //Communicate data
1084 std::map<int, std::map<long, bool> > dataBrec;
1085
1086 communicateBoolField(field, dataBrec);
1087
1088 for (long id : *targetCells){
1089
1090 if (mappingInfo[id].type == mapping::Type::TYPE_RENUMBERING){
1091 bool dataB;
1092 int rank = mappingInfo[id].ranks[0];
1093 if (rank != m_rank){
1094 dataB = dataBrec[rank][mappingInfo[id].ids[0]];
1095 }
1096 else{
1097 dataB = field.at(mappingInfo[id].ids[0]);
1098 }//end if other rank
1099 mappedField.set(id, dataB);
1100 }
1101 else if (mappingInfo[id].type == mapping::Type::TYPE_COARSENING){
1102 mappedField.set(id, true);
1103 bool dataB, dataMappedB = true;
1104 int ii = 0;
1105 for (long idd : mappingInfo[id].ids){
1106 int rank = mappingInfo[id].ranks[ii];
1107 if (rank != m_rank){
1108 dataB = dataBrec[rank][idd];
1109 }
1110 else{
1111 dataB = field.at(idd);
1112 }//end if other rank
1113 dataMappedB &= dataB;
1114 ii++;
1115 }
1116 mappedField.set(id, dataMappedB);
1117
1118 }
1119 else if (mappingInfo[id].type == mapping::Type::TYPE_REFINEMENT){
1120 bool dataB;
1121 int rank = mappingInfo[id].ranks[0];
1122 if (rank != m_rank){
1123 dataB = dataBrec[rank][mappingInfo[id].ids[0]];
1124 }
1125 else{
1126 dataB = field.at(mappingInfo[id].ids[0]);
1127 }//end if other rank
1128 mappedField.set(id, dataB);
1129 }
1130 }
1131
1132 }//end if partitioned
1133#endif
1134
1135}
1136
1143std::unordered_set<long> PODVolOctree::mapCellsToPOD(const std::unordered_set<long> * targetCells)
1144{
1145 std::unordered_set<long> mappedCells;
1146 const PiercedStorage<mapping::Info> &inverseMapping = getMapper()->getInverseMapping();
1147
1148#if BITPIT_ENABLE_MPI
1149 if ( getMapper()->checkPartition()){
1150#endif
1151
1152 for (const long & id : *targetCells){
1153 if (inverseMapping[id].type == mapping::Type::TYPE_RENUMBERING){
1154 mappedCells.insert(inverseMapping[id].ids[0]);
1155 }
1156 if (inverseMapping[id].type == mapping::Type::TYPE_REFINEMENT){
1157 mappedCells.insert(inverseMapping[id].ids[0]);
1158 }
1159 if (inverseMapping[id].type == mapping::Type::TYPE_COARSENING){
1160 for (long idd : inverseMapping[id].ids)
1161 mappedCells.insert(idd);
1162 }
1163 }//end targetcells
1164
1165#if BITPIT_ENABLE_MPI
1166 }
1167 else{
1168
1169 //Mapping between partitioned meshes
1170
1171 //Inset local POD cells and recover POD cells on other partitions by inverse mapping
1172 std::map<int, std::unordered_set<long> > sendPODcells;
1173
1174 for (const long & id : *targetCells){
1175 if (inverseMapping[id].type == mapping::Type::TYPE_RENUMBERING){
1176 int rank = inverseMapping[id].ranks[0];
1177 if (rank != m_rank){
1178 sendPODcells[rank].insert(inverseMapping[id].ids[0]);
1179 }
1180 else{
1181 mappedCells.insert(inverseMapping[id].ids[0]);
1182 }
1183 }
1184 if (inverseMapping[id].type == mapping::Type::TYPE_REFINEMENT){
1185 int rank = inverseMapping[id].ranks[0];
1186 if (rank != m_rank){
1187 sendPODcells[rank].insert(inverseMapping[id].ids[0]);
1188 }
1189 else{
1190 mappedCells.insert(inverseMapping[id].ids[0]);
1191 }
1192 }
1193 if (inverseMapping[id].type == mapping::Type::TYPE_COARSENING){
1194 int ii=0;
1195 for (long idd : inverseMapping[id].ids){
1196 int rank = inverseMapping[id].ranks[ii];
1197 if (rank != m_rank){
1198 sendPODcells[rank].insert(idd);
1199 }
1200 else{
1201 mappedCells.insert(idd);
1202 }
1203 ii++;
1204 }
1205 }
1206 }//end targetcells
1207
1208
1209 //Communicate cells to send
1210 //build send buffers
1211 DataCommunicator dataCommunicator(m_communicator);
1212 std::size_t bytes = sizeof(long);
1213 for (const auto &val : sendPODcells){
1214 std::size_t ncells = val.second.size();
1215 int rank = val.first;
1216 //set size
1217 std::size_t buffSize = ncells * bytes;
1218 dataCommunicator.setSend(rank,buffSize);
1219 //fill buffer with octants
1220 SendBuffer &sendBuffer = dataCommunicator.getSendBuffer(rank);
1221 for (long ID : val.second){
1222 sendBuffer << ID;
1223 }
1224 }
1225
1226 dataCommunicator.discoverRecvs();
1227 dataCommunicator.startAllRecvs();
1228 dataCommunicator.startAllSends();
1229
1230 std::vector<int> recvRanks = dataCommunicator.getRecvRanks();
1231 std::sort(recvRanks.begin(),recvRanks.end());
1232
1233 for (int rank : recvRanks){
1234 dataCommunicator.waitRecv(rank);
1235 RecvBuffer & recvBuffer = dataCommunicator.getRecvBuffer(rank);
1236 long nof = recvBuffer.getSize() / bytes;
1237 for (long ii = 0; ii < nof; ii++){
1238 long ID;
1239 recvBuffer >> ID;
1240 mappedCells.insert(ID);
1241 }
1242 }
1243
1244 dataCommunicator.waitAllSends();
1245
1246 }//endif partitioned
1247#endif
1248
1249 return mappedCells;
1250
1251}
1252
1253void PODVolOctree::adaptMeshToMesh(const VolumeKernel* meshToAdapt, const VolumeKernel * meshReference)
1254{
1255 BITPIT_UNUSED(meshToAdapt);
1256
1257 computeMapper(meshReference, false);
1258
1259 bool adapt = true;
1260
1261 while(adapt){
1262
1263 adapt = false;
1264
1265 const PiercedStorage<mapping::Info> &mappingInfo = getMapper()->getMapping();
1266
1267 for (Cell & cell : getMesh()->getCells()){
1268 long id = cell.getId();
1269 if (mappingInfo[id].type == mapping::Type::TYPE_COARSENING){
1270 getMesh()->markCellForRefinement(id);
1271 adapt = true;
1272 }
1273 }
1274
1275#if BITPIT_ENABLE_MPI==1
1276 MPI_Allreduce(MPI_IN_PLACE, &adapt, 1, MPI_C_BOOL, MPI_LOR, m_communicator);
1277#endif
1278
1279 std::vector<adaption::Info> infoAdapt = getMesh()->adaptionPrepare(true);
1280
1281 getMapper()->adaptionPrepare(infoAdapt);
1282
1283 infoAdapt = getMesh()->adaptionAlter(true);
1284
1285 getMapper()->adaptionAlter(infoAdapt, true, false);
1286
1287 getMesh()->adaptionCleanup();
1288
1289 getMapper()->adaptionCleanup();
1290
1291 }
1292
1293}
1294
1295# if BITPIT_ENABLE_MPI
1296
1297void PODVolOctree::communicatePODField(const pod::PODField & field, std::map<int, std::map<long, bool> > & dataBrec, std::map<int, std::map<long, std::vector<double> > > & dataSrec, std::map<int, std::map<long, std::vector<std::array<double,3> > > > & dataVrec, std::map<int, std::map<long, double> > & volrec)
1298{
1299
1300 {
1301 // Size of fields
1302 std::size_t nsf = field.scalar->getFieldCount();
1303 std::size_t nvf = field.vector->getFieldCount();
1304
1305 std::unordered_map<int, std::vector<long> > rankIDrec = m_meshmap->getReceivedMappedIds();
1306 std::unordered_map<int, std::vector<long> > rankIDsend = m_meshmap->getSentMappedIds();
1307
1308 //build send buffers
1309 DataCommunicator dataCommunicator(m_communicator);
1310 std::size_t bytes = sizeof(bool) + (nsf+(3*nvf)+1)*sizeof(double);
1311 for (const auto &val : rankIDsend){
1312 int rank = val.first;
1313 //set size
1314 std::size_t buffSize = val.second.size() * bytes;
1315 dataCommunicator.setSend(rank,buffSize);
1316 //fill buffer with octants
1317 SendBuffer &sendBuffer = dataCommunicator.getSendBuffer(rank);
1318 for (long ID : val.second){
1319 bool mask = field.mask->at(ID);
1320 sendBuffer << mask;
1321 for (std::size_t i=0; i<nsf; i++)
1322 sendBuffer << field.scalar->at(ID,i);
1323 for (std::size_t i=0; i<nvf; i++){
1324 for (std::size_t j=0; j<3; j++){
1325 sendBuffer << field.vector->at(ID,i)[j];
1326 }
1327 }
1328 sendBuffer << field.mesh->evalCellVolume(ID);
1329 }
1330 }
1331
1332 dataCommunicator.discoverRecvs();
1333 dataCommunicator.startAllRecvs();
1334 dataCommunicator.startAllSends();
1335
1336 for (const auto &val : rankIDrec){
1337 int rank = val.first;
1338 dataCommunicator.waitRecv(rank);
1339 RecvBuffer & recvBuffer = dataCommunicator.getRecvBuffer(rank);
1340 for (long ID : val.second){
1341 recvBuffer >> dataBrec[rank][ID];
1342 for (std::size_t i=0; i<nsf; i++){
1343 double data;
1344 recvBuffer >> data;
1345 dataSrec[rank][ID].push_back(data);
1346 }
1347 for (std::size_t i=0; i<nvf; i++){
1348 std::array<double,3> data;
1349 recvBuffer >> data[0];
1350 recvBuffer >> data[1];
1351 recvBuffer >> data[2];
1352 dataVrec[rank][ID].push_back(data);
1353 }
1354 recvBuffer >> volrec[rank][ID];
1355 }
1356 }
1357 dataCommunicator.waitAllSends();
1358
1359 }
1360
1361}
1362
1363void PODVolOctree::communicatePODFieldFromPOD(const pod::PODField & field, std::map<int, std::map<long, bool> > & dataBrec, std::map<int, std::map<long, std::vector<double> > > & dataSrec, std::map<int, std::map<long, std::vector<std::array<double,3> > > > & dataVrec, std::map<int, std::map<long, double> > & volrec)
1364{
1365
1366 {
1367 // Size of fields
1368 std::size_t nsf = field.scalar->getFieldCount();
1369 std::size_t nvf = field.vector->getFieldCount();
1370
1371 std::unordered_map<int, std::vector<long> > rankIDsend = m_meshmap->getSentReferenceIds();
1372
1373 //build send buffers
1374 DataCommunicator dataCommunicator(m_communicator);
1375 std::size_t bytes = sizeof(long) + sizeof(bool) + (nsf+3*nvf+1)*sizeof(double);
1376 for (const auto &val : rankIDsend){
1377 int rank = val.first;
1378 //set size
1379 std::size_t buffSize = val.second.size() * bytes;
1380 dataCommunicator.setSend(rank,buffSize);
1381 //fill buffer with octants
1382 SendBuffer &sendBuffer = dataCommunicator.getSendBuffer(rank);
1383 for (long ID : val.second){
1384 sendBuffer << ID;
1385 bool mask = field.mask->at(ID);
1386 sendBuffer << mask;
1387 for (std::size_t i=0; i<nsf; i++)
1388 sendBuffer << field.scalar->at(ID,i);
1389 for (std::size_t i=0; i<nvf; i++){
1390 for (std::size_t j=0; j<3; j++){
1391 sendBuffer << field.vector->at(ID,i)[j];
1392 }
1393 }
1394 sendBuffer << field.mesh->evalCellVolume(ID);
1395 }
1396 }
1397
1398 dataCommunicator.discoverRecvs();
1399 dataCommunicator.startAllRecvs();
1400 dataCommunicator.startAllSends();
1401
1402 std::vector<int> recvRanks = dataCommunicator.getRecvRanks();
1403 std::sort(recvRanks.begin(),recvRanks.end());
1404
1405 for (int rank : recvRanks){
1406 dataCommunicator.waitRecv(rank);
1407 RecvBuffer & recvBuffer = dataCommunicator.getRecvBuffer(rank);
1408 long nof = recvBuffer.getSize() / bytes;
1409 for (long ii = 0; ii < nof; ii++){
1410 long ID;
1411 recvBuffer >> ID;
1412 recvBuffer >> dataBrec[rank][ID];
1413 for (std::size_t i=0; i<nsf; i++){
1414 double data;
1415 recvBuffer >> data;
1416 dataSrec[rank][ID].push_back(data);
1417 }
1418 for (std::size_t i=0; i<nvf; i++){
1419 std::array<double,3> data;
1420 recvBuffer >> data[0];
1421 recvBuffer >> data[1];
1422 recvBuffer >> data[2];
1423 dataVrec[rank][ID].push_back(data);
1424 }
1425 recvBuffer >> volrec[rank][ID];
1426 }
1427 }
1428
1429 dataCommunicator.waitAllSends();
1430
1431 }
1432
1433}
1434
1435void PODVolOctree::communicateBoolField(const PiercedStorage<bool> & field, std::map<int, std::map<long, bool> > & dataBrec)
1436{
1437
1438 std::unordered_map<int, std::vector<long> > rankIDrec = m_meshmap->getReceivedMappedIds();
1439 std::unordered_map<int, std::vector<long> > rankIDsend = m_meshmap->getSentMappedIds();
1440
1441 //build send buffers
1442 DataCommunicator dataCommunicator(m_communicator);
1443 std::size_t bytes = sizeof(bool);
1444 for (const auto &val : rankIDsend){
1445 int rank = val.first;
1446 //set size
1447 std::size_t buffSize = val.second.size() * bytes;
1448 dataCommunicator.setSend(rank,buffSize);
1449 //fill buffer with octants
1450 SendBuffer &sendBuffer = dataCommunicator.getSendBuffer(rank);
1451 for (long ID : val.second){
1452 sendBuffer << field.at(ID);
1453 }
1454 }
1455
1456 dataCommunicator.discoverRecvs();
1457 dataCommunicator.startAllRecvs();
1458 dataCommunicator.startAllSends();
1459
1460 for (const auto &val : rankIDrec){
1461 int rank = val.first;
1462 dataCommunicator.waitRecv(rank);
1463 RecvBuffer & recvBuffer = dataCommunicator.getRecvBuffer(rank);
1464 for (long ID : val.second){
1465 recvBuffer >> dataBrec[rank][ID];
1466 }
1467 }
1468
1469 dataCommunicator.waitAllSends();
1470
1471}
1472
1473void PODVolOctree::communicateField(const PiercedStorage<double> & field, const VolumeKernel * mesh, std::map<int, std::map<long, std::vector<double> > > & datarec, std::map<int, std::map<long, double> > & volrec)
1474{
1475
1476 std::size_t nf = field.getFieldCount();
1477
1478 std::unordered_map<int, std::vector<long> > rankIDrec = m_meshmap->getReceivedMappedIds();
1479 std::unordered_map<int, std::vector<long> > rankIDsend = m_meshmap->getSentMappedIds();
1480
1481 //build send buffers
1482 DataCommunicator dataCommunicator(m_communicator);
1483 std::size_t bytes = (nf+1)*sizeof(double);
1484 for (const auto &val : rankIDsend){
1485 int rank = val.first;
1486 //set size
1487 std::size_t buffSize = val.second.size() * bytes;
1488 dataCommunicator.setSend(rank,buffSize);
1489 //fill buffer with octants
1490 SendBuffer &sendBuffer = dataCommunicator.getSendBuffer(rank);
1491 for (long ID : val.second){
1492 for (std::size_t ifield=0; ifield<nf; ifield++)
1493 sendBuffer << field.at(ID, ifield);
1494 sendBuffer << mesh->evalCellVolume(ID);
1495 }
1496 }
1497
1498 dataCommunicator.discoverRecvs();
1499 dataCommunicator.startAllRecvs();
1500 dataCommunicator.startAllSends();
1501
1502 for (const auto &val : rankIDrec){
1503 int rank = val.first;
1504 dataCommunicator.waitRecv(rank);
1505 RecvBuffer & recvBuffer = dataCommunicator.getRecvBuffer(rank);
1506 for (long ID : val.second){
1507 datarec[rank][ID].resize(nf);
1508 for (std::size_t ifield=0; ifield<nf; ifield++)
1509 recvBuffer >> datarec[rank][ID][ifield];
1510 recvBuffer >> volrec[rank][ID];
1511 }
1512 }
1513
1514 dataCommunicator.waitAllSends();
1515
1516}
1517
1518void PODVolOctree::communicateFieldFromPOD(const PiercedStorage<double> & field, const VolumeKernel * mesh, std::map<int, std::map<long, std::vector<double> > > & datarec, std::map<int, std::map<long, double> > & volrec)
1519{
1520
1521 std::size_t nf = field.getFieldCount();
1522
1523 std::unordered_map<int, std::vector<long> > rankIDsend = m_meshmap->getSentReferenceIds();
1524
1525 //build send buffers
1526 DataCommunicator dataCommunicator(m_communicator);
1527 std::size_t bytes = (nf+1)*sizeof(double)+sizeof(long);
1528 for (const auto &val : rankIDsend){
1529 int rank = val.first;
1530 //set size
1531 std::size_t buffSize = val.second.size() * bytes;
1532 dataCommunicator.setSend(rank,buffSize);
1533 //fill buffer with octants
1534 SendBuffer &sendBuffer = dataCommunicator.getSendBuffer(rank);
1535 for (long ID : val.second){
1536 sendBuffer << ID;
1537 for (std::size_t ifield=0; ifield<nf; ifield++)
1538 sendBuffer << field.at(ID, ifield);
1539 sendBuffer << mesh->evalCellVolume(ID);
1540 }
1541 }
1542
1543 dataCommunicator.discoverRecvs();
1544 dataCommunicator.startAllRecvs();
1545 dataCommunicator.startAllSends();
1546
1547 std::vector<int> recvRanks = dataCommunicator.getRecvRanks();
1548 std::sort(recvRanks.begin(),recvRanks.end());
1549
1550 for (int rank : recvRanks){
1551 dataCommunicator.waitRecv(rank);
1552 RecvBuffer & recvBuffer = dataCommunicator.getRecvBuffer(rank);
1553 long nof = recvBuffer.getSize() / bytes;
1554 for (long ii = 0; ii < nof; ii++){
1555 long ID;
1556 recvBuffer >> ID;
1557 datarec[rank][ID].resize(nf);
1558 for (std::size_t ifield=0; ifield<nf; ifield++)
1559 recvBuffer >> datarec[rank][ID][ifield];
1560 recvBuffer >> volrec[rank][ID];
1561 }
1562 }
1563
1564 dataCommunicator.waitAllSends();
1565
1566}
1567
1568#endif
1569
1570}
bitpit::PODKernel
The PODKernel class provides an interface to manage the mesh dependent members and functions of a POD...
Definition pod_kernel.hpp:41
bitpit::PODKernel::m_meshmap
std::unique_ptr< VolumeMapper > m_meshmap
Definition pod_kernel.hpp:88
bitpit::PODKernel::computeMapper
void computeMapper(const VolumeKernel *mesh, bool fillInv=true)
Definition pod_kernel.cpp:194
bitpit::PODKernel::m_communicator
MPI_Comm m_communicator
Definition pod_kernel.hpp:83
bitpit::PODKernel::getMapper
VolumeMapper * getMapper()
Definition pod_kernel.cpp:246
bitpit::PODKernel::getMesh
VolumeKernel * getMesh()
Definition pod_kernel.cpp:120
bitpit::PODKernel::getCellVolume
double getCellVolume(long id)
Definition pod_kernel.cpp:174
bitpit::PODKernel::getRawCellVolume
double getRawCellVolume(long rawIndex)
Definition pod_kernel.cpp:184
bitpit::PODKernel::m_meshPOD
std::unique_ptr< VolumeKernel > m_meshPOD
Definition pod_kernel.hpp:79
bitpit::PODKernel::getCommunicator
MPI_Comm getCommunicator() const
Definition pod_kernel.cpp:307
bitpit::PODVolOctree
The PODVolOctree is the specialized class of PODKernel for VolOctree meshes.
Definition pod_voloctree.hpp:40
bitpit::PODVolOctree::PODVolOctree
PODVolOctree(MPI_Comm comm=MPI_COMM_WORLD)
Definition pod_voloctree.cpp:45
bitpit::PatchKernel::markCellForRefinement
void markCellForRefinement(long id)
Definition patch_kernel.cpp:890
bitpit::PatchKernel::adaptionPrepare
std::vector< adaption::Info > adaptionPrepare(bool trackAdaption=true)
Definition patch_kernel.cpp:719
bitpit::PatchKernel::adaptionAlter
std::vector< adaption::Info > adaptionAlter(bool trackAdaption=true, bool squeezeStorage=false)
Definition patch_kernel.cpp:760
bitpit::VolOctree
The VolOctree defines a Octree patch.
Definition voloctree.hpp:37
bitpit::VolumeMapper::getInverseMapping
const bitpit::PiercedStorage< mapping::Info > & getInverseMapping() const
Definition volume_mapper.cpp:151
bitpit::VolumeMapper::getMapping
const bitpit::PiercedStorage< mapping::Info > & getMapping() const
Definition volume_mapper.cpp:143
BITPIT_UNUSED
#define BITPIT_UNUSED(variable)
Definition compiler.hpp:63
--- layout: doxygen_footer ---
Generated on Wed Jun 19 2024 10:21:05 for bitpit by   doxygen 1.10.0