FreeFOAM The Cross-Platform CFD Toolkit
regionSplit.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 "regionSplit.H"
29 #include <OpenFOAM/globalIndex.H>
30 #include <OpenFOAM/syncTools.H>
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 
37 defineTypeNameAndDebug(regionSplit, 0);
38 
39 }
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 // Handle (non-processor) coupled faces.
44 void Foam::regionSplit::transferCoupledFaceRegion
45 (
46  const label faceI,
47  const label otherFaceI,
48 
49  labelList& faceRegion,
50  DynamicList<label>& newChangedFaces
51 ) const
52 {
53  if (faceRegion[faceI] >= 0)
54  {
55  if (faceRegion[otherFaceI] == -1)
56  {
57  faceRegion[otherFaceI] = faceRegion[faceI];
58  newChangedFaces.append(otherFaceI);
59  }
60  else if (faceRegion[otherFaceI] == -2)
61  {
62  // otherFaceI blocked but faceI is not. Is illegal for coupled
63  // faces, not for explicit connections.
64  }
65  else if (faceRegion[otherFaceI] != faceRegion[faceI])
66  {
68  (
69  "regionSplit::transferCoupledFaceRegion"
70  "(const label, const label, labelList&, labelList&) const"
71  ) << "Problem : coupled face " << faceI
72  << " on patch " << mesh_.boundaryMesh().whichPatch(faceI)
73  << " has region " << faceRegion[faceI]
74  << " but coupled face " << otherFaceI
75  << " has region " << faceRegion[otherFaceI]
76  << endl
77  << "Is your blocked faces specification"
78  << " synchronized across coupled boundaries?"
79  << abort(FatalError);
80  }
81  }
82  else if (faceRegion[faceI] == -1)
83  {
84  if (faceRegion[otherFaceI] >= 0)
85  {
86  faceRegion[faceI] = faceRegion[otherFaceI];
87  newChangedFaces.append(faceI);
88  }
89  else if (faceRegion[otherFaceI] == -2)
90  {
91  // otherFaceI blocked but faceI is not. Is illegal for coupled
92  // faces, not for explicit connections.
93  }
94  }
95 }
96 
97 
98 void Foam::regionSplit::fillSeedMask
99 (
100  const List<labelPair>& explicitConnections,
101  labelList& cellRegion,
102  labelList& faceRegion,
103  const label seedCellID,
104  const label markValue
105 ) const
106 {
107  // Do seed cell
108  cellRegion[seedCellID] = markValue;
109 
110 
111  // Collect faces on seed cell
112  const cell& cFaces = mesh_.cells()[seedCellID];
113 
114  label nFaces = 0;
115 
116  labelList changedFaces(cFaces.size());
117 
118  forAll(cFaces, i)
119  {
120  label faceI = cFaces[i];
121 
122  if (faceRegion[faceI] == -1)
123  {
124  faceRegion[faceI] = markValue;
125  changedFaces[nFaces++] = faceI;
126  }
127  }
128  changedFaces.setSize(nFaces);
129 
130 
131  // Loop over changed faces. MeshWave in small.
132 
133  while (changedFaces.size())
134  {
135  //if (debug)
136  //{
137  // Pout<< "regionSplit::fillSeedMask : changedFaces:"
138  // << changedFaces.size() << endl;
139  //}
140 
141  DynamicList<label> changedCells(changedFaces.size());
142 
143  forAll(changedFaces, i)
144  {
145  label faceI = changedFaces[i];
146 
147  label own = mesh_.faceOwner()[faceI];
148 
149  if (cellRegion[own] == -1)
150  {
151  cellRegion[own] = markValue;
152  changedCells.append(own);
153  }
154 
155  if (mesh_.isInternalFace(faceI))
156  {
157  label nei = mesh_.faceNeighbour()[faceI];
158 
159  if (cellRegion[nei] == -1)
160  {
161  cellRegion[nei] = markValue;
162  changedCells.append(nei);
163  }
164  }
165  }
166 
167 
168  //if (debug)
169  //{
170  // Pout<< "regionSplit::fillSeedMask : changedCells:"
171  // << changedCells.size() << endl;
172  //}
173 
174  // Loop over changedCells and collect faces
175  DynamicList<label> newChangedFaces(changedCells.size());
176 
177  forAll(changedCells, i)
178  {
179  label cellI = changedCells[i];
180 
181  const cell& cFaces = mesh_.cells()[cellI];
182 
183  forAll(cFaces, cFaceI)
184  {
185  label faceI = cFaces[cFaceI];
186 
187  if (faceRegion[faceI] == -1)
188  {
189  faceRegion[faceI] = markValue;
190  newChangedFaces.append(faceI);
191  }
192  }
193  }
194 
195 
196  //if (debug)
197  //{
198  // Pout<< "regionSplit::fillSeedMask : changedFaces before sync:"
199  // << changedFaces.size() << endl;
200  //}
201 
202 
203  // Check for changes to any locally coupled face.
204  // Global connections are done later.
205 
206  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
207 
208  forAll(patches, patchI)
209  {
210  const polyPatch& pp = patches[patchI];
211 
212  if (isA<cyclicPolyPatch>(pp))
213  {
214  label faceI = pp.start();
215 
216  label halfSz = pp.size()/2;
217 
218  for (label i = 0; i < halfSz; i++)
219  {
220  label otherFaceI = refCast<const cyclicPolyPatch>(pp)
221  .transformGlobalFace(faceI);
222 
223  transferCoupledFaceRegion
224  (
225  faceI,
226  otherFaceI,
227  faceRegion,
228  newChangedFaces
229  );
230 
231  faceI++;
232  }
233  }
234  }
235  forAll(explicitConnections, i)
236  {
237  transferCoupledFaceRegion
238  (
239  explicitConnections[i][0],
240  explicitConnections[i][1],
241  faceRegion,
242  newChangedFaces
243  );
244  }
245 
246  //if (debug)
247  //{
248  // Pout<< "regionSplit::fillSeedMask : changedFaces after sync:"
249  // << newChangedFaces.size() << endl;
250  //}
251 
252  changedFaces.transfer(newChangedFaces);
253  }
254 }
255 
256 
257 Foam::label Foam::regionSplit::calcRegionSplit
258 (
259  const boolList& blockedFace,
260  const List<labelPair>& explicitConnections,
261 
262  labelList& cellRegion
263 ) const
264 {
265  if (debug)
266  {
267  if (blockedFace.size())
268  {
269  // Check that blockedFace is synced.
270  boolList syncBlockedFace(blockedFace);
271  syncTools::swapFaceList(mesh_, syncBlockedFace, false);
272 
273  forAll(syncBlockedFace, faceI)
274  {
275  if (syncBlockedFace[faceI] != blockedFace[faceI])
276  {
278  (
279  "regionSplit::calcRegionSplit(..)"
280  ) << "Face " << faceI << " not synchronised. My value:"
281  << blockedFace[faceI] << " coupled value:"
282  << syncBlockedFace[faceI]
283  << abort(FatalError);
284  }
285  }
286  }
287  }
288 
289  // Region per face.
290  // -1 unassigned
291  // -2 blocked
292  labelList faceRegion(mesh_.nFaces(), -1);
293 
294  if (blockedFace.size())
295  {
296  forAll(blockedFace, faceI)
297  {
298  if (blockedFace[faceI])
299  {
300  faceRegion[faceI] = -2;
301  }
302  }
303  }
304 
305 
306  // Assign local regions
307  // ~~~~~~~~~~~~~~~~~~~~
308 
309  // Start with region 0
310  label nRegions = 0;
311 
312  label unsetCellI = 0;
313 
314  do
315  {
316  // Find first unset cell
317 
318  for (; unsetCellI < mesh_.nCells(); unsetCellI++)
319  {
320  if (cellRegion[unsetCellI] == -1)
321  {
322  break;
323  }
324  }
325 
326  if (unsetCellI >= mesh_.nCells())
327  {
328  break;
329  }
330 
331  fillSeedMask
332  (
333  explicitConnections,
334  cellRegion,
335  faceRegion,
336  unsetCellI,
337  nRegions
338  );
339 
340  // Current unsetCell has now been handled. Go to next region.
341  nRegions++;
342  unsetCellI++;
343  }
344  while(true);
345 
346 
347  if (debug)
348  {
349  forAll(cellRegion, cellI)
350  {
351  if (cellRegion[cellI] < 0)
352  {
353  FatalErrorIn("regionSplit::calcRegionSplit(..)")
354  << "cell:" << cellI << " region:" << cellRegion[cellI]
355  << abort(FatalError);
356  }
357  }
358 
359  forAll(faceRegion, faceI)
360  {
361  if (faceRegion[faceI] == -1)
362  {
363  FatalErrorIn("regionSplit::calcRegionSplit(..)")
364  << "face:" << faceI << " region:" << faceRegion[faceI]
365  << abort(FatalError);
366  }
367  }
368  }
369 
370 
371 
372  // Assign global regions
373  // ~~~~~~~~~~~~~~~~~~~~~
374  // Offset local regions to create unique global regions.
375 
376  globalIndex globalRegions(nRegions);
377 
378 
379  // Merge global regions
380  // ~~~~~~~~~~~~~~~~~~~~
381  // Regions across non-blocked proc patches get merged.
382  // This will set merged global regions to be the min of both.
383  // (this will create gaps in the global region list so they will get
384  // merged later on)
385 
386  // Map from global to merged global
387  labelList mergedGlobal(identity(globalRegions.size()));
388 
389 
390  // See if any regions get merged. Only nessecary for parallel
391  while (Pstream::parRun())
392  {
393  if (debug)
394  {
395  Pout<< nl << "-- Starting Iteration --" << endl;
396  }
397 
398  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
399 
400  // Send global regions across (or -2 if blocked face)
401  forAll(patches, patchI)
402  {
403  const polyPatch& pp = patches[patchI];
404 
405  if (isA<processorPolyPatch>(pp))
406  {
407  labelList myGlobalRegions(pp.size());
408 
409  label faceI = pp.start();
410 
411  forAll(pp, i)
412  {
413  if (faceRegion[faceI] < 0)
414  {
415  myGlobalRegions[i] = faceRegion[faceI];
416  }
417  else
418  {
419  myGlobalRegions[i] = mergedGlobal
420  [globalRegions.toGlobal(faceRegion[faceI])];
421  }
422 
423  faceI++;
424  }
425 
426  OPstream toProcNbr
427  (
429  refCast<const processorPolyPatch>(pp).neighbProcNo()
430  );
431 
432  toProcNbr << myGlobalRegions;
433  }
434  }
435 
436 
437  // Receive global regions
438 
439  label nMerged = 0;
440 
441  forAll(patches, patchI)
442  {
443  const polyPatch& pp = patches[patchI];
444 
445  if (isA<processorPolyPatch>(pp))
446  {
447  const processorPolyPatch& procPp =
448  refCast<const processorPolyPatch>(pp);
449 
450  IPstream fromProcNbr(Pstream::blocking, procPp.neighbProcNo());
451 
452  labelList nbrRegions(fromProcNbr);
453 
454 
455  // Compare with my regions to see which get merged.
456 
457  label faceI = pp.start();
458 
459  forAll(pp, i)
460  {
461  if
462  (
463  faceRegion[faceI] < 0
464  || nbrRegions[i] < 0
465  )
466  {
467  if (faceRegion[faceI] != nbrRegions[i])
468  {
469  FatalErrorIn("regionSplit::calcRegionSplit(..)")
470  << "On patch:" << pp.name()
471  << " face:" << faceI
472  << " my local region:" << faceRegion[faceI]
473  << " neighbouring region:"
474  << nbrRegions[i] << nl
475  << "Maybe your blockedFaces are not"
476  << " synchronized across coupled faces?"
477  << abort(FatalError);
478  }
479  }
480  else
481  {
482  label uncompactGlobal =
483  globalRegions.toGlobal(faceRegion[faceI]);
484 
485  label myGlobal = mergedGlobal[uncompactGlobal];
486 
487  if (myGlobal != nbrRegions[i])
488  {
489  label minRegion = min(myGlobal, nbrRegions[i]);
490 
491  if (debug)
492  {
493  Pout<< "Merging region " << myGlobal
494  << " (on proc " << Pstream::myProcNo()
495  << ") and region " << nbrRegions[i]
496  << " (on proc " << procPp.neighbProcNo()
497  << ") into region " << minRegion << endl;
498  }
499 
500  mergedGlobal[uncompactGlobal] = minRegion;
501  mergedGlobal[myGlobal] = minRegion;
502  mergedGlobal[nbrRegions[i]] = minRegion;
503 
504  nMerged++;
505  }
506  }
507 
508  faceI++;
509  }
510  }
511  }
512 
513 
514  reduce(nMerged, sumOp<label>());
515 
516  if (debug)
517  {
518  Pout<< "nMerged:" << nMerged << endl;
519  }
520 
521  if (nMerged == 0)
522  {
523  break;
524  }
525 
526  // Merge the compacted regions.
527  Pstream::listCombineGather(mergedGlobal, minEqOp<label>());
528  Pstream::listCombineScatter(mergedGlobal);
529  }
530 
531 
532  // Compact global regions
533  // ~~~~~~~~~~~~~~~~~~~~~~
534 
535  // All procs will have the same global mergedGlobal region.
536  // There might be gaps in it however so compact.
537 
538  labelList mergedToCompacted(globalRegions.size(), -1);
539 
540  label compactI = 0;
541 
542  forAll(mergedGlobal, i)
543  {
544  label merged = mergedGlobal[i];
545 
546  if (mergedToCompacted[merged] == -1)
547  {
548  mergedToCompacted[merged] = compactI++;
549  }
550  }
551 
552  if (debug)
553  {
554  Pout<< "Compacted down to " << compactI << " regions." << endl;
555  }
556 
557  // Renumber cellRegion to be global regions
558  forAll(cellRegion, cellI)
559  {
560  label region = cellRegion[cellI];
561 
562  if (region >= 0)
563  {
564  label merged = mergedGlobal[globalRegions.toGlobal(region)];
565 
566  cellRegion[cellI] = mergedToCompacted[merged];
567  }
568  }
569 
570  return compactI;
571 }
572 
573 
574 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
575 
577 :
578  labelList(mesh.nCells(), -1),
579  mesh_(mesh),
580  nRegions_(calcRegionSplit(boolList(0, false), List<labelPair>(0), *this))
581 {}
582 
583 
585 (
586  const polyMesh& mesh,
587  const boolList& blockedFace
588 )
589 :
590  labelList(mesh.nCells(), -1),
591  mesh_(mesh),
592  nRegions_(calcRegionSplit(blockedFace, List<labelPair>(0), *this))
593 {}
594 
595 
597 (
598  const polyMesh& mesh,
599  const boolList& blockedFace,
600  const List<labelPair>& explicitConnections
601 )
602 :
603  labelList(mesh.nCells(), -1),
604  mesh_(mesh),
605  nRegions_(calcRegionSplit(blockedFace, explicitConnections, *this))
606 {}
607 
608 
609 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
610 
611 
612 // ************************ vim: set sw=4 sts=4 et: ************************ //