FreeFOAM The Cross-Platform CFD Toolkit
domainDecomposition.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 "domainDecomposition.H"
27 #include <OpenFOAM/dictionary.H>
28 #include <OpenFOAM/labelIOList.H>
30 #include <finiteVolume/fvMesh.H>
31 #include <OpenFOAM/OSspecific.H>
32 #include <OpenFOAM/Map.H>
34 #include <OpenFOAM/DynamicList.H>
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 void domainDecomposition::mark
39 (
40  const labelList& zoneElems,
41  const label zoneI,
42  labelList& elementToZone
43 )
44 {
45  forAll(zoneElems, i)
46  {
47  label pointi = zoneElems[i];
48 
49  if (elementToZone[pointi] == -1)
50  {
51  // First occurrence
52  elementToZone[pointi] = zoneI;
53  }
54  else if (elementToZone[pointi] >= 0)
55  {
56  // Multiple zones
57  elementToZone[pointi] = -2;
58  }
59  }
60 }
61 
62 
63 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
64 
65 // from components
67 :
68  fvMesh(io),
69  decompositionDict_
70  (
71  IOobject
72  (
73  "decomposeParDict",
74  time().system(),
75  *this,
76  IOobject::MUST_READ,
77  IOobject::NO_WRITE
78  )
79  ),
80  nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
81  distributed_(false),
82  cellToProc_(nCells()),
83  procPointAddressing_(nProcs_),
84  procFaceAddressing_(nProcs_),
85  procCellAddressing_(nProcs_),
86  procBoundaryAddressing_(nProcs_),
87  procPatchSize_(nProcs_),
88  procPatchStartIndex_(nProcs_),
89  procNeighbourProcessors_(nProcs_),
90  procProcessorPatchSize_(nProcs_),
91  procProcessorPatchStartIndex_(nProcs_),
92  globallySharedPoints_(0),
93  cyclicParallel_(false)
94 {
95  if (decompositionDict_.found("distributed"))
96  {
97  Switch distributed(decompositionDict_.lookup("distributed"));
98  distributed_ = distributed;
99  }
100 }
101 
102 
103 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
104 
106 {}
107 
108 
109 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
110 
112 {
113  Info<< "\nConstructing processor meshes" << endl;
114 
115  // Make a lookup map for globally shared points
116  Map<label> sharedPointLookup(2*globallySharedPoints_.size());
117 
118  forAll (globallySharedPoints_, pointi)
119  {
120  sharedPointLookup.insert(globallySharedPoints_[pointi], pointi);
121  }
122 
123 
124  // Mark point/faces/cells that are in zones.
125  // -1 : not in zone
126  // -2 : in multiple zones
127  // >= 0 : in single given zone
128  // This will give direct lookup of elements that are in a single zone
129  // and we'll only have to revert back to searching through all zones
130  // for the duplicate elements
131 
132  // Point zones
133  labelList pointToZone(points().size(), -1);
134 
135  forAll(pointZones(), zoneI)
136  {
137  mark(pointZones()[zoneI], zoneI, pointToZone);
138  }
139 
140  // Face zones
141  labelList faceToZone(faces().size(), -1);
142 
143  forAll(faceZones(), zoneI)
144  {
145  mark(faceZones()[zoneI], zoneI, faceToZone);
146  }
147 
148  // Cell zones
149  labelList cellToZone(nCells(), -1);
150 
151  forAll(cellZones(), zoneI)
152  {
153  mark(cellZones()[zoneI], zoneI, cellToZone);
154  }
155 
156 
157  label totProcFaces = 0;
158  label maxProcPatches = 0;
159  label maxProcFaces = 0;
160 
161 
162  // Write out the meshes
163  for (label procI = 0; procI < nProcs_; procI++)
164  {
165  // Create processor points
166  const labelList& curPointLabels = procPointAddressing_[procI];
167 
168  const pointField& meshPoints = points();
169 
170  labelList pointLookup(nPoints(), -1);
171 
172  pointField procPoints(curPointLabels.size());
173 
174  forAll (curPointLabels, pointi)
175  {
176  procPoints[pointi] = meshPoints[curPointLabels[pointi]];
177 
178  pointLookup[curPointLabels[pointi]] = pointi;
179  }
180 
181  // Create processor faces
182  const labelList& curFaceLabels = procFaceAddressing_[procI];
183 
184  const faceList& meshFaces = faces();
185 
186  labelList faceLookup(nFaces(), -1);
187 
188  faceList procFaces(curFaceLabels.size());
189 
190  forAll (curFaceLabels, facei)
191  {
192  // Mark the original face as used
193  // Remember to decrement the index by one (turning index)
194  //
195  label curF = mag(curFaceLabels[facei]) - 1;
196 
197  faceLookup[curF] = facei;
198 
199  // get the original face
200  labelList origFaceLabels;
201 
202  if (curFaceLabels[facei] >= 0)
203  {
204  // face not turned
205  origFaceLabels = meshFaces[curF];
206  }
207  else
208  {
209  origFaceLabels = meshFaces[curF].reverseFace();
210  }
211 
212  // translate face labels into local point list
213  face& procFaceLabels = procFaces[facei];
214 
215  procFaceLabels.setSize(origFaceLabels.size());
216 
217  forAll (origFaceLabels, pointi)
218  {
219  procFaceLabels[pointi] = pointLookup[origFaceLabels[pointi]];
220  }
221  }
222 
223  // Create processor cells
224  const labelList& curCellLabels = procCellAddressing_[procI];
225 
226  const cellList& meshCells = cells();
227 
228  cellList procCells(curCellLabels.size());
229 
230  forAll (curCellLabels, celli)
231  {
232  const labelList& origCellLabels = meshCells[curCellLabels[celli]];
233 
234  cell& curCell = procCells[celli];
235 
236  curCell.setSize(origCellLabels.size());
237 
238  forAll (origCellLabels, cellFaceI)
239  {
240  curCell[cellFaceI] = faceLookup[origCellLabels[cellFaceI]];
241  }
242  }
243 
244  // Create processor mesh without a boundary
245 
246  fileName processorCasePath
247  (
248  time().caseName()/fileName(word("processor") + Foam::name(procI))
249  );
250 
251  // make the processor directory
252  mkDir(time().rootPath()/processorCasePath);
253 
254  // create a database
255  Time processorDb
256  (
258  time().rootPath(),
259  processorCasePath,
260  "system",
261  "constant"
262  );
263 
264  // create the mesh
265  polyMesh procMesh
266  (
267  IOobject
268  (
269  this->polyMesh::name(), // region name of undecomposed mesh
270  pointsInstance(),
271  processorDb
272  ),
273  xferMove(procPoints),
274  xferMove(procFaces),
275  xferMove(procCells)
276  );
277 
278  // Create processor boundary patches
279  const labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
280 
281  const labelList& curPatchSizes = procPatchSize_[procI];
282 
283  const labelList& curPatchStarts = procPatchStartIndex_[procI];
284 
285  const labelList& curNeighbourProcessors =
286  procNeighbourProcessors_[procI];
287 
288  const labelList& curProcessorPatchSizes =
289  procProcessorPatchSize_[procI];
290 
291  const labelList& curProcessorPatchStarts =
292  procProcessorPatchStartIndex_[procI];
293 
294  const polyPatchList& meshPatches = boundaryMesh();
295 
296  List<polyPatch*> procPatches
297  (
298  curPatchSizes.size()
299  + curProcessorPatchSizes.size(),
300  reinterpret_cast<polyPatch*>(0)
301  );
302 
303  label nPatches = 0;
304 
305  forAll (curPatchSizes, patchi)
306  {
307  procPatches[nPatches] =
308  meshPatches[curBoundaryAddressing[patchi]].clone
309  (
310  procMesh.boundaryMesh(),
311  nPatches,
312  curPatchSizes[patchi],
313  curPatchStarts[patchi]
314  ).ptr();
315 
316  nPatches++;
317  }
318 
319  forAll (curProcessorPatchSizes, procPatchI)
320  {
321  procPatches[nPatches] =
323  (
324  word("procBoundary") + Foam::name(procI)
325  + word("to")
326  + Foam::name(curNeighbourProcessors[procPatchI]),
327  curProcessorPatchSizes[procPatchI],
328  curProcessorPatchStarts[procPatchI],
329  nPatches,
330  procMesh.boundaryMesh(),
331  procI,
332  curNeighbourProcessors[procPatchI]
333  );
334 
335  nPatches++;
336  }
337 
338  // Add boundary patches
339  procMesh.addPatches(procPatches);
340 
341  // Create and add zones
342 
343  // Point zones
344  {
345  const pointZoneMesh& pz = pointZones();
346 
347  // Go through all the zoned points and find out if they
348  // belong to a zone. If so, add it to the zone as
349  // necessary
350  List<DynamicList<label> > zonePoints(pz.size());
351 
352  // Estimate size
353  forAll(zonePoints, zoneI)
354  {
355  zonePoints[zoneI].setCapacity(pz[zoneI].size() / nProcs_);
356  }
357 
358  // Use the pointToZone map to find out the single zone (if any),
359  // use slow search only for shared points.
360  forAll (curPointLabels, pointi)
361  {
362  label curPoint = curPointLabels[pointi];
363 
364  label zoneI = pointToZone[curPoint];
365 
366  if (zoneI >= 0)
367  {
368  // Single zone.
369  zonePoints[zoneI].append(pointi);
370  }
371  else if (zoneI == -2)
372  {
373  // Multiple zones. Lookup.
374  forAll(pz, zoneI)
375  {
376  label index = pz[zoneI].whichPoint(curPoint);
377 
378  if (index != -1)
379  {
380  zonePoints[zoneI].append(pointi);
381  }
382  }
383  }
384  }
385 
386  procMesh.pointZones().clearAddressing();
387  procMesh.pointZones().setSize(zonePoints.size());
388  forAll(zonePoints, zoneI)
389  {
390  procMesh.pointZones().set
391  (
392  zoneI,
393  pz[zoneI].clone
394  (
395  procMesh.pointZones(),
396  zoneI,
397  zonePoints[zoneI].shrink()
398  )
399  );
400  }
401 
402  if (pz.size())
403  {
404  // Force writing on all processors
405  procMesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
406  }
407  }
408 
409  // Face zones
410  {
411  const faceZoneMesh& fz = faceZones();
412 
413  // Go through all the zoned face and find out if they
414  // belong to a zone. If so, add it to the zone as
415  // necessary
416  List<DynamicList<label> > zoneFaces(fz.size());
417  List<DynamicList<bool> > zoneFaceFlips(fz.size());
418 
419  // Estimate size
420  forAll(zoneFaces, zoneI)
421  {
422  label procSize = fz[zoneI].size() / nProcs_;
423 
424  zoneFaces[zoneI].setCapacity(procSize);
425  zoneFaceFlips[zoneI].setCapacity(procSize);
426  }
427 
428  // Go through all the zoned faces and find out if they
429  // belong to a zone. If so, add it to the zone as
430  // necessary
431  forAll (curFaceLabels, facei)
432  {
433  // Remember to decrement the index by one (turning index)
434  //
435  label curF = mag(curFaceLabels[facei]) - 1;
436 
437  label zoneI = faceToZone[curF];
438 
439  if (zoneI >= 0)
440  {
441  // Single zone. Add the face
442  zoneFaces[zoneI].append(facei);
443 
444  label index = fz[zoneI].whichFace(curF);
445 
446  bool flip = fz[zoneI].flipMap()[index];
447 
448  if (curFaceLabels[facei] < 0)
449  {
450  flip = !flip;
451  }
452 
453  zoneFaceFlips[zoneI].append(flip);
454  }
455  else if (zoneI == -2)
456  {
457  // Multiple zones. Lookup.
458  forAll(fz, zoneI)
459  {
460  label index = fz[zoneI].whichFace(curF);
461 
462  if (index != -1)
463  {
464  zoneFaces[zoneI].append(facei);
465 
466  bool flip = fz[zoneI].flipMap()[index];
467 
468  if (curFaceLabels[facei] < 0)
469  {
470  flip = !flip;
471  }
472 
473  zoneFaceFlips[zoneI].append(flip);
474  }
475  }
476  }
477  }
478 
479  procMesh.faceZones().clearAddressing();
480  procMesh.faceZones().setSize(zoneFaces.size());
481  forAll(zoneFaces, zoneI)
482  {
483  procMesh.faceZones().set
484  (
485  zoneI,
486  fz[zoneI].clone
487  (
488  zoneFaces[zoneI].shrink(), // addressing
489  zoneFaceFlips[zoneI].shrink(), // flipmap
490  zoneI,
491  procMesh.faceZones()
492  )
493  );
494  }
495 
496  if (fz.size())
497  {
498  // Force writing on all processors
499  procMesh.faceZones().writeOpt() = IOobject::AUTO_WRITE;
500  }
501  }
502 
503  // Cell zones
504  {
505  const cellZoneMesh& cz = cellZones();
506 
507  // Go through all the zoned cells and find out if they
508  // belong to a zone. If so, add it to the zone as
509  // necessary
510  List<DynamicList<label> > zoneCells(cz.size());
511 
512  // Estimate size
513  forAll(zoneCells, zoneI)
514  {
515  zoneCells[zoneI].setCapacity(cz[zoneI].size() / nProcs_);
516  }
517 
518  forAll (curCellLabels, celli)
519  {
520  label curCellI = curCellLabels[celli];
521 
522  label zoneI = cellToZone[curCellI];
523 
524  if (zoneI >= 0)
525  {
526  // Single zone.
527  zoneCells[zoneI].append(celli);
528  }
529  else if (zoneI == -2)
530  {
531  // Multiple zones. Lookup.
532  forAll(cz, zoneI)
533  {
534  label index = cz[zoneI].whichCell(curCellI);
535 
536  if (index != -1)
537  {
538  zoneCells[zoneI].append(celli);
539  }
540  }
541  }
542  }
543 
544  procMesh.cellZones().clearAddressing();
545  procMesh.cellZones().setSize(zoneCells.size());
546  forAll(zoneCells, zoneI)
547  {
548  procMesh.cellZones().set
549  (
550  zoneI,
551  cz[zoneI].clone
552  (
553  zoneCells[zoneI].shrink(),
554  zoneI,
555  procMesh.cellZones()
556  )
557  );
558  }
559 
560  if (cz.size())
561  {
562  // Force writing on all processors
563  procMesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
564  }
565  }
566 
567  // Set the precision of the points data to 10
569 
570  procMesh.write();
571 
572  Info<< endl
573  << "Processor " << procI << nl
574  << " Number of cells = " << procMesh.nCells()
575  << endl;
576 
577  label nBoundaryFaces = 0;
578  label nProcPatches = 0;
579  label nProcFaces = 0;
580 
581  forAll (procMesh.boundaryMesh(), patchi)
582  {
583  if
584  (
585  procMesh.boundaryMesh()[patchi].type()
586  == processorPolyPatch::typeName
587  )
588  {
589  const processorPolyPatch& ppp =
590  refCast<const processorPolyPatch>
591  (
592  procMesh.boundaryMesh()[patchi]
593  );
594 
595  Info<< " Number of faces shared with processor "
596  << ppp.neighbProcNo() << " = " << ppp.size() << endl;
597 
598  nProcPatches++;
599  nProcFaces += ppp.size();
600  }
601  else
602  {
603  nBoundaryFaces += procMesh.boundaryMesh()[patchi].size();
604  }
605  }
606 
607  Info<< " Number of processor patches = " << nProcPatches << nl
608  << " Number of processor faces = " << nProcFaces << nl
609  << " Number of boundary faces = " << nBoundaryFaces << endl;
610 
611  totProcFaces += nProcFaces;
612  maxProcPatches = max(maxProcPatches, nProcPatches);
613  maxProcFaces = max(maxProcFaces, nProcFaces);
614 
615  // create and write the addressing information
616  labelIOList pointProcAddressing
617  (
618  IOobject
619  (
620  "pointProcAddressing",
621  procMesh.facesInstance(),
622  procMesh.meshSubDir,
623  procMesh,
626  ),
627  procPointAddressing_[procI]
628  );
629  pointProcAddressing.write();
630 
631  labelIOList faceProcAddressing
632  (
633  IOobject
634  (
635  "faceProcAddressing",
636  procMesh.facesInstance(),
637  procMesh.meshSubDir,
638  procMesh,
641  ),
642  procFaceAddressing_[procI]
643  );
644  faceProcAddressing.write();
645 
646  labelIOList cellProcAddressing
647  (
648  IOobject
649  (
650  "cellProcAddressing",
651  procMesh.facesInstance(),
652  procMesh.meshSubDir,
653  procMesh,
656  ),
657  procCellAddressing_[procI]
658  );
659  cellProcAddressing.write();
660 
661  labelIOList boundaryProcAddressing
662  (
663  IOobject
664  (
665  "boundaryProcAddressing",
666  procMesh.facesInstance(),
667  procMesh.meshSubDir,
668  procMesh,
671  ),
672  procBoundaryAddressing_[procI]
673  );
674  boundaryProcAddressing.write();
675  }
676 
677  Info<< nl
678  << "Number of processor faces = " << totProcFaces/2 << nl
679  << "Max number of processor patches = " << maxProcPatches << nl
680  << "Max number of faces between processors = " << maxProcFaces
681  << endl;
682 
683  return true;
684 }
685 
686 
687 // ************************ vim: set sw=4 sts=4 et: ************************ //