lduPrimitiveMesh.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2013-2017 OpenFOAM Foundation
9 Copyright (C) 2019 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
29#include "lduPrimitiveMesh.H"
31#include "edgeHashes.H"
32#include "labelPair.H"
34
35// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36
37namespace Foam
38{
40
41 //- Less operator for pairs of <processor><index>
43 {
44 const labelPairList& list_;
45
46 public:
47
49 :
50 list_(list)
51 {}
52
53 bool operator()(const label a, const label b)
54 {
55 return list_[a].first() < list_[b].first();
56 }
57 };
58
59}
60
61
62// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
63
65(
66 const label size,
67 const labelUList& l,
68 const labelUList& u
69)
70{
71 forAll(l, facei)
72 {
73 if (u[facei] < l[facei])
74 {
76 << "Reversed face. Problem at face " << facei
77 << " l:" << l[facei] << " u:" << u[facei]
78 << abort(FatalError);
79 }
80 if (l[facei] < 0 || u[facei] < 0 || u[facei] >= size)
81 {
83 << "Illegal cell label. Problem at face " << facei
84 << " l:" << l[facei] << " u:" << u[facei]
85 << abort(FatalError);
86 }
87 }
88
89 for (label facei=1; facei < l.size(); ++facei)
90 {
91 if (l[facei-1] > l[facei])
92 {
94 << "Lower not in incremental cell order."
95 << " Problem at face " << facei
96 << " l:" << l[facei] << " u:" << u[facei]
97 << " previous l:" << l[facei-1]
98 << abort(FatalError);
99 }
100 else if (l[facei-1] == l[facei])
101 {
102 // Same cell.
103 if (u[facei-1] > u[facei])
104 {
106 << "Upper not in incremental cell order."
107 << " Problem at face " << facei
108 << " l:" << l[facei] << " u:" << u[facei]
109 << " previous u:" << u[facei-1]
110 << abort(FatalError);
111 }
112 }
113 }
114}
115
116
118(
119 const lduMesh& mesh,
120 const globalIndex& globalNumbering
121)
122{
123 const lduAddressing& addr = mesh.lduAddr();
124 lduInterfacePtrsList interfaces = mesh.interfaces();
125
126 const labelList globalIndices
127 (
129 (
130 addr.size(),
131 globalNumbering.localStart(UPstream::myProcNo(mesh.comm()))
132 )
133 );
134
135 // Get the interface cells
136 PtrList<labelList> nbrGlobalCells(interfaces.size());
137 {
138 // Initialise transfer of restrict addressing on the interface
139 forAll(interfaces, inti)
140 {
141 if (interfaces.set(inti))
142 {
143 interfaces[inti].initInternalFieldTransfer
144 (
146 globalIndices
147 );
148 }
149 }
150
151 if (Pstream::parRun())
152 {
154 }
155
156 forAll(interfaces, inti)
157 {
158 if (interfaces.set(inti))
159 {
160 nbrGlobalCells.set
161 (
162 inti,
163 new labelList
164 (
165 interfaces[inti].internalFieldTransfer
166 (
168 globalIndices
169 )
170 )
171 );
172 }
173 }
174 }
175
176
177 // Scan the neighbour list to find out how many times the cell
178 // appears as a neighbour of the face. Done this way to avoid guessing
179 // and resizing list
180 labelList nNbrs(addr.size(), Zero);
181
182 const labelUList& nbr = addr.upperAddr();
183 const labelUList& own = addr.lowerAddr();
184
185 {
186 forAll(nbr, facei)
187 {
188 nNbrs[nbr[facei]]++;
189 nNbrs[own[facei]]++;
190 }
191
192 forAll(interfaces, inti)
193 {
194 if (interfaces.set(inti))
195 {
196 for (const label celli : interfaces[inti].faceCells())
197 {
198 nNbrs[celli]++;
199 }
200 }
201 }
202 }
203
204
205 // Create cell-cells addressing
206 labelListList cellCells(addr.size());
207
208 forAll(cellCells, celli)
209 {
210 cellCells[celli].setSize(nNbrs[celli], -1);
211 }
212
213 // Reset the list of number of neighbours to zero
214 nNbrs = 0;
215
216 // Scatter the neighbour faces
217 forAll(nbr, facei)
218 {
219 const label c0 = own[facei];
220 const label c1 = nbr[facei];
221
222 cellCells[c0][nNbrs[c0]++] = globalIndices[c1];
223 cellCells[c1][nNbrs[c1]++] = globalIndices[c0];
224 }
225 forAll(interfaces, inti)
226 {
227 if (interfaces.set(inti))
228 {
229 const labelUList& faceCells = interfaces[inti].faceCells();
230
231 forAll(faceCells, facei)
232 {
233 const label c0 = faceCells[facei];
234 cellCells[c0][nNbrs[c0]++] = nbrGlobalCells[inti][facei];
235 }
236 }
237 }
238
239 return cellCells;
240}
241
242
244(
246)
247{
248 label size = 0;
249
250 for (const lduPrimitiveMesh& msh : meshes)
251 {
252 size += msh.lduAddr().size();
253 }
254 return size;
255}
256
257
259(
260 const label nCells,
261 const labelUList& lower,
262 const labelUList& upper
263)
264{
265 labelList nNbrs(nCells, Zero);
266
267 // Count number of upper neighbours
268 forAll(lower, facei)
269 {
270 if (upper[facei] < lower[facei])
271 {
273 << "Problem at face:" << facei
274 << " lower:" << lower[facei]
275 << " upper:" << upper[facei]
276 << exit(FatalError);
277 }
278 nNbrs[lower[facei]]++;
279 }
280
281 // Construct cell-upper cell addressing
282 labelList offsets(nCells+1);
283 offsets[0] = 0;
284 forAll(nNbrs, celli)
285 {
286 offsets[celli+1] = offsets[celli]+nNbrs[celli];
287 }
288
289 nNbrs = offsets;
290
291 labelList cellToFaces(offsets.last());
292 forAll(upper, facei)
293 {
294 const label celli = lower[facei];
295 cellToFaces[nNbrs[celli]++] = facei;
296 }
297
298 // Sort
299
300 labelList oldToNew(lower.size());
301
302 DynamicList<label> order;
304
305 label newFacei = 0;
306
307 for (label celli = 0; celli < nCells; ++celli)
308 {
309 const label startOfCell = offsets[celli];
310 const label nNbr = offsets[celli+1] - startOfCell;
311
312 nbr.resize(nNbr);
313 order.resize(nNbr);
314
315 forAll(nbr, i)
316 {
317 nbr[i] = upper[cellToFaces[offsets[celli]+i]];
318 }
319 sortedOrder(nbr, order);
320
321 for (const label index : order)
322 {
323 oldToNew[cellToFaces[startOfCell + index]] = newFacei++;
324 }
325 }
326
327 return oldToNew;
328}
329
330
331// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
332
334(
335 const label nCells,
336 labelList& l,
337 labelList& u,
338 const label comm,
339 bool reuse
340)
341:
342 lduAddressing(nCells),
343 lowerAddr_(l, reuse),
344 upperAddr_(u, reuse),
345 comm_(comm)
346{}
347
349(
350 lduInterfacePtrsList& interfaces,
351 const lduSchedule& ps
352)
353{
354 interfaces_ = interfaces;
355 patchSchedule_ = ps;
356
357 // Create interfaces
358 primitiveInterfaces_.setSize(interfaces_.size());
359 forAll(interfaces_, i)
360 {
361 if (interfaces_.set(i))
362 {
363 primitiveInterfaces_.set(i, &interfaces_[i]);
364 }
365 }
366}
367
368
370(
371 const label nCells
372)
373:
374 lduAddressing(nCells),
375 comm_(0)
376{}
377
378
380(
381 const label nCells,
382 labelList& l,
383 labelList& u,
384 PtrList<const lduInterface>& primitiveInterfaces,
385 const lduSchedule& ps,
386 const label comm
387)
388:
389 lduAddressing(nCells),
390 lowerAddr_(l, true),
391 upperAddr_(u, true),
392 primitiveInterfaces_(),
393 patchSchedule_(ps),
394 comm_(comm)
395{
396 primitiveInterfaces_.transfer(primitiveInterfaces);
397
398 // Create interfaces
399 interfaces_.setSize(primitiveInterfaces_.size());
400 forAll(primitiveInterfaces_, i)
401 {
402 if (primitiveInterfaces_.set(i))
403 {
404 interfaces_.set(i, &primitiveInterfaces_[i]);
405 }
406 }
407}
408
409
411(
412 const label comm,
413 const labelList& procAgglomMap,
414
415 const labelList& procIDs,
416 const lduMesh& myMesh,
417 const PtrList<lduPrimitiveMesh>& otherMeshes,
418
419 labelList& cellOffsets,
420 labelList& faceOffsets,
422 labelListList& boundaryMap,
423 labelListListList& boundaryFaceMap
424)
425:
426 lduAddressing(myMesh.lduAddr().size() + totalSize(otherMeshes)),
427 comm_(comm)
428{
429 const label currentComm = myMesh.comm();
430
431 forAll(otherMeshes, i)
432 {
433 if (otherMeshes[i].comm() != currentComm)
434 {
436 << "Communicator " << otherMeshes[i].comm()
437 << " at index " << i
438 << " differs from that of predecessor "
439 << currentComm
440 << endl; //exit(FatalError);
441 }
442 }
443
444 const label nMeshes = otherMeshes.size()+1;
445
446 const label myAgglom = procAgglomMap[UPstream::myProcNo(currentComm)];
447
448 if (lduPrimitiveMesh::debug)
449 {
450 Pout<< "I am " << UPstream::myProcNo(currentComm)
451 << " agglomerating into " << myAgglom
452 << " as are " << findIndices(procAgglomMap, myAgglom)
453 << endl;
454 }
455
456
457 forAll(procIDs, i)
458 {
459 if (procAgglomMap[procIDs[i]] != procAgglomMap[procIDs[0]])
460 {
462 << "Processor " << procIDs[i]
463 << " agglomerates to " << procAgglomMap[procIDs[i]]
464 << " whereas other processors " << procIDs
465 << " agglomerate to "
466 << labelUIndList(procAgglomMap, procIDs)
467 << exit(FatalError);
468 }
469 }
470
471
472 // Cells get added in order.
473 cellOffsets.setSize(nMeshes+1);
474 cellOffsets[0] = 0;
475 for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
476 {
477 const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
478
479 cellOffsets[procMeshI+1] =
480 cellOffsets[procMeshI]
481 + procMesh.lduAddr().size();
482 }
483
484
485 // Faces initially get added in order, sorted later
486 labelList internalFaceOffsets(nMeshes+1);
487 internalFaceOffsets[0] = 0;
488 for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
489 {
490 const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
491
492 internalFaceOffsets[procMeshI+1] =
493 internalFaceOffsets[procMeshI]
494 + procMesh.lduAddr().lowerAddr().size();
495 }
496
497 // Count how faces get added. Interfaces inbetween get merged.
498
499 // Merged interfaces: map from two coarse processors back to
500 // - procMeshes
501 // - interface in procMesh
502 // (estimate size from number of patches of mesh0)
503 EdgeMap<labelPairList> mergedMap(2*myMesh.interfaces().size());
504
505 // Unmerged interfaces: map from two coarse processors back to
506 // - procMeshes
507 // - interface in procMesh
508 EdgeMap<labelPairList> unmergedMap(mergedMap.size());
509
510 boundaryMap.setSize(nMeshes);
511 boundaryFaceMap.setSize(nMeshes);
512
513
514 label nOtherInterfaces = 0;
515 labelList nCoupledFaces(nMeshes, Zero);
516
517 for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
518 {
520 mesh(myMesh, otherMeshes, procMeshI).interfaces();
521
522 // Initialise all boundaries as merged
523 boundaryMap[procMeshI].setSize(interfaces.size(), -1);
524 boundaryFaceMap[procMeshI].setSize(interfaces.size());
525
526 // Get sorted order of processors
527 forAll(interfaces, intI)
528 {
529 if (interfaces.set(intI))
530 {
531 const lduInterface& ldui = interfaces[intI];
532
533 if (isA<processorLduInterface>(ldui))
534 {
535 const processorLduInterface& pldui =
536 refCast<const processorLduInterface>(ldui);
537
538 label agglom0 = procAgglomMap[pldui.myProcNo()];
539 label agglom1 = procAgglomMap[pldui.neighbProcNo()];
540
541 const edge procEdge(agglom0, agglom1);
542
543 if (agglom0 != myAgglom && agglom1 != myAgglom)
544 {
546 << "At mesh from processor " << procIDs[procMeshI]
547 << " have interface " << intI
548 << " with myProcNo:" << pldui.myProcNo()
549 << " with neighbProcNo:" << pldui.neighbProcNo()
550 << exit(FatalError);
551 }
552 else if (agglom0 == myAgglom && agglom1 == myAgglom)
553 {
554 // Merged interface
555 if (debug)
556 {
557 Pout<< "merged interface: myProcNo:"
558 << pldui.myProcNo()
559 << " nbr:" << pldui.neighbProcNo()
560 << " size:" << ldui.faceCells().size()
561 << endl;
562 }
563
564 const label nbrProcMeshI =
565 procIDs.find(pldui.neighbProcNo());
566
567 if (procMeshI < nbrProcMeshI)
568 {
569 // I am 'master' since get lowest numbered cells
570 nCoupledFaces[procMeshI] += ldui.faceCells().size();
571 }
572
573 mergedMap(procEdge).append
574 (
575 labelPair(procMeshI, intI)
576 );
577 }
578 else
579 {
580 if (debug)
581 {
582 Pout<< "external interface: myProcNo:"
583 << pldui.myProcNo()
584 << " nbr:" << pldui.neighbProcNo()
585 << " size:" << ldui.faceCells().size()
586 << endl;
587 }
588
589 unmergedMap(procEdge).append
590 (
591 labelPair(procMeshI, intI)
592 );
593 }
594 }
595 else
596 {
597 // Still external (non proc) interface
599 << "At mesh from processor " << procIDs[procMeshI]
600 << " have interface " << intI
601 << " of unhandled type " << interfaces[intI].type()
602 << exit(FatalError);
603
604 ++nOtherInterfaces;
605 }
606 }
607 }
608 }
609
610
611
612 if (debug)
613 {
614 Pout<< "Remaining interfaces:" << endl;
615 forAllConstIters(unmergedMap, iter)
616 {
617 Pout<< " agglom procEdge:" << iter.key() << endl;
618 const labelPairList& elems = iter.val();
619 forAll(elems, i)
620 {
621 label procMeshI = elems[i][0];
622 label interfacei = elems[i][1];
624 mesh(myMesh, otherMeshes, procMeshI).interfaces();
625
626 const processorLduInterface& pldui =
627 refCast<const processorLduInterface>
628 (
629 interfaces[interfacei]
630 );
631
632 Pout<< " proc:" << procIDs[procMeshI]
633 << " interfacei:" << interfacei
634 << " between:" << pldui.myProcNo()
635 << " and:" << pldui.neighbProcNo()
636 << endl;
637 }
638 Pout<< endl;
639 }
640 }
641 if (debug)
642 {
643 Pout<< "Merged interfaces:" << endl;
644 forAllConstIters(mergedMap, iter)
645 {
646 Pout<< " agglom procEdge:" << iter.key() << endl;
647 const labelPairList& elems = iter.val();
648
649 forAll(elems, i)
650 {
651 label procMeshI = elems[i][0];
652 label interfacei = elems[i][1];
654 mesh(myMesh, otherMeshes, procMeshI).interfaces();
655 const processorLduInterface& pldui =
656 refCast<const processorLduInterface>
657 (
658 interfaces[interfacei]
659 );
660
661 Pout<< " proc:" << procIDs[procMeshI]
662 << " interfacei:" << interfacei
663 << " between:" << pldui.myProcNo()
664 << " and:" << pldui.neighbProcNo()
665 << endl;
666 }
667 Pout<< endl;
668 }
669 }
670
671
672 // Adapt faceOffsets for internal interfaces
673 faceOffsets.setSize(nMeshes+1);
674 faceOffsets[0] = 0;
675 faceMap.setSize(nMeshes);
676 for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
677 {
678 const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
679 label nInternal = procMesh.lduAddr().lowerAddr().size();
680
681 faceOffsets[procMeshI+1] =
682 faceOffsets[procMeshI]
683 + nInternal
684 + nCoupledFaces[procMeshI];
685
686 labelList& map = faceMap[procMeshI];
687 map.setSize(nInternal);
688 forAll(map, i)
689 {
690 map[i] = faceOffsets[procMeshI] + i;
691 }
692 }
693
694
695 // Combine upper and lower
696 lowerAddr_.setSize(faceOffsets.last(), -1);
697 upperAddr_.setSize(lowerAddr_.size(), -1);
698
699
700 // Old internal faces and resolved coupled interfaces
701 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
702
703 for (label procMeshI = 0; procMeshI < nMeshes; ++procMeshI)
704 {
705 const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
706
707 const labelUList& l = procMesh.lduAddr().lowerAddr();
708 const labelUList& u = procMesh.lduAddr().upperAddr();
709
710 // Add internal faces
711 label allFacei = faceOffsets[procMeshI];
712
713 forAll(l, facei)
714 {
715 lowerAddr_[allFacei] = cellOffsets[procMeshI]+l[facei];
716 upperAddr_[allFacei] = cellOffsets[procMeshI]+u[facei];
717 allFacei++;
718 }
719
720
721 // Add merged interfaces
722 const lduInterfacePtrsList interfaces = procMesh.interfaces();
723
724 forAll(interfaces, intI)
725 {
726 if (interfaces.set(intI))
727 {
728 const lduInterface& ldui = interfaces[intI];
729
730 if (isA<processorLduInterface>(ldui))
731 {
732 const processorLduInterface& pldui =
733 refCast<const processorLduInterface>(ldui);
734
735 // Look up corresponding interfaces
736 label myP = pldui.myProcNo();
737 label nbrP = pldui.neighbProcNo();
738 label nbrProcMeshI = procIDs.find(nbrP);
739
740 if (procMeshI < nbrProcMeshI)
741 {
742 // I am 'master' since my cell numbers will be lower
743 // since cells get added in procMeshI order.
744
745 const label agglom0 = procAgglomMap[myP];
746 const label agglom1 = procAgglomMap[nbrP];
747
748 const auto fnd =
749 mergedMap.cfind(edge(agglom0, agglom1));
750
751 if (fnd.found())
752 {
753 const labelPairList& elems = fnd();
754
755 // Find nbrP in elems
756 label nbrIntI = -1;
757 forAll(elems, i)
758 {
759 label proci = elems[i][0];
760 label interfacei = elems[i][1];
762 mesh
763 (
764 myMesh,
765 otherMeshes,
766 proci
767 ).interfaces();
768 const processorLduInterface& pldui =
769 refCast<const processorLduInterface>
770 (
771 interfaces[interfacei]
772 );
773
774 if
775 (
776 elems[i][0] == nbrProcMeshI
777 && pldui.neighbProcNo() == procIDs[procMeshI]
778 )
779 {
780 nbrIntI = elems[i][1];
781 break;
782 }
783 }
784
785
786 if (nbrIntI == -1)
787 {
789 << "elems:" << elems << abort(FatalError);
790 }
791
792
793 const lduInterfacePtrsList nbrInterfaces = mesh
794 (
795 myMesh,
796 otherMeshes,
797 nbrProcMeshI
798 ).interfaces();
799
800
801 const labelUList& faceCells =
802 ldui.faceCells();
803 const labelUList& nbrFaceCells =
804 nbrInterfaces[nbrIntI].faceCells();
805
806 if (faceCells.size() != nbrFaceCells.size())
807 {
809 << "faceCells:" << faceCells
810 << " nbrFaceCells:" << nbrFaceCells
811 << abort(FatalError);
812 }
813
814
815 labelList& bfMap =
816 boundaryFaceMap[procMeshI][intI];
817 labelList& nbrBfMap =
818 boundaryFaceMap[nbrProcMeshI][nbrIntI];
819
820 bfMap.setSize(faceCells.size());
821 nbrBfMap.setSize(faceCells.size());
822
823 forAll(faceCells, pfI)
824 {
825 lowerAddr_[allFacei] =
826 cellOffsets[procMeshI]+faceCells[pfI];
827 bfMap[pfI] = allFacei;
828 upperAddr_[allFacei] =
829 cellOffsets[nbrProcMeshI]+nbrFaceCells[pfI];
830 nbrBfMap[pfI] = (-allFacei-1);
831 allFacei++;
832 }
833 }
834 }
835 }
836 }
837 }
838 }
839
840
841 // Sort upper-tri order
842 {
843 labelList oldToNew
844 (
846 (
847 cellOffsets.last(), //nCells
848 lowerAddr_,
849 upperAddr_
850 )
851 );
852
853 forAll(faceMap, procMeshI)
854 {
855 labelList& map = faceMap[procMeshI];
856 forAll(map, i)
857 {
858 if (map[i] >= 0)
859 {
860 map[i] = oldToNew[map[i]];
861 }
862 else
863 {
864 label allFacei = -map[i]-1;
865 map[i] = -oldToNew[allFacei]-1;
866 }
867 }
868 }
869
870
871 inplaceReorder(oldToNew, lowerAddr_);
872 inplaceReorder(oldToNew, upperAddr_);
873
874 forAll(boundaryFaceMap, proci)
875 {
876 const labelList& bMap = boundaryMap[proci];
877 forAll(bMap, intI)
878 {
879 if (bMap[intI] == -1)
880 {
881 // Merged interface
882 labelList& bfMap = boundaryFaceMap[proci][intI];
883
884 forAll(bfMap, i)
885 {
886 if (bfMap[i] >= 0)
887 {
888 bfMap[i] = oldToNew[bfMap[i]];
889 }
890 else
891 {
892 label allFacei = -bfMap[i]-1;
893 bfMap[i] = (-oldToNew[allFacei]-1);
894 }
895 }
896 }
897 }
898 }
899 }
900
901
902 // Kept interfaces
903 // ~~~~~~~~~~~~~~~
904
905 interfaces_.setSize(unmergedMap.size() + nOtherInterfaces);
906 primitiveInterfaces_.setSize(interfaces_.size());
907
908 label allInterfacei = 0;
909
910 forAllConstIters(unmergedMap, iter)
911 {
912 const labelPairList& elems = iter.val();
913
914 // Sort processors in increasing order so both sides walk through in
915 // same order.
916 labelPairList procPairs(elems.size());
917 forAll(elems, i)
918 {
919 const labelPair& elem = elems[i];
920 label procMeshI = elem[0];
921 label interfacei = elem[1];
923 (
924 myMesh,
925 otherMeshes,
926 procMeshI
927 ).interfaces();
928
929 const processorLduInterface& pldui =
930 refCast<const processorLduInterface>
931 (
932 interfaces[interfacei]
933 );
934 label myProcNo = pldui.myProcNo();
935 label nbrProcNo = pldui.neighbProcNo();
936 procPairs[i] = labelPair
937 (
938 min(myProcNo, nbrProcNo),
939 max(myProcNo, nbrProcNo)
940 );
941 }
942
943 labelList order(sortedOrder(procPairs));
944
945 // Count
946 label n = 0;
947
948 forAll(order, i)
949 {
950 const labelPair& elem = elems[order[i]];
951 label procMeshI = elem[0];
952 label interfacei = elem[1];
954 (
955 myMesh,
956 otherMeshes,
957 procMeshI
958 ).interfaces();
959
960 n += interfaces[interfacei].faceCells().size();
961 }
962
963 // Size
964 labelField allFaceCells(n);
965 labelField allFaceRestrictAddressing(n);
966 n = 0;
967
968 // Fill
969 forAll(order, i)
970 {
971 const labelPair& elem = elems[order[i]];
972 label procMeshI = elem[0];
973 label interfacei = elem[1];
975 (
976 myMesh,
977 otherMeshes,
978 procMeshI
979 ).interfaces();
980
981 boundaryMap[procMeshI][interfacei] = allInterfacei;
982 labelList& bfMap = boundaryFaceMap[procMeshI][interfacei];
983
984 const labelUList& l = interfaces[interfacei].faceCells();
985 bfMap.setSize(l.size());
986
987 forAll(l, facei)
988 {
989 allFaceCells[n] = cellOffsets[procMeshI]+l[facei];
990 allFaceRestrictAddressing[n] = n;
991 bfMap[facei] = n;
992 n++;
993 }
994 }
995
996
997 // Find out local and remote processor in new communicator
998
999 label neighbProcNo = -1;
1000
1001 // See what the two processors map onto
1002
1003 if (iter.key()[0] == myAgglom)
1004 {
1005 if (iter.key()[1] == myAgglom)
1006 {
1008 << "problem procEdge:" << iter.key()
1009 << exit(FatalError);
1010 }
1011
1012 neighbProcNo = iter.key()[1];
1013 }
1014 else
1015 {
1016 if (iter.key()[1] != myAgglom)
1017 {
1019 << "problem procEdge:" << iter.key()
1020 << exit(FatalError);
1021 }
1022
1023 neighbProcNo = iter.key()[0];
1024 }
1025
1026 primitiveInterfaces_.set
1027 (
1028 allInterfacei,
1030 (
1031 allInterfacei,
1032 interfaces_,
1033 allFaceCells,
1034 allFaceRestrictAddressing,
1035 comm_,
1036 myAgglom,
1037 neighbProcNo,
1038 tensorField(), // forwardT
1039 Pstream::msgType() // tag
1040 )
1041 );
1042 interfaces_.set(allInterfacei, &primitiveInterfaces_[allInterfacei]);
1043
1044 if (debug)
1045 {
1046 Pout<< "Created " << interfaces_[allInterfacei].type()
1047 << " interface at " << allInterfacei
1048 << " comm:" << comm_
1049 << " myProcNo:" << myAgglom
1050 << " neighbProcNo:" << neighbProcNo
1051 << " nFaces:" << allFaceCells.size()
1052 << endl;
1053 }
1054
1055
1056 allInterfacei++;
1057 }
1058
1059
1060 patchSchedule_ = nonBlockingSchedule<processorGAMGInterface>(interfaces_);
1061
1062 if (debug)
1063 {
1064 checkUpperTriangular(cellOffsets.last(), lowerAddr_, upperAddr_);
1065 }
1066}
1067
1068
1069// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1070
1072(
1073 const lduMesh& myMesh,
1074 const PtrList<lduPrimitiveMesh>& otherMeshes,
1075 const label meshI
1076)
1077{
1078 return (meshI == 0 ? myMesh : otherMeshes[meshI-1]);
1079}
1080
1081
1083(
1084 const label comm,
1085 const lduMesh& mesh,
1086 const labelList& procIDs,
1087 PtrList<lduPrimitiveMesh>& otherMeshes
1088)
1089{
1090 // Force calculation of schedule (since does parallel comms)
1091 (void)mesh.lduAddr().patchSchedule();
1092
1093 if (Pstream::myProcNo(comm) == procIDs[0])
1094 {
1095 otherMeshes.setSize(procIDs.size()-1);
1096
1097 // Slave meshes
1098 for (label i = 1; i < procIDs.size(); ++i)
1099 {
1100 //Pout<< "on master :"
1101 // << " receiving from slave " << procIDs[i] << endl;
1102
1103 IPstream fromSlave
1104 (
1106 procIDs[i],
1107 0, // bufSize
1109 comm
1110 );
1111
1112 label nCells = readLabel(fromSlave);
1113 labelList lowerAddr(fromSlave);
1114 labelList upperAddr(fromSlave);
1115 boolList validInterface(fromSlave);
1116
1117
1118 // Construct mesh without interfaces
1119 otherMeshes.set
1120 (
1121 i-1,
1123 (
1124 nCells,
1125 lowerAddr,
1126 upperAddr,
1127 comm,
1128 true // reuse
1129 )
1130 );
1131
1132 // Construct GAMGInterfaces
1133 lduInterfacePtrsList newInterfaces(validInterface.size());
1134 forAll(validInterface, intI)
1135 {
1136 if (validInterface[intI])
1137 {
1138 word coupleType(fromSlave);
1139
1140 newInterfaces.set
1141 (
1142 intI,
1144 (
1145 coupleType,
1146 intI,
1147 otherMeshes[i-1].rawInterfaces(),
1148 fromSlave
1149 ).ptr()
1150 );
1151 }
1152 }
1153
1154 otherMeshes[i-1].addInterfaces
1155 (
1156 newInterfaces,
1157 nonBlockingSchedule<processorGAMGInterface>
1158 (
1159 newInterfaces
1160 )
1161 );
1162 }
1163 }
1164 else if (procIDs.found(Pstream::myProcNo(comm)))
1165 {
1166 // Send to master
1167
1168 const lduAddressing& addressing = mesh.lduAddr();
1169 lduInterfacePtrsList interfaces(mesh.interfaces());
1170 boolList validInterface(interfaces.size());
1171 forAll(interfaces, intI)
1172 {
1173 validInterface[intI] = interfaces.set(intI);
1174 }
1175
1176 OPstream toMaster
1177 (
1179 procIDs[0],
1180 0,
1182 comm
1183 );
1184
1185 toMaster
1186 << addressing.size()
1187 << addressing.lowerAddr()
1188 << addressing.upperAddr()
1189 << validInterface;
1190
1191 forAll(interfaces, intI)
1192 {
1193 if (interfaces.set(intI))
1194 {
1195 const GAMGInterface& interface = refCast<const GAMGInterface>
1196 (
1197 interfaces[intI]
1198 );
1199
1200 toMaster << interface.type();
1201 interface.write(toMaster);
1202 }
1203 }
1204 }
1205}
1206
1207
1208// ************************************************************************* //
label n
label totalSize() const
The total addressed size.
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void resize(const label len)
Definition: DynamicListI.H:353
Map from edge (expressed as its endpoints) to value. For easier forward declaration it is currently i...
Definition: EdgeMap.H:54
Abstract base class for GAMG agglomerated interfaces.
Definition: GAMGInterface.H:57
const_iterator cfind(const Key &key) const
Find and return an const_iterator set at the hashed entry.
Definition: HashTableI.H:141
label size() const noexcept
The number of elements in table.
Definition: HashTableI.H:52
Input inter-processor communications stream.
Definition: IPstream.H:57
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
void setSize(const label n)
Alias for resize()
Definition: List.H:218
Output inter-processor communications stream.
Definition: OPstream.H:57
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers....
Definition: PtrList.H:73
const T * set(const label i) const
Definition: PtrList.H:138
void setSize(const label newLen)
Same as resize()
Definition: PtrList.H:151
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
T & first()
Return the first element of the list.
Definition: UListI.H:202
bool found(const T &val, label pos=0) const
True if the value if found in the list.
Definition: UListI.H:265
label find(const T &val, label pos=0) const
Find index of the first occurrence of the value.
Definition: UList.C:212
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
@ nonBlocking
"nonBlocking"
static int & msgType() noexcept
Message tag of standard messages.
Definition: UPstream.H:556
static void waitRequests(const label start=0)
Wait until all requests (from start onwards) have finished.
Definition: UPstream.C:100
static bool & parRun() noexcept
Test if this a parallel run.
Definition: UPstream.H:433
const T * set(const label i) const
Definition: UPtrList.H:248
void setSize(const label n)
Alias for resize()
Definition: UPtrList.H:261
label size() const noexcept
The number of elements in the list.
Definition: UPtrListI.H:106
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:66
Smooth ATC in cells next to a set of patches supplied by type.
Definition: faceCells.H:59
virtual lduInterfacePtrsList interfaces() const
Return a list of pointers for each patch.
Definition: fvMesh.C:735
virtual const lduAddressing & lduAddr() const
Return ldu addressing.
Definition: fvMesh.C:718
virtual label comm() const
Return communicator used for parallel communication.
Definition: fvMesh.H:326
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:68
label localStart() const
My local start.
Definition: globalIndexI.H:195
The class contains the addressing required by the lduMatrix: upper, lower and losort.
virtual const lduSchedule & patchSchedule() const =0
Return patch field evaluation schedule.
virtual const labelUList & upperAddr() const =0
Return upper addressing.
label size() const
Return number of equations.
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
An abstract base class for implicitly-coupled interfaces e.g. processor and cyclic patches.
Definition: lduInterface.H:58
virtual const labelUList & faceCells() const =0
Return faceCell addressing.
Abstract base class for meshes which provide LDU addressing for the construction of lduMatrix and LDU...
Definition: lduMesh.H:63
virtual label comm() const =0
Return communicator used for parallel communication.
virtual const lduAddressing & lduAddr() const =0
Return ldu addressing.
virtual lduInterfacePtrsList interfaces() const =0
Simplest concrete lduMesh that stores the addressing needed by lduMatrix.
static void checkUpperTriangular(const label size, const labelUList &l, const labelUList &u)
Check if in upper-triangular ordering.
static labelListList globalCellCells(const lduMesh &mesh, const globalIndex &globalNumbering)
Calculate global cell-cells.
void addInterfaces(lduInterfacePtrsList &interfaces, const lduSchedule &ps)
Add interfaces to a mesh constructed without.
PtrList< const lduInterface > & primitiveInterfaces()
Return a non-const list of primitive interfaces.
static labelList upperTriOrder(const label nCells, const labelUList &lower, const labelUList &upper)
Calculate upper-triangular order.
virtual label comm() const
Return communicator used for parallel communication.
virtual lduInterfacePtrsList interfaces() const
static void gather(const label comm, const lduMesh &mesh, const labelList &procIDs, PtrList< lduPrimitiveMesh > &otherMeshes)
Gather meshes from other processors onto procIDs[0].
Less operator for pairs of <processor><index>
procLess(const labelPairList &list)
bool operator()(const label a, const label b)
int myProcNo() const noexcept
Return processor number.
GAMG agglomerated processor interface.
An abstract base class for processor coupled interfaces.
virtual int neighbProcNo() const =0
Return neighbour processor number (rank in communicator)
virtual int myProcNo() const =0
Return processor number (rank in communicator)
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
dynamicFvMesh & mesh
Foam::PtrList< Foam::fvMesh > meshes(regionNames.size())
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
Pair< label > labelPair
A pair of labels.
Definition: Pair.H:57
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
To & refCast(From &r)
Reference type cast template function.
Definition: typeInfo.H:131
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
label readLabel(const char *buf)
Parse entire buffer as a label, skipping leading/trailing whitespace.
Definition: label.H:66
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
void inplaceReorder(const labelUList &oldToNew, ListType &input, const bool prune=false)
Inplace reorder the elements of a list.
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:144
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
error FatalError
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: IndirectList.H:68
volScalarField & b
Definition: createFields.H:27
interfaceProperties interface(alpha1, U, thermo->transportPropertiesDict())
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278