FreeFOAM The Cross-Platform CFD Toolkit
parMetisDecomp.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 | Copyright (C) 1991-2010 OpenCFD Ltd.
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 "parMetisDecomp.H"
29 #include <OpenFOAM/syncTools.H>
31 #include <OpenFOAM/floatScalar.H>
32 #include <OpenFOAM/polyMesh.H>
33 #include <OpenFOAM/Time.H>
34 #include <OpenFOAM/labelIOField.H>
35 #include <OpenFOAM/globalIndex.H>
37 
38 #include <mpi.h>
39 
40 extern "C"
41 {
42 # include <parmetis.h>
43 }
44 
45 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
46 
47 namespace Foam
48 {
49  defineTypeNameAndDebug(parMetisDecomp, 0);
50 
52  (
53  decompositionMethod,
54  parMetisDecomp,
55  dictionaryMesh
56  );
57 }
58 
59 
60 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
61 
62 //- Does prevention of 0 cell domains and calls parmetis.
63 Foam::label Foam::parMetisDecomp::decompose
64 (
65  Field<int>& xadj,
66  Field<int>& adjncy,
67  const pointField& cellCentres,
68  Field<int>& cellWeights,
69  Field<int>& faceWeights,
70  const List<int>& options,
71 
72  List<int>& finalDecomp
73 )
74 {
75  // C style numbering
76  int numFlag = 0;
77 
78  // Number of dimensions
79  int nDims = 3;
80 
81 
82  if (cellCentres.size() != xadj.size()-1)
83  {
84  FatalErrorIn("parMetisDecomp::decompose(..)")
85  << "cellCentres:" << cellCentres.size()
86  << " xadj:" << xadj.size()
87  << abort(FatalError);
88  }
89 
90 
91  // Get number of cells on all processors
92  List<int> nLocalCells(Pstream::nProcs());
93  nLocalCells[Pstream::myProcNo()] = xadj.size()-1;
94  Pstream::gatherList(nLocalCells);
95  Pstream::scatterList(nLocalCells);
96 
97  // Get cell offsets.
98  List<int> cellOffsets(Pstream::nProcs()+1);
99  int nGlobalCells = 0;
100  forAll(nLocalCells, procI)
101  {
102  cellOffsets[procI] = nGlobalCells;
103  nGlobalCells += nLocalCells[procI];
104  }
105  cellOffsets[Pstream::nProcs()] = nGlobalCells;
106 
107  // Convert pointField into float
108  Field<floatScalar> xyz(3*cellCentres.size());
109  int compI = 0;
110  forAll(cellCentres, cellI)
111  {
112  const point& cc = cellCentres[cellI];
113  xyz[compI++] = float(cc.x());
114  xyz[compI++] = float(cc.y());
115  xyz[compI++] = float(cc.z());
116  }
117 
118  // Make sure every domain has at least one cell
119  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
120  // (Metis falls over with zero sized domains)
121  // Trickle cells from processors that have them up to those that
122  // don't.
123 
124 
125  // Number of cells to send to the next processor
126  // (is same as number of cells next processor has to receive)
127  List<int> nSendCells(Pstream::nProcs(), 0);
128 
129  for (label procI = nLocalCells.size()-1; procI >=1; procI--)
130  {
131  if (nLocalCells[procI]-nSendCells[procI] < 1)
132  {
133  nSendCells[procI-1] = nSendCells[procI]-nLocalCells[procI]+1;
134  }
135  }
136 
137  // First receive (so increasing the sizes of all arrays)
138 
139  if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
140  {
141  // Receive cells from previous processor
142  IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
143 
144  Field<int> prevXadj(fromPrevProc);
145  Field<int> prevAdjncy(fromPrevProc);
146  Field<floatScalar> prevXyz(fromPrevProc);
147  Field<int> prevCellWeights(fromPrevProc);
148  Field<int> prevFaceWeights(fromPrevProc);
149 
150  if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1])
151  {
152  FatalErrorIn("parMetisDecomp::decompose(..)")
153  << "Expected from processor " << Pstream::myProcNo()-1
154  << " connectivity for " << nSendCells[Pstream::myProcNo()-1]
155  << " nCells but only received " << prevXadj.size()
156  << abort(FatalError);
157  }
158 
159  // Insert adjncy
160  prepend(prevAdjncy, adjncy);
161  // Adapt offsets and prepend xadj
162  xadj += prevAdjncy.size();
163  prepend(prevXadj, xadj);
164  // Coords
165  prepend(prevXyz, xyz);
166  // Weights
167  prepend(prevCellWeights, cellWeights);
168  prepend(prevFaceWeights, faceWeights);
169  }
170 
171 
172  // Send to my next processor
173 
174  if (nSendCells[Pstream::myProcNo()] > 0)
175  {
176  // Send cells to next processor
177  OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1);
178 
179  int nCells = nSendCells[Pstream::myProcNo()];
180  int startCell = xadj.size()-1 - nCells;
181  int startFace = xadj[startCell];
182  int nFaces = adjncy.size()-startFace;
183 
184  // Send for all cell data: last nCells elements
185  // Send for all face data: last nFaces elements
186  toNextProc
187  << Field<int>::subField(xadj, nCells, startCell)-startFace
188  << Field<int>::subField(adjncy, nFaces, startFace)
189  << SubField<floatScalar>(xyz, nDims*nCells, nDims*startCell)
190  <<
191  (
192  cellWeights.size()
193  ? static_cast<const Field<int>&>
194  (
195  Field<int>::subField(cellWeights, nCells, startCell)
196  )
197  : Field<int>(0)
198  )
199  <<
200  (
201  faceWeights.size()
202  ? static_cast<const Field<int>&>
203  (
204  Field<int>::subField(faceWeights, nFaces, startFace)
205  )
206  : Field<int>(0)
207  );
208 
209  // Remove data that has been sent
210  if (faceWeights.size())
211  {
212  faceWeights.setSize(faceWeights.size()-nFaces);
213  }
214  if (cellWeights.size())
215  {
216  cellWeights.setSize(cellWeights.size()-nCells);
217  }
218  xyz.setSize(xyz.size()-nDims*nCells);
219  adjncy.setSize(adjncy.size()-nFaces);
220  xadj.setSize(xadj.size() - nCells);
221  }
222 
223 
224 
225  // Adapt number of cells
226  forAll(nSendCells, procI)
227  {
228  // Sent cells
229  nLocalCells[procI] -= nSendCells[procI];
230 
231  if (procI >= 1)
232  {
233  // Received cells
234  nLocalCells[procI] += nSendCells[procI-1];
235  }
236  }
237  // Adapt cellOffsets
238  nGlobalCells = 0;
239  forAll(nLocalCells, procI)
240  {
241  cellOffsets[procI] = nGlobalCells;
242  nGlobalCells += nLocalCells[procI];
243  }
244 
245 
246  if (nLocalCells[Pstream::myProcNo()] != (xadj.size()-1))
247  {
248  FatalErrorIn("parMetisDecomp::decompose(..)")
249  << "Have connectivity for " << xadj.size()-1
250  << " cells but nLocalCells:" << nLocalCells[Pstream::myProcNo()]
251  << abort(FatalError);
252  }
253 
254  // Weight info
255  int wgtFlag = 0;
256  int* vwgtPtr = NULL;
257  int* adjwgtPtr = NULL;
258 
259  if (cellWeights.size())
260  {
261  vwgtPtr = cellWeights.begin();
262  wgtFlag += 2; // Weights on vertices
263  }
264  if (faceWeights.size())
265  {
266  adjwgtPtr = faceWeights.begin();
267  wgtFlag += 1; // Weights on edges
268  }
269 
270 
271  // Number of weights or balance constraints
272  int nCon = 1;
273  // Per processor, per constraint the weight
274  Field<floatScalar> tpwgts(nCon*nProcessors_, 1./nProcessors_);
275  // Imbalance tolerance
276  Field<floatScalar> ubvec(nCon, 1.02);
277  if (nProcessors_ == 1)
278  {
279  // If only one processor there is no imbalance.
280  ubvec[0] = 1;
281  }
282 
283  MPI_Comm comm = MPI_COMM_WORLD;
284 
285  // output: cell -> processor addressing
286  finalDecomp.setSize(nLocalCells[Pstream::myProcNo()]);
287 
288  // output: number of cut edges
289  int edgeCut = 0;
290 
291 
292  ParMETIS_V3_PartGeomKway
293  (
294  cellOffsets.begin(), // vtxDist
295  xadj.begin(),
296  adjncy.begin(),
297  vwgtPtr, // vertexweights
298  adjwgtPtr, // edgeweights
299  &wgtFlag,
300  &numFlag,
301  &nDims,
302  xyz.begin(),
303  &nCon,
304  &nProcessors_, // nParts
305  tpwgts.begin(),
306  ubvec.begin(),
307  const_cast<List<int>&>(options).begin(),
308  &edgeCut,
309  finalDecomp.begin(),
310  &comm
311  );
312 
313 
314  // If we sent cells across make sure we undo it
315  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
316 
317  // Receive back from next processor if I sent something
318  if (nSendCells[Pstream::myProcNo()] > 0)
319  {
320  IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1);
321 
322  List<int> nextFinalDecomp(fromNextProc);
323 
324  if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()])
325  {
326  FatalErrorIn("parMetisDecomp::decompose(..)")
327  << "Expected from processor " << Pstream::myProcNo()+1
328  << " decomposition for " << nSendCells[Pstream::myProcNo()]
329  << " nCells but only received " << nextFinalDecomp.size()
330  << abort(FatalError);
331  }
332 
333  append(nextFinalDecomp, finalDecomp);
334  }
335 
336  // Send back to previous processor.
337  if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0)
338  {
339  OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1);
340 
341  int nToPrevious = nSendCells[Pstream::myProcNo()-1];
342 
343  toPrevProc <<
344  SubList<int>
345  (
346  finalDecomp,
347  nToPrevious,
348  finalDecomp.size()-nToPrevious
349  );
350 
351  // Remove locally what has been sent
352  finalDecomp.setSize(finalDecomp.size()-nToPrevious);
353  }
354 
355  return edgeCut;
356 }
357 
358 
359 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
360 
361 Foam::parMetisDecomp::parMetisDecomp
362 (
363  const dictionary& decompositionDict,
364  const polyMesh& mesh
365 )
366 :
367  decompositionMethod(decompositionDict),
368  mesh_(mesh)
369 {
370  if (Pstream::nProcs() <= 1)
371  {
372  if(!dlLibraryTable::open("metisDecomp.so"))
373  {
375  (
376  "parMetisDecomp::parMetisDecomp()"
377  ) << "Failed to load the library 'metisDecomp.so'" << endl
379  }
380  }
381 }
382 
383 
384 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
385 
386 Foam::labelList Foam::parMetisDecomp::decompose
387 (
388  const pointField& cc,
389  const scalarField& cWeights
390 )
391 {
392  if (cc.size() != mesh_.nCells())
393  {
395  (
396  "parMetisDecomp::decompose"
397  "(const pointField&, const scalarField&)"
398  ) << "Can use this decomposition method only for the whole mesh"
399  << endl
400  << "and supply one coordinate (cellCentre) for every cell." << endl
401  << "The number of coordinates " << cc.size() << endl
402  << "The number of cells in the mesh " << mesh_.nCells()
403  << exit(FatalError);
404  }
405 
406  // For running sequential ...
407  if (Pstream::nProcs() <= 1)
408  {
409  return metisDecomp(decompositionDict_, mesh_).decompose
410  (
411  cc,
412  cWeights
413  );
414  }
415 
416 
417  // Connections
418  Field<int> adjncy;
419  // Offsets into adjncy
420  Field<int> xadj;
421  calcMetisDistributedCSR
422  (
423  mesh_,
424  adjncy,
425  xadj
426  );
427 
428 
429  // decomposition options. 0 = use defaults
430  List<int> options(3, 0);
431  //options[0] = 1; // don't use defaults but use values below
432  //options[1] = -1; // full debug info
433  //options[2] = 15; // random number seed
434 
435  // cell weights (so on the vertices of the dual)
436  Field<int> cellWeights;
437 
438  // face weights (so on the edges of the dual)
439  Field<int> faceWeights;
440 
441 
442  // Check for externally provided cellweights and if so initialise weights
443  scalar minWeights = gMin(cWeights);
444  if (cWeights.size() > 0)
445  {
446  if (minWeights <= 0)
447  {
448  WarningIn
449  (
450  "metisDecomp::decompose"
451  "(const pointField&, const scalarField&)"
452  ) << "Illegal minimum weight " << minWeights
453  << endl;
454  }
455 
456  if (cWeights.size() != mesh_.nCells())
457  {
459  (
460  "parMetisDecomp::decompose"
461  "(const pointField&, const scalarField&)"
462  ) << "Number of cell weights " << cWeights.size()
463  << " does not equal number of cells " << mesh_.nCells()
464  << exit(FatalError);
465  }
466 
467  // Convert to integers.
468  cellWeights.setSize(cWeights.size());
469  forAll(cellWeights, i)
470  {
471  cellWeights[i] = int(cWeights[i]/minWeights);
472  }
473  }
474 
475 
476  // Check for user supplied weights and decomp options
477  if (decompositionDict_.found("metisCoeffs"))
478  {
479  const dictionary& metisCoeffs =
480  decompositionDict_.subDict("metisCoeffs");
481  word weightsFile;
482 
483  if (metisCoeffs.readIfPresent("cellWeightsFile", weightsFile))
484  {
485  Info<< "parMetisDecomp : Using cell-based weights read from "
486  << weightsFile << endl;
487 
488  labelIOField cellIOWeights
489  (
490  IOobject
491  (
492  weightsFile,
493  mesh_.time().timeName(),
494  mesh_,
497  )
498  );
499  cellWeights.transfer(cellIOWeights);
500 
501  if (cellWeights.size() != mesh_.nCells())
502  {
504  (
505  "parMetisDecomp::decompose"
506  "(const pointField&, const scalarField&)"
507  ) << "Number of cell weights " << cellWeights.size()
508  << " read from " << cellIOWeights.objectPath()
509  << " does not equal number of cells " << mesh_.nCells()
510  << exit(FatalError);
511  }
512  }
513 
514  if (metisCoeffs.readIfPresent("faceWeightsFile", weightsFile))
515  {
516  Info<< "parMetisDecomp : Using face-based weights read from "
517  << weightsFile << endl;
518 
519  labelIOField weights
520  (
521  IOobject
522  (
523  weightsFile,
524  mesh_.time().timeName(),
525  mesh_,
528  )
529  );
530 
531  if (weights.size() != mesh_.nFaces())
532  {
534  (
535  "parMetisDecomp::decompose"
536  "(const pointField&, const scalarField&)"
537  ) << "Number of face weights " << weights.size()
538  << " does not equal number of internal and boundary faces "
539  << mesh_.nFaces()
540  << exit(FatalError);
541  }
542 
543  faceWeights.setSize(adjncy.size());
544 
545  // Assume symmetric weights. Keep same ordering as adjncy.
546  List<int> nFacesPerCell(mesh_.nCells(), 0);
547 
548  // Handle internal faces
549  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
550  {
551  label w = weights[faceI];
552 
553  label own = mesh_.faceOwner()[faceI];
554  label nei = mesh_.faceNeighbour()[faceI];
555 
556  faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
557  faceWeights[xadj[nei] + nFacesPerCell[nei]++] = w;
558  }
559  // Coupled boundary faces
560  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
561 
562  forAll(patches, patchI)
563  {
564  const polyPatch& pp = patches[patchI];
565 
566  if (pp.coupled())
567  {
568  label faceI = pp.start();
569 
570  forAll(pp, i)
571  {
572  label w = weights[faceI];
573  label own = mesh_.faceOwner()[faceI];
574  faceWeights[xadj[own] + nFacesPerCell[own]++] = w;
575  faceI++;
576  }
577  }
578  }
579  }
580 
581  if (metisCoeffs.readIfPresent("options", options))
582  {
583  Info<< "Using Metis options " << options
584  << nl << endl;
585 
586  if (options.size() != 3)
587  {
589  (
590  "parMetisDecomp::decompose"
591  "(const pointField&, const scalarField&)"
592  ) << "Number of options " << options.size()
593  << " should be three." << exit(FatalError);
594  }
595  }
596  }
597 
598 
599  // Do actual decomposition
600  List<int> finalDecomp;
601  decompose
602  (
603  xadj,
604  adjncy,
605  cc,
606  cellWeights,
607  faceWeights,
608  options,
609 
610  finalDecomp
611  );
612 
613  // Copy back to labelList
614  labelList decomp(finalDecomp.size());
615  forAll(decomp, i)
616  {
617  decomp[i] = finalDecomp[i];
618  }
619  return decomp;
620 }
621 
622 
623 Foam::labelList Foam::parMetisDecomp::decompose
624 (
625  const labelList& cellToRegion,
626  const pointField& regionPoints,
627  const scalarField& regionWeights
628 )
629 {
630  const labelList& faceOwner = mesh_.faceOwner();
631  const labelList& faceNeighbour = mesh_.faceNeighbour();
632  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
633 
634  if (cellToRegion.size() != mesh_.nCells())
635  {
637  (
638  "parMetisDecomp::decompose(const labelList&, const pointField&)"
639  ) << "Size of cell-to-coarse map " << cellToRegion.size()
640  << " differs from number of cells in mesh " << mesh_.nCells()
641  << exit(FatalError);
642  }
643 
644 
645  // Global region numbering engine
646  globalIndex globalRegions(regionPoints.size());
647 
648 
649  // Get renumbered owner region on other side of coupled faces
650  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
651 
652  List<int> globalNeighbour(mesh_.nFaces()-mesh_.nInternalFaces());
653 
654  forAll(patches, patchI)
655  {
656  const polyPatch& pp = patches[patchI];
657 
658  if (pp.coupled())
659  {
660  label faceI = pp.start();
661  label bFaceI = pp.start() - mesh_.nInternalFaces();
662 
663  forAll(pp, i)
664  {
665  label ownRegion = cellToRegion[faceOwner[faceI]];
666  globalNeighbour[bFaceI++] = globalRegions.toGlobal(ownRegion);
667  faceI++;
668  }
669  }
670  }
671 
672  // Get the cell on the other side of coupled patches
673  syncTools::swapBoundaryFaceList(mesh_, globalNeighbour, false);
674 
675 
676  // Get globalCellCells on coarse mesh
677  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
678 
679  labelListList globalRegionRegions;
680  {
681  List<DynamicList<label> > dynRegionRegions(regionPoints.size());
682 
683  // Internal faces first
684  forAll(faceNeighbour, faceI)
685  {
686  label ownRegion = cellToRegion[faceOwner[faceI]];
687  label neiRegion = cellToRegion[faceNeighbour[faceI]];
688 
689  if (ownRegion != neiRegion)
690  {
691  label globalOwn = globalRegions.toGlobal(ownRegion);
692  label globalNei = globalRegions.toGlobal(neiRegion);
693 
694  if (findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
695  {
696  dynRegionRegions[ownRegion].append(globalNei);
697  }
698  if (findIndex(dynRegionRegions[neiRegion], globalOwn) == -1)
699  {
700  dynRegionRegions[neiRegion].append(globalOwn);
701  }
702  }
703  }
704 
705  // Coupled boundary faces
706  forAll(patches, patchI)
707  {
708  const polyPatch& pp = patches[patchI];
709 
710  if (pp.coupled())
711  {
712  label faceI = pp.start();
713  label bFaceI = pp.start() - mesh_.nInternalFaces();
714 
715  forAll(pp, i)
716  {
717  label ownRegion = cellToRegion[faceOwner[faceI]];
718  label globalNei = globalNeighbour[bFaceI++];
719  faceI++;
720 
721  if (findIndex(dynRegionRegions[ownRegion], globalNei) == -1)
722  {
723  dynRegionRegions[ownRegion].append(globalNei);
724  }
725  }
726  }
727  }
728 
729  globalRegionRegions.setSize(dynRegionRegions.size());
730  forAll(dynRegionRegions, i)
731  {
732  globalRegionRegions[i].transfer(dynRegionRegions[i]);
733  }
734  }
735 
736  labelList regionDecomp
737  (
738  decompose
739  (
740  globalRegionRegions,
741  regionPoints,
742  regionWeights
743  )
744  );
745 
746  // Rework back into decomposition for original mesh_
747  labelList cellDistribution(cellToRegion.size());
748 
749  forAll(cellDistribution, cellI)
750  {
751  cellDistribution[cellI] = regionDecomp[cellToRegion[cellI]];
752  }
753  return cellDistribution;
754 }
755 
756 
757 Foam::labelList Foam::parMetisDecomp::decompose
758 (
759  const labelListList& globalCellCells,
760  const pointField& cellCentres,
761  const scalarField& cWeights
762 )
763 {
764  if (cellCentres.size() != globalCellCells.size())
765  {
767  (
768  "parMetisDecomp::decompose(const labelListList&"
769  ", const pointField&, const scalarField&)"
770  ) << "Inconsistent number of cells (" << globalCellCells.size()
771  << ") and number of cell centres (" << cellCentres.size()
772  << ") or weights (" << cWeights.size() << ")." << exit(FatalError);
773  }
774 
775  // For running sequential ...
776  if (Pstream::nProcs() <= 1)
777  {
778  return metisDecomp(decompositionDict_, mesh_)
779  .decompose(globalCellCells, cellCentres, cWeights);
780  }
781 
782 
783  // Make Metis Distributed CSR (Compressed Storage Format) storage
784 
785  // Connections
786  Field<int> adjncy;
787  // Offsets into adjncy
788  Field<int> xadj;
789  scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
790 
791  // decomposition options. 0 = use defaults
792  List<int> options(3, 0);
793  //options[0] = 1; // don't use defaults but use values below
794  //options[1] = -1; // full debug info
795  //options[2] = 15; // random number seed
796 
797  // cell weights (so on the vertices of the dual)
798  Field<int> cellWeights;
799 
800  // face weights (so on the edges of the dual)
801  Field<int> faceWeights;
802 
803 
804  // Check for externally provided cellweights and if so initialise weights
805  scalar minWeights = gMin(cWeights);
806  if (cWeights.size() > 0)
807  {
808  if (minWeights <= 0)
809  {
810  WarningIn
811  (
812  "parMetisDecomp::decompose(const labelListList&"
813  ", const pointField&, const scalarField&)"
814  ) << "Illegal minimum weight " << minWeights
815  << endl;
816  }
817 
818  if (cWeights.size() != globalCellCells.size())
819  {
821  (
822  "parMetisDecomp::decompose(const labelListList&"
823  ", const pointField&, const scalarField&)"
824  ) << "Number of cell weights " << cWeights.size()
825  << " does not equal number of cells " << globalCellCells.size()
826  << exit(FatalError);
827  }
828 
829  // Convert to integers.
830  cellWeights.setSize(cWeights.size());
831  forAll(cellWeights, i)
832  {
833  cellWeights[i] = int(cWeights[i]/minWeights);
834  }
835  }
836 
837 
838  // Check for user supplied weights and decomp options
839  if (decompositionDict_.found("metisCoeffs"))
840  {
841  const dictionary& metisCoeffs =
842  decompositionDict_.subDict("metisCoeffs");
843 
844  if (metisCoeffs.readIfPresent("options", options))
845  {
846  Info<< "Using Metis options " << options
847  << nl << endl;
848 
849  if (options.size() != 3)
850  {
852  (
853  "parMetisDecomp::decompose(const labelListList&"
854  ", const pointField&, const scalarField&)"
855  ) << "Number of options " << options.size()
856  << " should be three." << exit(FatalError);
857  }
858  }
859  }
860 
861 
862  // Do actual decomposition
863  List<int> finalDecomp;
864  decompose
865  (
866  xadj,
867  adjncy,
868  cellCentres,
869  cellWeights,
870  faceWeights,
871  options,
872 
873  finalDecomp
874  );
875 
876  // Copy back to labelList
877  labelList decomp(finalDecomp.size());
878  forAll(decomp, i)
879  {
880  decomp[i] = finalDecomp[i];
881  }
882  return decomp;
883 }
884 
885 
886 void Foam::parMetisDecomp::calcMetisDistributedCSR
887 (
888  const polyMesh& mesh,
889  List<int>& adjncy,
890  List<int>& xadj
891 )
892 {
893  // Create global cell numbers
894  // ~~~~~~~~~~~~~~~~~~~~~~~~~~
895 
896  globalIndex globalCells(mesh.nCells());
897 
898 
899  //
900  // Make Metis Distributed CSR (Compressed Storage Format) storage
901  // adjncy : contains cellCells (= edges in graph)
902  // xadj(celli) : start of information in adjncy for celli
903  //
904 
905 
906  const labelList& faceOwner = mesh.faceOwner();
907  const labelList& faceNeighbour = mesh.faceNeighbour();
908  const polyBoundaryMesh& patches = mesh.boundaryMesh();
909 
910 
911  // Get renumbered owner on other side of coupled faces
912  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
913 
914  List<int> globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
915 
916  forAll(patches, patchI)
917  {
918  const polyPatch& pp = patches[patchI];
919 
920  if (pp.coupled())
921  {
922  label faceI = pp.start();
923  label bFaceI = pp.start() - mesh.nInternalFaces();
924 
925  forAll(pp, i)
926  {
927  globalNeighbour[bFaceI++] = globalCells.toGlobal
928  (
929  faceOwner[faceI++]
930  );
931  }
932  }
933  }
934 
935  // Get the cell on the other side of coupled patches
936  syncTools::swapBoundaryFaceList(mesh, globalNeighbour, false);
937 
938 
939  // Count number of faces (internal + coupled)
940  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
941 
942  // Number of faces per cell
943  List<int> nFacesPerCell(mesh.nCells(), 0);
944 
945  // Number of coupled faces
946  label nCoupledFaces = 0;
947 
948  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
949  {
950  nFacesPerCell[faceOwner[faceI]]++;
951  nFacesPerCell[faceNeighbour[faceI]]++;
952  }
953  // Handle coupled faces
954  HashSet<edge, Hash<edge> > cellPair(mesh.nFaces()-mesh.nInternalFaces());
955 
956  forAll(patches, patchI)
957  {
958  const polyPatch& pp = patches[patchI];
959 
960  if (pp.coupled())
961  {
962  label faceI = pp.start();
963 
964  forAll(pp, i)
965  {
966  label own = faceOwner[faceI];
967  label globalNei = globalNeighbour[faceI-mesh.nInternalFaces()];
968 
969  if (cellPair.insert(edge(own, globalNei)))
970  {
971  nCoupledFaces++;
972  nFacesPerCell[faceOwner[faceI]]++;
973  }
974  faceI++;
975  }
976  }
977  }
978 
979 
980  // Fill in xadj
981  // ~~~~~~~~~~~~
982 
983  xadj.setSize(mesh.nCells()+1);
984 
985  int freeAdj = 0;
986 
987  for (label cellI = 0; cellI < mesh.nCells(); cellI++)
988  {
989  xadj[cellI] = freeAdj;
990 
991  freeAdj += nFacesPerCell[cellI];
992  }
993  xadj[mesh.nCells()] = freeAdj;
994 
995 
996 
997  // Fill in adjncy
998  // ~~~~~~~~~~~~~~
999 
1000  adjncy.setSize(2*mesh.nInternalFaces() + nCoupledFaces);
1001 
1002  nFacesPerCell = 0;
1003 
1004  // For internal faces is just offsetted owner and neighbour
1005  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
1006  {
1007  label own = faceOwner[faceI];
1008  label nei = faceNeighbour[faceI];
1009 
1010  adjncy[xadj[own] + nFacesPerCell[own]++] = globalCells.toGlobal(nei);
1011  adjncy[xadj[nei] + nFacesPerCell[nei]++] = globalCells.toGlobal(own);
1012  }
1013  // For boundary faces is offsetted coupled neighbour
1014  cellPair.clear();
1015  forAll(patches, patchI)
1016  {
1017  const polyPatch& pp = patches[patchI];
1018 
1019  if (pp.coupled())
1020  {
1021  label faceI = pp.start();
1022  label bFaceI = pp.start()-mesh.nInternalFaces();
1023 
1024  forAll(pp, i)
1025  {
1026  label own = faceOwner[faceI];
1027  label globalNei = globalNeighbour[bFaceI];
1028 
1029  if (cellPair.insert(edge(own, globalNei)))
1030  {
1031  adjncy[xadj[own] + nFacesPerCell[own]++] =
1032  globalNeighbour[bFaceI];
1033  }
1034 
1035  faceI++;
1036  bFaceI++;
1037  }
1038  }
1039  }
1040 }
1041 
1042 
1043 // ************************ vim: set sw=4 sts=4 et: ************************ //