FreeFOAM The Cross-Platform CFD Toolkit
decomposeMesh.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 InClass
25  domainDecomposition
26 
27 Description
28  Private member of domainDecomposition.
29  Decomposes the mesh into bits
30 
31 \*---------------------------------------------------------------------------*/
32 
33 #include "domainDecomposition.H"
34 #include <OpenFOAM/IOstreams.H>
35 #include <OpenFOAM/SLPtrList.H>
36 #include <OpenFOAM/boolList.H>
37 #include <OpenFOAM/primitiveMesh.H>
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void domainDecomposition::decomposeMesh(const bool filterEmptyPatches)
43 {
44  // Decide which cell goes to which processor
45  distributeCells();
46 
47  // Distribute the cells according to the given processor label
48 
49  // calculate the addressing information for the original mesh
50  Info<< "\nCalculating original mesh data" << endl;
51 
52  // set references to the original mesh
54  const faceList& fcs = faces();
55  const labelList& owner = faceOwner();
56  const labelList& neighbour = faceNeighbour();
57 
58  // loop through the list of processor labels for the cell and add the
59  // cell shape to the list of cells for the appropriate processor
60 
61  Info<< "\nDistributing cells to processors" << endl;
62 
63  // Memory management
64  {
65  List<SLList<label> > procCellList(nProcs_);
66 
67  forAll (cellToProc_, celli)
68  {
69  if (cellToProc_[celli] >= nProcs_)
70  {
71  FatalErrorIn("domainDecomposition::decomposeMesh()")
72  << "Impossible processor label " << cellToProc_[celli]
73  << "for cell " << celli
74  << abort(FatalError);
75  }
76  else
77  {
78  procCellList[cellToProc_[celli]].append(celli);
79  }
80  }
81 
82  // Convert linked lists into normal lists
83  forAll (procCellList, procI)
84  {
85  procCellAddressing_[procI] = procCellList[procI];
86  }
87  }
88 
89  Info << "\nDistributing faces to processors" << endl;
90 
91  // Loop through all internal faces and decide which processor they belong to
92  // First visit all internal faces. If cells at both sides belong to the
93  // same processor, the face is an internal face. If they are different,
94  // it belongs to both processors.
95 
96  // Memory management
97  {
98  List<SLList<label> > procFaceList(nProcs_);
99 
100  forAll (neighbour, facei)
101  {
102  if (cellToProc_[owner[facei]] == cellToProc_[neighbour[facei]])
103  {
104  // Face internal to processor
105  procFaceList[cellToProc_[owner[facei]]].append(facei);
106  }
107  }
108 
109  // Detect inter-processor boundaries
110 
111  // Neighbour processor for each subdomain
112  List<SLList<label> > interProcBoundaries(nProcs_);
113 
114  // Face labels belonging to each inter-processor boundary
115  List<SLList<SLList<label> > > interProcBFaces(nProcs_);
116 
117  List<SLList<label> > procPatchIndex(nProcs_);
118 
119  forAll (neighbour, facei)
120  {
121  if (cellToProc_[owner[facei]] != cellToProc_[neighbour[facei]])
122  {
123  // inter - processor patch face found. Go through the list of
124  // inside boundaries for the owner processor and try to find
125  // this inter-processor patch.
126 
127  label ownerProc = cellToProc_[owner[facei]];
128  label neighbourProc = cellToProc_[neighbour[facei]];
129 
130  SLList<label>::iterator curInterProcBdrsOwnIter =
131  interProcBoundaries[ownerProc].begin();
132 
133  SLList<SLList<label> >::iterator curInterProcBFacesOwnIter =
134  interProcBFaces[ownerProc].begin();
135 
136  bool interProcBouFound = false;
137 
138  // WARNING: Synchronous SLList iterators
139 
140  for
141  (
142  ;
143  curInterProcBdrsOwnIter
144  != interProcBoundaries[ownerProc].end()
145  && curInterProcBFacesOwnIter
146  != interProcBFaces[ownerProc].end();
147  ++curInterProcBdrsOwnIter, ++curInterProcBFacesOwnIter
148  )
149  {
150  if (curInterProcBdrsOwnIter() == neighbourProc)
151  {
152  // the inter - processor boundary exists. Add the face
153  interProcBouFound = true;
154 
155  curInterProcBFacesOwnIter().append(facei);
156 
157  SLList<label>::iterator curInterProcBdrsNeiIter =
158  interProcBoundaries[neighbourProc].begin();
159 
160  SLList<SLList<label> >::iterator
161  curInterProcBFacesNeiIter =
162  interProcBFaces[neighbourProc].begin();
163 
164  bool neighbourFound = false;
165 
166  // WARNING: Synchronous SLList iterators
167 
168  for
169  (
170  ;
171  curInterProcBdrsNeiIter !=
172  interProcBoundaries[neighbourProc].end()
173  && curInterProcBFacesNeiIter !=
174  interProcBFaces[neighbourProc].end();
175  ++curInterProcBdrsNeiIter,
176  ++curInterProcBFacesNeiIter
177  )
178  {
179  if (curInterProcBdrsNeiIter() == ownerProc)
180  {
181  // boundary found. Add the face
182  neighbourFound = true;
183 
184  curInterProcBFacesNeiIter().append(facei);
185  }
186 
187  if (neighbourFound) break;
188  }
189 
190  if (interProcBouFound && !neighbourFound)
191  {
192  FatalErrorIn("domainDecomposition::decomposeMesh()")
193  << "Inconsistency in inter - "
194  << "processor boundary lists for processors "
195  << ownerProc << " and " << neighbourProc
196  << abort(FatalError);
197  }
198  }
199 
200  if (interProcBouFound) break;
201  }
202 
203  if (!interProcBouFound)
204  {
205  // inter - processor boundaries do not exist and need to
206  // be created
207 
208  // set the new addressing information
209 
210  // owner
211  interProcBoundaries[ownerProc].append(neighbourProc);
212  interProcBFaces[ownerProc].append(SLList<label>(facei));
213 
214  // neighbour
215  interProcBoundaries[neighbourProc].append(ownerProc);
216  interProcBFaces[neighbourProc].append(SLList<label>(facei));
217  }
218  }
219  }
220 
221  // Loop through patches. For cyclic boundaries detect inter-processor
222  // faces; for all other, add faces to the face list and remember start
223  // and size of all patches.
224 
225  // for all processors, set the size of start index and patch size
226  // lists to the number of patches in the mesh
227  forAll (procPatchSize_, procI)
228  {
229  procPatchSize_[procI].setSize(patches.size());
230  procPatchStartIndex_[procI].setSize(patches.size());
231  }
232 
233  forAll (patches, patchi)
234  {
235  // Reset size and start index for all processors
236  forAll (procPatchSize_, procI)
237  {
238  procPatchSize_[procI][patchi] = 0;
239  procPatchStartIndex_[procI][patchi] =
240  procFaceList[procI].size();
241  }
242 
243  const label patchStart = patches[patchi].start();
244 
245  if (!isA<cyclicPolyPatch>(patches[patchi]))
246  {
247  // Normal patch. Add faces to processor where the cell
248  // next to the face lives
249 
250  const unallocLabelList& patchFaceCells =
251  patches[patchi].faceCells();
252 
253  forAll (patchFaceCells, facei)
254  {
255  const label curProc = cellToProc_[patchFaceCells[facei]];
256 
257  // add the face
258  procFaceList[curProc].append(patchStart + facei);
259 
260  // increment the number of faces for this patch
261  procPatchSize_[curProc][patchi]++;
262  }
263  }
264  else
265  {
266  // Cyclic patch special treatment
267 
268  const polyPatch& cPatch = patches[patchi];
269 
270  const label cycOffset = cPatch.size()/2;
271 
272  // Set reference to faceCells for both patches
273  const labelList::subList firstFaceCells
274  (
275  cPatch.faceCells(),
276  cycOffset
277  );
278 
279  const labelList::subList secondFaceCells
280  (
281  cPatch.faceCells(),
282  cycOffset,
283  cycOffset
284  );
285 
286  forAll (firstFaceCells, facei)
287  {
288  if
289  (
290  cellToProc_[firstFaceCells[facei]]
291  != cellToProc_[secondFaceCells[facei]]
292  )
293  {
294  // This face becomes an inter-processor boundary face
295  // inter - processor patch face found. Go through
296  // the list of inside boundaries for the owner
297  // processor and try to find this inter-processor
298  // patch.
299 
300  cyclicParallel_ = true;
301 
302  label ownerProc = cellToProc_[firstFaceCells[facei]];
303  label neighbourProc =
304  cellToProc_[secondFaceCells[facei]];
305 
306  SLList<label>::iterator curInterProcBdrsOwnIter =
307  interProcBoundaries[ownerProc].begin();
308 
309  SLList<SLList<label> >::iterator
310  curInterProcBFacesOwnIter =
311  interProcBFaces[ownerProc].begin();
312 
313  bool interProcBouFound = false;
314 
315  // WARNING: Synchronous SLList iterators
316 
317  for
318  (
319  ;
320  curInterProcBdrsOwnIter !=
321  interProcBoundaries[ownerProc].end()
322  && curInterProcBFacesOwnIter !=
323  interProcBFaces[ownerProc].end();
324  ++curInterProcBdrsOwnIter,
325  ++curInterProcBFacesOwnIter
326  )
327  {
328  if (curInterProcBdrsOwnIter() == neighbourProc)
329  {
330  // the inter - processor boundary exists.
331  // Add the face
332  interProcBouFound = true;
333 
334  curInterProcBFacesOwnIter().append
335  (patchStart + facei);
336 
337  SLList<label>::iterator curInterProcBdrsNeiIter
338  = interProcBoundaries[neighbourProc].begin();
339 
340  SLList<SLList<label> >::iterator
341  curInterProcBFacesNeiIter =
342  interProcBFaces[neighbourProc].begin();
343 
344  bool neighbourFound = false;
345 
346  // WARNING: Synchronous SLList iterators
347 
348  for
349  (
350  ;
351  curInterProcBdrsNeiIter
352  != interProcBoundaries[neighbourProc].end()
353  && curInterProcBFacesNeiIter
354  != interProcBFaces[neighbourProc].end();
355  ++curInterProcBdrsNeiIter,
356  ++curInterProcBFacesNeiIter
357  )
358  {
359  if (curInterProcBdrsNeiIter() == ownerProc)
360  {
361  // boundary found. Add the face
362  neighbourFound = true;
363 
364  curInterProcBFacesNeiIter()
365  .append
366  (
367  patchStart
368  + cycOffset
369  + facei
370  );
371  }
372 
373  if (neighbourFound) break;
374  }
375 
376  if (interProcBouFound && !neighbourFound)
377  {
379  (
380  "domainDecomposition::decomposeMesh()"
381  ) << "Inconsistency in inter-processor "
382  << "boundary lists for processors "
383  << ownerProc << " and " << neighbourProc
384  << " in cyclic boundary matching"
385  << abort(FatalError);
386  }
387  }
388 
389  if (interProcBouFound) break;
390  }
391 
392  if (!interProcBouFound)
393  {
394  // inter - processor boundaries do not exist
395  // and need to be created
396 
397  // set the new addressing information
398 
399  // owner
400  interProcBoundaries[ownerProc]
401  .append(neighbourProc);
402  interProcBFaces[ownerProc]
403  .append(SLList<label>(patchStart + facei));
404 
405  // neighbour
406  interProcBoundaries[neighbourProc]
407  .append(ownerProc);
408  interProcBFaces[neighbourProc]
409  .append
410  (
412  (
413  patchStart
414  + cycOffset
415  + facei
416  )
417  );
418  }
419  }
420  else
421  {
422  // This cyclic face remains on the processor
423  label ownerProc = cellToProc_[firstFaceCells[facei]];
424 
425  // add the face
426  procFaceList[ownerProc].append(patchStart + facei);
427 
428  // increment the number of faces for this patch
429  procPatchSize_[ownerProc][patchi]++;
430 
431  // Note: I cannot add the other side of the cyclic
432  // boundary here because this would violate the order.
433  // They will be added in a separate loop below
434  //
435  }
436  }
437 
438  // Ordering in cyclic boundaries is important.
439  // Add the other half of cyclic faces for cyclic boundaries
440  // that remain on the processor
441  forAll (secondFaceCells, facei)
442  {
443  if
444  (
445  cellToProc_[firstFaceCells[facei]]
446  == cellToProc_[secondFaceCells[facei]]
447  )
448  {
449  // This cyclic face remains on the processor
450  label ownerProc = cellToProc_[firstFaceCells[facei]];
451 
452  // add the second face
453  procFaceList[ownerProc].append
454  (patchStart + cycOffset + facei);
455 
456  // increment the number of faces for this patch
457  procPatchSize_[ownerProc][patchi]++;
458  }
459  }
460  }
461  }
462 
463  // Convert linked lists into normal lists
464  // Add inter-processor boundaries and remember start indices
465  forAll (procFaceList, procI)
466  {
467  // Get internal and regular boundary processor faces
468  SLList<label>& curProcFaces = procFaceList[procI];
469 
470  // Get reference to processor face addressing
471  labelList& curProcFaceAddressing = procFaceAddressing_[procI];
472 
473  labelList& curProcNeighbourProcessors =
474  procNeighbourProcessors_[procI];
475 
476  labelList& curProcProcessorPatchSize =
477  procProcessorPatchSize_[procI];
478 
479  labelList& curProcProcessorPatchStartIndex =
480  procProcessorPatchStartIndex_[procI];
481 
482  // calculate the size
483  label nFacesOnProcessor = curProcFaces.size();
484 
485  for
486  (
487  SLList<SLList<label> >::iterator curInterProcBFacesIter =
488  interProcBFaces[procI].begin();
489  curInterProcBFacesIter != interProcBFaces[procI].end();
490  ++curInterProcBFacesIter
491  )
492  {
493  nFacesOnProcessor += curInterProcBFacesIter().size();
494  }
495 
496  curProcFaceAddressing.setSize(nFacesOnProcessor);
497 
498  // Fill in the list. Calculate turning index.
499  // Turning index will be -1 only for some faces on processor
500  // boundaries, i.e. the ones where the current processor ID
501  // is in the cell which is a face neighbour.
502  // Turning index is stored as the sign of the face addressing list
503 
504  label nFaces = 0;
505 
506  // Add internal and boundary faces
507  // Remember to increment the index by one such that the
508  // turning index works properly.
509  for
510  (
511  SLList<label>::iterator curProcFacesIter = curProcFaces.begin();
512  curProcFacesIter != curProcFaces.end();
513  ++curProcFacesIter
514  )
515  {
516  curProcFaceAddressing[nFaces] = curProcFacesIter() + 1;
517  nFaces++;
518  }
519 
520  // Add inter-processor boundary faces. At the beginning of each
521  // patch, grab the patch start index and size
522 
523  curProcNeighbourProcessors.setSize
524  (
525  interProcBoundaries[procI].size()
526  );
527 
528  curProcProcessorPatchSize.setSize
529  (
530  interProcBoundaries[procI].size()
531  );
532 
533  curProcProcessorPatchStartIndex.setSize
534  (
535  interProcBoundaries[procI].size()
536  );
537 
538  label nProcPatches = 0;
539 
540  SLList<label>::iterator curInterProcBdrsIter =
541  interProcBoundaries[procI].begin();
542 
543  SLList<SLList<label> >::iterator curInterProcBFacesIter =
544  interProcBFaces[procI].begin();
545 
546  for
547  (
548  ;
549  curInterProcBdrsIter != interProcBoundaries[procI].end()
550  && curInterProcBFacesIter != interProcBFaces[procI].end();
551  ++curInterProcBdrsIter, ++curInterProcBFacesIter
552  )
553  {
554  curProcNeighbourProcessors[nProcPatches] =
555  curInterProcBdrsIter();
556 
557  // Get start index for processor patch
558  curProcProcessorPatchStartIndex[nProcPatches] = nFaces;
559 
560  label& curSize =
561  curProcProcessorPatchSize[nProcPatches];
562 
563  curSize = 0;
564 
565  // add faces for this processor boundary
566 
567  for
568  (
569  SLList<label>::iterator curFacesIter =
570  curInterProcBFacesIter().begin();
571  curFacesIter != curInterProcBFacesIter().end();
572  ++curFacesIter
573  )
574  {
575  // add the face
576 
577  // Remember to increment the index by one such that the
578  // turning index works properly.
579  if (cellToProc_[owner[curFacesIter()]] == procI)
580  {
581  curProcFaceAddressing[nFaces] = curFacesIter() + 1;
582  }
583  else
584  {
585  // turning face
586  curProcFaceAddressing[nFaces] = -(curFacesIter() + 1);
587  }
588 
589  // increment the size
590  curSize++;
591 
592  nFaces++;
593  }
594 
595  nProcPatches++;
596  }
597  }
598  }
599 
600  Info << "\nCalculating processor boundary addressing" << endl;
601  // For every patch of processor boundary, find the index of the original
602  // patch. Mis-alignment is caused by the fact that patches with zero size
603  // are omitted. For processor patches, set index to -1.
604  // At the same time, filter the procPatchSize_ and procPatchStartIndex_
605  // lists to exclude zero-size patches
606  forAll (procPatchSize_, procI)
607  {
608  // Make a local copy of old lists
609  const labelList oldPatchSizes = procPatchSize_[procI];
610 
611  const labelList oldPatchStarts = procPatchStartIndex_[procI];
612 
613  labelList& curPatchSizes = procPatchSize_[procI];
614 
615  labelList& curPatchStarts = procPatchStartIndex_[procI];
616 
617  const labelList& curProcessorPatchSizes =
618  procProcessorPatchSize_[procI];
619 
620  labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
621 
622  curBoundaryAddressing.setSize
623  (
624  oldPatchSizes.size()
625  + curProcessorPatchSizes.size()
626  );
627 
628  label nPatches = 0;
629 
630  forAll (oldPatchSizes, patchi)
631  {
632  if (!filterEmptyPatches || oldPatchSizes[patchi] > 0)
633  {
634  curBoundaryAddressing[nPatches] = patchi;
635 
636  curPatchSizes[nPatches] = oldPatchSizes[patchi];
637 
638  curPatchStarts[nPatches] = oldPatchStarts[patchi];
639 
640  nPatches++;
641  }
642  }
643 
644  // reset to the size of live patches
645  curPatchSizes.setSize(nPatches);
646  curPatchStarts.setSize(nPatches);
647 
648  forAll (curProcessorPatchSizes, procPatchI)
649  {
650  curBoundaryAddressing[nPatches] = -1;
651 
652  nPatches++;
653  }
654 
655  curBoundaryAddressing.setSize(nPatches);
656  }
657 
658  Info << "\nDistributing points to processors" << endl;
659  // For every processor, loop through the list of faces for the processor.
660  // For every face, loop through the list of points and mark the point as
661  // used for the processor. Collect the list of used points for the
662  // processor.
663 
664  forAll (procPointAddressing_, procI)
665  {
666  boolList pointLabels(nPoints(), false);
667 
668  // Get reference to list of used faces
669  const labelList& procFaceLabels = procFaceAddressing_[procI];
670 
671  forAll (procFaceLabels, facei)
672  {
673  // Because of the turning index, some labels may be negative
674  const labelList& facePoints = fcs[mag(procFaceLabels[facei]) - 1];
675 
676  forAll (facePoints, pointi)
677  {
678  // Mark the point as used
679  pointLabels[facePoints[pointi]] = true;
680  }
681  }
682 
683  // Collect the used points
684  labelList& procPointLabels = procPointAddressing_[procI];
685 
686  procPointLabels.setSize(pointLabels.size());
687 
688  label nUsedPoints = 0;
689 
690  forAll (pointLabels, pointi)
691  {
692  if (pointLabels[pointi])
693  {
694  procPointLabels[nUsedPoints] = pointi;
695 
696  nUsedPoints++;
697  }
698  }
699 
700  // Reset the size of used points
701  procPointLabels.setSize(nUsedPoints);
702  }
703 
704  // Gather data about globally shared points
705 
706  // Memory management
707  {
708  labelList pointsUsage(nPoints(), 0);
709 
710  // Globally shared points are the ones used by more than 2 processors
711  // Size the list approximately and gather the points
712  labelHashSet gSharedPoints
713  (
714  min(100, nPoints()/1000)
715  );
716 
717  // Loop through all the processors and mark up points used by
718  // processor boundaries. When a point is used twice, it is a
719  // globally shared point
720 
721  for (label procI = 0; procI < nProcs_; procI++)
722  {
723  // Get list of face labels
724  const labelList& curFaceLabels = procFaceAddressing_[procI];
725 
726  // Get start of processor faces
727  const labelList& curProcessorPatchStarts =
728  procProcessorPatchStartIndex_[procI];
729 
730  const labelList& curProcessorPatchSizes =
731  procProcessorPatchSize_[procI];
732 
733  // Reset the lookup list
734  pointsUsage = 0;
735 
736  forAll (curProcessorPatchStarts, patchi)
737  {
738  const label curStart = curProcessorPatchStarts[patchi];
739  const label curEnd = curStart + curProcessorPatchSizes[patchi];
740 
741  for
742  (
743  label facei = curStart;
744  facei < curEnd;
745  facei++
746  )
747  {
748  // Mark the original face as used
749  // Remember to decrement the index by one (turning index)
750  //
751  const label curF = mag(curFaceLabels[facei]) - 1;
752 
753  const face& f = fcs[curF];
754 
755  forAll (f, pointi)
756  {
757  if (pointsUsage[f[pointi]] == 0)
758  {
759  // Point not previously used
760  pointsUsage[f[pointi]] = patchi + 1;
761  }
762  else if (pointsUsage[f[pointi]] != patchi + 1)
763  {
764  // Point used by some other patch = global point!
765  gSharedPoints.insert(f[pointi]);
766  }
767  }
768  }
769  }
770  }
771 
772  // Grab the result from the hash list
773  globallySharedPoints_ = gSharedPoints.toc();
774  sort(globallySharedPoints_);
775  }
776 }