FreeFOAM The Cross-Platform CFD Toolkit
meshRefinementBaffles.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 "meshRefinement.H"
27 #include <finiteVolume/fvMesh.H>
28 #include <OpenFOAM/syncTools.H>
29 #include <OpenFOAM/Time.H>
31 #include <meshTools/pointSet.H>
32 #include <meshTools/faceSet.H>
35 #include <meshTools/meshTools.H>
43 #include <OpenFOAM/OFstream.H>
44 #include <meshTools/regionSplit.H>
46 
47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48 
49 // Repatches external face or creates baffle for internal face
50 // with user specified patches (might be different for both sides).
51 // Returns label of added face.
52 Foam::label Foam::meshRefinement::createBaffle
53 (
54  const label faceI,
55  const label ownPatch,
56  const label neiPatch,
57  polyTopoChange& meshMod
58 ) const
59 {
60  const face& f = mesh_.faces()[faceI];
61  label zoneID = mesh_.faceZones().whichZone(faceI);
62  bool zoneFlip = false;
63 
64  if (zoneID >= 0)
65  {
66  const faceZone& fZone = mesh_.faceZones()[zoneID];
67  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
68  }
69 
70  meshMod.setAction
71  (
72  polyModifyFace
73  (
74  f, // modified face
75  faceI, // label of face
76  mesh_.faceOwner()[faceI], // owner
77  -1, // neighbour
78  false, // face flip
79  ownPatch, // patch for face
80  false, // remove from zone
81  zoneID, // zone for face
82  zoneFlip // face flip in zone
83  )
84  );
85 
86 
87  label dupFaceI = -1;
88 
89  if (mesh_.isInternalFace(faceI))
90  {
91  if (neiPatch == -1)
92  {
94  (
95  "meshRefinement::createBaffle"
96  "(const label, const label, const label, polyTopoChange&)"
97  ) << "No neighbour patch for internal face " << faceI
98  << " fc:" << mesh_.faceCentres()[faceI]
99  << " ownPatch:" << ownPatch << abort(FatalError);
100  }
101 
102  bool reverseFlip = false;
103  if (zoneID >= 0)
104  {
105  reverseFlip = !zoneFlip;
106  }
107 
108  dupFaceI = meshMod.setAction
109  (
110  polyAddFace
111  (
112  f.reverseFace(), // modified face
113  mesh_.faceNeighbour()[faceI],// owner
114  -1, // neighbour
115  -1, // masterPointID
116  -1, // masterEdgeID
117  faceI, // masterFaceID,
118  true, // face flip
119  neiPatch, // patch for face
120  zoneID, // zone for face
121  reverseFlip // face flip in zone
122  )
123  );
124  }
125  return dupFaceI;
126 }
127 
128 
129 // Get an estimate for the patch the internal face should get. Bit heuristic.
130 Foam::label Foam::meshRefinement::getBafflePatch
131 (
132  const labelList& facePatch,
133  const label faceI
134 ) const
135 {
136  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
137 
138  // Loop over face points
139  // for each point check all faces patch IDs
140  // as soon as an ID >= 0 is found, break and assign that ID
141  // to the current face.
142  // Check first for real patch (so proper surface intersection and then
143  // in facePatch array for patches to block off faces
144 
145  forAll(mesh_.faces()[faceI], fp)
146  {
147  label pointI = mesh_.faces()[faceI][fp];
148 
149  forAll(mesh_.pointFaces()[pointI], pf)
150  {
151  label pFaceI = mesh_.pointFaces()[pointI][pf];
152 
153  label patchI = patches.whichPatch(pFaceI);
154 
155  if (patchI != -1 && !patches[patchI].coupled())
156  {
157  return patchI;
158  }
159  else if (facePatch[pFaceI] != -1)
160  {
161  return facePatch[pFaceI];
162  }
163  }
164  }
165 
166  // Loop over owner and neighbour cells, looking for the first face with a
167  // valid patch number
168  const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
169 
170  forAll(ownFaces, i)
171  {
172  label cFaceI = ownFaces[i];
173 
174  label patchI = patches.whichPatch(cFaceI);
175 
176  if (patchI != -1 && !patches[patchI].coupled())
177  {
178  return patchI;
179  }
180  else if (facePatch[cFaceI] != -1)
181  {
182  return facePatch[cFaceI];
183  }
184  }
185 
186  if (mesh_.isInternalFace(faceI))
187  {
188  const cell& neiFaces = mesh_.cells()[mesh_.faceNeighbour()[faceI]];
189 
190  forAll(neiFaces, i)
191  {
192  label cFaceI = neiFaces[i];
193 
194  label patchI = patches.whichPatch(cFaceI);
195 
196  if (patchI != -1 && !patches[patchI].coupled())
197  {
198  return patchI;
199  }
200  else if (facePatch[cFaceI] != -1)
201  {
202  return facePatch[cFaceI];
203  }
204  }
205  }
206 
207  WarningIn
208  (
209  "meshRefinement::getBafflePatch(const labelList&, const label)"
210  ) << "Could not find boundary face neighbouring internal face "
211  << faceI << " with face centre " << mesh_.faceCentres()[faceI]
212  << nl
213  << "Using arbitrary patch " << 0 << " instead." << endl;
214 
215  return 0;
216 }
217 
218 
219 // Determine patches for baffles.
220 void Foam::meshRefinement::getBafflePatches
221 (
222  const labelList& globalToPatch,
223  const labelList& neiLevel,
224  const pointField& neiCc,
225 
226  labelList& ownPatch,
227  labelList& neiPatch
228 ) const
229 {
230  autoPtr<OFstream> str;
231  label vertI = 0;
232  if (debug&OBJINTERSECTIONS)
233  {
234  str.reset
235  (
236  new OFstream
237  (
238  mesh_.time().path()/timeName()/"intersections.obj"
239  )
240  );
241 
242  Pout<< "getBafflePatches : Writing surface intersections to file "
243  << str().name() << nl << endl;
244  }
245 
246  const pointField& cellCentres = mesh_.cellCentres();
247 
248  // Build list of surfaces that are not to be baffled.
249  const wordList& cellZoneNames = surfaces_.cellZoneNames();
250 
251  labelList surfacesToBaffle(cellZoneNames.size());
252  label baffleI = 0;
253  forAll(cellZoneNames, surfI)
254  {
255  if (cellZoneNames[surfI].size())
256  {
257  if (debug)
258  {
259  Pout<< "getBafflePatches : Not baffling surface "
260  << surfaces_.names()[surfI] << endl;
261  }
262  }
263  else
264  {
265  surfacesToBaffle[baffleI++] = surfI;
266  }
267  }
268  surfacesToBaffle.setSize(baffleI);
269 
270  ownPatch.setSize(mesh_.nFaces());
271  ownPatch = -1;
272  neiPatch.setSize(mesh_.nFaces());
273  neiPatch = -1;
274 
275 
276  // Collect candidate faces
277  // ~~~~~~~~~~~~~~~~~~~~~~~
278 
279  labelList testFaces(intersectedFaces());
280 
281  // Collect segments
282  // ~~~~~~~~~~~~~~~~
283 
284  pointField start(testFaces.size());
285  pointField end(testFaces.size());
286 
287  forAll(testFaces, i)
288  {
289  label faceI = testFaces[i];
290 
291  label own = mesh_.faceOwner()[faceI];
292 
293  if (mesh_.isInternalFace(faceI))
294  {
295  start[i] = cellCentres[own];
296  end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
297  }
298  else
299  {
300  start[i] = cellCentres[own];
301  end[i] = neiCc[faceI-mesh_.nInternalFaces()];
302  }
303  }
304 
305 
306  // Do test for intersections
307  // ~~~~~~~~~~~~~~~~~~~~~~~~~
308  labelList surface1;
309  List<pointIndexHit> hit1;
310  labelList region1;
311  labelList surface2;
312  List<pointIndexHit> hit2;
313  labelList region2;
314  surfaces_.findNearestIntersection
315  (
316  surfacesToBaffle,
317  start,
318  end,
319 
320  surface1,
321  hit1,
322  region1,
323 
324  surface2,
325  hit2,
326  region2
327  );
328 
329  forAll(testFaces, i)
330  {
331  label faceI = testFaces[i];
332 
333  if (hit1[i].hit() && hit2[i].hit())
334  {
335  if (str.valid())
336  {
337  meshTools::writeOBJ(str(), start[i]);
338  vertI++;
339  meshTools::writeOBJ(str(), hit1[i].rawPoint());
340  vertI++;
341  meshTools::writeOBJ(str(), hit2[i].rawPoint());
342  vertI++;
343  meshTools::writeOBJ(str(), end[i]);
344  vertI++;
345  str()<< "l " << vertI-3 << ' ' << vertI-2 << nl;
346  str()<< "l " << vertI-2 << ' ' << vertI-1 << nl;
347  str()<< "l " << vertI-1 << ' ' << vertI << nl;
348  }
349 
350  // Pick up the patches
351  ownPatch[faceI] = globalToPatch
352  [
353  surfaces_.globalRegion(surface1[i], region1[i])
354  ];
355  neiPatch[faceI] = globalToPatch
356  [
357  surfaces_.globalRegion(surface2[i], region2[i])
358  ];
359 
360  if (ownPatch[faceI] == -1 || neiPatch[faceI] == -1)
361  {
362  FatalErrorIn("getBafflePatches(..)")
363  << "problem." << abort(FatalError);
364  }
365  }
366  }
367 
368  // No need to parallel sync since intersection data (surfaceIndex_ etc.)
369  // already guaranteed to be synced...
370  // However:
371  // - owncc and neicc are reversed on different procs so might pick
372  // up different regions reversed? No problem. Neighbour on one processor
373  // might not be owner on the other processor but the neighbour is
374  // not used when creating baffles from proc faces.
375  // - tolerances issues occasionally crop up.
376  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
377  syncTools::syncFaceList(mesh_, neiPatch, maxEqOp<label>(), false);
378 }
379 
380 
381 // Create baffle for every face where ownPatch != -1
383 (
384  const labelList& ownPatch,
385  const labelList& neiPatch
386 )
387 {
388  if
389  (
390  ownPatch.size() != mesh_.nFaces()
391  || neiPatch.size() != mesh_.nFaces()
392  )
393  {
395  (
396  "meshRefinement::createBaffles"
397  "(const labelList&, const labelList&)"
398  ) << "Illegal size :"
399  << " ownPatch:" << ownPatch.size()
400  << " neiPatch:" << neiPatch.size()
401  << ". Should be number of faces:" << mesh_.nFaces()
402  << abort(FatalError);
403  }
404 
405  if (debug)
406  {
407  labelList syncedOwnPatch(ownPatch);
408  syncTools::syncFaceList(mesh_, syncedOwnPatch, maxEqOp<label>(), false);
409  labelList syncedNeiPatch(neiPatch);
410  syncTools::syncFaceList(mesh_, syncedNeiPatch, maxEqOp<label>(), false);
411 
412  forAll(syncedOwnPatch, faceI)
413  {
414  if
415  (
416  (ownPatch[faceI] == -1 && syncedOwnPatch[faceI] != -1)
417  || (neiPatch[faceI] == -1 && syncedNeiPatch[faceI] != -1)
418  )
419  {
421  (
422  "meshRefinement::createBaffles"
423  "(const labelList&, const labelList&)"
424  ) << "Non synchronised at face:" << faceI
425  << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
426  << " fc:" << mesh_.faceCentres()[faceI] << endl
427  << "ownPatch:" << ownPatch[faceI]
428  << " syncedOwnPatch:" << syncedOwnPatch[faceI]
429  << " neiPatch:" << neiPatch[faceI]
430  << " syncedNeiPatch:" << syncedNeiPatch[faceI]
431  << abort(FatalError);
432  }
433  }
434  }
435 
436  // Topochange container
437  polyTopoChange meshMod(mesh_);
438 
439  label nBaffles = 0;
440 
441  forAll(ownPatch, faceI)
442  {
443  if (ownPatch[faceI] != -1)
444  {
445  // Create baffle or repatch face. Return label of inserted baffle
446  // face.
447  createBaffle
448  (
449  faceI,
450  ownPatch[faceI], // owner side patch
451  neiPatch[faceI], // neighbour side patch
452  meshMod
453  );
454  nBaffles++;
455  }
456  }
457 
458  // Change the mesh (no inflation, parallel sync)
459  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
460 
461  // Update fields
462  mesh_.updateMesh(map);
463 
464  // Move mesh if in inflation mode
465  if (map().hasMotionPoints())
466  {
467  mesh_.movePoints(map().preMotionPoints());
468  }
469  else
470  {
471  // Delete mesh volumes.
472  mesh_.clearOut();
473  }
474 
475  if (overwrite())
476  {
477  mesh_.setInstance(oldInstance());
478  }
479 
480  //- Redo the intersections on the newly create baffle faces. Note that
481  // this changes also the cell centre positions.
482  faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
483 
484  const labelList& reverseFaceMap = map().reverseFaceMap();
485  const labelList& faceMap = map().faceMap();
486 
487  // Pick up owner side of baffle
488  forAll(ownPatch, oldFaceI)
489  {
490  label faceI = reverseFaceMap[oldFaceI];
491 
492  if (ownPatch[oldFaceI] != -1 && faceI >= 0)
493  {
494  const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
495 
496  forAll(ownFaces, i)
497  {
498  baffledFacesSet.insert(ownFaces[i]);
499  }
500  }
501  }
502  // Pick up neighbour side of baffle (added faces)
503  forAll(faceMap, faceI)
504  {
505  label oldFaceI = faceMap[faceI];
506 
507  if (oldFaceI >= 0 && reverseFaceMap[oldFaceI] != faceI)
508  {
509  const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
510 
511  forAll(ownFaces, i)
512  {
513  baffledFacesSet.insert(ownFaces[i]);
514  }
515  }
516  }
517  baffledFacesSet.sync(mesh_);
518 
519  updateMesh(map, baffledFacesSet.toc());
520 
521  return map;
522 }
523 
524 
525 // Return a list of coupled face pairs, i.e. faces that use the same vertices.
526 // (this information is recalculated instead of maintained since would be too
527 // hard across splitMeshRegions).
529 (
530  const labelList& testFaces
531 ) const
532 {
533  labelList duplicateFace
534  (
536  (
537  mesh_,
538  testFaces
539  )
540  );
541 
542  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
543 
544  // Convert into list of coupled face pairs (mesh face labels).
545  List<labelPair> duplicateFaces(testFaces.size());
546  label dupI = 0;
547 
548  forAll(duplicateFace, i)
549  {
550  label otherFaceI = duplicateFace[i];
551 
552  if (otherFaceI != -1 && i < otherFaceI)
553  {
554  label meshFace0 = testFaces[i];
555  label patch0 = patches.whichPatch(meshFace0);
556  label meshFace1 = testFaces[otherFaceI];
557  label patch1 = patches.whichPatch(meshFace1);
558 
559  if
560  (
561  (patch0 != -1 && isA<processorPolyPatch>(patches[patch0]))
562  || (patch1 != -1 && isA<processorPolyPatch>(patches[patch1]))
563  )
564  {
566  (
567  "meshRefinement::getDuplicateFaces"
568  "(const bool, const labelList&)"
569  ) << "One of two duplicate faces is on"
570  << " processorPolyPatch."
571  << "This is not allowed." << nl
572  << "Face:" << meshFace0
573  << " is on patch:" << patches[patch0].name()
574  << nl
575  << "Face:" << meshFace1
576  << " is on patch:" << patches[patch1].name()
577  << abort(FatalError);
578  }
579 
580  duplicateFaces[dupI++] = labelPair(meshFace0, meshFace1);
581  }
582  }
583  duplicateFaces.setSize(dupI);
584 
585  Info<< "getDuplicateFaces : found " << returnReduce(dupI, sumOp<label>())
586  << " pairs of duplicate faces." << nl << endl;
587 
588 
589  if (debug)
590  {
591  faceSet duplicateFaceSet(mesh_, "duplicateFaces", 2*dupI);
592 
593  forAll(duplicateFaces, i)
594  {
595  duplicateFaceSet.insert(duplicateFaces[i][0]);
596  duplicateFaceSet.insert(duplicateFaces[i][1]);
597  }
598  Pout<< "Writing duplicate faces (baffles) to faceSet "
599  << duplicateFaceSet.name() << nl << endl;
600  duplicateFaceSet.write();
601  }
602 
603  return duplicateFaces;
604 }
605 
606 
607 // Extract those baffles (duplicate) faces that are on the edge of a baffle
608 // region. These are candidates for merging.
609 // Done by counting the number of baffles faces per mesh edge. If edge
610 // has 2 boundary faces and both are baffle faces it is the edge of a baffle
611 // region.
612 Foam::List<Foam::labelPair> Foam::meshRefinement::filterDuplicateFaces
613 (
614  const List<labelPair>& couples
615 ) const
616 {
617  // Construct addressing engine for all duplicate faces (only one
618  // for each pair)
619 
620  // Duplicate faces in mesh labels (first face of each pair only)
621  // (reused later on to mark off filtered couples. see below)
622  labelList duplicateFaces(couples.size());
623 
624 
625  // All duplicate faces on edge of the patch are to be merged.
626  // So we count for all edges of duplicate faces how many duplicate
627  // faces use them.
628  labelList nBafflesPerEdge(mesh_.nEdges(), 0);
629 
630 
631 
632  // Count number of boundary faces per edge
633  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
634 
635  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
636 
637  forAll(patches, patchI)
638  {
639  const polyPatch& pp = patches[patchI];
640 
641  // Count number of boundary faces. Discard coupled boundary faces.
642  if (!pp.coupled())
643  {
644  label faceI = pp.start();
645 
646  forAll(pp, i)
647  {
648  const labelList& fEdges = mesh_.faceEdges(faceI);
649 
650  forAll(fEdges, fEdgeI)
651  {
652  nBafflesPerEdge[fEdges[fEdgeI]]++;
653  }
654  faceI++;
655  }
656  }
657  }
658 
659 
660  DynamicList<label> fe0;
661  DynamicList<label> fe1;
662 
663 
664  // Count number of duplicate boundary faces per edge
665  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
666 
667  forAll(couples, i)
668  {
669  const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0);
670 
671  forAll(fEdges0, fEdgeI)
672  {
673  nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000;
674  }
675 
676  const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1);
677 
678  forAll(fEdges1, fEdgeI)
679  {
680  nBafflesPerEdge[fEdges1[fEdgeI]] += 1000000;
681  }
682  }
683 
684  // Add nBaffles on shared edges
686  (
687  mesh_,
688  nBafflesPerEdge,
689  plusEqOp<label>(), // in-place add
690  0, // initial value
691  false // no separation
692  );
693 
694 
695  // Baffles which are not next to other boundaries and baffles will have
696  // value 2*1000000+2*1
697 
698  List<labelPair> filteredCouples(couples.size());
699  label filterI = 0;
700 
701  forAll(couples, i)
702  {
703  const labelPair& couple = couples[i];
704 
705  if
706  (
707  patches.whichPatch(couple.first())
708  == patches.whichPatch(couple.second())
709  )
710  {
711  const labelList& fEdges = mesh_.faceEdges(couples[i].first());
712 
713  forAll(fEdges, fEdgeI)
714  {
715  label edgeI = fEdges[fEdgeI];
716 
717  if (nBafflesPerEdge[edgeI] == 2*1000000+2*1)
718  {
719  filteredCouples[filterI++] = couples[i];
720  break;
721  }
722  }
723  }
724  }
725  filteredCouples.setSize(filterI);
726 
727  //Info<< "filterDuplicateFaces : from "
728  // << returnReduce(couples.size(), sumOp<label>())
729  // << " down to "
730  // << returnReduce(filteredCouples.size(), sumOp<label>())
731  // << " baffles." << nl << endl;
732 
733  return filteredCouples;
734 }
735 
736 
737 // Merge baffles. Gets pairs of faces.
739 (
740  const List<labelPair>& couples
741 )
742 {
743  // Mesh change engine
744  polyTopoChange meshMod(mesh_);
745 
746  const faceList& faces = mesh_.faces();
747  const labelList& faceOwner = mesh_.faceOwner();
748  const faceZoneMesh& faceZones = mesh_.faceZones();
749 
750  forAll(couples, i)
751  {
752  label face0 = couples[i].first();
753  label face1 = couples[i].second();
754 
755  // face1 < 0 signals a coupled face that has been converted to baffle.
756 
757  label own0 = faceOwner[face0];
758  label own1 = faceOwner[face1];
759 
760  if (face1 < 0 || own0 < own1)
761  {
762  // Use face0 as the new internal face.
763  label zoneID = faceZones.whichZone(face0);
764  bool zoneFlip = false;
765 
766  if (zoneID >= 0)
767  {
768  const faceZone& fZone = faceZones[zoneID];
769  zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
770  }
771 
772  label nei = (face1 < 0 ? -1 : own1);
773 
774  meshMod.setAction(polyRemoveFace(face1));
775  meshMod.setAction
776  (
778  (
779  faces[face0], // modified face
780  face0, // label of face being modified
781  own0, // owner
782  nei, // neighbour
783  false, // face flip
784  -1, // patch for face
785  false, // remove from zone
786  zoneID, // zone for face
787  zoneFlip // face flip in zone
788  )
789  );
790  }
791  else
792  {
793  // Use face1 as the new internal face.
794  label zoneID = faceZones.whichZone(face1);
795  bool zoneFlip = false;
796 
797  if (zoneID >= 0)
798  {
799  const faceZone& fZone = faceZones[zoneID];
800  zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
801  }
802 
803  meshMod.setAction(polyRemoveFace(face0));
804  meshMod.setAction
805  (
807  (
808  faces[face1], // modified face
809  face1, // label of face being modified
810  own1, // owner
811  own0, // neighbour
812  false, // face flip
813  -1, // patch for face
814  false, // remove from zone
815  zoneID, // zone for face
816  zoneFlip // face flip in zone
817  )
818  );
819  }
820  }
821 
822  // Change the mesh (no inflation)
823  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
824 
825  // Update fields
826  mesh_.updateMesh(map);
827 
828  // Move mesh (since morphing does not do this)
829  if (map().hasMotionPoints())
830  {
831  mesh_.movePoints(map().preMotionPoints());
832  }
833  else
834  {
835  // Delete mesh volumes.
836  mesh_.clearOut();
837  }
838 
839  if (overwrite())
840  {
841  mesh_.setInstance(oldInstance());
842  }
843 
844  // Update intersections. Recalculate intersections on merged faces since
845  // this seems to give problems? Note: should not be nessecary since
846  // baffles preserve intersections from when they were created.
847  labelList newExposedFaces(2*couples.size());
848  label newI = 0;
849 
850  forAll(couples, i)
851  {
852  label newFace0 = map().reverseFaceMap()[couples[i].first()];
853  if (newFace0 != -1)
854  {
855  newExposedFaces[newI++] = newFace0;
856  }
857 
858  label newFace1 = map().reverseFaceMap()[couples[i].second()];
859  if (newFace1 != -1)
860  {
861  newExposedFaces[newI++] = newFace1;
862  }
863  }
864  newExposedFaces.setSize(newI);
865  updateMesh(map, newExposedFaces);
866 
867  return map;
868 }
869 
870 
871 // Finds region per cell for cells inside closed named surfaces
872 void Foam::meshRefinement::findCellZoneGeometric
873 (
874  const labelList& closedNamedSurfaces, // indices of closed surfaces
875  labelList& namedSurfaceIndex, // per face index of named surface
876  const labelList& surfaceToCellZone, // cell zone index per surface
877 
878  labelList& cellToZone
879 ) const
880 {
881  const pointField& cellCentres = mesh_.cellCentres();
882  const labelList& faceOwner = mesh_.faceOwner();
883  const labelList& faceNeighbour = mesh_.faceNeighbour();
884 
885  // Check if cell centre is inside
886  labelList insideSurfaces;
887  surfaces_.findInside
888  (
889  closedNamedSurfaces,
890  cellCentres,
891  insideSurfaces
892  );
893 
894  forAll(insideSurfaces, cellI)
895  {
896  if (cellToZone[cellI] == -2)
897  {
898  label surfI = insideSurfaces[cellI];
899 
900  if (surfI != -1)
901  {
902  cellToZone[cellI] = surfaceToCellZone[surfI];
903  }
904  }
905  }
906 
907 
908  // Some cells with cell centres close to surface might have
909  // had been put into wrong surface. Recheck with perturbed cell centre.
910 
911 
912  // 1. Collect points
913 
914  // Count points to test.
915  label nCandidates = 0;
916  forAll(namedSurfaceIndex, faceI)
917  {
918  label surfI = namedSurfaceIndex[faceI];
919 
920  if (surfI != -1)
921  {
922  if (mesh_.isInternalFace(faceI))
923  {
924  nCandidates += 2;
925  }
926  else
927  {
928  nCandidates += 1;
929  }
930  }
931  }
932 
933  // Collect points.
934  pointField candidatePoints(nCandidates);
935  nCandidates = 0;
936  forAll(namedSurfaceIndex, faceI)
937  {
938  label surfI = namedSurfaceIndex[faceI];
939 
940  if (surfI != -1)
941  {
942  label own = faceOwner[faceI];
943  const point& ownCc = cellCentres[own];
944 
945  if (mesh_.isInternalFace(faceI))
946  {
947  label nei = faceNeighbour[faceI];
948  const point& neiCc = cellCentres[nei];
949  // Perturbed cc
950  const vector d = 1E-4*(neiCc - ownCc);
951  candidatePoints[nCandidates++] = ownCc-d;
952  candidatePoints[nCandidates++] = neiCc+d;
953  }
954  else
955  {
956  const point& neiFc = mesh_.faceCentres()[faceI];
957  // Perturbed cc
958  const vector d = 1E-4*(neiFc - ownCc);
959  candidatePoints[nCandidates++] = ownCc-d;
960  }
961  }
962  }
963 
964 
965  // 2. Test points for inside
966 
967  surfaces_.findInside
968  (
969  closedNamedSurfaces,
970  candidatePoints,
971  insideSurfaces
972  );
973 
974 
975  // 3. Update zone information
976 
977  nCandidates = 0;
978  forAll(namedSurfaceIndex, faceI)
979  {
980  label surfI = namedSurfaceIndex[faceI];
981 
982  if (surfI != -1)
983  {
984  label own = faceOwner[faceI];
985 
986  if (mesh_.isInternalFace(faceI))
987  {
988  label ownSurfI = insideSurfaces[nCandidates++];
989  if (ownSurfI != -1)
990  {
991  cellToZone[own] = surfaceToCellZone[ownSurfI];
992  }
993 
994  label neiSurfI = insideSurfaces[nCandidates++];
995  if (neiSurfI != -1)
996  {
997  label nei = faceNeighbour[faceI];
998 
999  cellToZone[nei] = surfaceToCellZone[neiSurfI];
1000  }
1001  }
1002  else
1003  {
1004  label ownSurfI = insideSurfaces[nCandidates++];
1005  if (ownSurfI != -1)
1006  {
1007  cellToZone[own] = surfaceToCellZone[ownSurfI];
1008  }
1009  }
1010  }
1011  }
1012 
1013 
1014  // Adapt the namedSurfaceIndex
1015  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1016  // for if any cells were not completely covered.
1017 
1018  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1019  {
1020  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1021  label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
1022 
1023  if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
1024  {
1025  // Give face the zone of max cell zone
1026  namedSurfaceIndex[faceI] = findIndex
1027  (
1028  surfaceToCellZone,
1029  max(ownZone, neiZone)
1030  );
1031  }
1032  }
1033 
1034  labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
1035  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1036 
1037  forAll(patches, patchI)
1038  {
1039  const polyPatch& pp = patches[patchI];
1040 
1041  if (pp.coupled())
1042  {
1043  forAll(pp, i)
1044  {
1045  label faceI = pp.start()+i;
1046  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1047  neiCellZone[faceI-mesh_.nInternalFaces()] = ownZone;
1048  }
1049  }
1050  }
1051  syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
1052 
1053  forAll(patches, patchI)
1054  {
1055  const polyPatch& pp = patches[patchI];
1056 
1057  if (pp.coupled())
1058  {
1059  forAll(pp, i)
1060  {
1061  label faceI = pp.start()+i;
1062  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1063  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
1064 
1065  if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
1066  {
1067  // Give face the max cell zone
1068  namedSurfaceIndex[faceI] = findIndex
1069  (
1070  surfaceToCellZone,
1071  max(ownZone, neiZone)
1072  );
1073  }
1074  }
1075  }
1076  }
1077 
1078  // Sync
1080  (
1081  mesh_,
1082  namedSurfaceIndex,
1083  maxEqOp<label>(),
1084  false
1085  );
1086 }
1087 
1088 
1089 bool Foam::meshRefinement::calcRegionToZone
1090 (
1091  const label surfZoneI,
1092  const label ownRegion,
1093  const label neiRegion,
1094 
1095  labelList& regionToCellZone
1096 ) const
1097 {
1098  bool changed = false;
1099 
1100  // Check whether inbetween different regions
1101  if (ownRegion != neiRegion)
1102  {
1103  // Jump. Change one of the sides to my type.
1104 
1105  // 1. Interface between my type and unset region.
1106  // Set region to keepRegion
1107 
1108  if (regionToCellZone[ownRegion] == -2)
1109  {
1110  if (regionToCellZone[neiRegion] == surfZoneI)
1111  {
1112  // Face between unset and my region. Put unset
1113  // region into keepRegion
1114  regionToCellZone[ownRegion] = -1;
1115  changed = true;
1116  }
1117  else if (regionToCellZone[neiRegion] != -2)
1118  {
1119  // Face between unset and other region.
1120  // Put unset region into my region
1121  regionToCellZone[ownRegion] = surfZoneI;
1122  changed = true;
1123  }
1124  }
1125  else if (regionToCellZone[neiRegion] == -2)
1126  {
1127  if (regionToCellZone[ownRegion] == surfZoneI)
1128  {
1129  // Face between unset and my region. Put unset
1130  // region into keepRegion
1131  regionToCellZone[neiRegion] = -1;
1132  changed = true;
1133  }
1134  else if (regionToCellZone[ownRegion] != -2)
1135  {
1136  // Face between unset and other region.
1137  // Put unset region into my region
1138  regionToCellZone[neiRegion] = surfZoneI;
1139  changed = true;
1140  }
1141  }
1142  }
1143  return changed;
1144 }
1145 
1146 
1147 // Finds region per cell. Assumes:
1148 // - region containing keepPoint does not go into a cellZone
1149 // - all other regions can be found by crossing faces marked in
1150 // namedSurfaceIndex.
1151 void Foam::meshRefinement::findCellZoneTopo
1152 (
1153  const point& keepPoint,
1154  const labelList& namedSurfaceIndex,
1155  const labelList& surfaceToCellZone,
1156  labelList& cellToZone
1157 ) const
1158 {
1159  // Analyse regions. Reuse regionsplit
1160  boolList blockedFace(mesh_.nFaces());
1161 
1162  forAll(namedSurfaceIndex, faceI)
1163  {
1164  if (namedSurfaceIndex[faceI] == -1)
1165  {
1166  blockedFace[faceI] = false;
1167  }
1168  else
1169  {
1170  blockedFace[faceI] = true;
1171  }
1172  }
1173  // No need to sync since namedSurfaceIndex already is synced
1174 
1175  // Set region per cell based on walking
1176  regionSplit cellRegion(mesh_, blockedFace);
1177  blockedFace.clear();
1178 
1179  // Per mesh region the zone the cell should be put in.
1180  // -2 : not analysed yet
1181  // -1 : keepPoint region. Not put into any cellzone.
1182  // >= 0 : index of cellZone
1183  labelList regionToCellZone(cellRegion.nRegions(), -2);
1184 
1185  // See which cells already are set in the cellToZone (from geometric
1186  // searching) and use these to take over their zones.
1187  // Note: could be improved to count number of cells per region.
1188  forAll(cellToZone, cellI)
1189  {
1190  if (cellToZone[cellI] != -2)
1191  {
1192  regionToCellZone[cellRegion[cellI]] = cellToZone[cellI];
1193  }
1194  }
1195 
1196 
1197 
1198  // Find the region containing the keepPoint
1199  label keepRegionI = -1;
1200 
1201  label cellI = mesh_.findCell(keepPoint);
1202 
1203  if (cellI != -1)
1204  {
1205  keepRegionI = cellRegion[cellI];
1206  }
1207  reduce(keepRegionI, maxOp<label>());
1208 
1209  Info<< "Found point " << keepPoint << " in cell " << cellI
1210  << " in global region " << keepRegionI
1211  << " out of " << cellRegion.nRegions() << " regions." << endl;
1212 
1213  if (keepRegionI == -1)
1214  {
1215  FatalErrorIn
1216  (
1217  "meshRefinement::findCellZoneTopo"
1218  "(const point&, const labelList&, const labelList&, labelList&)"
1219  ) << "Point " << keepPoint
1220  << " is not inside the mesh." << nl
1221  << "Bounding box of the mesh:" << mesh_.globalData().bb()
1222  << exit(FatalError);
1223  }
1224 
1225  // Mark default region with zone -1.
1226  if (regionToCellZone[keepRegionI] == -2)
1227  {
1228  regionToCellZone[keepRegionI] = -1;
1229  }
1230 
1231 
1232  // Find correspondence between cell zone and surface
1233  // by changing cell zone every time we cross a surface.
1234  while (true)
1235  {
1236  // Synchronise regionToCellZone.
1237  // Note:
1238  // - region numbers are identical on all processors
1239  // - keepRegion is identical ,,
1240  // - cellZones are identical ,,
1241  // This done at top of loop to account for geometric matching
1242  // not being synchronised.
1243  Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
1244  Pstream::listCombineScatter(regionToCellZone);
1245 
1246 
1247  bool changed = false;
1248 
1249  // Internal faces
1250 
1251  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1252  {
1253  label surfI = namedSurfaceIndex[faceI];
1254 
1255  if (surfI != -1)
1256  {
1257  // Calculate region to zone from cellRegions on either side
1258  // of internal face.
1259  bool changedCell = calcRegionToZone
1260  (
1261  surfaceToCellZone[surfI],
1262  cellRegion[mesh_.faceOwner()[faceI]],
1263  cellRegion[mesh_.faceNeighbour()[faceI]],
1264  regionToCellZone
1265  );
1266 
1267  changed = changed | changedCell;
1268  }
1269  }
1270 
1271  // Boundary faces
1272 
1273  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1274 
1275  // Get coupled neighbour cellRegion
1276  labelList neiCellRegion(mesh_.nFaces()-mesh_.nInternalFaces());
1277  forAll(patches, patchI)
1278  {
1279  const polyPatch& pp = patches[patchI];
1280 
1281  if (pp.coupled())
1282  {
1283  forAll(pp, i)
1284  {
1285  label faceI = pp.start()+i;
1286  neiCellRegion[faceI-mesh_.nInternalFaces()] =
1287  cellRegion[mesh_.faceOwner()[faceI]];
1288  }
1289  }
1290  }
1291  syncTools::swapBoundaryFaceList(mesh_, neiCellRegion, false);
1292 
1293  // Calculate region to zone from cellRegions on either side of coupled
1294  // face.
1295  forAll(patches, patchI)
1296  {
1297  const polyPatch& pp = patches[patchI];
1298 
1299  if (pp.coupled())
1300  {
1301  forAll(pp, i)
1302  {
1303  label faceI = pp.start()+i;
1304 
1305  label surfI = namedSurfaceIndex[faceI];
1306 
1307  if (surfI != -1)
1308  {
1309  bool changedCell = calcRegionToZone
1310  (
1311  surfaceToCellZone[surfI],
1312  cellRegion[mesh_.faceOwner()[faceI]],
1313  neiCellRegion[faceI-mesh_.nInternalFaces()],
1314  regionToCellZone
1315  );
1316 
1317  changed = changed | changedCell;
1318  }
1319  }
1320  }
1321  }
1322 
1323  if (!returnReduce(changed, orOp<bool>()))
1324  {
1325  break;
1326  }
1327  }
1328 
1329 
1330  forAll(regionToCellZone, regionI)
1331  {
1332  label zoneI = regionToCellZone[regionI];
1333 
1334  if (zoneI == -2)
1335  {
1336  FatalErrorIn
1337  (
1338  "meshRefinement::findCellZoneTopo"
1339  "(const point&, const labelList&, const labelList&, labelList&)"
1340  ) << "For region " << regionI << " haven't set cell zone."
1341  << exit(FatalError);
1342  }
1343  }
1344 
1345  if (debug)
1346  {
1347  forAll(regionToCellZone, regionI)
1348  {
1349  Pout<< "Region " << regionI
1350  << " becomes cellZone:" << regionToCellZone[regionI]
1351  << endl;
1352  }
1353  }
1354 
1355  // Rework into cellToZone
1356  forAll(cellToZone, cellI)
1357  {
1358  cellToZone[cellI] = regionToCellZone[cellRegion[cellI]];
1359  }
1360 }
1361 
1362 
1363 // Make namedSurfaceIndex consistent with cellToZone - clear out any blocked
1364 // faces inbetween same cell zone.
1365 void Foam::meshRefinement::makeConsistentFaceIndex
1366 (
1367  const labelList& cellToZone,
1368  labelList& namedSurfaceIndex
1369 ) const
1370 {
1371  const labelList& faceOwner = mesh_.faceOwner();
1372  const labelList& faceNeighbour = mesh_.faceNeighbour();
1373 
1374  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1375  {
1376  label ownZone = cellToZone[faceOwner[faceI]];
1377  label neiZone = cellToZone[faceNeighbour[faceI]];
1378 
1379  if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
1380  {
1381  namedSurfaceIndex[faceI] = -1;
1382  }
1383  else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
1384  {
1385  FatalErrorIn("meshRefinement::zonify()")
1386  << "Different cell zones on either side of face " << faceI
1387  << " at " << mesh_.faceCentres()[faceI]
1388  << " but face not marked with a surface."
1389  << abort(FatalError);
1390  }
1391  }
1392 
1393  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1394 
1395  // Get coupled neighbour cellZone
1396  labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
1397  forAll(patches, patchI)
1398  {
1399  const polyPatch& pp = patches[patchI];
1400 
1401  if (pp.coupled())
1402  {
1403  forAll(pp, i)
1404  {
1405  label faceI = pp.start()+i;
1406  neiCellZone[faceI-mesh_.nInternalFaces()] =
1407  cellToZone[mesh_.faceOwner()[faceI]];
1408  }
1409  }
1410  }
1411  syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
1412 
1413  // Use coupled cellZone to do check
1414  forAll(patches, patchI)
1415  {
1416  const polyPatch& pp = patches[patchI];
1417 
1418  if (pp.coupled())
1419  {
1420  forAll(pp, i)
1421  {
1422  label faceI = pp.start()+i;
1423 
1424  label ownZone = cellToZone[faceOwner[faceI]];
1425  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
1426 
1427  if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
1428  {
1429  namedSurfaceIndex[faceI] = -1;
1430  }
1431  else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
1432  {
1433  FatalErrorIn("meshRefinement::zonify()")
1434  << "Different cell zones on either side of face "
1435  << faceI << " at " << mesh_.faceCentres()[faceI]
1436  << " but face not marked with a surface."
1437  << abort(FatalError);
1438  }
1439  }
1440  }
1441  }
1442 }
1443 
1444 
1445 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1446 
1447 // Split off unreachable areas of mesh.
1450  const bool handleSnapProblems,
1451  const bool removeEdgeConnectedCells,
1452  const scalarField& perpendicularAngle,
1453  const bool mergeFreeStanding,
1454  const dictionary& motionDict,
1455  Time& runTime,
1456  const labelList& globalToPatch,
1457  const point& keepPoint
1458 )
1459 {
1460  // Introduce baffles
1461  // ~~~~~~~~~~~~~~~~~
1462 
1463  // Split the mesh along internal faces wherever there is a pierce between
1464  // two cell centres.
1465 
1466  Info<< "Introducing baffles for "
1467  << returnReduce(countHits(), sumOp<label>())
1468  << " faces that are intersected by the surface." << nl << endl;
1469 
1470  // Swap neighbouring cell centres and cell level
1471  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1472  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
1473  calcNeighbourData(neiLevel, neiCc);
1474 
1475  labelList ownPatch, neiPatch;
1476  getBafflePatches
1477  (
1478  globalToPatch,
1479  neiLevel,
1480  neiCc,
1481 
1482  ownPatch,
1483  neiPatch
1484  );
1485 
1486  if (debug)
1487  {
1488  runTime++;
1489  }
1490 
1491  createBaffles(ownPatch, neiPatch);
1492 
1493  if (debug)
1494  {
1495  // Debug:test all is still synced across proc patches
1496  checkData();
1497  }
1498 
1499  Info<< "Created baffles in = "
1500  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1501 
1502  printMeshInfo(debug, "After introducing baffles");
1503 
1504  if (debug)
1505  {
1506  Pout<< "Writing baffled mesh to time " << timeName()
1507  << endl;
1508  write(debug, runTime.path()/"baffles");
1509  Pout<< "Dumped debug data in = "
1510  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1511  }
1512 
1513 
1514  // Introduce baffles to delete problem cells
1515  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1516 
1517  // Create some additional baffles where we want surface cells removed.
1518 
1519  if (handleSnapProblems)
1520  {
1521  Info<< nl
1522  << "Introducing baffles to block off problem cells" << nl
1523  << "----------------------------------------------" << nl
1524  << endl;
1525 
1526  labelList facePatch
1527  (
1528  markFacesOnProblemCells
1529  (
1530  motionDict,
1531  removeEdgeConnectedCells,
1532  perpendicularAngle,
1533  globalToPatch
1534  )
1535  //markFacesOnProblemCellsGeometric(motionDict)
1536  );
1537  Info<< "Analyzed problem cells in = "
1538  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1539 
1540  if (debug)
1541  {
1542  faceSet problemTopo(mesh_, "problemFacesTopo", 100);
1543 
1544  const labelList facePatchTopo
1545  (
1546  markFacesOnProblemCells
1547  (
1548  motionDict,
1549  removeEdgeConnectedCells,
1550  perpendicularAngle,
1551  globalToPatch
1552  )
1553  );
1554  forAll(facePatchTopo, faceI)
1555  {
1556  if (facePatchTopo[faceI] != -1)
1557  {
1558  problemTopo.insert(faceI);
1559  }
1560  }
1561  Pout<< "Dumping " << problemTopo.size()
1562  << " problem faces to " << problemTopo.objectPath() << endl;
1563  problemTopo.write();
1564  }
1565 
1566  Info<< "Introducing baffles to delete problem cells." << nl << endl;
1567 
1568  if (debug)
1569  {
1570  runTime++;
1571  }
1572 
1573  // Create baffles with same owner and neighbour for now.
1574  createBaffles(facePatch, facePatch);
1575 
1576  if (debug)
1577  {
1578  // Debug:test all is still synced across proc patches
1579  checkData();
1580  }
1581  Info<< "Created baffles in = "
1582  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1583 
1584  printMeshInfo(debug, "After introducing baffles");
1585 
1586  if (debug)
1587  {
1588  Pout<< "Writing extra baffled mesh to time "
1589  << timeName() << endl;
1590  write(debug, runTime.path()/"extraBaffles");
1591  Pout<< "Dumped debug data in = "
1592  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1593  }
1594  }
1595 
1596 
1597  // Select part of mesh
1598  // ~~~~~~~~~~~~~~~~~~~
1599 
1600  Info<< nl
1601  << "Remove unreachable sections of mesh" << nl
1602  << "-----------------------------------" << nl
1603  << endl;
1604 
1605  if (debug)
1606  {
1607  runTime++;
1608  }
1609 
1610  splitMeshRegions(keepPoint);
1611 
1612  if (debug)
1613  {
1614  // Debug:test all is still synced across proc patches
1615  checkData();
1616  }
1617  Info<< "Split mesh in = "
1618  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1619 
1620  printMeshInfo(debug, "After subsetting");
1621 
1622  if (debug)
1623  {
1624  Pout<< "Writing subsetted mesh to time " << timeName()
1625  << endl;
1626  write(debug, runTime.path()/timeName());
1627  Pout<< "Dumped debug data in = "
1628  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1629  }
1630 
1631 
1632  // Merge baffles
1633  // ~~~~~~~~~~~~~
1634 
1635  if (mergeFreeStanding)
1636  {
1637  Info<< nl
1638  << "Merge free-standing baffles" << nl
1639  << "---------------------------" << nl
1640  << endl;
1641 
1642  if (debug)
1643  {
1644  runTime++;
1645  }
1646 
1647  // List of pairs of freestanding baffle faces.
1648  List<labelPair> couples
1649  (
1650  filterDuplicateFaces // filter out freestanding baffles
1651  (
1652  getDuplicateFaces // get all baffles
1653  (
1654  identity(mesh_.nFaces()-mesh_.nInternalFaces())
1655  +mesh_.nInternalFaces()
1656  )
1657  )
1658  );
1659 
1660  label nCouples = couples.size();
1661  reduce(nCouples, sumOp<label>());
1662 
1663  Info<< "Detected free-standing baffles : " << nCouples << endl;
1664 
1665  if (nCouples > 0)
1666  {
1667  // Actually merge baffles. Note: not exactly parallellized. Should
1668  // convert baffle faces into processor faces if they resulted
1669  // from them.
1670  mergeBaffles(couples);
1671 
1672  if (debug)
1673  {
1674  // Debug:test all is still synced across proc patches
1675  checkData();
1676  }
1677  }
1678  Info<< "Merged free-standing baffles in = "
1679  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1680  }
1681 }
1682 
1683 
1684 // Split off (with optional buffer layers) unreachable areas of mesh.
1687  const label nBufferLayers,
1688  const labelList& globalToPatch,
1689  const point& keepPoint
1690 )
1691 {
1692  // Determine patches to put intersections into
1693  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1694 
1695  // Swap neighbouring cell centres and cell level
1696  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1697  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
1698  calcNeighbourData(neiLevel, neiCc);
1699 
1700  labelList ownPatch, neiPatch;
1701  getBafflePatches
1702  (
1703  globalToPatch,
1704  neiLevel,
1705  neiCc,
1706 
1707  ownPatch,
1708  neiPatch
1709  );
1710 
1711  // Analyse regions. Reuse regionsplit
1712  boolList blockedFace(mesh_.nFaces(), false);
1713 
1714  forAll(ownPatch, faceI)
1715  {
1716  if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
1717  {
1718  blockedFace[faceI] = true;
1719  }
1720  }
1721  syncTools::syncFaceList(mesh_, blockedFace, orEqOp<bool>(), false);
1722 
1723  // Set region per cell based on walking
1724  regionSplit cellRegion(mesh_, blockedFace);
1725  blockedFace.clear();
1726 
1727  // Find the region containing the keepPoint
1728  label keepRegionI = -1;
1729 
1730  label cellI = mesh_.findCell(keepPoint);
1731 
1732  if (cellI != -1)
1733  {
1734  keepRegionI = cellRegion[cellI];
1735  }
1736  reduce(keepRegionI, maxOp<label>());
1737 
1738  Info<< "Found point " << keepPoint << " in cell " << cellI
1739  << " in global region " << keepRegionI
1740  << " out of " << cellRegion.nRegions() << " regions." << endl;
1741 
1742  if (keepRegionI == -1)
1743  {
1744  FatalErrorIn
1745  (
1746  "meshRefinement::splitMesh"
1747  "(const label, const labelList&, const point&)"
1748  ) << "Point " << keepPoint
1749  << " is not inside the mesh." << nl
1750  << "Bounding box of the mesh:" << mesh_.globalData().bb()
1751  << exit(FatalError);
1752  }
1753 
1754 
1755  // Walk out nBufferlayers from region split
1756  // (modifies cellRegion, ownPatch)
1757  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1758  // Takes over face patch onto points and then back to faces and cells
1759  // (so cell-face-point walk)
1760 
1761  const labelList& faceOwner = mesh_.faceOwner();
1762  const labelList& faceNeighbour = mesh_.faceNeighbour();
1763 
1764  // Patch for exposed faces for lack of anything sensible.
1765  label defaultPatch = 0;
1766  if (globalToPatch.size())
1767  {
1768  defaultPatch = globalToPatch[0];
1769  }
1770 
1771  for (label i = 0; i < nBufferLayers; i++)
1772  {
1773  // 1. From cells (via faces) to points
1774 
1775  labelList pointBaffle(mesh_.nPoints(), -1);
1776 
1777  forAll(faceNeighbour, faceI)
1778  {
1779  const face& f = mesh_.faces()[faceI];
1780 
1781  label ownRegion = cellRegion[faceOwner[faceI]];
1782  label neiRegion = cellRegion[faceNeighbour[faceI]];
1783 
1784  if (ownRegion == keepRegionI && neiRegion != keepRegionI)
1785  {
1786  // Note max(..) since possibly regionSplit might have split
1787  // off extra unreachable parts of mesh. Note: or can this only
1788  // happen for boundary faces?
1789  forAll(f, fp)
1790  {
1791  pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
1792  }
1793  }
1794  else if (ownRegion != keepRegionI && neiRegion == keepRegionI)
1795  {
1796  label newPatchI = neiPatch[faceI];
1797  if (newPatchI == -1)
1798  {
1799  newPatchI = max(defaultPatch, ownPatch[faceI]);
1800  }
1801  forAll(f, fp)
1802  {
1803  pointBaffle[f[fp]] = newPatchI;
1804  }
1805  }
1806  }
1807  for
1808  (
1809  label faceI = mesh_.nInternalFaces();
1810  faceI < mesh_.nFaces();
1811  faceI++
1812  )
1813  {
1814  const face& f = mesh_.faces()[faceI];
1815 
1816  label ownRegion = cellRegion[faceOwner[faceI]];
1817 
1818  if (ownRegion == keepRegionI)
1819  {
1820  forAll(f, fp)
1821  {
1822  pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
1823  }
1824  }
1825  }
1826 
1827  // Sync
1829  (
1830  mesh_,
1831  pointBaffle,
1832  maxEqOp<label>(),
1833  -1, // null value
1834  false // no separation
1835  );
1836 
1837 
1838  // 2. From points back to faces
1839 
1840  const labelListList& pointFaces = mesh_.pointFaces();
1841 
1842  forAll(pointFaces, pointI)
1843  {
1844  if (pointBaffle[pointI] != -1)
1845  {
1846  const labelList& pFaces = pointFaces[pointI];
1847 
1848  forAll(pFaces, pFaceI)
1849  {
1850  label faceI = pFaces[pFaceI];
1851 
1852  if (ownPatch[faceI] == -1)
1853  {
1854  ownPatch[faceI] = pointBaffle[pointI];
1855  }
1856  }
1857  }
1858  }
1859  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
1860 
1861 
1862  // 3. From faces to cells (cellRegion) and back to faces (ownPatch)
1863 
1864  labelList newOwnPatch(ownPatch);
1865 
1866  forAll(ownPatch, faceI)
1867  {
1868  if (ownPatch[faceI] != -1)
1869  {
1870  label own = faceOwner[faceI];
1871 
1872  if (cellRegion[own] != keepRegionI)
1873  {
1874  cellRegion[own] = keepRegionI;
1875 
1876  const cell& ownFaces = mesh_.cells()[own];
1877  forAll(ownFaces, j)
1878  {
1879  if (ownPatch[ownFaces[j]] == -1)
1880  {
1881  newOwnPatch[ownFaces[j]] = ownPatch[faceI];
1882  }
1883  }
1884  }
1885  if (mesh_.isInternalFace(faceI))
1886  {
1887  label nei = faceNeighbour[faceI];
1888 
1889  if (cellRegion[nei] != keepRegionI)
1890  {
1891  cellRegion[nei] = keepRegionI;
1892 
1893  const cell& neiFaces = mesh_.cells()[nei];
1894  forAll(neiFaces, j)
1895  {
1896  if (ownPatch[neiFaces[j]] == -1)
1897  {
1898  newOwnPatch[neiFaces[j]] = ownPatch[faceI];
1899  }
1900  }
1901  }
1902  }
1903  }
1904  }
1905 
1906  ownPatch.transfer(newOwnPatch);
1907 
1908  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>(), false);
1909  }
1910 
1911 
1912 
1913  // Subset
1914  // ~~~~~~
1915 
1916  // Get cells to remove
1917  DynamicList<label> cellsToRemove(mesh_.nCells());
1918  forAll(cellRegion, cellI)
1919  {
1920  if (cellRegion[cellI] != keepRegionI)
1921  {
1922  cellsToRemove.append(cellI);
1923  }
1924  }
1925  cellsToRemove.shrink();
1926 
1927  label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
1928  reduce(nCellsToKeep, sumOp<label>());
1929 
1930  Info<< "Keeping all cells in region " << keepRegionI
1931  << " containing point " << keepPoint << endl
1932  << "Selected for keeping : " << nCellsToKeep
1933  << " cells." << endl;
1934 
1935 
1936  // Remove cells
1937  removeCells cellRemover(mesh_);
1938 
1939  // Pick up patches for exposed faces
1940  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
1941  labelList exposedPatches(exposedFaces.size());
1942 
1943  forAll(exposedFaces, i)
1944  {
1945  label faceI = exposedFaces[i];
1946 
1947  if (ownPatch[faceI] != -1)
1948  {
1949  exposedPatches[i] = ownPatch[faceI];
1950  }
1951  else
1952  {
1953  WarningIn("meshRefinement::splitMesh(..)")
1954  << "For exposed face " << faceI
1955  << " fc:" << mesh_.faceCentres()[faceI]
1956  << " found no patch." << endl
1957  << " Taking patch " << defaultPatch
1958  << " instead." << endl;
1959  exposedPatches[i] = defaultPatch;
1960  }
1961  }
1962 
1963  return doRemoveCells
1964  (
1965  cellsToRemove,
1966  exposedFaces,
1967  exposedPatches,
1968  cellRemover
1969  );
1970 }
1971 
1972 
1973 // Find boundary points that connect to more than one cell region and
1974 // split them.
1976 {
1977  // Topochange container
1978  polyTopoChange meshMod(mesh_);
1979 
1980 
1981  // Analyse which points need to be duplicated
1983 
1984  label nNonManifPoints = returnReduce
1985  (
1986  regionSide.meshPointMap().size(),
1987  sumOp<label>()
1988  );
1989 
1990  Info<< "dupNonManifoldPoints : Found : " << nNonManifPoints
1991  << " non-manifold points (out of "
1992  << mesh_.globalData().nTotalPoints()
1993  << ')' << endl;
1994 
1995  // Topo change engine
1996  duplicatePoints pointDuplicator(mesh_);
1997 
1998  // Insert changes into meshMod
1999  pointDuplicator.setRefinement(regionSide, meshMod);
2000 
2001  // Change the mesh (no inflation, parallel sync)
2002  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
2003 
2004  // Update fields
2005  mesh_.updateMesh(map);
2006 
2007  // Move mesh if in inflation mode
2008  if (map().hasMotionPoints())
2009  {
2010  mesh_.movePoints(map().preMotionPoints());
2011  }
2012  else
2013  {
2014  // Delete mesh volumes.
2015  mesh_.clearOut();
2016  }
2017 
2018  if (overwrite())
2019  {
2020  mesh_.setInstance(oldInstance());
2021  }
2022 
2023  // Update intersections. Is mapping only (no faces created, positions stay
2024  // same) so no need to recalculate intersections.
2025  updateMesh(map, labelList(0));
2026 
2027  return map;
2028 }
2029 
2030 
2031 // Zoning
2034  const point& keepPoint,
2035  const bool allowFreeStandingZoneFaces
2036 )
2037 {
2038  const wordList& cellZoneNames = surfaces_.cellZoneNames();
2039  const wordList& faceZoneNames = surfaces_.faceZoneNames();
2040 
2041  labelList namedSurfaces(surfaces_.getNamedSurfaces());
2042 
2043  boolList isNamedSurface(cellZoneNames.size(), false);
2044 
2045  forAll(namedSurfaces, i)
2046  {
2047  label surfI = namedSurfaces[i];
2048 
2049  isNamedSurface[surfI] = true;
2050 
2051  Info<< "Surface : " << surfaces_.names()[surfI] << nl
2052  << " faceZone : " << faceZoneNames[surfI] << nl
2053  << " cellZone : " << cellZoneNames[surfI] << endl;
2054  }
2055 
2056 
2057  // Add zones to mesh
2058 
2059  labelList surfaceToFaceZone(faceZoneNames.size(), -1);
2060  {
2061  faceZoneMesh& faceZones = mesh_.faceZones();
2062 
2063  forAll(namedSurfaces, i)
2064  {
2065  label surfI = namedSurfaces[i];
2066 
2067  label zoneI = faceZones.findZoneID(faceZoneNames[surfI]);
2068 
2069  if (zoneI == -1)
2070  {
2071  zoneI = faceZones.size();
2072  faceZones.setSize(zoneI+1);
2073  faceZones.set
2074  (
2075  zoneI,
2076  new faceZone
2077  (
2078  faceZoneNames[surfI], //name
2079  labelList(0), //addressing
2080  boolList(0), //flipmap
2081  zoneI, //index
2082  faceZones //faceZoneMesh
2083  )
2084  );
2085  }
2086 
2087  if (debug)
2088  {
2089  Pout<< "Faces on " << surfaces_.names()[surfI]
2090  << " will go into faceZone " << zoneI << endl;
2091  }
2092  surfaceToFaceZone[surfI] = zoneI;
2093  }
2094 
2095  // Check they are synced
2096  List<wordList> allFaceZones(Pstream::nProcs());
2097  allFaceZones[Pstream::myProcNo()] = faceZones.names();
2098  Pstream::gatherList(allFaceZones);
2099  Pstream::scatterList(allFaceZones);
2100 
2101  for (label procI = 1; procI < allFaceZones.size(); procI++)
2102  {
2103  if (allFaceZones[procI] != allFaceZones[0])
2104  {
2105  FatalErrorIn
2106  (
2107  "meshRefinement::zonify"
2108  "(const label, const point&)"
2109  ) << "Zones not synchronised among processors." << nl
2110  << " Processor0 has faceZones:" << allFaceZones[0]
2111  << " , processor" << procI
2112  << " has faceZones:" << allFaceZones[procI]
2113  << exit(FatalError);
2114  }
2115  }
2116  }
2117 
2118  labelList surfaceToCellZone(cellZoneNames.size(), -1);
2119  {
2120  cellZoneMesh& cellZones = mesh_.cellZones();
2121 
2122  forAll(namedSurfaces, i)
2123  {
2124  label surfI = namedSurfaces[i];
2125 
2126  label zoneI = cellZones.findZoneID(cellZoneNames[surfI]);
2127 
2128  if (zoneI == -1)
2129  {
2130  zoneI = cellZones.size();
2131  cellZones.setSize(zoneI+1);
2132  cellZones.set
2133  (
2134  zoneI,
2135  new cellZone
2136  (
2137  cellZoneNames[surfI], //name
2138  labelList(0), //addressing
2139  zoneI, //index
2140  cellZones //cellZoneMesh
2141  )
2142  );
2143  }
2144 
2145  if (debug)
2146  {
2147  Pout<< "Cells inside " << surfaces_.names()[surfI]
2148  << " will go into cellZone " << zoneI << endl;
2149  }
2150  surfaceToCellZone[surfI] = zoneI;
2151  }
2152 
2153  // Check they are synced
2154  List<wordList> allCellZones(Pstream::nProcs());
2155  allCellZones[Pstream::myProcNo()] = cellZones.names();
2156  Pstream::gatherList(allCellZones);
2157  Pstream::scatterList(allCellZones);
2158 
2159  for (label procI = 1; procI < allCellZones.size(); procI++)
2160  {
2161  if (allCellZones[procI] != allCellZones[0])
2162  {
2163  FatalErrorIn
2164  (
2165  "meshRefinement::zonify"
2166  "(const label, const point&)"
2167  ) << "Zones not synchronised among processors." << nl
2168  << " Processor0 has cellZones:" << allCellZones[0]
2169  << " , processor" << procI
2170  << " has cellZones:" << allCellZones[procI]
2171  << exit(FatalError);
2172  }
2173  }
2174  }
2175 
2176 
2177 
2178  const pointField& cellCentres = mesh_.cellCentres();
2179  const labelList& faceOwner = mesh_.faceOwner();
2180  const labelList& faceNeighbour = mesh_.faceNeighbour();
2181  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2182 
2183 
2184  // Mark faces intersecting zoned surfaces
2185  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2186 
2187 
2188  // Like surfaceIndex_ but only for named surfaces.
2189  labelList namedSurfaceIndex(mesh_.nFaces(), -1);
2190 
2191  {
2192  // Statistics: number of faces per faceZone
2193  labelList nSurfFaces(faceZoneNames.size(), 0);
2194 
2195  // Swap neighbouring cell centres and cell level
2196  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2197  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2198  calcNeighbourData(neiLevel, neiCc);
2199 
2200  // Note: for all internal faces? internal + coupled?
2201  // Since zonify is run after baffling the surfaceIndex_ on baffles is
2202  // not synchronised across both baffle faces. Fortunately we don't
2203  // do zonify baffle faces anyway (they are normal boundary faces).
2204 
2205  // Collect candidate faces
2206  // ~~~~~~~~~~~~~~~~~~~~~~~
2207 
2208  labelList testFaces(intersectedFaces());
2209 
2210  // Collect segments
2211  // ~~~~~~~~~~~~~~~~
2212 
2213  pointField start(testFaces.size());
2214  pointField end(testFaces.size());
2215 
2216  forAll(testFaces, i)
2217  {
2218  label faceI = testFaces[i];
2219 
2220  if (mesh_.isInternalFace(faceI))
2221  {
2222  start[i] = cellCentres[faceOwner[faceI]];
2223  end[i] = cellCentres[faceNeighbour[faceI]];
2224  }
2225  else
2226  {
2227  start[i] = cellCentres[faceOwner[faceI]];
2228  end[i] = neiCc[faceI-mesh_.nInternalFaces()];
2229  }
2230  }
2231 
2232 
2233  // Do test for intersections
2234  // ~~~~~~~~~~~~~~~~~~~~~~~~~
2235  // Note that we intersect all intersected faces again. Could reuse
2236  // the information already in surfaceIndex_.
2237 
2238  labelList surface1;
2239  labelList surface2;
2240  {
2241  List<pointIndexHit> hit1;
2242  labelList region1;
2243  List<pointIndexHit> hit2;
2244  labelList region2;
2245  surfaces_.findNearestIntersection
2246  (
2247  namedSurfaces,
2248  start,
2249  end,
2250 
2251  surface1,
2252  hit1,
2253  region1,
2254  surface2,
2255  hit2,
2256  region2
2257  );
2258  }
2259 
2260  forAll(testFaces, i)
2261  {
2262  label faceI = testFaces[i];
2263 
2264  if (surface1[i] != -1)
2265  {
2266  // If both hit should probably choose nearest. For later.
2267  namedSurfaceIndex[faceI] = surface1[i];
2268  nSurfFaces[surface1[i]]++;
2269  }
2270  else if (surface2[i] != -1)
2271  {
2272  namedSurfaceIndex[faceI] = surface2[i];
2273  nSurfFaces[surface2[i]]++;
2274  }
2275  }
2276 
2277 
2278  // surfaceIndex migh have different surfaces on both sides if
2279  // there happen to be a (obviously thin) surface with different
2280  // regions between the cell centres. If one is on a named surface
2281  // and the other is not this might give problems so sync.
2283  (
2284  mesh_,
2285  namedSurfaceIndex,
2286  maxEqOp<label>(),
2287  false
2288  );
2289 
2290  // Print a bit
2291  if (debug)
2292  {
2293  forAll(nSurfFaces, surfI)
2294  {
2295  Pout<< "Surface:"
2296  << surfaces_.names()[surfI]
2297  << " nZoneFaces:" << nSurfFaces[surfI] << nl;
2298  }
2299  Pout<< endl;
2300  }
2301  }
2302 
2303 
2304  // Put the cells into the correct zone
2305  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2306 
2307  // Closed surfaces with cellZone specified.
2308  labelList closedNamedSurfaces(surfaces_.getClosedNamedSurfaces());
2309 
2310  // Zone per cell:
2311  // -2 : unset
2312  // -1 : not in any zone
2313  // >=0: zoneID
2314  labelList cellToZone(mesh_.nCells(), -2);
2315 
2316 
2317  // Set using geometric test
2318  // ~~~~~~~~~~~~~~~~~~~~~~~~
2319 
2320  if (closedNamedSurfaces.size())
2321  {
2322  findCellZoneGeometric
2323  (
2324  closedNamedSurfaces, // indices of closed surfaces
2325  namedSurfaceIndex, // per face index of named surface
2326  surfaceToCellZone, // cell zone index per surface
2327  cellToZone
2328  );
2329  }
2330 
2331  // Set using walking
2332  // ~~~~~~~~~~~~~~~~~
2333 
2334  //if (returnReduce(nSet, sumOp<label>()) < mesh_.globalData().nTotalCells())
2335  {
2336  // Topological walk
2337  findCellZoneTopo
2338  (
2339  keepPoint,
2340  namedSurfaceIndex,
2341  surfaceToCellZone,
2342  cellToZone
2343  );
2344  }
2345 
2346 
2347  // Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
2348  if (!allowFreeStandingZoneFaces)
2349  {
2350  makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
2351  }
2352 
2353  // Topochange container
2354  polyTopoChange meshMod(mesh_);
2355 
2356 
2357  // Put the faces into the correct zone
2358  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2359 
2360  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
2361  {
2362  label surfI = namedSurfaceIndex[faceI];
2363 
2364  if (surfI != -1)
2365  {
2366  // Orient face zone to have slave cells in max cell zone.
2367  label ownZone = cellToZone[faceOwner[faceI]];
2368  label neiZone = cellToZone[faceNeighbour[faceI]];
2369 
2370  bool flip;
2371  if (ownZone == max(ownZone, neiZone))
2372  {
2373  flip = false;
2374  }
2375  else
2376  {
2377  flip = true;
2378  }
2379 
2380  meshMod.setAction
2381  (
2383  (
2384  mesh_.faces()[faceI], // modified face
2385  faceI, // label of face
2386  faceOwner[faceI], // owner
2387  faceNeighbour[faceI], // neighbour
2388  false, // face flip
2389  -1, // patch for face
2390  false, // remove from zone
2391  surfaceToFaceZone[surfI], // zone for face
2392  flip // face flip in zone
2393  )
2394  );
2395  }
2396  }
2397 
2398  // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
2399  labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
2400  forAll(patches, patchI)
2401  {
2402  const polyPatch& pp = patches[patchI];
2403 
2404  if (pp.coupled())
2405  {
2406  forAll(pp, i)
2407  {
2408  label faceI = pp.start()+i;
2409  neiCellZone[faceI-mesh_.nInternalFaces()] =
2410  cellToZone[mesh_.faceOwner()[faceI]];
2411  }
2412  }
2413  }
2414  syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
2415 
2416  // Get per face whether is it master (of a coupled set of faces)
2417  PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
2418 
2419  // Set owner as no-flip
2420  forAll(patches, patchI)
2421  {
2422  const polyPatch& pp = patches[patchI];
2423 
2424  label faceI = pp.start();
2425 
2426  forAll(pp, i)
2427  {
2428  label surfI = namedSurfaceIndex[faceI];
2429 
2430  if (surfI != -1)
2431  {
2432  label ownZone = cellToZone[faceOwner[faceI]];
2433  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
2434 
2435  bool flip;
2436 
2437  label maxZone = max(ownZone, neiZone);
2438 
2439  if (maxZone == -1)
2440  {
2441  flip = false;
2442  }
2443  else if (ownZone == neiZone)
2444  {
2445  // Free-standing zone face or coupled boundary. Keep master
2446  // face unflipped.
2447  flip = !isMasterFace[faceI];
2448  }
2449  else if (neiZone == maxZone)
2450  {
2451  flip = true;
2452  }
2453  else
2454  {
2455  flip = false;
2456  }
2457 
2458  meshMod.setAction
2459  (
2461  (
2462  mesh_.faces()[faceI], // modified face
2463  faceI, // label of face
2464  faceOwner[faceI], // owner
2465  -1, // neighbour
2466  false, // face flip
2467  patchI, // patch for face
2468  false, // remove from zone
2469  surfaceToFaceZone[surfI], // zone for face
2470  flip // face flip in zone
2471  )
2472  );
2473  }
2474  faceI++;
2475  }
2476  }
2477 
2478 
2479  // Put the cells into the correct zone
2480  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2481 
2482  forAll(cellToZone, cellI)
2483  {
2484  label zoneI = cellToZone[cellI];
2485 
2486  if (zoneI >= 0)
2487  {
2488  meshMod.setAction
2489  (
2491  (
2492  cellI,
2493  false, // removeFromZone
2494  zoneI
2495  )
2496  );
2497  }
2498  }
2499 
2500  // Change the mesh (no inflation, parallel sync)
2501  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
2502 
2503  // Update fields
2504  mesh_.updateMesh(map);
2505 
2506  // Move mesh if in inflation mode
2507  if (map().hasMotionPoints())
2508  {
2509  mesh_.movePoints(map().preMotionPoints());
2510  }
2511  else
2512  {
2513  // Delete mesh volumes.
2514  mesh_.clearOut();
2515  }
2516 
2517  if (overwrite())
2518  {
2519  mesh_.setInstance(oldInstance());
2520  }
2521 
2522  // None of the faces has changed, only the zones. Still...
2523  updateMesh(map, labelList());
2524 
2525  return map;
2526 }
2527 
2528 
2529 // ************************ vim: set sw=4 sts=4 et: ************************ //