decompositionMethod.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) 2011-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 InClass
25  decompositionMethod
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "decompositionMethod.H"
30 #include "globalIndex.H"
31 #include "syncTools.H"
32 #include "Tuple2.H"
33 #include "faceSet.H"
34 #include "regionSplit.H"
35 #include "localPointRegion.H"
36 #include "minData.H"
37 #include "FaceCellWave.H"
38 
43 
44 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
45 
46 namespace Foam
47 {
48  defineTypeNameAndDebug(decompositionMethod, 0);
49  defineRunTimeSelectionTable(decompositionMethod, dictionary);
50 }
51 
52 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
53 
55 (
56  const dictionary& decompositionDict
57 )
58 :
59  decompositionDict_(decompositionDict),
60  nProcessors_
61  (
62  decompositionDict.lookup<label>("numberOfSubdomains")
63  )
64 {
65  // Read any constraints
66  wordList constraintTypes_;
67  if (decompositionDict_.found("constraints"))
68  {
69  // PtrList<dictionary> constraintsList
70  //(
71  // decompositionDict_.lookup("constraints")
72  //);
73  // forAll(constraintsList, i)
74  //{
75  // const dictionary& dict = constraintsList[i];
76  const dictionary& constraintsList = decompositionDict_.subDict
77  (
78  "constraints"
79  );
80  forAllConstIter(dictionary, constraintsList, iter)
81  {
82  const dictionary& dict = iter().dict();
83 
84  constraintTypes_.append(dict.lookup("type"));
85 
86  constraints_.append
87  (
89  (
90  dict,
91  constraintTypes_.last()
92  )
93  );
94  }
95  }
96 
97  // Backwards compatibility
98  if
99  (
100  decompositionDict_.found("preserveBaffles")
101  && findIndex
102  (
103  constraintTypes_,
104  decompositionConstraints::preserveBafflesConstraint::typeName
105  ) == -1
106  )
107  {
108  constraints_.append
109  (
111  );
112  }
113 
114  if
115  (
116  decompositionDict_.found("preservePatches")
117  && findIndex
118  (
119  constraintTypes_,
120  decompositionConstraints::preservePatchesConstraint::typeName
121  ) == -1
122  )
123  {
124  const wordReList pNames(decompositionDict_.lookup("preservePatches"));
125 
126  constraints_.append
127  (
129  );
130  }
131 
132  if
133  (
134  decompositionDict_.found("preserveFaceZones")
135  && findIndex
136  (
137  constraintTypes_,
138  decompositionConstraints::preserveFaceZonesConstraint::typeName
139  ) == -1
140  )
141  {
142  const wordReList zNames(decompositionDict_.lookup("preserveFaceZones"));
143 
144  constraints_.append
145  (
147  );
148  }
149 
150  if
151  (
152  decompositionDict_.found("singleProcessorFaceSets")
153  && findIndex
154  (
155  constraintTypes_,
156  decompositionConstraints::preserveFaceZonesConstraint::typeName
157  ) == -1
158  )
159  {
160  const List<Tuple2<word, label>> zNameAndProcs
161  (
162  decompositionDict_.lookup("singleProcessorFaceSets")
163  );
164 
165  constraints_.append
166  (
168  (
169  zNameAndProcs
170  )
171  );
172  }
173 }
174 
175 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
176 
178 (
179  const dictionary& decompositionDict
180 )
181 {
182  const word methodType(decompositionDict.lookup("method"));
183 
184  Pout<< "Selecting decompositionMethod " << methodType << endl;
185 
186  dictionaryConstructorTable::iterator cstrIter =
187  dictionaryConstructorTablePtr_->find(methodType);
188 
189  if (cstrIter == dictionaryConstructorTablePtr_->end())
190  {
192  << "Unknown decompositionMethod "
193  << methodType << nl << nl
194  << "Valid decompositionMethods are : " << endl
195  << dictionaryConstructorTablePtr_->sortedToc()
196  << exit(FatalError);
197  }
198 
199  return autoPtr<decompositionMethod>(cstrIter()(decompositionDict));
200 }
201 
202 
204 (
205  const polyMesh& mesh,
206  const pointField& points
207 )
208 {
209  scalarField weights(points.size(), 1.0);
210 
211  return decompose(mesh, points, weights);
212 }
213 
214 
216 (
217  const polyMesh& mesh,
218  const labelList& fineToCoarse,
219  const pointField& coarsePoints,
220  const scalarField& coarseWeights
221 )
222 {
223  CompactListList<label> coarseCellCells;
224  calcCellCells
225  (
226  mesh,
227  fineToCoarse,
228  coarsePoints.size(),
229  true, // use global cell labels
230  coarseCellCells
231  );
232 
233  // Decompose based on agglomerated points
234  labelList coarseDistribution
235  (
236  decompose
237  (
238  coarseCellCells(),
239  coarsePoints,
240  coarseWeights
241  )
242  );
243 
244  // Rework back into decomposition for original mesh_
245  labelList fineDistribution(fineToCoarse.size());
246 
247  forAll(fineDistribution, i)
248  {
249  fineDistribution[i] = coarseDistribution[fineToCoarse[i]];
250  }
251 
252  return fineDistribution;
253 }
254 
255 
257 (
258  const polyMesh& mesh,
259  const labelList& fineToCoarse,
260  const pointField& coarsePoints
261 )
262 {
263  scalarField cWeights(coarsePoints.size(), 1.0);
264 
265  return decompose
266  (
267  mesh,
268  fineToCoarse,
269  coarsePoints,
270  cWeights
271  );
272 }
273 
274 
276 (
277  const labelListList& globalCellCells,
278  const pointField& cc
279 )
280 {
281  scalarField cWeights(cc.size(), 1.0);
282 
283  return decompose(globalCellCells, cc, cWeights);
284 }
285 
286 
288 (
289  const polyMesh& mesh,
290  const labelList& agglom,
291  const label nLocalCoarse,
292  const bool parallel,
293  CompactListList<label>& cellCells
294 )
295 {
296  const labelList& faceOwner = mesh.faceOwner();
297  const labelList& faceNeighbour = mesh.faceNeighbour();
298  const polyBoundaryMesh& patches = mesh.boundaryMesh();
299 
300 
301  // Create global cell numbers
302  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
303 
304  globalIndex globalAgglom
305  (
306  nLocalCoarse,
309  parallel
310  );
311 
312 
313  // Get agglomerate owner on other side of coupled faces
314  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
315 
316  labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
317 
318  forAll(patches, patchi)
319  {
320  const polyPatch& pp = patches[patchi];
321 
322  if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
323  {
324  label facei = pp.start();
325  label bFacei = pp.start() - mesh.nInternalFaces();
326 
327  forAll(pp, i)
328  {
329  globalNeighbour[bFacei] = globalAgglom.toGlobal
330  (
331  agglom[faceOwner[facei]]
332  );
333 
334  bFacei++;
335  facei++;
336  }
337  }
338  }
339 
340  // Get the cell on the other side of coupled patches
341  syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
342 
343 
344  // Count number of faces (internal + coupled)
345  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
346 
347  // Number of faces per coarse cell
348  labelList nFacesPerCell(nLocalCoarse, 0);
349 
350  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
351  {
352  label own = agglom[faceOwner[facei]];
353  label nei = agglom[faceNeighbour[facei]];
354 
355  nFacesPerCell[own]++;
356  nFacesPerCell[nei]++;
357  }
358 
359  forAll(patches, patchi)
360  {
361  const polyPatch& pp = patches[patchi];
362 
363  if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
364  {
365  label facei = pp.start();
366  label bFacei = pp.start()-mesh.nInternalFaces();
367 
368  forAll(pp, i)
369  {
370  label own = agglom[faceOwner[facei]];
371 
372  label globalNei = globalNeighbour[bFacei];
373  if
374  (
375  !globalAgglom.isLocal(globalNei)
376  || globalAgglom.toLocal(globalNei) != own
377  )
378  {
379  nFacesPerCell[own]++;
380  }
381 
382  facei++;
383  bFacei++;
384  }
385  }
386  }
387 
388 
389  // Fill in offset and data
390  // ~~~~~~~~~~~~~~~~~~~~~~~
391 
392  cellCells.setSize(nFacesPerCell);
393 
394  nFacesPerCell = 0;
395 
396  labelList& m = cellCells.m();
397  const labelList& offsets = cellCells.offsets();
398 
399  // For internal faces is just offsetted owner and neighbour
400  for (label facei = 0; facei < mesh.nInternalFaces(); facei++)
401  {
402  label own = agglom[faceOwner[facei]];
403  label nei = agglom[faceNeighbour[facei]];
404 
405  m[offsets[own] + nFacesPerCell[own]++] = globalAgglom.toGlobal(nei);
406  m[offsets[nei] + nFacesPerCell[nei]++] = globalAgglom.toGlobal(own);
407  }
408 
409  // For boundary faces is offsetted coupled neighbour
410  forAll(patches, patchi)
411  {
412  const polyPatch& pp = patches[patchi];
413 
414  if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
415  {
416  label facei = pp.start();
417  label bFacei = pp.start()-mesh.nInternalFaces();
418 
419  forAll(pp, i)
420  {
421  label own = agglom[faceOwner[facei]];
422 
423  label globalNei = globalNeighbour[bFacei];
424 
425  if
426  (
427  !globalAgglom.isLocal(globalNei)
428  || globalAgglom.toLocal(globalNei) != own
429  )
430  {
431  m[offsets[own] + nFacesPerCell[own]++] = globalNei;
432  }
433 
434  facei++;
435  bFacei++;
436  }
437  }
438  }
439 
440 
441  // Check for duplicates connections between cells
442  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
443  // Done as postprocessing step since we now have cellCells.
444  label newIndex = 0;
445  labelHashSet nbrCells;
446 
447 
448  if (cellCells.size() == 0)
449  {
450  return;
451  }
452 
453  label startIndex = cellCells.offsets()[0];
454 
455  forAll(cellCells, celli)
456  {
457  nbrCells.clear();
458  nbrCells.insert(globalAgglom.toGlobal(celli));
459 
460  label endIndex = cellCells.offsets()[celli+1];
461 
462  for (label i = startIndex; i < endIndex; i++)
463  {
464  if (nbrCells.insert(cellCells.m()[i]))
465  {
466  cellCells.m()[newIndex++] = cellCells.m()[i];
467  }
468  }
469  startIndex = endIndex;
470  cellCells.offsets()[celli+1] = newIndex;
471  }
472 
473  cellCells.m().setSize(newIndex);
474 }
475 
476 
478 (
479  const polyMesh& mesh,
480  const labelList& agglom,
481  const label nLocalCoarse,
482  const bool parallel,
483  CompactListList<label>& cellCells,
484  CompactListList<scalar>& cellCellWeights
485 )
486 {
487  const labelList& faceOwner = mesh.faceOwner();
488  const labelList& faceNeighbour = mesh.faceNeighbour();
489  const polyBoundaryMesh& patches = mesh.boundaryMesh();
490 
491 
492  // Create global cell numbers
493  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
494 
495  globalIndex globalAgglom
496  (
497  nLocalCoarse,
500  parallel
501  );
502 
503 
504  // Get agglomerate owner on other side of coupled faces
505  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
506 
507  labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
508 
509  forAll(patches, patchi)
510  {
511  const polyPatch& pp = patches[patchi];
512 
513  if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
514  {
515  label faceI = pp.start();
516  label bFaceI = pp.start() - mesh.nInternalFaces();
517 
518  forAll(pp, i)
519  {
520  globalNeighbour[bFaceI] = globalAgglom.toGlobal
521  (
522  agglom[faceOwner[faceI]]
523  );
524 
525  bFaceI++;
526  faceI++;
527  }
528  }
529  }
530 
531  // Get the cell on the other side of coupled patches
532  syncTools::swapBoundaryFaceList(mesh, globalNeighbour);
533 
534 
535  // Count number of faces (internal + coupled)
536  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
537 
538  // Number of faces per coarse cell
539  labelList nFacesPerCell(nLocalCoarse, 0);
540 
541  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
542  {
543  label own = agglom[faceOwner[faceI]];
544  label nei = agglom[faceNeighbour[faceI]];
545 
546  nFacesPerCell[own]++;
547  nFacesPerCell[nei]++;
548  }
549 
550  forAll(patches, patchi)
551  {
552  const polyPatch& pp = patches[patchi];
553 
554  if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
555  {
556  label faceI = pp.start();
557  label bFaceI = pp.start()-mesh.nInternalFaces();
558 
559  forAll(pp, i)
560  {
561  label own = agglom[faceOwner[faceI]];
562 
563  label globalNei = globalNeighbour[bFaceI];
564  if
565  (
566  !globalAgglom.isLocal(globalNei)
567  || globalAgglom.toLocal(globalNei) != own
568  )
569  {
570  nFacesPerCell[own]++;
571  }
572 
573  faceI++;
574  bFaceI++;
575  }
576  }
577  }
578 
579 
580  // Fill in offset and data
581  // ~~~~~~~~~~~~~~~~~~~~~~~
582 
583  cellCells.setSize(nFacesPerCell);
584  cellCellWeights.setSize(nFacesPerCell);
585 
586  nFacesPerCell = 0;
587 
588  labelList& m = cellCells.m();
589  scalarList& w = cellCellWeights.m();
590  const labelList& offsets = cellCells.offsets();
591 
592  // For internal faces is just offsetted owner and neighbour
593  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
594  {
595  label own = agglom[faceOwner[faceI]];
596  label nei = agglom[faceNeighbour[faceI]];
597 
598  label ownIndex = offsets[own] + nFacesPerCell[own]++;
599  label neiIndex = offsets[nei] + nFacesPerCell[nei]++;
600 
601  m[ownIndex] = globalAgglom.toGlobal(nei);
602  w[ownIndex] = mesh.magFaceAreas()[faceI];
603  m[neiIndex] = globalAgglom.toGlobal(own);
604  w[ownIndex] = mesh.magFaceAreas()[faceI];
605  }
606 
607  // For boundary faces is offsetted coupled neighbour
608  forAll(patches, patchi)
609  {
610  const polyPatch& pp = patches[patchi];
611 
612  if (pp.coupled() && (parallel || !isA<processorPolyPatch>(pp)))
613  {
614  label faceI = pp.start();
615  label bFaceI = pp.start()-mesh.nInternalFaces();
616 
617  forAll(pp, i)
618  {
619  label own = agglom[faceOwner[faceI]];
620 
621  label globalNei = globalNeighbour[bFaceI];
622 
623  if
624  (
625  !globalAgglom.isLocal(globalNei)
626  || globalAgglom.toLocal(globalNei) != own
627  )
628  {
629  label ownIndex = offsets[own] + nFacesPerCell[own]++;
630  m[ownIndex] = globalNei;
631  w[ownIndex] = mesh.magFaceAreas()[faceI];
632  }
633 
634  faceI++;
635  bFaceI++;
636  }
637  }
638  }
639 
640 
641  // Check for duplicates connections between cells
642  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
643  // Done as postprocessing step since we now have cellCells.
644  label newIndex = 0;
645  labelHashSet nbrCells;
646 
647 
648  if (cellCells.size() == 0)
649  {
650  return;
651  }
652 
653  label startIndex = cellCells.offsets()[0];
654 
655  forAll(cellCells, cellI)
656  {
657  nbrCells.clear();
658  nbrCells.insert(globalAgglom.toGlobal(cellI));
659 
660  label endIndex = cellCells.offsets()[cellI+1];
661 
662  for (label i = startIndex; i < endIndex; i++)
663  {
664  if (nbrCells.insert(cellCells.m()[i]))
665  {
666  cellCells.m()[newIndex] = cellCells.m()[i];
667  cellCellWeights.m()[newIndex] = cellCellWeights.m()[i];
668  newIndex++;
669  }
670  }
671  startIndex = endIndex;
672  cellCells.offsets()[cellI+1] = newIndex;
673  cellCellWeights.offsets()[cellI+1] = newIndex;
674  }
675 
676  cellCells.m().setSize(newIndex);
677  cellCellWeights.m().setSize(newIndex);
678 }
679 
680 
682 (
683  const polyMesh& mesh,
684  const scalarField& cellWeights,
685 
686  //- Whether owner and neighbour should be on same processor
687  // (takes priority over explicitConnections)
688  const boolList& blockedFace,
689 
690  //- Whether whole sets of faces (and point neighbours) need to be kept
691  // on single processor
692  const PtrList<labelList>& specifiedProcessorFaces,
693  const labelList& specifiedProcessor,
694 
695  //- Additional connections between boundary faces
696  const List<labelPair>& explicitConnections
697 )
698 {
699  // Any weights specified?
700  label nWeights = returnReduce(cellWeights.size(), sumOp<label>());
701 
702  if (nWeights > 0 && cellWeights.size() != mesh.nCells())
703  {
705  << "Number of weights " << cellWeights.size()
706  << " differs from number of cells " << mesh.nCells()
707  << exit(FatalError);
708  }
709 
710 
711  // Any processor sets?
712  label nProcSets = 0;
713  forAll(specifiedProcessorFaces, setI)
714  {
715  nProcSets += specifiedProcessorFaces[setI].size();
716  }
717  reduce(nProcSets, sumOp<label>());
718 
719  // Any non-mesh connections?
720  label nConnections = returnReduce
721  (
722  explicitConnections.size(),
723  sumOp<label>()
724  );
725 
726  // Any faces not blocked?
727  label nUnblocked = 0;
728  forAll(blockedFace, facei)
729  {
730  if (!blockedFace[facei])
731  {
732  nUnblocked++;
733  }
734  }
735  reduce(nUnblocked, sumOp<label>());
736 
737 
738 
739  // Either do decomposition on cell centres or on agglomeration
740 
741  labelList finalDecomp;
742 
743 
744  if (nProcSets+nConnections+nUnblocked == 0)
745  {
746  // No constraints, possibly weights
747 
748  if (nWeights > 0)
749  {
750  finalDecomp = decompose
751  (
752  mesh,
753  mesh.cellCentres(),
754  cellWeights
755  );
756  }
757  else
758  {
759  finalDecomp = decompose(mesh, mesh.cellCentres());
760  }
761  }
762  else
763  {
764  if (debug)
765  {
766  Info<< "Constrained decomposition:" << endl
767  << " faces with same owner and neighbour processor : "
768  << nUnblocked << endl
769  << " baffle faces with same owner processor : "
770  << nConnections << endl
771  << " faces all on same processor : "
772  << nProcSets << endl << endl;
773  }
774 
775  // Determine local regions, separated by blockedFaces
776  regionSplit localRegion(mesh, blockedFace, explicitConnections, false);
777 
778 
779  if (debug)
780  {
781  Info<< "Constrained decomposition:" << endl
782  << " split into " << localRegion.nLocalRegions()
783  << " regions."
784  << endl;
785  }
786 
787  // Determine region cell centres
788  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
789 
790  // This just takes the first cell in the region. Otherwise the problem
791  // is with cyclics - if we'd average the region centre might be
792  // somewhere in the middle of the domain which might not be anywhere
793  // near any of the cells.
794 
795  pointField regionCentres(localRegion.nLocalRegions(), point::max);
796 
797  forAll(localRegion, celli)
798  {
799  label regionI = localRegion[celli];
800 
801  if (regionCentres[regionI] == point::max)
802  {
803  regionCentres[regionI] = mesh.cellCentres()[celli];
804  }
805  }
806 
807  // Do decomposition on agglomeration
808  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
809 
810  scalarField regionWeights(localRegion.nLocalRegions(), 0);
811 
812  if (nWeights > 0)
813  {
814  forAll(localRegion, celli)
815  {
816  label regionI = localRegion[celli];
817 
818  regionWeights[regionI] += cellWeights[celli];
819  }
820  }
821  else
822  {
823  forAll(localRegion, celli)
824  {
825  label regionI = localRegion[celli];
826 
827  regionWeights[regionI] += 1.0;
828  }
829  }
830 
831  finalDecomp = decompose
832  (
833  mesh,
834  localRegion,
835  regionCentres,
836  regionWeights
837  );
838 
839 
840 
841  // Implement the explicitConnections since above decompose
842  // does not know about them
843  forAll(explicitConnections, i)
844  {
845  const labelPair& baffle = explicitConnections[i];
846  label f0 = baffle.first();
847  label f1 = baffle.second();
848 
849  if (!blockedFace[f0] && !blockedFace[f1])
850  {
851  // Note: what if internal faces and owner and neighbour on
852  // different processor? So for now just push owner side
853  // proc
854 
855  const label proci = finalDecomp[mesh.faceOwner()[f0]];
856 
857  finalDecomp[mesh.faceOwner()[f1]] = proci;
858  if (mesh.isInternalFace(f1))
859  {
860  finalDecomp[mesh.faceNeighbour()[f1]] = proci;
861  }
862  }
863  else if (blockedFace[f0] != blockedFace[f1])
864  {
866  << "On explicit connection between faces " << f0
867  << " and " << f1
868  << " the two blockedFace status are not equal : "
869  << blockedFace[f0] << " and " << blockedFace[f1]
870  << exit(FatalError);
871  }
872  }
873 
874 
875  // blockedFaces corresponding to processor faces need to be handled
876  // separately since not handled by local regionSplit. We need to
877  // walk now across coupled faces and make sure to move a whole
878  // global region across
879  if (Pstream::parRun())
880  {
881  // Re-do regionSplit
882 
883  // Field on cells and faces.
884  List<minData> cellData(mesh.nCells());
885  List<minData> faceData(mesh.nFaces());
886 
887  // Take over blockedFaces by seeding a negative number
888  // (so is always less than the decomposition)
889  label nUnblocked = 0;
890  forAll(blockedFace, facei)
891  {
892  if (blockedFace[facei])
893  {
894  faceData[facei] = minData(-123);
895  }
896  else
897  {
898  nUnblocked++;
899  }
900  }
901 
902  // Seed unblocked faces with destination processor
903  labelList seedFaces(nUnblocked);
904  List<minData> seedData(nUnblocked);
905  nUnblocked = 0;
906 
907  forAll(blockedFace, facei)
908  {
909  if (!blockedFace[facei])
910  {
911  label own = mesh.faceOwner()[facei];
912  seedFaces[nUnblocked] = facei;
913  seedData[nUnblocked] = minData(finalDecomp[own]);
914  nUnblocked++;
915  }
916  }
917 
918 
919  // Propagate information inwards
920  FaceCellWave<minData> deltaCalc
921  (
922  mesh,
923  seedFaces,
924  seedData,
925  faceData,
926  cellData,
927  mesh.globalData().nTotalCells()+1
928  );
929 
930  // And extract
931  forAll(finalDecomp, celli)
932  {
933  if (cellData[celli].valid(deltaCalc.data()))
934  {
935  finalDecomp[celli] = cellData[celli].data();
936  }
937  }
938  }
939 
940 
941  // For specifiedProcessorFaces rework the cellToProc to enforce
942  // all on one processor since we can't guarantee that the input
943  // to regionSplit was a single region.
944  // E.g. faceSet 'a' with the cells split into two regions
945  // by a notch formed by two walls
946  //
947  // \ /
948  // \ /
949  // ---a----+-----a-----
950  //
951  //
952  // Note that reworking the cellToProc might make the decomposition
953  // unbalanced.
954  forAll(specifiedProcessorFaces, setI)
955  {
956  const labelList& set = specifiedProcessorFaces[setI];
957 
958  label proci = specifiedProcessor[setI];
959  if (proci == -1)
960  {
961  // If no processor specified use the one from the
962  // 0th element
963  proci = finalDecomp[mesh.faceOwner()[set[0]]];
964  }
965 
966  forAll(set, fI)
967  {
968  const face& f = mesh.faces()[set[fI]];
969  forAll(f, fp)
970  {
971  const labelList& pFaces = mesh.pointFaces()[f[fp]];
972  forAll(pFaces, i)
973  {
974  label facei = pFaces[i];
975 
976  finalDecomp[mesh.faceOwner()[facei]] = proci;
977  if (mesh.isInternalFace(facei))
978  {
979  finalDecomp[mesh.faceNeighbour()[facei]] = proci;
980  }
981  }
982  }
983  }
984  }
985 
986 
987  if (debug && Pstream::parRun())
988  {
989  labelList nbrDecomp;
990  syncTools::swapBoundaryCellList(mesh, finalDecomp, nbrDecomp);
991 
992  const polyBoundaryMesh& patches = mesh.boundaryMesh();
993  forAll(patches, patchi)
994  {
995  const polyPatch& pp = patches[patchi];
996  if (pp.coupled())
997  {
998  forAll(pp, i)
999  {
1000  label facei = pp.start()+i;
1001  label own = mesh.faceOwner()[facei];
1002  label bFacei = facei-mesh.nInternalFaces();
1003 
1004  if (!blockedFace[facei])
1005  {
1006  label ownProc = finalDecomp[own];
1007  label nbrProc = nbrDecomp[bFacei];
1008  if (ownProc != nbrProc)
1009  {
1011  << "patch:" << pp.name()
1012  << " face:" << facei
1013  << " at:" << mesh.faceCentres()[facei]
1014  << " ownProc:" << ownProc
1015  << " nbrProc:" << nbrProc
1016  << exit(FatalError);
1017  }
1018  }
1019  }
1020  }
1021  }
1022  }
1023  }
1024 
1025  return finalDecomp;
1026 }
1027 
1028 
1031  const polyMesh& mesh,
1032  boolList& blockedFace,
1033  PtrList<labelList>& specifiedProcessorFaces,
1034  labelList& specifiedProcessor,
1035  List<labelPair>& explicitConnections
1036 )
1037 {
1038  blockedFace.setSize(mesh.nFaces());
1039  blockedFace = true;
1040 
1041  specifiedProcessorFaces.clear();
1042  explicitConnections.clear();
1043 
1044  forAll(constraints_, constraintI)
1045  {
1046  constraints_[constraintI].add
1047  (
1048  mesh,
1049  blockedFace,
1050  specifiedProcessorFaces,
1051  specifiedProcessor,
1052  explicitConnections
1053  );
1054  }
1055 }
1056 
1057 
1060  const polyMesh& mesh,
1061  const boolList& blockedFace,
1062  const PtrList<labelList>& specifiedProcessorFaces,
1063  const labelList& specifiedProcessor,
1064  const List<labelPair>& explicitConnections,
1065  labelList& decomposition
1066 )
1067 {
1068  forAll(constraints_, constraintI)
1069  {
1070  constraints_[constraintI].apply
1071  (
1072  mesh,
1073  blockedFace,
1074  specifiedProcessorFaces,
1075  specifiedProcessor,
1076  explicitConnections,
1077  decomposition
1078  );
1079  }
1080 }
1081 
1082 
1085  const polyMesh& mesh,
1086  const scalarField& cellWeights
1087 )
1088 {
1089  // Collect all constraints
1090 
1091  boolList blockedFace;
1092  PtrList<labelList> specifiedProcessorFaces;
1093  labelList specifiedProcessor;
1094  List<labelPair> explicitConnections;
1095  setConstraints
1096  (
1097  mesh,
1098  blockedFace,
1099  specifiedProcessorFaces,
1100  specifiedProcessor,
1101  explicitConnections
1102  );
1103 
1104 
1105  // Construct decomposition method and either do decomposition on
1106  // cell centres or on agglomeration
1107 
1108  labelList finalDecomp = decompose
1109  (
1110  mesh,
1111  cellWeights, // optional weights
1112  blockedFace, // any cells to be combined
1113  specifiedProcessorFaces,// any whole cluster of cells to be kept
1114  specifiedProcessor,
1115  explicitConnections // baffles
1116  );
1117 
1118 
1119  // Give any constraint the option of modifying the decomposition
1120 
1121  applyConstraints
1122  (
1123  mesh,
1124  blockedFace,
1125  specifiedProcessorFaces,
1126  specifiedProcessor,
1127  explicitConnections,
1128  finalDecomp
1129  );
1130 
1131  return finalDecomp;
1132 }
1133 
1134 
1135 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
bool isLocal(const label i) const
Is on local processor.
Definition: globalIndexI.H:95
virtual labelList decompose(const pointField &points, const scalarField &pointWeights)
Return for every coordinate the wanted processor number.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition: regionSplit.H:119
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
const word & name() const
Return name.
label nLocalRegions() const
Return local number of regions.
Definition: regionSplit.H:195
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
label nInternalFaces() const
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1175
label nFaces() const
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
const List< T > & m() const
Return the packed matrix of data.
label nCells() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Wave propagation of information through grid. Every iteration information goes through one layer of c...
Definition: FaceCellWave.H:76
label nTotalCells() const
Return total number of cells in decomposed mesh.
static int & msgType()
Message tag of standard messages.
Definition: UPstream.H:476
Constraint to keep/move owner and neighbour of faceZone onto same processor.
patches[0]
static label worldComm
Default communicator (all processors)
Definition: UPstream.H:278
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:111
scalar f1
Definition: createFields.H:15
static void calcCellCells(const polyMesh &mesh, const labelList &agglom, const label nLocalCoarse, const bool global, CompactListList< label > &cellCells)
Helper: determine (local or global) cellCells from mesh.
const dictionary & subDict(const word &) const
Find and return a sub-dictionary.
Definition: dictionary.C:982
const scalarField & magFaceAreas() const
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
Detects baffles and keeps owner and neighbour on same processor.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:319
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:125
A class for handling words, derived from string.
Definition: word.H:59
void clear()
Clear all entries from table.
Definition: HashTable.C:468
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:178
void setConstraints(const polyMesh &mesh, boolList &blockedFace, PtrList< labelList > &specifiedProcessorFaces, labelList &specifiedProcessor, List< labelPair > &explicitConnections)
Helper: extract constraints:
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
label size() const
Return the primary size, i.e. the number of rows.
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1394
label toLocal(const label i) const
From global to local on current processor.
Definition: globalIndexI.H:117
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
const vectorField & cellCentres() const
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
decompositionMethod(const dictionary &decompositionDict)
Construct given the decomposition dictionary.
static autoPtr< decompositionMethod > New(const dictionary &decompositionDict)
Return a reference to the selected decomposition method.
Foam::polyBoundaryMesh.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints, -1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){ pointMap[i]=i;}for(label i=0;i< nPoints;i++){ if(f[i] > 0.0) { hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei) { if(edges[ei].mag(points)< small) { label start=pointMap[edges[ei].start()];while(start !=pointMap[start]) { start=pointMap[start];} label end=pointMap[edges[ei].end()];while(end !=pointMap[end]) { end=pointMap[end];} label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;} } cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){ cellShape &cs=cellShapes[celli];forAll(cs, i) { cs[i]=pointMap[cs[i]];} cs.collapse();}label bcIDs[11]={-1, 0, 2, 4, -1, 5, -1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={ &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &wallPolyPatch::typeName, &symmetryPolyPatch::typeName, &wedgePolyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, &mergedCyclicPolyPatch::typeName};enum patchTypeNames{ PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={ "piston", "valve", "liner", "cylinderHead", "axis", "wedge", "inflow", "outflow", "presin", "presout", "symmetryPlane", "cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
A packed storage unstructured matrix of objects of type <T> using an offset table for access...
void setSize(const label mRows)
Reset size of CompactListList.
static const char nl
Definition: Ostream.H:260
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
defineTypeNameAndDebug(combustionModel, 0)
labelList f(nPoints)
static autoPtr< decompositionConstraint > New(const dictionary &constraintsDict, const word &type)
Return a reference to the selected decompositionConstraint.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
const Type & second() const
Return second.
Definition: Pair.H:110
Constraint to keep owner and neighbour of (cyclic) patch on same processor.
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
const vectorField & faceCentres() const
void setSize(const label)
Reset size of List.
Definition: List.C:281
label toGlobal(const label i) const
From local to global.
Definition: globalIndexI.H:82
T * data()
Return a pointer to the first data element,.
Definition: UListI.H:149
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
void applyConstraints(const polyMesh &mesh, const boolList &blockedFace, const PtrList< labelList > &specifiedProcessorFaces, const labelList &specifiedProcessor, const List< labelPair > &explicitConnections, labelList &finalDecomp)
Helper: apply constraints to a decomposition. This gives.
label patchi
const List< label > & offsets() const
Return the offset table (= size()+1)
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
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:309
void clear()
Clear the PtrList, i.e. set size to zero deleting all the.
Definition: PtrList.C:174
messageStream Info
const labelListList & pointFaces() const
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Constraint to keep all cells connected to face or point of faceSet on a single processor.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:430
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
T & last()
Return the last element of the list.
Definition: UListI.H:128
For use with FaceCellWave. Transports minimum passive data.
Definition: minData.H:61
const Type & first() const
Return first.
Definition: Pair.H:98
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844