lduPrimitiveMesh.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2013-2020 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "lduPrimitiveMesh.H"
27 #include "processorLduInterface.H"
28 #include "EdgeMap.H"
29 #include "labelPair.H"
30 #include "processorGAMGInterface.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(lduPrimitiveMesh, 0);
37 
38  //- Less operator for pairs of <processor><index>
39  class procLess
40  {
41  const labelPairList& lst_;
42 
43  public:
44 
45  procLess(const labelPairList& lst)
46  :
47  lst_(lst)
48  {}
49 
50  bool operator()(const label a, const label b)
51  {
52  return lst_[a].first() < lst_[b].first();
53  }
54  };
55 
56 }
57 
58 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
59 
60 void Foam::lduPrimitiveMesh::checkUpperTriangular
61 (
62  const label size,
63  const labelUList& l,
64  const labelUList& u
65 )
66 {
67  forAll(l, facei)
68  {
69  if (u[facei] < l[facei])
70  {
72  << "Reversed face. Problem at face " << facei
73  << " l:" << l[facei] << " u:" << u[facei]
74  << abort(FatalError);
75  }
76  if (l[facei] < 0 || u[facei] < 0 || u[facei] >= size)
77  {
79  << "Illegal cell label. Problem at face " << facei
80  << " l:" << l[facei] << " u:" << u[facei]
81  << abort(FatalError);
82  }
83  }
84 
85  for (label facei=1; facei < l.size(); facei++)
86  {
87  if (l[facei-1] > l[facei])
88  {
90  << "Lower not in incremental cell order."
91  << " Problem at face " << facei
92  << " l:" << l[facei] << " u:" << u[facei]
93  << " previous l:" << l[facei-1]
94  << abort(FatalError);
95  }
96  else if (l[facei-1] == l[facei])
97  {
98  // Same cell.
99  if (u[facei-1] > u[facei])
100  {
102  << "Upper not in incremental cell order."
103  << " Problem at face " << facei
104  << " l:" << l[facei] << " u:" << u[facei]
105  << " previous u:" << u[facei-1]
106  << abort(FatalError);
107  }
108  }
109  }
110 }
111 
112 
113 Foam::label Foam::lduPrimitiveMesh::totalSize
114 (
115  const PtrList<lduPrimitiveMesh>& meshes
116 )
117 {
118  label size = 0;
119 
120  forAll(meshes, i)
121  {
122  size += meshes[i].lduAddr().size();
123  }
124  return size;
125 }
126 
127 
128 Foam::labelList Foam::lduPrimitiveMesh::upperTriOrder
129 (
130  const label nCells,
131  const labelUList& lower,
132  const labelUList& upper
133 )
134 {
135  labelList nNbrs(nCells, 0);
136 
137  // Count number of upper neighbours
138  forAll(lower, facei)
139  {
140  if (upper[facei] < lower[facei])
141  {
143  << "Problem at face:" << facei
144  << " lower:" << lower[facei]
145  << " upper:" << upper[facei]
146  << exit(FatalError);
147  }
148  nNbrs[lower[facei]]++;
149  }
150 
151  // Construct cell-upper cell addressing
152  labelList offsets(nCells+1);
153  offsets[0] = 0;
154  forAll(nNbrs, celli)
155  {
156  offsets[celli+1] = offsets[celli]+nNbrs[celli];
157  }
158 
159  nNbrs = offsets;
160 
161  labelList cellToFaces(offsets.last());
162  forAll(upper, facei)
163  {
164  label celli = lower[facei];
165  cellToFaces[nNbrs[celli]++] = facei;
166  }
167 
168  // Sort
169 
170  labelList oldToNew(lower.size());
171 
172  labelList order;
173  labelList nbr;
174 
175  label newFacei = 0;
176 
177  for (label celli = 0; celli < nCells; celli++)
178  {
179  label startOfCell = offsets[celli];
180  label nNbr = offsets[celli+1] - startOfCell;
181 
182  nbr.setSize(nNbr);
183  order.setSize(nNbr);
184  forAll(order, i)
185  {
186  nbr[i] = upper[cellToFaces[offsets[celli]+i]];
187  }
188  sortedOrder(nbr, order);
189 
190  forAll(order, i)
191  {
192  label index = order[i];
193  oldToNew[cellToFaces[startOfCell + index]] = newFacei++;
194  }
195  }
196 
197  return oldToNew;
198 }
199 
200 
201 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
202 
204 (
205  const label nCells,
206  labelList& l,
207  labelList& u,
208  const label comm,
209  bool reuse
210 )
211 :
212  lduAddressing(nCells),
213  lowerAddr_(l, reuse),
214  upperAddr_(u, reuse),
215  comm_(comm)
216 {}
217 
218 
220 (
221  lduInterfacePtrsList& interfaces,
222  const lduSchedule& ps
223 )
224 {
225  interfaces_ = interfaces;
226  patchSchedule_ = ps;
227 
228  // Create interfaces
229  primitiveInterfaces_.setSize(interfaces_.size());
230  forAll(interfaces_, i)
231  {
232  if (interfaces_.set(i))
233  {
234  primitiveInterfaces_.set(i, &interfaces_[i]);
235  }
236  }
237 }
238 
239 
241 (
242  const label nCells,
243  labelList& l,
244  labelList& u,
245  PtrList<const lduInterface>& primitiveInterfaces,
246  const lduSchedule& ps,
247  const label comm
248 )
249 :
250  lduAddressing(nCells),
251  lowerAddr_(l, true),
252  upperAddr_(u, true),
253  primitiveInterfaces_(0),
254  patchSchedule_(ps),
255  comm_(comm)
256 {
257  primitiveInterfaces_.transfer(primitiveInterfaces);
258 
259  // Create interfaces
260  interfaces_.setSize(primitiveInterfaces_.size());
261  forAll(primitiveInterfaces_, i)
262  {
263  if (primitiveInterfaces_.set(i))
264  {
265  interfaces_.set(i, &primitiveInterfaces_[i]);
266  }
267  }
268 }
269 
270 
272 (
273  const label comm,
274  const labelList& procAgglomMap,
275 
276  const labelList& procIDs,
277  const lduMesh& myMesh,
278  const PtrList<lduPrimitiveMesh>& otherMeshes,
279 
280  labelList& cellOffsets,
281  labelList& faceOffsets,
282  labelListList& faceMap,
283  labelListList& boundaryMap,
284  labelListListList& boundaryFaceMap
285 )
286 :
287  lduAddressing(myMesh.lduAddr().size() + totalSize(otherMeshes)),
288  lowerAddr_(0),
289  upperAddr_(0),
290  interfaces_(0),
291  patchSchedule_(0),
292  comm_(comm)
293 {
294  const label currentComm = myMesh.comm();
295 
296  forAll(otherMeshes, i)
297  {
298  if (otherMeshes[i].comm() != currentComm)
299  {
301  << "Communicator " << otherMeshes[i].comm()
302  << " at index " << i
303  << " differs from that of predecessor "
304  << currentComm
305  << endl; // exit(FatalError);
306  }
307  }
308 
309  const label nMeshes = otherMeshes.size()+1;
310 
311  const label myAgglom = procAgglomMap[UPstream::myProcNo(currentComm)];
312 
313  if (lduPrimitiveMesh::debug)
314  {
315  Pout<< "I am " << UPstream::myProcNo(currentComm)
316  << " agglomerating into " << myAgglom
317  << " as are " << findIndices(procAgglomMap, myAgglom)
318  << endl;
319  }
320 
321 
322  forAll(procIDs, i)
323  {
324  if (procAgglomMap[procIDs[i]] != procAgglomMap[procIDs[0]])
325  {
327  << "Processor " << procIDs[i]
328  << " agglomerates to " << procAgglomMap[procIDs[i]]
329  << " whereas other processors " << procIDs
330  << " agglomerate to "
331  << UIndirectList<label>(procAgglomMap, procIDs)
332  << exit(FatalError);
333  }
334  }
335 
336 
337  // Cells get added in order.
338  cellOffsets.setSize(nMeshes+1);
339  cellOffsets[0] = 0;
340  for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
341  {
342  const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
343 
344  cellOffsets[procMeshI+1] =
345  cellOffsets[procMeshI]
346  + procMesh.lduAddr().size();
347  }
348 
349 
350  // Faces initially get added in order, sorted later
351  labelList internalFaceOffsets(nMeshes+1);
352  internalFaceOffsets[0] = 0;
353  for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
354  {
355  const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
356 
357  internalFaceOffsets[procMeshI+1] =
358  internalFaceOffsets[procMeshI]
359  + procMesh.lduAddr().lowerAddr().size();
360  }
361 
362  // Count how faces get added. Interfaces in between get merged.
363 
364  // Merged interfaces: map from two coarse processors back to
365  // - procMeshes
366  // - interface in procMesh
367  // (estimate size from number of patches of mesh0)
368  EdgeMap<labelPairList> mergedMap(2*myMesh.interfaces().size());
369 
370  // Unmerged interfaces: map from two coarse processors back to
371  // - procMeshes
372  // - interface in procMesh
373  EdgeMap<labelPairList> unmergedMap(mergedMap.size());
374 
375  boundaryMap.setSize(nMeshes);
376  boundaryFaceMap.setSize(nMeshes);
377 
378 
379  label nOtherInterfaces = 0;
380  labelList nCoupledFaces(nMeshes, 0);
381 
382  for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
383  {
384  const lduInterfacePtrsList interfaces =
385  mesh(myMesh, otherMeshes, procMeshI).interfaces();
386 
387  // Initialise all boundaries as merged
388  boundaryMap[procMeshI].setSize(interfaces.size(), -1);
389  boundaryFaceMap[procMeshI].setSize(interfaces.size());
390 
391  // Get sorted order of processors
392  forAll(interfaces, intI)
393  {
394  if (interfaces.set(intI))
395  {
396  const lduInterface& ldui = interfaces[intI];
397 
398  if (isA<processorLduInterface>(ldui))
399  {
400  const processorLduInterface& pldui =
401  refCast<const processorLduInterface>(ldui);
402 
403  label agglom0 = procAgglomMap[pldui.myProcNo()];
404  label agglom1 = procAgglomMap[pldui.neighbProcNo()];
405 
406  const edge procEdge(agglom0, agglom1);
407 
408  if (agglom0 != myAgglom && agglom1 != myAgglom)
409  {
411  << "At mesh from processor " << procIDs[procMeshI]
412  << " have interface " << intI
413  << " with myProcNo:" << pldui.myProcNo()
414  << " with neighbProcNo:" << pldui.neighbProcNo()
415  << exit(FatalError);
416  }
417  else if (agglom0 == myAgglom && agglom1 == myAgglom)
418  {
419  // Merged interface
420  if (debug)
421  {
422  Pout<< "merged interface: myProcNo:"
423  << pldui.myProcNo()
424  << " nbr:" << pldui.neighbProcNo()
425  << " size:" << ldui.faceCells().size()
426  << endl;
427  }
428 
429  label nbrProcMeshI = findIndex
430  (
431  procIDs,
432  pldui.neighbProcNo()
433  );
434 
435  if (procMeshI < nbrProcMeshI)
436  {
437  // I am 'master' since get lowest numbered cells
438  nCoupledFaces[procMeshI] += ldui.faceCells().size();
439  }
440 
442  mergedMap.find(procEdge);
443 
444  if (iter != mergedMap.end())
445  {
446  iter().append(labelPair(procMeshI, intI));
447  }
448  else
449  {
450  mergedMap.insert
451  (
452  procEdge,
453  labelPairList(1, labelPair(procMeshI, intI))
454  );
455  }
456  }
457  else
458  {
459  if (debug)
460  {
461  Pout<< "external interface: myProcNo:"
462  << pldui.myProcNo()
463  << " nbr:" << pldui.neighbProcNo()
464  << " size:" << ldui.faceCells().size()
465  << endl;
466  }
467 
469  unmergedMap.find(procEdge);
470 
471  if (iter != unmergedMap.end())
472  {
473  iter().append(labelPair(procMeshI, intI));
474  }
475  else
476  {
477  unmergedMap.insert
478  (
479  procEdge,
480  labelPairList(1, labelPair(procMeshI, intI))
481  );
482  }
483  }
484  }
485  else
486  {
487  // Still external (non proc) interface
489  << "At mesh from processor " << procIDs[procMeshI]
490  << " have interface " << intI
491  << " of unhandled type " << interfaces[intI].type()
492  << exit(FatalError);
493 
494  nOtherInterfaces++;
495  }
496  }
497  }
498  }
499 
500 
501 
502  if (debug)
503  {
504  Pout<< "Remaining interfaces:" << endl;
505  forAllConstIter(EdgeMap<labelPairList>, unmergedMap, iter)
506  {
507  Pout<< " agglom procEdge:" << iter.key() << endl;
508  const labelPairList& elems = iter();
509  forAll(elems, i)
510  {
511  label procMeshI = elems[i][0];
512  label interfacei = elems[i][1];
513  const lduInterfacePtrsList interfaces =
514  mesh(myMesh, otherMeshes, procMeshI).interfaces();
515 
516  const processorLduInterface& pldui =
517  refCast<const processorLduInterface>
518  (
519  interfaces[interfacei]
520  );
521 
522  Pout<< " proc:" << procIDs[procMeshI]
523  << " interfacei:" << interfacei
524  << " between:" << pldui.myProcNo()
525  << " and:" << pldui.neighbProcNo()
526  << endl;
527  }
528  Pout<< endl;
529  }
530  }
531  if (debug)
532  {
533  Pout<< "Merged interfaces:" << endl;
534  forAllConstIter(EdgeMap<labelPairList>, mergedMap, iter)
535  {
536  Pout<< " agglom procEdge:" << iter.key() << endl;
537  const labelPairList& elems = iter();
538 
539  forAll(elems, i)
540  {
541  label procMeshI = elems[i][0];
542  label interfacei = elems[i][1];
543  const lduInterfacePtrsList interfaces =
544  mesh(myMesh, otherMeshes, procMeshI).interfaces();
545  const processorLduInterface& pldui =
546  refCast<const processorLduInterface>
547  (
548  interfaces[interfacei]
549  );
550 
551  Pout<< " proc:" << procIDs[procMeshI]
552  << " interfacei:" << interfacei
553  << " between:" << pldui.myProcNo()
554  << " and:" << pldui.neighbProcNo()
555  << endl;
556  }
557  Pout<< endl;
558  }
559  }
560 
561 
562  // Adapt faceOffsets for internal interfaces
563  faceOffsets.setSize(nMeshes+1);
564  faceOffsets[0] = 0;
565  faceMap.setSize(nMeshes);
566  for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
567  {
568  const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
569  label nInternal = procMesh.lduAddr().lowerAddr().size();
570 
571  faceOffsets[procMeshI+1] =
572  faceOffsets[procMeshI]
573  + nInternal
574  + nCoupledFaces[procMeshI];
575 
576  labelList& map = faceMap[procMeshI];
577  map.setSize(nInternal);
578  forAll(map, i)
579  {
580  map[i] = faceOffsets[procMeshI] + i;
581  }
582  }
583 
584 
585  // Combine upper and lower
586  lowerAddr_.setSize(faceOffsets.last(), -1);
587  upperAddr_.setSize(lowerAddr_.size(), -1);
588 
589 
590  // Old internal faces and resolved coupled interfaces
591  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
592 
593  for (label procMeshI = 0; procMeshI < nMeshes; procMeshI++)
594  {
595  const lduMesh& procMesh = mesh(myMesh, otherMeshes, procMeshI);
596 
597  const labelUList& l = procMesh.lduAddr().lowerAddr();
598  const labelUList& u = procMesh.lduAddr().upperAddr();
599 
600  // Add internal faces
601  label allFacei = faceOffsets[procMeshI];
602 
603  forAll(l, facei)
604  {
605  lowerAddr_[allFacei] = cellOffsets[procMeshI]+l[facei];
606  upperAddr_[allFacei] = cellOffsets[procMeshI]+u[facei];
607  allFacei++;
608  }
609 
610 
611  // Add merged interfaces
612  const lduInterfacePtrsList interfaces = procMesh.interfaces();
613 
614  forAll(interfaces, intI)
615  {
616  if (interfaces.set(intI))
617  {
618  const lduInterface& ldui = interfaces[intI];
619 
620  if (isA<processorLduInterface>(ldui))
621  {
622  const processorLduInterface& pldui =
623  refCast<const processorLduInterface>(ldui);
624 
625  // Look up corresponding interfaces
626  label myP = pldui.myProcNo();
627  label nbrP = pldui.neighbProcNo();
628  label nbrProcMeshI = findIndex(procIDs, nbrP);
629 
630  if (procMeshI < nbrProcMeshI)
631  {
632  // I am 'master' since my cell numbers will be lower
633  // since cells get added in procMeshI order.
634 
635  label agglom0 = procAgglomMap[myP];
636  label agglom1 = procAgglomMap[nbrP];
637 
639  mergedMap.find(edge(agglom0, agglom1));
640 
641  if (fnd != mergedMap.end())
642  {
643  const labelPairList& elems = fnd();
644 
645  // Find nbrP in elems
646  label nbrIntI = -1;
647  forAll(elems, i)
648  {
649  label proci = elems[i][0];
650  label interfacei = elems[i][1];
651  const lduInterfacePtrsList interfaces =
652  mesh
653  (
654  myMesh,
655  otherMeshes,
656  proci
657  ).interfaces();
658  const processorLduInterface& pldui =
659  refCast<const processorLduInterface>
660  (
661  interfaces[interfacei]
662  );
663 
664  if
665  (
666  elems[i][0] == nbrProcMeshI
667  && pldui.neighbProcNo() == procIDs[procMeshI]
668  )
669  {
670  nbrIntI = elems[i][1];
671  break;
672  }
673  }
674 
675 
676  if (nbrIntI == -1)
677  {
679  << "elems:" << elems << abort(FatalError);
680  }
681 
682 
683  const lduInterfacePtrsList nbrInterfaces = mesh
684  (
685  myMesh,
686  otherMeshes,
687  nbrProcMeshI
688  ).interfaces();
689 
690 
691  const labelUList& faceCells =
692  ldui.faceCells();
693  const labelUList& nbrFaceCells =
694  nbrInterfaces[nbrIntI].faceCells();
695 
696  if (faceCells.size() != nbrFaceCells.size())
697  {
699  << "faceCells:" << faceCells
700  << " nbrFaceCells:" << nbrFaceCells
701  << abort(FatalError);
702  }
703 
704 
705  labelList& bfMap =
706  boundaryFaceMap[procMeshI][intI];
707  labelList& nbrBfMap =
708  boundaryFaceMap[nbrProcMeshI][nbrIntI];
709 
710  bfMap.setSize(faceCells.size());
711  nbrBfMap.setSize(faceCells.size());
712 
713  forAll(faceCells, pfI)
714  {
715  lowerAddr_[allFacei] =
716  cellOffsets[procMeshI]+faceCells[pfI];
717  bfMap[pfI] = allFacei;
718  upperAddr_[allFacei] =
719  cellOffsets[nbrProcMeshI]+nbrFaceCells[pfI];
720  nbrBfMap[pfI] = (-allFacei-1);
721  allFacei++;
722  }
723  }
724  }
725  }
726  }
727  }
728  }
729 
730 
731  // Sort upper-tri order
732  {
733  labelList oldToNew
734  (
735  upperTriOrder
736  (
737  cellOffsets.last(), // nCells
738  lowerAddr_,
739  upperAddr_
740  )
741  );
742 
743  forAll(faceMap, procMeshI)
744  {
745  labelList& map = faceMap[procMeshI];
746  forAll(map, i)
747  {
748  if (map[i] >= 0)
749  {
750  map[i] = oldToNew[map[i]];
751  }
752  else
753  {
754  label allFacei = -map[i]-1;
755  map[i] = -oldToNew[allFacei]-1;
756  }
757  }
758  }
759 
760 
761  inplaceReorder(oldToNew, lowerAddr_);
762  inplaceReorder(oldToNew, upperAddr_);
763 
764  forAll(boundaryFaceMap, proci)
765  {
766  const labelList& bMap = boundaryMap[proci];
767  forAll(bMap, intI)
768  {
769  if (bMap[intI] == -1)
770  {
771  // Merged interface
772  labelList& bfMap = boundaryFaceMap[proci][intI];
773 
774  forAll(bfMap, i)
775  {
776  if (bfMap[i] >= 0)
777  {
778  bfMap[i] = oldToNew[bfMap[i]];
779  }
780  else
781  {
782  label allFacei = -bfMap[i]-1;
783  bfMap[i] = (-oldToNew[allFacei]-1);
784  }
785  }
786  }
787  }
788  }
789  }
790 
791 
792  // Kept interfaces
793  // ~~~~~~~~~~~~~~~
794 
795  interfaces_.setSize(unmergedMap.size() + nOtherInterfaces);
796  primitiveInterfaces_.setSize(interfaces_.size());
797 
798  label allInterfacei = 0;
799 
800  forAllConstIter(EdgeMap<labelPairList>, unmergedMap, iter)
801  {
802  const labelPairList& elems = iter();
803 
804  // Sort processors in increasing order so both sides walk through in
805  // same order.
806  labelPairList procPairs(elems.size());
807  forAll(elems, i)
808  {
809  const labelPair& elem = elems[i];
810  label procMeshI = elem[0];
811  label interfacei = elem[1];
812  const lduInterfacePtrsList interfaces = mesh
813  (
814  myMesh,
815  otherMeshes,
816  procMeshI
817  ).interfaces();
818 
819  const processorLduInterface& pldui =
820  refCast<const processorLduInterface>
821  (
822  interfaces[interfacei]
823  );
824  label myProcNo = pldui.myProcNo();
825  label nbrProcNo = pldui.neighbProcNo();
826  procPairs[i] = labelPair
827  (
828  min(myProcNo, nbrProcNo),
829  max(myProcNo, nbrProcNo)
830  );
831  }
832 
833  labelList order;
834  sortedOrder(procPairs, order);
835 
836  // Count
837  label n = 0;
838 
839  forAll(order, i)
840  {
841  const labelPair& elem = elems[order[i]];
842  label procMeshI = elem[0];
843  label interfacei = elem[1];
844  const lduInterfacePtrsList interfaces = mesh
845  (
846  myMesh,
847  otherMeshes,
848  procMeshI
849  ).interfaces();
850 
851  n += interfaces[interfacei].faceCells().size();
852  }
853 
854  // Size
855  labelField allFaceCells(n);
856  labelField allFaceRestrictAddressing(n);
857  n = 0;
858 
859  // Fill
860  forAll(order, i)
861  {
862  const labelPair& elem = elems[order[i]];
863  label procMeshI = elem[0];
864  label interfacei = elem[1];
865  const lduInterfacePtrsList interfaces = mesh
866  (
867  myMesh,
868  otherMeshes,
869  procMeshI
870  ).interfaces();
871 
872  boundaryMap[procMeshI][interfacei] = allInterfacei;
873  labelList& bfMap = boundaryFaceMap[procMeshI][interfacei];
874 
875  const labelUList& l = interfaces[interfacei].faceCells();
876  bfMap.setSize(l.size());
877 
878  forAll(l, facei)
879  {
880  allFaceCells[n] = cellOffsets[procMeshI]+l[facei];
881  allFaceRestrictAddressing[n] = n;
882  bfMap[facei] = n;
883  n++;
884  }
885  }
886 
887 
888  // Find out local and remote processor in new communicator
889 
890  label neighbProcNo = -1;
891 
892  // See what the two processors map onto
893 
894  if (iter.key()[0] == myAgglom)
895  {
896  if (iter.key()[1] == myAgglom)
897  {
899  << "problem procEdge:" << iter.key()
900  << exit(FatalError);
901  }
902 
903  neighbProcNo = iter.key()[1];
904  }
905  else
906  {
907  if (iter.key()[1] != myAgglom)
908  {
910  << "problem procEdge:" << iter.key()
911  << exit(FatalError);
912  }
913 
914  neighbProcNo = iter.key()[0];
915  }
916 
917  primitiveInterfaces_.set
918  (
919  allInterfacei,
921  (
922  allInterfacei,
923  interfaces_,
924  allFaceCells,
925  allFaceRestrictAddressing,
926  comm_,
927  myAgglom,
928  neighbProcNo,
929  transformer(), // forwardT
930  Pstream::msgType() // tag
931  )
932  );
933  interfaces_.set(allInterfacei, &primitiveInterfaces_[allInterfacei]);
934 
935  if (debug)
936  {
937  Pout<< "Created " << interfaces_[allInterfacei].type()
938  << " interface at " << allInterfacei
939  << " comm:" << comm_
940  << " myProcNo:" << myAgglom
941  << " neighbProcNo:" << neighbProcNo
942  << " nFaces:" << allFaceCells.size()
943  << endl;
944  }
945 
946 
947  allInterfacei++;
948  }
949 
950 
951  patchSchedule_ = nonBlockingSchedule<processorGAMGInterface>(interfaces_);
952 
953  if (debug)
954  {
955  checkUpperTriangular(cellOffsets.last(), lowerAddr_, upperAddr_);
956  }
957 }
958 
959 
960 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
961 
963 (
964  const lduMesh& myMesh,
965  const PtrList<lduPrimitiveMesh>& otherMeshes,
966  const label meshI
967 )
968 {
969  return (meshI == 0 ? myMesh : otherMeshes[meshI-1]);
970 }
971 
972 
974 (
975  const label comm,
976  const lduMesh& mesh,
977  const labelList& procIDs,
978  PtrList<lduPrimitiveMesh>& otherMeshes
979 )
980 {
981  // Force calculation of schedule (since does parallel comms)
982  (void)mesh.lduAddr().patchSchedule();
983 
984  if (Pstream::myProcNo(comm) == procIDs[0])
985  {
986  otherMeshes.setSize(procIDs.size()-1);
987 
988  // Slave meshes
989  for (label i = 1; i < procIDs.size(); i++)
990  {
991  // Pout<< "on master :"
992  // << " receiving from slave " << procIDs[i] << endl;
993 
994  IPstream fromSlave
995  (
997  procIDs[i],
998  0, // bufSize
1000  comm
1001  );
1002 
1003  label nCells = readLabel(fromSlave);
1004  labelList lowerAddr(fromSlave);
1005  labelList upperAddr(fromSlave);
1006  boolList validInterface(fromSlave);
1007 
1008 
1009  // Construct mesh without interfaces
1010  otherMeshes.set
1011  (
1012  i-1,
1013  new lduPrimitiveMesh
1014  (
1015  nCells,
1016  lowerAddr,
1017  upperAddr,
1018  comm,
1019  true // reuse
1020  )
1021  );
1022 
1023  // Construct GAMGInterfaces
1024  lduInterfacePtrsList newInterfaces(validInterface.size());
1025  forAll(validInterface, intI)
1026  {
1027  if (validInterface[intI])
1028  {
1029  word coupleType(fromSlave);
1030 
1031  newInterfaces.set
1032  (
1033  intI,
1035  (
1036  coupleType,
1037  intI,
1038  otherMeshes[i-1].rawInterfaces(),
1039  fromSlave
1040  ).ptr()
1041  );
1042  }
1043  }
1044 
1045  otherMeshes[i-1].addInterfaces
1046  (
1047  newInterfaces,
1048  nonBlockingSchedule<processorGAMGInterface>
1049  (
1050  newInterfaces
1051  )
1052  );
1053  }
1054  }
1055  else if (findIndex(procIDs, Pstream::myProcNo(comm)) != -1)
1056  {
1057  // Send to master
1058 
1059  const lduAddressing& addressing = mesh.lduAddr();
1060  lduInterfacePtrsList interfaces(mesh.interfaces());
1061  boolList validInterface(interfaces.size());
1062  forAll(interfaces, intI)
1063  {
1064  validInterface[intI] = interfaces.set(intI);
1065  }
1066 
1067  OPstream toMaster
1068  (
1070  procIDs[0],
1071  0,
1072  Pstream::msgType(),
1073  comm
1074  );
1075 
1076  toMaster
1077  << addressing.size()
1078  << addressing.lowerAddr()
1079  << addressing.upperAddr()
1080  << validInterface;
1081 
1082  forAll(interfaces, intI)
1083  {
1084  if (interfaces.set(intI))
1085  {
1086  const GAMGInterface& interface = refCast<const GAMGInterface>
1087  (
1088  interfaces[intI]
1089  );
1090 
1091  toMaster << interface.type();
1092  interface.write(toMaster);
1093  }
1094  }
1095  }
1096 }
1097 
1098 
1099 // ************************************************************************* //
virtual lduInterfacePtrsList interfaces() const =0
Return a list of pointers for each patch.
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
static autoPtr< GAMGInterface > New(const label index, const lduInterfacePtrsList &coarseInterfaces, const lduInterface &fineInterface, const labelField &localRestrictAddressing, const labelField &neighbourRestrictAddressing, const label fineLevelIndex, const label coarseComm)
Return a pointer to a new interface created on freestore given.
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
Vector-tensor class used to perform translations, rotations and scaling operations in 3D space...
Definition: transformer.H:83
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
virtual const labelUList & faceCells() const =0
Return faceCell addressing.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
void addInterfaces(lduInterfacePtrsList &interfaces, const lduSchedule &ps)
Add interfaces to a mesh constructed without.
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
#define forAllConstIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:477
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
To & refCast(From &r)
Reference type cast template function.
Definition: typeInfo.H:106
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Abstract base class for meshes which provide LDU addressing for the construction of lduMatrix and LDU...
Definition: lduMesh.H:59
static int & msgType()
Message tag of standard messages.
Definition: UPstream.H:476
virtual label comm() const =0
Return communicator used for parallel communication.
labelList findIndices(const ListType &, typename ListType::const_reference, const label start=0)
Find all occurrences of given element. Linear search.
fvMesh & mesh
T & first()
Return the first element of the list.
Definition: UListI.H:114
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
Less operator for pairs of <processor><index>
Input inter-processor communications stream.
Definition: IPstream.H:50
bool insert(const edge &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:142
lduPrimitiveMesh(const label nCells, labelList &l, labelList &u, const label comm, bool reuse)
Construct from components but without interfaces. Add interfaces.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
A class for handling words, derived from string.
Definition: word.H:59
virtual const labelUList & upperAddr() const =0
Return upper addressing.
Simplest contrete lduMesh which stores the addressing needed by lduMatrix.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
virtual const lduAddressing & lduAddr() const =0
Return ldu addressing.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
List< labelPair > labelPairList
List of labelPairs.
Definition: labelPair.H:49
static const lduMesh & mesh(const lduMesh &mesh0, const PtrList< lduPrimitiveMesh > &otherMeshes, const label meshI)
Select either mesh0 (meshI is 0) or otherMeshes[meshI-1].
errorManip< error > abort(error &err)
Definition: errorManip.H:131
label readLabel(Istream &is)
Definition: label.H:64
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
bool set(const label) const
Is element set.
Definition: UPtrListI.H:78
An abstract base class for processor coupled interfaces.
Map from edge (expressed as its endpoints) to value.
Definition: EdgeMap.H:47
defineTypeNameAndDebug(combustionModel, 0)
static bool write(const commsTypes commsType, const int toProcNo, const char *buf, const std::streamsize bufSize, const int tag=UPstream::msgType(), const label communicator=0)
Write given buffer to given processor.
Definition: UOPwrite.C:34
GAMG agglomerated processor interface.
Output inter-processor communications stream.
Definition: OPstream.H:50
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
Abstract base class for GAMG agglomerated interfaces.
Definition: GAMGInterface.H:51
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
void setSize(const label)
Reset size of List.
Definition: List.C:281
#define WarningInFunction
Report a warning using Foam::Warning.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
A List with indirect addressing.
Definition: fvMatrix.H:106
virtual int neighbProcNo() const =0
Return neighbour processor number (rank in communicator)
An abstract base class for implicitly-coupled interfaces e.g. processor and cyclic patches...
Definition: lduInterface.H:53
label n
bool operator()(const label a, const label b)
static void gather(const label comm, const lduMesh &mesh, const labelList &procIDs, PtrList< lduPrimitiveMesh > &otherMeshes)
Gather meshes from other processors onto procIDs[0].
The class contains the addressing required by the lduMatrix: upper, lower and losort.
virtual int myProcNo() const =0
Return processor number (rank in communicator)
T & last()
Return the last element of the list.
Definition: UListI.H:128
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
procLess(const labelPairList &lst)
Namespace for OpenFOAM.
label size() const
Return number of equations.
virtual const lduSchedule & patchSchedule() const =0