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) 2011-2015 OpenFOAM Foundation
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 "fvMesh.H"
28 #include "syncTools.H"
29 #include "Time.H"
30 #include "refinementSurfaces.H"
31 #include "pointSet.H"
32 #include "faceSet.H"
33 #include "indirectPrimitivePatch.H"
34 #include "polyTopoChange.H"
35 #include "meshTools.H"
36 #include "polyModifyFace.H"
37 #include "polyModifyCell.H"
38 #include "polyAddFace.H"
39 #include "polyRemoveFace.H"
40 #include "polyAddPoint.H"
41 #include "localPointRegion.H"
42 #include "duplicatePoints.H"
43 #include "OFstream.H"
44 #include "regionSplit.H"
45 #include "removeCells.H"
46 #include "unitConversion.H"
47 #include "OBJstream.H"
48 #include "patchFaceOrientation.H"
49 #include "PatchEdgeFaceWave.H"
50 #include "patchEdgeFaceRegion.H"
51 
52 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
53 
54 Foam::label Foam::meshRefinement::createBaffle
55 (
56  const label faceI,
57  const label ownPatch,
58  const label neiPatch,
59  polyTopoChange& meshMod
60 ) const
61 {
62  const face& f = mesh_.faces()[faceI];
63  label zoneID = mesh_.faceZones().whichZone(faceI);
64  bool zoneFlip = false;
65 
66  if (zoneID >= 0)
67  {
68  const faceZone& fZone = mesh_.faceZones()[zoneID];
69  zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
70  }
71 
72  meshMod.setAction
73  (
74  polyModifyFace
75  (
76  f, // modified face
77  faceI, // label of face
78  mesh_.faceOwner()[faceI], // owner
79  -1, // neighbour
80  false, // face flip
81  ownPatch, // patch for face
82  false, // remove from zone
83  zoneID, // zone for face
84  zoneFlip // face flip in zone
85  )
86  );
87 
88 
89  label dupFaceI = -1;
90 
91  if (mesh_.isInternalFace(faceI))
92  {
93  if (neiPatch == -1)
94  {
96  << "No neighbour patch for internal face " << faceI
97  << " fc:" << mesh_.faceCentres()[faceI]
98  << " ownPatch:" << ownPatch << abort(FatalError);
99  }
100 
101  bool reverseFlip = false;
102  if (zoneID >= 0)
103  {
104  reverseFlip = !zoneFlip;
105  }
106 
107  dupFaceI = meshMod.setAction
108  (
109  polyAddFace
110  (
111  f.reverseFace(), // modified face
112  mesh_.faceNeighbour()[faceI],// owner
113  -1, // neighbour
114  -1, // masterPointID
115  -1, // masterEdgeID
116  faceI, // masterFaceID,
117  true, // face flip
118  neiPatch, // patch for face
119  zoneID, // zone for face
120  reverseFlip // face flip in zone
121  )
122  );
123  }
124  return dupFaceI;
125 }
126 
127 
128 void Foam::meshRefinement::getBafflePatches
129 (
130  const labelList& globalToMasterPatch,
131  const labelList& neiLevel,
132  const pointField& neiCc,
133  labelList& ownPatch,
134  labelList& neiPatch
135 ) const
136 {
137  autoPtr<OFstream> str;
138  label vertI = 0;
139  if (debug&OBJINTERSECTIONS)
140  {
141  mkDir(mesh_.time().path()/timeName());
142  str.reset
143  (
144  new OFstream
145  (
146  mesh_.time().path()/timeName()/"intersections.obj"
147  )
148  );
149 
150  Pout<< "getBafflePatches : Writing surface intersections to file "
151  << str().name() << nl << endl;
152  }
153 
154  const pointField& cellCentres = mesh_.cellCentres();
155 
156  // Surfaces that need to be baffled
157  const labelList surfacesToBaffle
158  (
160  );
161 
162  ownPatch.setSize(mesh_.nFaces());
163  ownPatch = -1;
164  neiPatch.setSize(mesh_.nFaces());
165  neiPatch = -1;
166 
167 
168  // Collect candidate faces
169  // ~~~~~~~~~~~~~~~~~~~~~~~
170 
171  labelList testFaces(intersectedFaces());
172 
173  // Collect segments
174  // ~~~~~~~~~~~~~~~~
175 
176  pointField start(testFaces.size());
177  pointField end(testFaces.size());
178 
179  forAll(testFaces, i)
180  {
181  label faceI = testFaces[i];
182 
183  label own = mesh_.faceOwner()[faceI];
184 
185  if (mesh_.isInternalFace(faceI))
186  {
187  start[i] = cellCentres[own];
188  end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
189  }
190  else
191  {
192  start[i] = cellCentres[own];
193  end[i] = neiCc[faceI-mesh_.nInternalFaces()];
194  }
195  }
196 
197  // Extend segments a bit
198  {
199  const vectorField smallVec(ROOTSMALL*(end-start));
200  start -= smallVec;
201  end += smallVec;
202  }
203 
204 
205  // Do test for intersections
206  // ~~~~~~~~~~~~~~~~~~~~~~~~~
207  labelList surface1;
208  List<pointIndexHit> hit1;
209  labelList region1;
210  labelList surface2;
211  List<pointIndexHit> hit2;
212  labelList region2;
213  surfaces_.findNearestIntersection
214  (
215  surfacesToBaffle,
216  start,
217  end,
218 
219  surface1,
220  hit1,
221  region1,
222 
223  surface2,
224  hit2,
225  region2
226  );
227 
228  forAll(testFaces, i)
229  {
230  label faceI = testFaces[i];
231 
232  if (hit1[i].hit() && hit2[i].hit())
233  {
234  if (str.valid())
235  {
236  meshTools::writeOBJ(str(), start[i]);
237  vertI++;
238  meshTools::writeOBJ(str(), hit1[i].rawPoint());
239  vertI++;
240  meshTools::writeOBJ(str(), hit2[i].rawPoint());
241  vertI++;
242  meshTools::writeOBJ(str(), end[i]);
243  vertI++;
244  str()<< "l " << vertI-3 << ' ' << vertI-2 << nl;
245  str()<< "l " << vertI-2 << ' ' << vertI-1 << nl;
246  str()<< "l " << vertI-1 << ' ' << vertI << nl;
247  }
248 
249  // Pick up the patches
250  ownPatch[faceI] = globalToMasterPatch
251  [
252  surfaces_.globalRegion(surface1[i], region1[i])
253  ];
254  neiPatch[faceI] = globalToMasterPatch
255  [
256  surfaces_.globalRegion(surface2[i], region2[i])
257  ];
258 
259  if (ownPatch[faceI] == -1 || neiPatch[faceI] == -1)
260  {
262  << "problem." << abort(FatalError);
263  }
264  }
265  }
266 
267  // No need to parallel sync since intersection data (surfaceIndex_ etc.)
268  // already guaranteed to be synced...
269  // However:
270  // - owncc and neicc are reversed on different procs so might pick
271  // up different regions reversed? No problem. Neighbour on one processor
272  // might not be owner on the other processor but the neighbour is
273  // not used when creating baffles from proc faces.
274  // - tolerances issues occasionally crop up.
275  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
276  syncTools::syncFaceList(mesh_, neiPatch, maxEqOp<label>());
277 }
278 
279 
280 Foam::Map<Foam::labelPair> Foam::meshRefinement::getZoneBafflePatches
281 (
282  const bool allowBoundary,
283  const labelList& globalToMasterPatch,
284  const labelList& globalToSlavePatch
285 ) const
286 {
287  Map<labelPair> bafflePatch(mesh_.nFaces()/1000);
288 
289  const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
290  const faceZoneMesh& fZones = mesh_.faceZones();
291 
292  forAll(surfZones, surfI)
293  {
294  const word& faceZoneName = surfZones[surfI].faceZoneName();
295 
296  if (faceZoneName.size())
297  {
298  // Get zone
299  label zoneI = fZones.findZoneID(faceZoneName);
300 
301  const faceZone& fZone = fZones[zoneI];
302 
303  // Get patch allocated for zone
304  label globalRegionI = surfaces_.globalRegion(surfI, 0);
305  labelPair zPatches
306  (
307  globalToMasterPatch[globalRegionI],
308  globalToSlavePatch[globalRegionI]
309  );
310 
311  Info<< "For zone " << fZone.name() << " found patches "
312  << mesh_.boundaryMesh()[zPatches[0]].name() << " and "
313  << mesh_.boundaryMesh()[zPatches[1]].name()
314  << endl;
315 
316  forAll(fZone, i)
317  {
318  label faceI = fZone[i];
319 
320  if (allowBoundary || mesh_.isInternalFace(faceI))
321  {
322  labelPair patches = zPatches;
323  if (fZone.flipMap()[i])
324  {
325  patches = reverse(patches);
326  }
327 
328  if (!bafflePatch.insert(faceI, patches))
329  {
331  << "Face " << faceI
332  << " fc:" << mesh_.faceCentres()[faceI]
333  << " in zone " << fZone.name()
334  << " is in multiple zones!"
335  << abort(FatalError);
336  }
337  }
338  }
339  }
340  }
341  return bafflePatch;
342 }
343 
344 
346 (
347  const labelList& ownPatch,
348  const labelList& neiPatch
349 )
350 {
351  if
352  (
353  ownPatch.size() != mesh_.nFaces()
354  || neiPatch.size() != mesh_.nFaces()
355  )
356  {
358  << "Illegal size :"
359  << " ownPatch:" << ownPatch.size()
360  << " neiPatch:" << neiPatch.size()
361  << ". Should be number of faces:" << mesh_.nFaces()
362  << abort(FatalError);
363  }
364 
365  if (debug)
366  {
367  labelList syncedOwnPatch(ownPatch);
368  syncTools::syncFaceList(mesh_, syncedOwnPatch, maxEqOp<label>());
369  labelList syncedNeiPatch(neiPatch);
370  syncTools::syncFaceList(mesh_, syncedNeiPatch, maxEqOp<label>());
371 
372  forAll(syncedOwnPatch, faceI)
373  {
374  if
375  (
376  (ownPatch[faceI] == -1 && syncedOwnPatch[faceI] != -1)
377  || (neiPatch[faceI] == -1 && syncedNeiPatch[faceI] != -1)
378  )
379  {
381  << "Non synchronised at face:" << faceI
382  << " on patch:" << mesh_.boundaryMesh().whichPatch(faceI)
383  << " fc:" << mesh_.faceCentres()[faceI] << endl
384  << "ownPatch:" << ownPatch[faceI]
385  << " syncedOwnPatch:" << syncedOwnPatch[faceI]
386  << " neiPatch:" << neiPatch[faceI]
387  << " syncedNeiPatch:" << syncedNeiPatch[faceI]
388  << abort(FatalError);
389  }
390  }
391  }
392 
393  // Topochange container
394  polyTopoChange meshMod(mesh_);
395 
396  label nBaffles = 0;
397 
398  forAll(ownPatch, faceI)
399  {
400  if (ownPatch[faceI] != -1)
401  {
402  // Create baffle or repatch face. Return label of inserted baffle
403  // face.
404  createBaffle
405  (
406  faceI,
407  ownPatch[faceI], // owner side patch
408  neiPatch[faceI], // neighbour side patch
409  meshMod
410  );
411  nBaffles++;
412  }
413  }
414 
415  // Change the mesh (no inflation, parallel sync)
416  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
417 
418  // Update fields
419  mesh_.updateMesh(map);
420 
421  // Move mesh if in inflation mode
422  if (map().hasMotionPoints())
423  {
424  mesh_.movePoints(map().preMotionPoints());
425  }
426  else
427  {
428  // Delete mesh volumes.
429  mesh_.clearOut();
430  }
431 
432 
433  // Reset the instance for if in overwrite mode
434  mesh_.setInstance(timeName());
435 
436  //- Redo the intersections on the newly create baffle faces. Note that
437  // this changes also the cell centre positions.
438  faceSet baffledFacesSet(mesh_, "baffledFacesSet", 2*nBaffles);
439 
440  const labelList& reverseFaceMap = map().reverseFaceMap();
441  const labelList& faceMap = map().faceMap();
442 
443  // Pick up owner side of baffle
444  forAll(ownPatch, oldFaceI)
445  {
446  label faceI = reverseFaceMap[oldFaceI];
447 
448  if (ownPatch[oldFaceI] != -1 && faceI >= 0)
449  {
450  const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
451 
452  forAll(ownFaces, i)
453  {
454  baffledFacesSet.insert(ownFaces[i]);
455  }
456  }
457  }
458  // Pick up neighbour side of baffle (added faces)
459  forAll(faceMap, faceI)
460  {
461  label oldFaceI = faceMap[faceI];
462 
463  if (oldFaceI >= 0 && reverseFaceMap[oldFaceI] != faceI)
464  {
465  const cell& ownFaces = mesh_.cells()[mesh_.faceOwner()[faceI]];
466 
467  forAll(ownFaces, i)
468  {
469  baffledFacesSet.insert(ownFaces[i]);
470  }
471  }
472  }
473  baffledFacesSet.sync(mesh_);
474 
475  updateMesh(map, baffledFacesSet.toc());
476 
477  return map;
478 }
479 
480 
482 {
483  const faceZoneMesh& fZones = mesh_.faceZones();
484 
485  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
486 
487  forAll(pbm, patchI)
488  {
489  const polyPatch& pp = pbm[patchI];
490 
491  if (isA<processorPolyPatch>(pp))
492  {
493  forAll(pp, i)
494  {
495  label faceI = pp.start()+i;
496  label zoneI = fZones.whichZone(faceI);
497 
498  if (zoneI != -1)
499  {
501  << "face:" << faceI << " on patch " << pp.name()
502  << " is in zone " << fZones[zoneI].name()
503  << exit(FatalError);
504  }
505  }
506  }
507  }
508 }
509 
510 
512 (
513  const labelList& globalToMasterPatch,
514  const labelList& globalToSlavePatch,
515  List<labelPair>& baffles
516 )
517 {
518  const labelList zonedSurfaces
519  (
521  );
522 
524 
525  // No need to sync; all processors will have all same zonedSurfaces.
526  if (zonedSurfaces.size())
527  {
528  // Split internal faces on interface surfaces
529  Info<< "Converting zoned faces into baffles ..." << endl;
530 
531  // Get faces (internal only) to be baffled. Map from face to patch
532  // label.
533  Map<labelPair> faceToPatch
534  (
535  getZoneBafflePatches
536  (
537  false,
538  globalToMasterPatch,
539  globalToSlavePatch
540  )
541  );
542 
543  label nZoneFaces = returnReduce(faceToPatch.size(), sumOp<label>());
544  if (nZoneFaces > 0)
545  {
546  // Convert into labelLists
547  labelList ownPatch(mesh_.nFaces(), -1);
548  labelList neiPatch(mesh_.nFaces(), -1);
549  forAllConstIter(Map<labelPair>, faceToPatch, iter)
550  {
551  ownPatch[iter.key()] = iter().first();
552  neiPatch[iter.key()] = iter().second();
553  }
554 
555  // Create baffles. both sides same patch.
556  map = createBaffles(ownPatch, neiPatch);
557 
558  // Get pairs of faces created.
559  // Just loop over faceMap and store baffle if we encounter a slave
560  // face.
561 
562  baffles.setSize(faceToPatch.size());
563  label baffleI = 0;
564 
565  const labelList& faceMap = map().faceMap();
566  const labelList& reverseFaceMap = map().reverseFaceMap();
567 
568  forAll(faceMap, faceI)
569  {
570  label oldFaceI = faceMap[faceI];
571 
572  // Does face originate from face-to-patch
573  Map<labelPair>::const_iterator iter = faceToPatch.find
574  (
575  oldFaceI
576  );
577 
578  if (iter != faceToPatch.end())
579  {
580  label masterFaceI = reverseFaceMap[oldFaceI];
581  if (faceI != masterFaceI)
582  {
583  baffles[baffleI++] = labelPair(masterFaceI, faceI);
584  }
585  }
586  }
587 
588  if (baffleI != faceToPatch.size())
589  {
591  << "Had " << faceToPatch.size() << " patches to create "
592  << " but encountered " << baffleI
593  << " slave faces originating from patcheable faces."
594  << abort(FatalError);
595  }
596 
597  if (debug&MESH)
598  {
599  const_cast<Time&>(mesh_.time())++;
600  Pout<< "Writing zone-baffled mesh to time " << timeName()
601  << endl;
602  write
603  (
604  debugType(debug),
606  mesh_.time().path()/"baffles"
607  );
608  }
609  }
610  Info<< "Created " << nZoneFaces << " baffles in = "
611  << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
612  }
613  return map;
614 }
615 
616 
617 Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
618 (
619  const List<labelPair>& couples,
620  const scalar planarAngle
621 ) const
622 {
623  // Done by counting the number of baffles faces per mesh edge. If edge
624  // has 2 boundary faces and both are baffle faces it is the edge of a baffle
625  // region.
626 
627  // All duplicate faces on edge of the patch are to be merged.
628  // So we count for all edges of duplicate faces how many duplicate
629  // faces use them.
630  labelList nBafflesPerEdge(mesh_.nEdges(), 0);
631 
632 
633  // This algorithm is quite tricky. We don't want to use edgeFaces and
634  // also want it to run in parallel so it is now an algorithm over
635  // all (boundary) faces instead.
636  // We want to pick up any edges that are only used by the baffle
637  // or internal faces but not by any other boundary faces. So
638  // - increment count on an edge by 1 if it is used by any (uncoupled)
639  // boundary face.
640  // - increment count on an edge by 1000000 if it is used by a baffle face
641  // - sum in parallel
642  //
643  // So now any edge that is used by baffle faces only will have the
644  // value 2*1000000+2*1.
645 
646 
647  const label baffleValue = 1000000;
648 
649 
650  // Count number of boundary faces per edge
651  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
652 
653  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
654 
655  forAll(patches, patchI)
656  {
657  const polyPatch& pp = patches[patchI];
658 
659  // Count number of boundary faces. Discard coupled boundary faces.
660  if (!pp.coupled())
661  {
662  label faceI = pp.start();
663 
664  forAll(pp, i)
665  {
666  const labelList& fEdges = mesh_.faceEdges(faceI);
667 
668  forAll(fEdges, fEdgeI)
669  {
670  nBafflesPerEdge[fEdges[fEdgeI]]++;
671  }
672  faceI++;
673  }
674  }
675  }
676 
677 
678  DynamicList<label> fe0;
679  DynamicList<label> fe1;
680 
681 
682  // Count number of duplicate boundary faces per edge
683  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
684 
685  forAll(couples, i)
686  {
687  {
688  label f0 = couples[i].first();
689  const labelList& fEdges0 = mesh_.faceEdges(f0, fe0);
690  forAll(fEdges0, fEdgeI)
691  {
692  nBafflesPerEdge[fEdges0[fEdgeI]] += baffleValue;
693  }
694  }
695 
696  {
697  label f1 = couples[i].second();
698  const labelList& fEdges1 = mesh_.faceEdges(f1, fe1);
699  forAll(fEdges1, fEdgeI)
700  {
701  nBafflesPerEdge[fEdges1[fEdgeI]] += baffleValue;
702  }
703  }
704  }
705 
706  // Add nBaffles on shared edges
708  (
709  mesh_,
710  nBafflesPerEdge,
711  plusEqOp<label>(), // in-place add
712  label(0) // initial value
713  );
714 
715 
716  // Baffles which are not next to other boundaries and baffles will have
717  // nBafflesPerEdge value 2*baffleValue+2*1 (from 2 boundary faces which
718  // are both baffle faces)
719 
720  List<labelPair> filteredCouples(couples.size());
721  label filterI = 0;
722 
723  forAll(couples, i)
724  {
725  const labelPair& couple = couples[i];
726 
727  if
728  (
729  patches.whichPatch(couple.first())
730  == patches.whichPatch(couple.second())
731  )
732  {
733  const labelList& fEdges = mesh_.faceEdges(couple.first());
734 
735  forAll(fEdges, fEdgeI)
736  {
737  label edgeI = fEdges[fEdgeI];
738 
739  if (nBafflesPerEdge[edgeI] == 2*baffleValue+2*1)
740  {
741  filteredCouples[filterI++] = couple;
742  break;
743  }
744  }
745  }
746  }
747  filteredCouples.setSize(filterI);
748 
749 
750  label nFiltered = returnReduce(filteredCouples.size(), sumOp<label>());
751 
752  Info<< "freeStandingBaffles : detected "
753  << nFiltered
754  << " free-standing baffles out of "
755  << returnReduce(couples.size(), sumOp<label>())
756  << nl << endl;
757 
758 
759  if (nFiltered > 0)
760  {
761  // Collect segments
762  // ~~~~~~~~~~~~~~~~
763 
764  pointField start(filteredCouples.size());
765  pointField end(filteredCouples.size());
766 
767  const pointField& cellCentres = mesh_.cellCentres();
768 
769  forAll(filteredCouples, i)
770  {
771  const labelPair& couple = filteredCouples[i];
772  start[i] = cellCentres[mesh_.faceOwner()[couple.first()]];
773  end[i] = cellCentres[mesh_.faceOwner()[couple.second()]];
774  }
775 
776  // Extend segments a bit
777  {
778  const vectorField smallVec(ROOTSMALL*(end-start));
779  start -= smallVec;
780  end += smallVec;
781  }
782 
783 
784  // Do test for intersections
785  // ~~~~~~~~~~~~~~~~~~~~~~~~~
786  labelList surface1;
787  List<pointIndexHit> hit1;
788  labelList region1;
789  vectorField normal1;
790 
791  labelList surface2;
792  List<pointIndexHit> hit2;
793  labelList region2;
794  vectorField normal2;
795 
796  surfaces_.findNearestIntersection
797  (
798  identity(surfaces_.surfaces().size()),
799  start,
800  end,
801 
802  surface1,
803  hit1,
804  region1,
805  normal1,
806 
807  surface2,
808  hit2,
809  region2,
810  normal2
811  );
812 
813  //mkDir(mesh_.time().path()/timeName());
814  //OBJstream str
815  //(
816  // mesh_.time().path()/timeName()/"flatBaffles.obj"
817  //);
818 
819  const scalar planarAngleCos = Foam::cos(degToRad(planarAngle));
820 
821  label filterI = 0;
822  forAll(filteredCouples, i)
823  {
824  const labelPair& couple = filteredCouples[i];
825 
826  if
827  (
828  hit1[i].hit()
829  && hit2[i].hit()
830  && (
831  surface1[i] != surface2[i]
832  || hit1[i].index() != hit2[i].index()
833  )
834  )
835  {
836  // Two different hits. Check angle.
837  //str.write
838  //(
839  // linePointRef(hit1[i].hitPoint(), hit2[i].hitPoint()),
840  // normal1[i],
841  // normal2[i]
842  //);
843 
844  if ((normal1[i]&normal2[i]) > planarAngleCos)
845  {
846  // Both normals aligned
847  vector n = end[i]-start[i];
848  scalar magN = mag(n);
849  if (magN > VSMALL)
850  {
851  filteredCouples[filterI++] = couple;
852  }
853  }
854  }
855  else if (hit1[i].hit() || hit2[i].hit())
856  {
857  // Single hit. Do not include in freestanding baffles.
858  }
859  }
860 
861  filteredCouples.setSize(filterI);
862 
863  Info<< "freeStandingBaffles : detected "
864  << returnReduce(filterI, sumOp<label>())
865  << " planar (within " << planarAngle
866  << " degrees) free-standing baffles out of "
867  << nFiltered
868  << nl << endl;
869  }
870 
871  return filteredCouples;
872 }
873 
874 
876 (
877  const List<labelPair>& couples
878 )
879 {
880  // Mesh change engine
881  polyTopoChange meshMod(mesh_);
882 
883  const faceList& faces = mesh_.faces();
884  const labelList& faceOwner = mesh_.faceOwner();
885  const faceZoneMesh& faceZones = mesh_.faceZones();
886 
887  forAll(couples, i)
888  {
889  label face0 = couples[i].first();
890  label face1 = couples[i].second();
891 
892  // face1 < 0 signals a coupled face that has been converted to baffle.
893 
894  label own0 = faceOwner[face0];
895  label own1 = faceOwner[face1];
896 
897  if (face1 < 0 || own0 < own1)
898  {
899  // Use face0 as the new internal face.
900  label zoneID = faceZones.whichZone(face0);
901  bool zoneFlip = false;
902 
903  if (zoneID >= 0)
904  {
905  const faceZone& fZone = faceZones[zoneID];
906  zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
907  }
908 
909  label nei = (face1 < 0 ? -1 : own1);
910 
911  meshMod.setAction(polyRemoveFace(face1));
912  meshMod.setAction
913  (
915  (
916  faces[face0], // modified face
917  face0, // label of face being modified
918  own0, // owner
919  nei, // neighbour
920  false, // face flip
921  -1, // patch for face
922  false, // remove from zone
923  zoneID, // zone for face
924  zoneFlip // face flip in zone
925  )
926  );
927  }
928  else
929  {
930  // Use face1 as the new internal face.
931  label zoneID = faceZones.whichZone(face1);
932  bool zoneFlip = false;
933 
934  if (zoneID >= 0)
935  {
936  const faceZone& fZone = faceZones[zoneID];
937  zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
938  }
939 
940  meshMod.setAction(polyRemoveFace(face0));
941  meshMod.setAction
942  (
944  (
945  faces[face1], // modified face
946  face1, // label of face being modified
947  own1, // owner
948  own0, // neighbour
949  false, // face flip
950  -1, // patch for face
951  false, // remove from zone
952  zoneID, // zone for face
953  zoneFlip // face flip in zone
954  )
955  );
956  }
957  }
958 
959  // Change the mesh (no inflation)
960  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
961 
962  // Update fields
963  mesh_.updateMesh(map);
964 
965  // Move mesh (since morphing does not do this)
966  if (map().hasMotionPoints())
967  {
968  mesh_.movePoints(map().preMotionPoints());
969  }
970  else
971  {
972  // Delete mesh volumes.
973  mesh_.clearOut();
974  }
975 
976  // Reset the instance for if in overwrite mode
977  mesh_.setInstance(timeName());
978 
979  // Update intersections. Recalculate intersections on merged faces since
980  // this seems to give problems? Note: should not be necessary since
981  // baffles preserve intersections from when they were created.
982  labelList newExposedFaces(2*couples.size());
983  label newI = 0;
984 
985  forAll(couples, i)
986  {
987  label newFace0 = map().reverseFaceMap()[couples[i].first()];
988  if (newFace0 != -1)
989  {
990  newExposedFaces[newI++] = newFace0;
991  }
992 
993  label newFace1 = map().reverseFaceMap()[couples[i].second()];
994  if (newFace1 != -1)
995  {
996  newExposedFaces[newI++] = newFace1;
997  }
998  }
999  newExposedFaces.setSize(newI);
1000  updateMesh(map, newExposedFaces);
1001 
1002  return map;
1003 }
1004 
1005 
1006 // Finds region per cell for cells inside closed named surfaces
1007 void Foam::meshRefinement::findCellZoneGeometric
1008 (
1009  const pointField& neiCc,
1010  const labelList& closedNamedSurfaces, // indices of closed surfaces
1011  labelList& namedSurfaceIndex, // per face index of named surface
1012  const labelList& surfaceToCellZone, // cell zone index per surface
1013 
1014  labelList& cellToZone
1015 ) const
1016 {
1017  const pointField& cellCentres = mesh_.cellCentres();
1018  const labelList& faceOwner = mesh_.faceOwner();
1019  const labelList& faceNeighbour = mesh_.faceNeighbour();
1020 
1021  // Check if cell centre is inside
1022  labelList insideSurfaces;
1023  surfaces_.findInside
1024  (
1025  closedNamedSurfaces,
1026  cellCentres,
1027  insideSurfaces
1028  );
1029 
1030  forAll(insideSurfaces, cellI)
1031  {
1032  if (cellToZone[cellI] == -2)
1033  {
1034  label surfI = insideSurfaces[cellI];
1035 
1036  if (surfI != -1)
1037  {
1038  cellToZone[cellI] = surfaceToCellZone[surfI];
1039  }
1040  }
1041  }
1042 
1043 
1044  // Some cells with cell centres close to surface might have
1045  // had been put into wrong surface. Recheck with perturbed cell centre.
1046 
1047 
1048  // 1. Collect points
1049 
1050  // Count points to test.
1051  label nCandidates = 0;
1052  forAll(namedSurfaceIndex, faceI)
1053  {
1054  label surfI = namedSurfaceIndex[faceI];
1055 
1056  if (surfI != -1)
1057  {
1058  if (mesh_.isInternalFace(faceI))
1059  {
1060  nCandidates += 2;
1061  }
1062  else
1063  {
1064  nCandidates += 1;
1065  }
1066  }
1067  }
1068 
1069  // Collect points.
1070  pointField candidatePoints(nCandidates);
1071  nCandidates = 0;
1072  forAll(namedSurfaceIndex, faceI)
1073  {
1074  label surfI = namedSurfaceIndex[faceI];
1075 
1076  if (surfI != -1)
1077  {
1078  label own = faceOwner[faceI];
1079  const point& ownCc = cellCentres[own];
1080 
1081  if (mesh_.isInternalFace(faceI))
1082  {
1083  label nei = faceNeighbour[faceI];
1084  const point& neiCc = cellCentres[nei];
1085  // Perturbed cc
1086  const vector d = 1e-4*(neiCc - ownCc);
1087  candidatePoints[nCandidates++] = ownCc-d;
1088  candidatePoints[nCandidates++] = neiCc+d;
1089  }
1090  else
1091  {
1092  //const point& neiFc = mesh_.faceCentres()[faceI];
1093  const point& neiFc = neiCc[faceI-mesh_.nInternalFaces()];
1094 
1095  // Perturbed cc
1096  const vector d = 1e-4*(neiFc - ownCc);
1097  candidatePoints[nCandidates++] = ownCc-d;
1098  }
1099  }
1100  }
1101 
1102 
1103  // 2. Test points for inside
1104 
1105  surfaces_.findInside
1106  (
1107  closedNamedSurfaces,
1108  candidatePoints,
1109  insideSurfaces
1110  );
1111 
1112 
1113  // 3. Update zone information
1114 
1115  nCandidates = 0;
1116  forAll(namedSurfaceIndex, faceI)
1117  {
1118  label surfI = namedSurfaceIndex[faceI];
1119 
1120  if (surfI != -1)
1121  {
1122  label own = faceOwner[faceI];
1123 
1124  if (mesh_.isInternalFace(faceI))
1125  {
1126  label ownSurfI = insideSurfaces[nCandidates++];
1127  if (ownSurfI != -1)
1128  {
1129  cellToZone[own] = surfaceToCellZone[ownSurfI];
1130  }
1131 
1132  label neiSurfI = insideSurfaces[nCandidates++];
1133  if (neiSurfI != -1)
1134  {
1135  label nei = faceNeighbour[faceI];
1136 
1137  cellToZone[nei] = surfaceToCellZone[neiSurfI];
1138  }
1139  }
1140  else
1141  {
1142  label ownSurfI = insideSurfaces[nCandidates++];
1143  if (ownSurfI != -1)
1144  {
1145  cellToZone[own] = surfaceToCellZone[ownSurfI];
1146  }
1147  }
1148  }
1149  }
1150 
1151 
1152  // Adapt the namedSurfaceIndex
1153  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
1154  // for if any cells were not completely covered.
1155 
1156  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1157  {
1158  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1159  label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
1160 
1161  if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
1162  {
1163  // Give face the zone of max cell zone
1164  namedSurfaceIndex[faceI] = findIndex
1165  (
1166  surfaceToCellZone,
1167  max(ownZone, neiZone)
1168  );
1169  }
1170  }
1171 
1172  labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
1173  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1174 
1175  forAll(patches, patchI)
1176  {
1177  const polyPatch& pp = patches[patchI];
1178 
1179  if (pp.coupled())
1180  {
1181  forAll(pp, i)
1182  {
1183  label faceI = pp.start()+i;
1184  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1185  neiCellZone[faceI-mesh_.nInternalFaces()] = ownZone;
1186  }
1187  }
1188  }
1189  syncTools::swapBoundaryFaceList(mesh_, neiCellZone);
1190 
1191  forAll(patches, patchI)
1192  {
1193  const polyPatch& pp = patches[patchI];
1194 
1195  if (pp.coupled())
1196  {
1197  forAll(pp, i)
1198  {
1199  label faceI = pp.start()+i;
1200  label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
1201  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
1202 
1203  if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
1204  {
1205  // Give face the max cell zone
1206  namedSurfaceIndex[faceI] = findIndex
1207  (
1208  surfaceToCellZone,
1209  max(ownZone, neiZone)
1210  );
1211  }
1212  }
1213  }
1214  }
1215 
1216  // Sync
1217  syncTools::syncFaceList(mesh_, namedSurfaceIndex, maxEqOp<label>());
1218 }
1219 
1220 
1221 void Foam::meshRefinement::findCellZoneInsideWalk
1222 (
1223  const labelList& locationSurfaces, // indices of surfaces with inside point
1224  const labelList& namedSurfaceIndex, // per face index of named surface
1225  const labelList& surfaceToCellZone, // cell zone index per surface
1226 
1227  labelList& cellToZone
1228 ) const
1229 {
1230  // Analyse regions. Reuse regionsplit
1231  boolList blockedFace(mesh_.nFaces());
1232  //selectSeparatedCoupledFaces(blockedFace);
1233 
1234  forAll(namedSurfaceIndex, faceI)
1235  {
1236  if (namedSurfaceIndex[faceI] == -1)
1237  {
1238  blockedFace[faceI] = false;
1239  }
1240  else
1241  {
1242  blockedFace[faceI] = true;
1243  }
1244  }
1245  // No need to sync since namedSurfaceIndex already is synced
1246 
1247  // Set region per cell based on walking
1248  regionSplit cellRegion(mesh_, blockedFace);
1249  blockedFace.clear();
1250 
1251 
1252  // Force calculation of face decomposition (used in findCell)
1253  (void)mesh_.tetBasePtIs();
1254 
1255  const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
1256 
1257  // For all locationSurface find the cell
1258  forAll(locationSurfaces, i)
1259  {
1260  label surfI = locationSurfaces[i];
1261 
1262  const point& insidePoint = surfZones[surfI].zoneInsidePoint();
1263 
1264  Info<< "For surface " << surfaces_.names()[surfI]
1265  << " finding inside point " << insidePoint
1266  << endl;
1267 
1268  // Find the region containing the insidePoint
1269  label keepRegionI = findRegion
1270  (
1271  mesh_,
1272  cellRegion,
1273  mergeDistance_*vector(1,1,1),
1274  insidePoint
1275  );
1276 
1277  Info<< "For surface " << surfaces_.names()[surfI]
1278  << " found point " << insidePoint
1279  << " in global region " << keepRegionI
1280  << " out of " << cellRegion.nRegions() << " regions." << endl;
1281 
1282  if (keepRegionI == -1)
1283  {
1285  << "Point " << insidePoint
1286  << " is not inside the mesh." << nl
1287  << "Bounding box of the mesh:" << mesh_.bounds()
1288  << exit(FatalError);
1289  }
1290 
1291  // Set all cells with this region
1292  forAll(cellRegion, cellI)
1293  {
1294  if (cellRegion[cellI] == keepRegionI)
1295  {
1296  if (cellToZone[cellI] == -2)
1297  {
1298  cellToZone[cellI] = surfaceToCellZone[surfI];
1299  }
1300  else if (cellToZone[cellI] != surfaceToCellZone[surfI])
1301  {
1303  << "Cell " << cellI
1304  << " at " << mesh_.cellCentres()[cellI]
1305  << " is inside surface " << surfaces_.names()[surfI]
1306  << " but already marked as being in zone "
1307  << cellToZone[cellI] << endl
1308  << "This can happen if your surfaces are not"
1309  << " (sufficiently) closed."
1310  << endl;
1311  }
1312  }
1313  }
1314  }
1315 }
1316 
1317 
1318 bool Foam::meshRefinement::calcRegionToZone
1319 (
1320  const label surfZoneI,
1321  const label ownRegion,
1322  const label neiRegion,
1323 
1324  labelList& regionToCellZone
1325 ) const
1326 {
1327  bool changed = false;
1328 
1329  // Check whether inbetween different regions
1330  if (ownRegion != neiRegion)
1331  {
1332  // Jump. Change one of the sides to my type.
1333 
1334  // 1. Interface between my type and unset region.
1335  // Set region to keepRegion
1336 
1337  if (regionToCellZone[ownRegion] == -2)
1338  {
1339  if (regionToCellZone[neiRegion] == surfZoneI)
1340  {
1341  // Face between unset and my region. Put unset
1342  // region into keepRegion
1343  regionToCellZone[ownRegion] = -1;
1344  changed = true;
1345  }
1346  else if (regionToCellZone[neiRegion] != -2)
1347  {
1348  // Face between unset and other region.
1349  // Put unset region into my region
1350  regionToCellZone[ownRegion] = surfZoneI;
1351  changed = true;
1352  }
1353  }
1354  else if (regionToCellZone[neiRegion] == -2)
1355  {
1356  if (regionToCellZone[ownRegion] == surfZoneI)
1357  {
1358  // Face between unset and my region. Put unset
1359  // region into keepRegion
1360  regionToCellZone[neiRegion] = -1;
1361  changed = true;
1362  }
1363  else if (regionToCellZone[ownRegion] != -2)
1364  {
1365  // Face between unset and other region.
1366  // Put unset region into my region
1367  regionToCellZone[neiRegion] = surfZoneI;
1368  changed = true;
1369  }
1370  }
1371  }
1372  return changed;
1373 }
1374 
1375 
1376 void Foam::meshRefinement::findCellZoneTopo
1377 (
1378  const point& keepPoint,
1379  const labelList& namedSurfaceIndex,
1380  const labelList& surfaceToCellZone,
1381  labelList& cellToZone
1382 ) const
1383 {
1384  // Assumes:
1385  // - region containing keepPoint does not go into a cellZone
1386  // - all other regions can be found by crossing faces marked in
1387  // namedSurfaceIndex.
1388 
1389  // Analyse regions. Reuse regionsplit
1390  boolList blockedFace(mesh_.nFaces());
1391 
1392  forAll(namedSurfaceIndex, faceI)
1393  {
1394  if (namedSurfaceIndex[faceI] == -1)
1395  {
1396  blockedFace[faceI] = false;
1397  }
1398  else
1399  {
1400  blockedFace[faceI] = true;
1401  }
1402  }
1403  // No need to sync since namedSurfaceIndex already is synced
1404 
1405  // Set region per cell based on walking
1406  regionSplit cellRegion(mesh_, blockedFace);
1407  blockedFace.clear();
1408 
1409  // Per mesh region the zone the cell should be put in.
1410  // -2 : not analysed yet
1411  // -1 : keepPoint region. Not put into any cellzone.
1412  // >= 0 : index of cellZone
1413  labelList regionToCellZone(cellRegion.nRegions(), -2);
1414 
1415  // See which cells already are set in the cellToZone (from geometric
1416  // searching) and use these to take over their zones.
1417  // Note: could be improved to count number of cells per region.
1418  forAll(cellToZone, cellI)
1419  {
1420  if (cellToZone[cellI] != -2)
1421  {
1422  regionToCellZone[cellRegion[cellI]] = cellToZone[cellI];
1423  }
1424  }
1425 
1426 
1427  // Find the region containing the keepPoint
1428  label keepRegionI = findRegion
1429  (
1430  mesh_,
1431  cellRegion,
1432  mergeDistance_*vector(1,1,1),
1433  keepPoint
1434  );
1435 
1436  Info<< "Found point " << keepPoint
1437  << " in global region " << keepRegionI
1438  << " out of " << cellRegion.nRegions() << " regions." << endl;
1439 
1440  if (keepRegionI == -1)
1441  {
1443  << "Point " << keepPoint
1444  << " is not inside the mesh." << nl
1445  << "Bounding box of the mesh:" << mesh_.bounds()
1446  << exit(FatalError);
1447  }
1448 
1449  // Mark default region with zone -1.
1450  if (regionToCellZone[keepRegionI] == -2)
1451  {
1452  regionToCellZone[keepRegionI] = -1;
1453  }
1454 
1455 
1456  // Find correspondence between cell zone and surface
1457  // by changing cell zone every time we cross a surface.
1458  while (true)
1459  {
1460  // Synchronise regionToCellZone.
1461  // Note:
1462  // - region numbers are identical on all processors
1463  // - keepRegion is identical ,,
1464  // - cellZones are identical ,,
1465  // This done at top of loop to account for geometric matching
1466  // not being synchronised.
1467  Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
1468  Pstream::listCombineScatter(regionToCellZone);
1469 
1470 
1471  bool changed = false;
1472 
1473  // Internal faces
1474 
1475  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1476  {
1477  label surfI = namedSurfaceIndex[faceI];
1478 
1479  // Connected even if no cellZone defined for surface
1480  if (surfI != -1)
1481  {
1482  // Calculate region to zone from cellRegions on either side
1483  // of internal face.
1484  bool changedCell = calcRegionToZone
1485  (
1486  surfaceToCellZone[surfI],
1487  cellRegion[mesh_.faceOwner()[faceI]],
1488  cellRegion[mesh_.faceNeighbour()[faceI]],
1489  regionToCellZone
1490  );
1491 
1492  changed = changed | changedCell;
1493  }
1494  }
1495 
1496  // Boundary faces
1497 
1498  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1499 
1500  // Get coupled neighbour cellRegion
1501  labelList neiCellRegion(mesh_.nFaces()-mesh_.nInternalFaces());
1502  forAll(patches, patchI)
1503  {
1504  const polyPatch& pp = patches[patchI];
1505 
1506  if (pp.coupled())
1507  {
1508  forAll(pp, i)
1509  {
1510  label faceI = pp.start()+i;
1511  neiCellRegion[faceI-mesh_.nInternalFaces()] =
1512  cellRegion[mesh_.faceOwner()[faceI]];
1513  }
1514  }
1515  }
1516  syncTools::swapBoundaryFaceList(mesh_, neiCellRegion);
1517 
1518  // Calculate region to zone from cellRegions on either side of coupled
1519  // face.
1520  forAll(patches, patchI)
1521  {
1522  const polyPatch& pp = patches[patchI];
1523 
1524  if (pp.coupled())
1525  {
1526  forAll(pp, i)
1527  {
1528  label faceI = pp.start()+i;
1529 
1530  label surfI = namedSurfaceIndex[faceI];
1531 
1532  // Connected even if no cellZone defined for surface
1533  if (surfI != -1)
1534  {
1535  bool changedCell = calcRegionToZone
1536  (
1537  surfaceToCellZone[surfI],
1538  cellRegion[mesh_.faceOwner()[faceI]],
1539  neiCellRegion[faceI-mesh_.nInternalFaces()],
1540  regionToCellZone
1541  );
1542 
1543  changed = changed | changedCell;
1544  }
1545  }
1546  }
1547  }
1548 
1549  if (!returnReduce(changed, orOp<bool>()))
1550  {
1551  break;
1552  }
1553  }
1554 
1555 
1556  forAll(regionToCellZone, regionI)
1557  {
1558  label zoneI = regionToCellZone[regionI];
1559 
1560  if (zoneI == -2)
1561  {
1563  << "For region " << regionI << " haven't set cell zone."
1564  << exit(FatalError);
1565  }
1566  }
1567 
1568  if (debug)
1569  {
1570  forAll(regionToCellZone, regionI)
1571  {
1572  Pout<< "Region " << regionI
1573  << " becomes cellZone:" << regionToCellZone[regionI]
1574  << endl;
1575  }
1576  }
1577 
1578  // Rework into cellToZone
1579  forAll(cellToZone, cellI)
1580  {
1581  cellToZone[cellI] = regionToCellZone[cellRegion[cellI]];
1582  }
1583 }
1584 
1585 
1586 void Foam::meshRefinement::makeConsistentFaceIndex
1587 (
1588  const labelList& cellToZone,
1589  labelList& namedSurfaceIndex
1590 ) const
1591 {
1592  const labelList& faceOwner = mesh_.faceOwner();
1593  const labelList& faceNeighbour = mesh_.faceNeighbour();
1594 
1595  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1596  {
1597  label ownZone = cellToZone[faceOwner[faceI]];
1598  label neiZone = cellToZone[faceNeighbour[faceI]];
1599 
1600  if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
1601  {
1602  namedSurfaceIndex[faceI] = -1;
1603  }
1604  else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
1605  {
1607  << "Different cell zones on either side of face " << faceI
1608  << " at " << mesh_.faceCentres()[faceI]
1609  << " but face not marked with a surface."
1610  << abort(FatalError);
1611  }
1612  }
1613 
1614  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1615 
1616  // Get coupled neighbour cellZone
1617  labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces());
1618  forAll(patches, patchI)
1619  {
1620  const polyPatch& pp = patches[patchI];
1621 
1622  if (pp.coupled())
1623  {
1624  forAll(pp, i)
1625  {
1626  label faceI = pp.start()+i;
1627  neiCellZone[faceI-mesh_.nInternalFaces()] =
1628  cellToZone[mesh_.faceOwner()[faceI]];
1629  }
1630  }
1631  }
1632  syncTools::swapBoundaryFaceList(mesh_, neiCellZone);
1633 
1634  // Use coupled cellZone to do check
1635  forAll(patches, patchI)
1636  {
1637  const polyPatch& pp = patches[patchI];
1638 
1639  if (pp.coupled())
1640  {
1641  forAll(pp, i)
1642  {
1643  label faceI = pp.start()+i;
1644 
1645  label ownZone = cellToZone[faceOwner[faceI]];
1646  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
1647 
1648  if (ownZone == neiZone && namedSurfaceIndex[faceI] != -1)
1649  {
1650  namedSurfaceIndex[faceI] = -1;
1651  }
1652  else if (ownZone != neiZone && namedSurfaceIndex[faceI] == -1)
1653  {
1655  << "Different cell zones on either side of face "
1656  << faceI << " at " << mesh_.faceCentres()[faceI]
1657  << " but face not marked with a surface."
1658  << abort(FatalError);
1659  }
1660  }
1661  }
1662  else
1663  {
1664  // Unzonify boundary faces
1665  forAll(pp, i)
1666  {
1667  label faceI = pp.start()+i;
1668  namedSurfaceIndex[faceI] = -1;
1669  }
1670  }
1671  }
1672 }
1673 
1674 
1675 void Foam::meshRefinement::handleSnapProblems
1676 (
1677  const snapParameters& snapParams,
1678  const bool useTopologicalSnapDetection,
1679  const bool removeEdgeConnectedCells,
1680  const scalarField& perpendicularAngle,
1681  const dictionary& motionDict,
1682  Time& runTime,
1683  const labelList& globalToMasterPatch,
1684  const labelList& globalToSlavePatch
1685 )
1686 {
1687  Info<< nl
1688  << "Introducing baffles to block off problem cells" << nl
1689  << "----------------------------------------------" << nl
1690  << endl;
1691 
1692  labelList facePatch;
1693  if (useTopologicalSnapDetection)
1694  {
1695  facePatch = markFacesOnProblemCells
1696  (
1697  motionDict,
1698  removeEdgeConnectedCells,
1699  perpendicularAngle,
1700  globalToMasterPatch
1701  );
1702  }
1703  else
1704  {
1705  facePatch = markFacesOnProblemCellsGeometric(snapParams, motionDict);
1706  }
1707  Info<< "Analyzed problem cells in = "
1708  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1709 
1710  if (debug&MESH)
1711  {
1712  faceSet problemFaces(mesh_, "problemFaces", mesh_.nFaces()/100);
1713 
1714  forAll(facePatch, faceI)
1715  {
1716  if (facePatch[faceI] != -1)
1717  {
1718  problemFaces.insert(faceI);
1719  }
1720  }
1721  problemFaces.instance() = timeName();
1722  Pout<< "Dumping " << problemFaces.size()
1723  << " problem faces to " << problemFaces.objectPath() << endl;
1724  problemFaces.write();
1725  }
1726 
1727  Info<< "Introducing baffles to delete problem cells." << nl << endl;
1728 
1729  if (debug)
1730  {
1731  runTime++;
1732  }
1733 
1734  // Create baffles with same owner and neighbour for now.
1735  createBaffles(facePatch, facePatch);
1736 
1737  if (debug)
1738  {
1739  // Debug:test all is still synced across proc patches
1740  checkData();
1741  }
1742  Info<< "Created baffles in = "
1743  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1744 
1745  printMeshInfo(debug, "After introducing baffles");
1746 
1747  if (debug&MESH)
1748  {
1749  Pout<< "Writing extra baffled mesh to time "
1750  << timeName() << endl;
1751  write
1752  (
1753  debugType(debug),
1755  runTime.path()/"extraBaffles"
1756  );
1757  Pout<< "Dumped debug data in = "
1758  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
1759  }
1760 }
1761 
1762 
1763 Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
1764 (
1765  const labelList& faceToZone,
1766  const labelList& cellToZone,
1767  const labelList& neiCellZone
1768 ) const
1769 {
1770  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
1771  const labelList& faceOwner = mesh_.faceOwner();
1772  const labelList& faceNeighbour = mesh_.faceNeighbour();
1773 
1774 
1775  // We want to pick up the faces to orient. These faces come in
1776  // two variants:
1777  // - faces originating from stand-alone faceZones
1778  // (these will most likely have no cellZone on either side so
1779  // ownZone and neiZone both -1)
1780  // - sticky-up faces originating from a 'bulge' in a outside of
1781  // a cellZone. These will have the same cellZone on either side.
1782  // How to orient these is not really clearly defined so do them
1783  // same as stand-alone faceZone faces for now. (Normally these will
1784  // already have been removed by the 'allowFreeStandingZoneFaces=false'
1785  // default setting)
1786 
1787  // Note that argument neiCellZone will have -1 on uncoupled boundaries.
1788 
1789  DynamicList<label> faceLabels(mesh_.nFaces()/100);
1790 
1791  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1792  {
1793  if (faceToZone[faceI] != -1)
1794  {
1795  // Free standing baffle?
1796  label ownZone = cellToZone[faceOwner[faceI]];
1797  label neiZone = cellToZone[faceNeighbour[faceI]];
1798  if (ownZone == neiZone)
1799  {
1800  faceLabels.append(faceI);
1801  }
1802  }
1803  }
1804  forAll(patches, patchI)
1805  {
1806  const polyPatch& pp = patches[patchI];
1807 
1808  forAll(pp, i)
1809  {
1810  label faceI = pp.start()+i;
1811  if (faceToZone[faceI] != -1)
1812  {
1813  // Free standing baffle?
1814  label ownZone = cellToZone[faceOwner[faceI]];
1815  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
1816  if (ownZone == neiZone)
1817  {
1818  faceLabels.append(faceI);
1819  }
1820  }
1821  }
1822  }
1823  return faceLabels.shrink();
1824 }
1825 
1826 
1827 void Foam::meshRefinement::calcPatchNumMasterFaces
1828 (
1829  const PackedBoolList& isMasterFace,
1830  const indirectPrimitivePatch& patch,
1831  labelList& nMasterFacesPerEdge
1832 ) const
1833 {
1834  // Number of (master)faces per edge
1835  nMasterFacesPerEdge.setSize(patch.nEdges());
1836  nMasterFacesPerEdge = 0;
1837 
1838  forAll(patch.addressing(), faceI)
1839  {
1840  const label meshFaceI = patch.addressing()[faceI];
1841 
1842  if (isMasterFace[meshFaceI])
1843  {
1844  const labelList& fEdges = patch.faceEdges()[faceI];
1845  forAll(fEdges, fEdgeI)
1846  {
1847  nMasterFacesPerEdge[fEdges[fEdgeI]]++;
1848  }
1849  }
1850  }
1851 
1853  (
1854  mesh_,
1855  patch.meshEdges(mesh_.edges(), mesh_.pointEdges()),
1856  nMasterFacesPerEdge,
1857  plusEqOp<label>(),
1858  label(0)
1859  );
1860 }
1861 
1862 
1863 Foam::label Foam::meshRefinement::markPatchZones
1864 (
1865  const indirectPrimitivePatch& patch,
1866  const labelList& nMasterFacesPerEdge,
1867  labelList& faceToZone
1868 ) const
1869 {
1870  List<patchEdgeFaceRegion> allEdgeInfo(patch.nEdges());
1871  List<patchEdgeFaceRegion> allFaceInfo(patch.size());
1872 
1873 
1874  // Protect all non-manifold edges
1875  {
1876  label nProtected = 0;
1877 
1878  forAll(nMasterFacesPerEdge, edgeI)
1879  {
1880  if (nMasterFacesPerEdge[edgeI] > 2)
1881  {
1882  allEdgeInfo[edgeI] = -2;
1883  nProtected++;
1884  }
1885  }
1886  //Info<< "Protected from visiting "
1887  // << returnReduce(nProtected, sumOp<label>())
1888  // << " non-manifold edges" << nl << endl;
1889  }
1890 
1891 
1892  // Hand out zones
1893 
1894  DynamicList<label> changedEdges;
1896 
1897  const scalar tol = PatchEdgeFaceWave
1898  <
1901  >::propagationTol();
1902 
1903  int dummyTrackData;
1904 
1905  const globalIndex globalFaces(patch.size());
1906 
1907  label faceI = 0;
1908 
1909  label currentZoneI = 0;
1910 
1911  while (true)
1912  {
1913  // Pick an unset face
1914  label globalSeed = labelMax;
1915  for (; faceI < allFaceInfo.size(); faceI++)
1916  {
1917  if (!allFaceInfo[faceI].valid(dummyTrackData))
1918  {
1919  globalSeed = globalFaces.toGlobal(faceI);
1920  break;
1921  }
1922  }
1923 
1924  reduce(globalSeed, minOp<label>());
1925 
1926  if (globalSeed == labelMax)
1927  {
1928  break;
1929  }
1930 
1931  label procI = globalFaces.whichProcID(globalSeed);
1932  label seedFaceI = globalFaces.toLocal(procI, globalSeed);
1933 
1934  //Info<< "Seeding zone " << currentZoneI
1935  // << " from processor " << procI << " face " << seedFaceI
1936  // << endl;
1937 
1938  if (procI == Pstream::myProcNo())
1939  {
1940  patchEdgeFaceRegion& faceInfo = allFaceInfo[seedFaceI];
1941 
1942 
1943  // Set face
1944  faceInfo = currentZoneI;
1945 
1946  // .. and seed its edges
1947  const labelList& fEdges = patch.faceEdges()[seedFaceI];
1948  forAll(fEdges, fEdgeI)
1949  {
1950  label edgeI = fEdges[fEdgeI];
1951 
1952  patchEdgeFaceRegion& edgeInfo = allEdgeInfo[edgeI];
1953 
1954  if
1955  (
1956  edgeInfo.updateEdge<int>
1957  (
1958  mesh_,
1959  patch,
1960  edgeI,
1961  seedFaceI,
1962  faceInfo,
1963  tol,
1964  dummyTrackData
1965  )
1966  )
1967  {
1968  changedEdges.append(edgeI);
1969  changedInfo.append(edgeInfo);
1970  }
1971  }
1972  }
1973 
1974 
1975  if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
1976  {
1977  break;
1978  }
1979 
1980 
1981  // Walk
1983  <
1986  > calc
1987  (
1988  mesh_,
1989  patch,
1990  changedEdges,
1991  changedInfo,
1992  allEdgeInfo,
1993  allFaceInfo,
1994  returnReduce(patch.nEdges(), sumOp<label>())
1995  );
1996 
1997  currentZoneI++;
1998  }
1999 
2000 
2001  faceToZone.setSize(patch.size());
2002  forAll(allFaceInfo, faceI)
2003  {
2004  if (!allFaceInfo[faceI].valid(dummyTrackData))
2005  {
2007  << "Problem: unvisited face " << faceI
2008  << " at " << patch.faceCentres()[faceI]
2009  << exit(FatalError);
2010  }
2011  faceToZone[faceI] = allFaceInfo[faceI].region();
2012  }
2013 
2014  return currentZoneI;
2015 }
2016 
2017 
2018 void Foam::meshRefinement::consistentOrientation
2019 (
2020  const PackedBoolList& isMasterFace,
2021  const indirectPrimitivePatch& patch,
2022  const labelList& nMasterFacesPerEdge,
2023  const labelList& faceToZone,
2024  const Map<label>& zoneToOrientation,
2025  boolList& meshFlipMap
2026 ) const
2027 {
2028  const polyBoundaryMesh& bm = mesh_.boundaryMesh();
2029 
2030  // Data on all edges and faces
2031  List<patchFaceOrientation> allEdgeInfo(patch.nEdges());
2032  List<patchFaceOrientation> allFaceInfo(patch.size());
2033 
2034  // Make sure we don't walk through
2035  // - slaves of coupled faces
2036  // - non-manifold edges
2037  {
2038  label nProtected = 0;
2039 
2040  forAll(patch.addressing(), faceI)
2041  {
2042  const label meshFaceI = patch.addressing()[faceI];
2043  const label patchI = bm.whichPatch(meshFaceI);
2044 
2045  if
2046  (
2047  patchI != -1
2048  && bm[patchI].coupled()
2049  && !isMasterFace[meshFaceI]
2050  )
2051  {
2052  // Slave side. Mark so doesn't get visited.
2053  allFaceInfo[faceI] = orientedSurface::NOFLIP;
2054  nProtected++;
2055  }
2056  }
2057  //Info<< "Protected from visiting "
2058  // << returnReduce(nProtected, sumOp<label>())
2059  // << " slaves of coupled faces" << nl << endl;
2060  }
2061  {
2062  label nProtected = 0;
2063 
2064  forAll(nMasterFacesPerEdge, edgeI)
2065  {
2066  if (nMasterFacesPerEdge[edgeI] > 2)
2067  {
2068  allEdgeInfo[edgeI] = orientedSurface::NOFLIP;
2069  nProtected++;
2070  }
2071  }
2072 
2073  Info<< "Protected from visiting "
2074  << returnReduce(nProtected, sumOp<label>())
2075  << " non-manifold edges" << nl << endl;
2076  }
2077 
2078 
2079 
2080  DynamicList<label> changedEdges;
2082 
2083  const scalar tol = PatchEdgeFaceWave
2084  <
2087  >::propagationTol();
2088 
2089  int dummyTrackData;
2090 
2091  globalIndex globalFaces(patch.size());
2092 
2093  while (true)
2094  {
2095  // Pick an unset face
2096  label globalSeed = labelMax;
2097  forAll(allFaceInfo, faceI)
2098  {
2099  if (allFaceInfo[faceI] == orientedSurface::UNVISITED)
2100  {
2101  globalSeed = globalFaces.toGlobal(faceI);
2102  break;
2103  }
2104  }
2105 
2106  reduce(globalSeed, minOp<label>());
2107 
2108  if (globalSeed == labelMax)
2109  {
2110  break;
2111  }
2112 
2113  label procI = globalFaces.whichProcID(globalSeed);
2114  label seedFaceI = globalFaces.toLocal(procI, globalSeed);
2115 
2116  //Info<< "Seeding from processor " << procI << " face " << seedFaceI
2117  // << endl;
2118 
2119  if (procI == Pstream::myProcNo())
2120  {
2121  // Determine orientation of seedFace
2122 
2123  patchFaceOrientation& faceInfo = allFaceInfo[seedFaceI];
2124 
2125  // Start off with correct orientation
2126  faceInfo = orientedSurface::NOFLIP;
2127 
2128  if (zoneToOrientation[faceToZone[seedFaceI]] < 0)
2129  {
2130  faceInfo.flip();
2131  }
2132 
2133 
2134  const labelList& fEdges = patch.faceEdges()[seedFaceI];
2135  forAll(fEdges, fEdgeI)
2136  {
2137  label edgeI = fEdges[fEdgeI];
2138 
2139  patchFaceOrientation& edgeInfo = allEdgeInfo[edgeI];
2140 
2141  if
2142  (
2143  edgeInfo.updateEdge<int>
2144  (
2145  mesh_,
2146  patch,
2147  edgeI,
2148  seedFaceI,
2149  faceInfo,
2150  tol,
2151  dummyTrackData
2152  )
2153  )
2154  {
2155  changedEdges.append(edgeI);
2156  changedInfo.append(edgeInfo);
2157  }
2158  }
2159  }
2160 
2161 
2162  if (returnReduce(changedEdges.size(), sumOp<label>()) == 0)
2163  {
2164  break;
2165  }
2166 
2167 
2168 
2169  // Walk
2171  <
2174  > calc
2175  (
2176  mesh_,
2177  patch,
2178  changedEdges,
2179  changedInfo,
2180  allEdgeInfo,
2181  allFaceInfo,
2182  returnReduce(patch.nEdges(), sumOp<label>())
2183  );
2184  }
2185 
2186 
2187  // Push master zone info over to slave (since slave faces never visited)
2188  {
2189  labelList neiStatus
2190  (
2191  mesh_.nFaces()-mesh_.nInternalFaces(),
2193  );
2194 
2195  forAll(patch.addressing(), i)
2196  {
2197  const label meshFaceI = patch.addressing()[i];
2198  if (!mesh_.isInternalFace(meshFaceI))
2199  {
2200  neiStatus[meshFaceI-mesh_.nInternalFaces()] =
2201  allFaceInfo[i].flipStatus();
2202  }
2203  }
2204  syncTools::swapBoundaryFaceList(mesh_, neiStatus);
2205 
2206  forAll(patch.addressing(), i)
2207  {
2208  const label meshFaceI = patch.addressing()[i];
2209  const label patchI = bm.whichPatch(meshFaceI);
2210 
2211  if
2212  (
2213  patchI != -1
2214  && bm[patchI].coupled()
2215  && !isMasterFace[meshFaceI]
2216  )
2217  {
2218  // Slave side. Take flipped from neighbour
2219  label bFaceI = meshFaceI-mesh_.nInternalFaces();
2220 
2221  if (neiStatus[bFaceI] == orientedSurface::NOFLIP)
2222  {
2223  allFaceInfo[i] = orientedSurface::FLIP;
2224  }
2225  else if (neiStatus[bFaceI] == orientedSurface::FLIP)
2226  {
2227  allFaceInfo[i] = orientedSurface::NOFLIP;
2228  }
2229  else
2230  {
2232  << "Incorrect status for face " << meshFaceI
2233  << abort(FatalError);
2234  }
2235  }
2236  }
2237  }
2238 
2239 
2240  // Convert to meshFlipMap and adapt faceZones
2241 
2242  meshFlipMap.setSize(mesh_.nFaces());
2243  meshFlipMap = false;
2244 
2245  forAll(allFaceInfo, faceI)
2246  {
2247  label meshFaceI = patch.addressing()[faceI];
2248 
2249  if (allFaceInfo[faceI] == orientedSurface::NOFLIP)
2250  {
2251  meshFlipMap[meshFaceI] = false;
2252  }
2253  else if (allFaceInfo[faceI] == orientedSurface::FLIP)
2254  {
2255  meshFlipMap[meshFaceI] = true;
2256  }
2257  else
2258  {
2260  << "Problem : unvisited face " << faceI
2261  << " centre:" << mesh_.faceCentres()[meshFaceI]
2262  << abort(FatalError);
2263  }
2264  }
2265 }
2266 
2267 
2268 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2269 
2272  const bool doHandleSnapProblems,
2273  const snapParameters& snapParams,
2274  const bool useTopologicalSnapDetection,
2275  const bool removeEdgeConnectedCells,
2276  const scalarField& perpendicularAngle,
2277  const bool mergeFreeStanding,
2278  const scalar planarAngle,
2279  const dictionary& motionDict,
2280  Time& runTime,
2281  const labelList& globalToMasterPatch,
2282  const labelList& globalToSlavePatch,
2283  const point& keepPoint
2284 )
2285 {
2286  // Introduce baffles
2287  // ~~~~~~~~~~~~~~~~~
2288 
2289  // Split the mesh along internal faces wherever there is a pierce between
2290  // two cell centres.
2291 
2292  Info<< "Introducing baffles for "
2294  << " faces that are intersected by the surface." << nl << endl;
2295 
2296  // Swap neighbouring cell centres and cell level
2297  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2298  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2299  calcNeighbourData(neiLevel, neiCc);
2300 
2301  labelList ownPatch, neiPatch;
2302  getBafflePatches
2303  (
2304  globalToMasterPatch,
2305  neiLevel,
2306  neiCc,
2307 
2308  ownPatch,
2309  neiPatch
2310  );
2311 
2312  createBaffles(ownPatch, neiPatch);
2313 
2314  if (debug)
2315  {
2316  // Debug:test all is still synced across proc patches
2317  checkData();
2318  }
2319 
2320  Info<< "Created baffles in = "
2321  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
2322 
2323  printMeshInfo(debug, "After introducing baffles");
2324 
2325  if (debug&MESH)
2326  {
2327  Pout<< "Writing baffled mesh to time " << timeName()
2328  << endl;
2329  write
2330  (
2331  debugType(debug),
2333  runTime.path()/"baffles"
2334  );
2335  Pout<< "Dumped debug data in = "
2336  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
2337  }
2338 
2339 
2340  // Introduce baffles to delete problem cells
2341  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2342 
2343  // Create some additional baffles where we want surface cells removed.
2344 
2345  if (doHandleSnapProblems)
2346  {
2347  handleSnapProblems
2348  (
2349  snapParams,
2350  useTopologicalSnapDetection,
2351  removeEdgeConnectedCells,
2352  perpendicularAngle,
2353  motionDict,
2354  runTime,
2355  globalToMasterPatch,
2356  globalToSlavePatch
2357  );
2358  }
2359 
2360 
2361  // Select part of mesh
2362  // ~~~~~~~~~~~~~~~~~~~
2363 
2364  Info<< nl
2365  << "Remove unreachable sections of mesh" << nl
2366  << "-----------------------------------" << nl
2367  << endl;
2368 
2369  if (debug)
2370  {
2371  runTime++;
2372  }
2373 
2374  splitMeshRegions(globalToMasterPatch, globalToSlavePatch, keepPoint);
2375 
2376  if (debug)
2377  {
2378  // Debug:test all is still synced across proc patches
2379  checkData();
2380  }
2381  Info<< "Split mesh in = "
2382  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
2383 
2384  printMeshInfo(debug, "After subsetting");
2385 
2386  if (debug&MESH)
2387  {
2388  Pout<< "Writing subsetted mesh to time " << timeName()
2389  << endl;
2390  write
2391  (
2392  debugType(debug),
2394  runTime.path()/timeName()
2395  );
2396  Pout<< "Dumped debug data in = "
2397  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
2398  }
2399 
2400 
2401  // Merge baffles
2402  // ~~~~~~~~~~~~~
2403 
2404  if (mergeFreeStanding)
2405  {
2406  Info<< nl
2407  << "Merge free-standing baffles" << nl
2408  << "---------------------------" << nl
2409  << endl;
2410 
2411 
2412  // List of pairs of freestanding baffle faces.
2413  List<labelPair> couples
2414  (
2415  freeStandingBaffles // filter out freestanding baffles
2416  (
2418  planarAngle
2419  )
2420  );
2421 
2422  label nCouples = couples.size();
2423  reduce(nCouples, sumOp<label>());
2424 
2425  Info<< "Detected free-standing baffles : " << nCouples << endl;
2426 
2427  if (nCouples > 0)
2428  {
2429  // Actually merge baffles. Note: not exactly parallellized. Should
2430  // convert baffle faces into processor faces if they resulted
2431  // from them.
2432  mergeBaffles(couples);
2433 
2434  // Detect any problem cells resulting from merging of baffles
2435  // and delete them
2436  handleSnapProblems
2437  (
2438  snapParams,
2439  useTopologicalSnapDetection,
2440  removeEdgeConnectedCells,
2441  perpendicularAngle,
2442  motionDict,
2443  runTime,
2444  globalToMasterPatch,
2445  globalToSlavePatch
2446  );
2447 
2448  if (debug)
2449  {
2450  // Debug:test all is still synced across proc patches
2451  checkData();
2452  }
2453  }
2454  Info<< "Merged free-standing baffles in = "
2455  << runTime.cpuTimeIncrement() << " s\n" << nl << endl;
2456  }
2457 }
2458 
2459 
2462  const label nBufferLayers,
2463  const labelList& globalToMasterPatch,
2464  const labelList& globalToSlavePatch,
2465  const point& keepPoint
2466 )
2467 {
2468  // Determine patches to put intersections into
2469  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2470 
2471  // Swap neighbouring cell centres and cell level
2472  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2473  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2474  calcNeighbourData(neiLevel, neiCc);
2475 
2476  labelList ownPatch, neiPatch;
2477  getBafflePatches
2478  (
2479  globalToMasterPatch,
2480  neiLevel,
2481  neiCc,
2482 
2483  ownPatch,
2484  neiPatch
2485  );
2486 
2487  // Analyse regions. Reuse regionsplit
2488  boolList blockedFace(mesh_.nFaces(), false);
2489 
2490  forAll(ownPatch, faceI)
2491  {
2492  if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
2493  {
2494  blockedFace[faceI] = true;
2495  }
2496  }
2497  syncTools::syncFaceList(mesh_, blockedFace, orEqOp<bool>());
2498 
2499  // Set region per cell based on walking
2500  regionSplit cellRegion(mesh_, blockedFace);
2501  blockedFace.clear();
2502 
2503  // Find the region containing the keepPoint
2504  const label keepRegionI = findRegion
2505  (
2506  mesh_,
2507  cellRegion,
2508  mergeDistance_*vector(1,1,1),
2509  keepPoint
2510  );
2511 
2512  Info<< "Found point " << keepPoint
2513  << " in global region " << keepRegionI
2514  << " out of " << cellRegion.nRegions() << " regions." << endl;
2515 
2516  if (keepRegionI == -1)
2517  {
2519  << "Point " << keepPoint
2520  << " is not inside the mesh." << nl
2521  << "Bounding box of the mesh:" << mesh_.bounds()
2522  << exit(FatalError);
2523  }
2524 
2525 
2526  // Walk out nBufferlayers from region split
2527  // (modifies cellRegion, ownPatch)
2528  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2529  // Takes over face patch onto points and then back to faces and cells
2530  // (so cell-face-point walk)
2531 
2532  const labelList& faceOwner = mesh_.faceOwner();
2533  const labelList& faceNeighbour = mesh_.faceNeighbour();
2534 
2535  // Patch for exposed faces for lack of anything sensible.
2536  label defaultPatch = 0;
2537  if (globalToMasterPatch.size())
2538  {
2539  defaultPatch = globalToMasterPatch[0];
2540  }
2541 
2542  for (label i = 0; i < nBufferLayers; i++)
2543  {
2544  // 1. From cells (via faces) to points
2545 
2546  labelList pointBaffle(mesh_.nPoints(), -1);
2547 
2548  forAll(faceNeighbour, faceI)
2549  {
2550  const face& f = mesh_.faces()[faceI];
2551 
2552  label ownRegion = cellRegion[faceOwner[faceI]];
2553  label neiRegion = cellRegion[faceNeighbour[faceI]];
2554 
2555  if (ownRegion == keepRegionI && neiRegion != keepRegionI)
2556  {
2557  // Note max(..) since possibly regionSplit might have split
2558  // off extra unreachable parts of mesh. Note: or can this only
2559  // happen for boundary faces?
2560  forAll(f, fp)
2561  {
2562  pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
2563  }
2564  }
2565  else if (ownRegion != keepRegionI && neiRegion == keepRegionI)
2566  {
2567  label newPatchI = neiPatch[faceI];
2568  if (newPatchI == -1)
2569  {
2570  newPatchI = max(defaultPatch, ownPatch[faceI]);
2571  }
2572  forAll(f, fp)
2573  {
2574  pointBaffle[f[fp]] = newPatchI;
2575  }
2576  }
2577  }
2578  for
2579  (
2580  label faceI = mesh_.nInternalFaces();
2581  faceI < mesh_.nFaces();
2582  faceI++
2583  )
2584  {
2585  const face& f = mesh_.faces()[faceI];
2586 
2587  label ownRegion = cellRegion[faceOwner[faceI]];
2588 
2589  if (ownRegion == keepRegionI)
2590  {
2591  forAll(f, fp)
2592  {
2593  pointBaffle[f[fp]] = max(defaultPatch, ownPatch[faceI]);
2594  }
2595  }
2596  }
2597 
2598  // Sync
2600  (
2601  mesh_,
2602  pointBaffle,
2603  maxEqOp<label>(),
2604  label(-1) // null value
2605  );
2606 
2607 
2608  // 2. From points back to faces
2609 
2610  const labelListList& pointFaces = mesh_.pointFaces();
2611 
2612  forAll(pointFaces, pointI)
2613  {
2614  if (pointBaffle[pointI] != -1)
2615  {
2616  const labelList& pFaces = pointFaces[pointI];
2617 
2618  forAll(pFaces, pFaceI)
2619  {
2620  label faceI = pFaces[pFaceI];
2621 
2622  if (ownPatch[faceI] == -1)
2623  {
2624  ownPatch[faceI] = pointBaffle[pointI];
2625  }
2626  }
2627  }
2628  }
2629  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
2630 
2631 
2632  // 3. From faces to cells (cellRegion) and back to faces (ownPatch)
2633 
2634  labelList newOwnPatch(ownPatch);
2635 
2636  forAll(ownPatch, faceI)
2637  {
2638  if (ownPatch[faceI] != -1)
2639  {
2640  label own = faceOwner[faceI];
2641 
2642  if (cellRegion[own] != keepRegionI)
2643  {
2644  cellRegion[own] = keepRegionI;
2645 
2646  const cell& ownFaces = mesh_.cells()[own];
2647  forAll(ownFaces, j)
2648  {
2649  if (ownPatch[ownFaces[j]] == -1)
2650  {
2651  newOwnPatch[ownFaces[j]] = ownPatch[faceI];
2652  }
2653  }
2654  }
2655  if (mesh_.isInternalFace(faceI))
2656  {
2657  label nei = faceNeighbour[faceI];
2658 
2659  if (cellRegion[nei] != keepRegionI)
2660  {
2661  cellRegion[nei] = keepRegionI;
2662 
2663  const cell& neiFaces = mesh_.cells()[nei];
2664  forAll(neiFaces, j)
2665  {
2666  if (ownPatch[neiFaces[j]] == -1)
2667  {
2668  newOwnPatch[neiFaces[j]] = ownPatch[faceI];
2669  }
2670  }
2671  }
2672  }
2673  }
2674  }
2675 
2676  ownPatch.transfer(newOwnPatch);
2677 
2678  syncTools::syncFaceList(mesh_, ownPatch, maxEqOp<label>());
2679  }
2680 
2681 
2682 
2683  // Subset
2684  // ~~~~~~
2685 
2686  // Get cells to remove
2687  DynamicList<label> cellsToRemove(mesh_.nCells());
2688  forAll(cellRegion, cellI)
2689  {
2690  if (cellRegion[cellI] != keepRegionI)
2691  {
2692  cellsToRemove.append(cellI);
2693  }
2694  }
2695  cellsToRemove.shrink();
2696 
2697  label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
2698  reduce(nCellsToKeep, sumOp<label>());
2699 
2700  Info<< "Keeping all cells in region " << keepRegionI
2701  << " containing point " << keepPoint << endl
2702  << "Selected for keeping : " << nCellsToKeep
2703  << " cells." << endl;
2704 
2705 
2706  // Remove cells
2707  removeCells cellRemover(mesh_);
2708 
2709  // Pick up patches for exposed faces
2710  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
2711  labelList exposedPatches(exposedFaces.size());
2712 
2713  forAll(exposedFaces, i)
2714  {
2715  label faceI = exposedFaces[i];
2716 
2717  if (ownPatch[faceI] != -1)
2718  {
2719  exposedPatches[i] = ownPatch[faceI];
2720  }
2721  else
2722  {
2724  << "For exposed face " << faceI
2725  << " fc:" << mesh_.faceCentres()[faceI]
2726  << " found no patch." << endl
2727  << " Taking patch " << defaultPatch
2728  << " instead." << endl;
2729  exposedPatches[i] = defaultPatch;
2730  }
2731  }
2732 
2733  return doRemoveCells
2734  (
2735  cellsToRemove,
2736  exposedFaces,
2737  exposedPatches,
2738  cellRemover
2739  );
2740 }
2741 
2742 
2745  const localPointRegion& regionSide
2746 )
2747 {
2748  // Topochange container
2749  polyTopoChange meshMod(mesh_);
2750 
2751  label nNonManifPoints = returnReduce
2752  (
2753  regionSide.meshPointMap().size(),
2754  sumOp<label>()
2755  );
2756 
2757  Info<< "dupNonManifoldPoints : Found : " << nNonManifPoints
2758  << " non-manifold points (out of "
2759  << mesh_.globalData().nTotalPoints()
2760  << ')' << endl;
2761 
2762  // Topo change engine
2763  duplicatePoints pointDuplicator(mesh_);
2764 
2765  // Insert changes into meshMod
2766  pointDuplicator.setRefinement(regionSide, meshMod);
2767 
2768  // Change the mesh (no inflation, parallel sync)
2769  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
2770 
2771  // Update fields
2772  mesh_.updateMesh(map);
2773 
2774  // Move mesh if in inflation mode
2775  if (map().hasMotionPoints())
2776  {
2777  mesh_.movePoints(map().preMotionPoints());
2778  }
2779  else
2780  {
2781  // Delete mesh volumes.
2782  mesh_.clearOut();
2783  }
2784 
2785  // Reset the instance for if in overwrite mode
2786  mesh_.setInstance(timeName());
2787 
2788  // Update intersections. Is mapping only (no faces created, positions stay
2789  // same) so no need to recalculate intersections.
2790  updateMesh(map, labelList(0));
2791 
2792  return map;
2793 }
2794 
2795 
2797 {
2798  // Analyse which points need to be duplicated
2799  localPointRegion regionSide(mesh_);
2800 
2801  return dupNonManifoldPoints(regionSide);
2802 }
2803 
2804 
2807  const point& keepPoint,
2808  const bool allowFreeStandingZoneFaces
2809 )
2810 {
2811  const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
2812 
2813  labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
2814 
2815  forAll(namedSurfaces, i)
2816  {
2817  label surfI = namedSurfaces[i];
2818 
2819  Info<< "Surface : " << surfaces_.names()[surfI] << nl
2820  << " faceZone : " << surfZones[surfI].faceZoneName() << nl
2821  << " cellZone : " << surfZones[surfI].cellZoneName() << endl;
2822  }
2823 
2824 
2825  // Add zones to mesh
2826  labelList surfaceToFaceZone =
2828  (
2829  surfZones,
2830  namedSurfaces,
2831  mesh_
2832  );
2833 
2834  labelList surfaceToCellZone =
2836  (
2837  surfZones,
2838  namedSurfaces,
2839  mesh_
2840  );
2841 
2842 
2843  const pointField& cellCentres = mesh_.cellCentres();
2844  const labelList& faceOwner = mesh_.faceOwner();
2845  const labelList& faceNeighbour = mesh_.faceNeighbour();
2846  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2847 
2848 
2849  // Swap neighbouring cell centres and cell level
2850  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2851  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2852  calcNeighbourData(neiLevel, neiCc);
2853 
2854 
2855  // Mark faces intersecting zoned surfaces
2856  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2857 
2858 
2859  // Like surfaceIndex_ but only for named surfaces.
2860  labelList namedSurfaceIndex(mesh_.nFaces(), -1);
2861  PackedBoolList posOrientation(mesh_.nFaces());
2862 
2863  {
2864  // Statistics: number of faces per faceZone
2865  labelList nSurfFaces(surfZones.size(), 0);
2866 
2867  // Note: for all internal faces? internal + coupled?
2868  // Since zonify is run after baffling the surfaceIndex_ on baffles is
2869  // not synchronised across both baffle faces. Fortunately we don't
2870  // do zonify baffle faces anyway (they are normal boundary faces).
2871 
2872  // Collect candidate faces
2873  // ~~~~~~~~~~~~~~~~~~~~~~~
2874 
2875  labelList testFaces(intersectedFaces());
2876 
2877  // Collect segments
2878  // ~~~~~~~~~~~~~~~~
2879 
2880  pointField start(testFaces.size());
2881  pointField end(testFaces.size());
2882 
2883  forAll(testFaces, i)
2884  {
2885  label faceI = testFaces[i];
2886 
2887  if (mesh_.isInternalFace(faceI))
2888  {
2889  start[i] = cellCentres[faceOwner[faceI]];
2890  end[i] = cellCentres[faceNeighbour[faceI]];
2891  }
2892  else
2893  {
2894  start[i] = cellCentres[faceOwner[faceI]];
2895  end[i] = neiCc[faceI-mesh_.nInternalFaces()];
2896  }
2897  }
2898 
2899  // Extend segments a bit
2900  {
2901  const vectorField smallVec(ROOTSMALL*(end-start));
2902  start -= smallVec;
2903  end += smallVec;
2904  }
2905 
2906 
2907  // Do test for intersections
2908  // ~~~~~~~~~~~~~~~~~~~~~~~~~
2909  // Note that we intersect all intersected faces again. Could reuse
2910  // the information already in surfaceIndex_.
2911 
2912  labelList surface1;
2913  List<pointIndexHit> hit1;
2914  vectorField normal1;
2915  labelList surface2;
2916  List<pointIndexHit> hit2;
2917  vectorField normal2;
2918  {
2919  labelList region1;
2920  labelList region2;
2921  surfaces_.findNearestIntersection
2922  (
2923  namedSurfaces,
2924  start,
2925  end,
2926 
2927  surface1,
2928  hit1,
2929  region1,
2930  normal1,
2931 
2932  surface2,
2933  hit2,
2934  region2,
2935  normal2
2936  );
2937  }
2938 
2939  forAll(testFaces, i)
2940  {
2941  label faceI = testFaces[i];
2942  const vector& area = mesh_.faceAreas()[faceI];
2943 
2944  if (surface1[i] != -1)
2945  {
2946  // If both hit should probably choose 'nearest'
2947  if
2948  (
2949  surface2[i] != -1
2950  && (
2951  magSqr(hit2[i].hitPoint())
2952  < magSqr(hit1[i].hitPoint())
2953  )
2954  )
2955  {
2956  namedSurfaceIndex[faceI] = surface2[i];
2957  posOrientation[faceI] = ((area&normal2[i]) > 0);
2958  nSurfFaces[surface2[i]]++;
2959  }
2960  else
2961  {
2962  namedSurfaceIndex[faceI] = surface1[i];
2963  posOrientation[faceI] = ((area&normal1[i]) > 0);
2964  nSurfFaces[surface1[i]]++;
2965  }
2966  }
2967  else if (surface2[i] != -1)
2968  {
2969  namedSurfaceIndex[faceI] = surface2[i];
2970  posOrientation[faceI] = ((area&normal2[i]) > 0);
2971  nSurfFaces[surface2[i]]++;
2972  }
2973  }
2974 
2975 
2976  // surfaceIndex might have different surfaces on both sides if
2977  // there happen to be a (obviously thin) surface with different
2978  // regions between the cell centres. If one is on a named surface
2979  // and the other is not this might give problems so sync.
2981  (
2982  mesh_,
2983  namedSurfaceIndex,
2984  maxEqOp<label>()
2985  );
2986 
2987  // Print a bit
2988  if (debug)
2989  {
2990  forAll(nSurfFaces, surfI)
2991  {
2992  Pout<< "Surface:"
2993  << surfaces_.names()[surfI]
2994  << " nZoneFaces:" << nSurfFaces[surfI] << nl;
2995  }
2996  Pout<< endl;
2997  }
2998  }
2999 
3000 
3001  // Put the cells into the correct zone
3002  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3003 
3004  // Zone per cell:
3005  // -2 : unset
3006  // -1 : not in any zone
3007  // >=0: zoneID
3008  labelList cellToZone(mesh_.nCells(), -2);
3009 
3010 
3011  // Set using geometric test
3012  // ~~~~~~~~~~~~~~~~~~~~~~~~
3013 
3014  // Closed surfaces with cellZone specified.
3015  labelList closedNamedSurfaces
3016  (
3018  (
3019  surfZones,
3020  surfaces_.geometry(),
3021  surfaces_.surfaces()
3022  )
3023  );
3024 
3025  if (closedNamedSurfaces.size())
3026  {
3027  Info<< "Found " << closedNamedSurfaces.size()
3028  << " closed, named surfaces. Assigning cells in/outside"
3029  << " these surfaces to the corresponding cellZone."
3030  << nl << endl;
3031 
3032  findCellZoneGeometric
3033  (
3034  neiCc,
3035  closedNamedSurfaces, // indices of closed surfaces
3036  namedSurfaceIndex, // per face index of named surface
3037  surfaceToCellZone, // cell zone index per surface
3038 
3039  cellToZone
3040  );
3041  }
3042 
3043 
3044  // Set using provided locations
3045  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3046 
3047  labelList locationSurfaces
3048  (
3050  );
3051 
3052  if (locationSurfaces.size())
3053  {
3054  Info<< "Found " << locationSurfaces.size()
3055  << " named surfaces with a provided inside point."
3056  << " Assigning cells inside these surfaces"
3057  << " to the corresponding cellZone."
3058  << nl << endl;
3059 
3060  findCellZoneInsideWalk
3061  (
3062  locationSurfaces, // indices of closed surfaces
3063  namedSurfaceIndex, // per face index of named surface
3064  surfaceToCellZone, // cell zone index per surface
3065 
3066  cellToZone
3067  );
3068  }
3069 
3070 
3071  // Set using walking
3072  // ~~~~~~~~~~~~~~~~~
3073 
3074  {
3075  Info<< "Walking from location-in-mesh " << keepPoint
3076  << " to assign cellZones "
3077  << "- crossing a faceZone face changes cellZone" << nl << endl;
3078 
3079  // Topological walk
3080  findCellZoneTopo
3081  (
3082  keepPoint,
3083  namedSurfaceIndex,
3084  surfaceToCellZone,
3085 
3086  cellToZone
3087  );
3088  }
3089 
3090 
3091  // Make sure namedSurfaceIndex is unset inbetween same cell cell zones.
3092  if (!allowFreeStandingZoneFaces)
3093  {
3094  Info<< "Only keeping zone faces inbetween different cellZones."
3095  << nl << endl;
3096 
3097  makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
3098  }
3099 
3100 
3101  //- Per face index of faceZone or -1
3102  labelList faceToZone(mesh_.nFaces(), -1);
3103 
3104  // Convert namedSurfaceIndex (index of named surfaces) to
3105  // actual faceZone index
3106 
3107  forAll(namedSurfaceIndex, faceI)
3108  {
3109  label surfI = namedSurfaceIndex[faceI];
3110  if (surfI != -1)
3111  {
3112  faceToZone[faceI] = surfaceToFaceZone[surfI];
3113  }
3114  }
3115 
3116 
3117  // Topochange container
3118  polyTopoChange meshMod(mesh_);
3119 
3120 
3121 
3122  // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
3123  labelList neiCellZone;
3124  syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
3125  forAll(patches, patchI)
3126  {
3127  const polyPatch& pp = patches[patchI];
3128 
3129  if (!pp.coupled())
3130  {
3131  label bFaceI = pp.start()-mesh_.nInternalFaces();
3132  forAll(pp, i)
3133  {
3134  neiCellZone[bFaceI++] = -1;
3135  }
3136  }
3137  }
3138 
3139 
3140 
3141  // Get per face whether is it master (of a coupled set of faces)
3142  const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
3143 
3144 
3145 
3146  // faceZones
3147  // ~~~~~~~~~
3148  // Faces on faceZones come in two variants:
3149  // - faces on the outside of a cellZone. They will be oriented to
3150  // point out of the maximum cellZone.
3151  // - free-standing faces. These will be oriented according to the
3152  // local surface normal. We do this in a two step algorithm:
3153  // - do a consistent orientation
3154  // - check number of faces with consistent orientation
3155  // - if <0 flip the whole patch
3156  boolList meshFlipMap(mesh_.nFaces(), false);
3157  {
3158  // Collect all data on zone faces without cellZones on either side.
3159  const indirectPrimitivePatch patch
3160  (
3162  (
3163  mesh_.faces(),
3164  freeStandingBaffleFaces
3165  (
3166  faceToZone,
3167  cellToZone,
3168  neiCellZone
3169  )
3170  ),
3171  mesh_.points()
3172  );
3173 
3174  label nFreeStanding = returnReduce(patch.size(), sumOp<label>());
3175  if (nFreeStanding > 0)
3176  {
3177  Info<< "Detected " << nFreeStanding << " free-standing zone faces"
3178  << endl;
3179 
3180  if (debug)
3181  {
3182  OBJstream str(mesh_.time().path()/"freeStanding.obj");
3183  str.write(patch.localFaces(), patch.localPoints(), false);
3184  }
3185 
3186 
3187  // Detect non-manifold edges
3188  labelList nMasterFacesPerEdge;
3189  calcPatchNumMasterFaces(isMasterFace, patch, nMasterFacesPerEdge);
3190 
3191  // Mark zones. Even a single original surface might create multiple
3192  // disconnected/non-manifold-connected zones
3193  labelList faceToConnectedZone;
3194  const label nZones = markPatchZones
3195  (
3196  patch,
3197  nMasterFacesPerEdge,
3198  faceToConnectedZone
3199  );
3200 
3201  Map<label> nPosOrientation(2*nZones);
3202  for (label zoneI = 0; zoneI < nZones; zoneI++)
3203  {
3204  nPosOrientation.insert(zoneI, 0);
3205  }
3206 
3207  // Make orientations consistent in a topological way. This just
3208  // checks the first face per zone for whether nPosOrientation
3209  // is negative (which it never is at this point)
3210  consistentOrientation
3211  (
3212  isMasterFace,
3213  patch,
3214  nMasterFacesPerEdge,
3215  faceToConnectedZone,
3216  nPosOrientation,
3217 
3218  meshFlipMap
3219  );
3220 
3221  // Count per region the number of orientations (taking the new
3222  // flipMap into account)
3223  forAll(patch.addressing(), faceI)
3224  {
3225  label meshFaceI = patch.addressing()[faceI];
3226 
3227  if (isMasterFace[meshFaceI])
3228  {
3229  label n = 1;
3230  if
3231  (
3232  bool(posOrientation[meshFaceI])
3233  == meshFlipMap[meshFaceI]
3234  )
3235  {
3236  n = -1;
3237  }
3238 
3239  nPosOrientation.find(faceToConnectedZone[faceI])() += n;
3240  }
3241  }
3242  Pstream::mapCombineGather(nPosOrientation, plusEqOp<label>());
3243  Pstream::mapCombineScatter(nPosOrientation);
3244 
3245 
3246  Info<< "Split " << nFreeStanding << " free-standing zone faces"
3247  << " into " << nZones << " disconnected regions with size"
3248  << " (negative denotes wrong orientation) :"
3249  << endl;
3250 
3251  for (label zoneI = 0; zoneI < nZones; zoneI++)
3252  {
3253  Info<< " " << zoneI << "\t" << nPosOrientation[zoneI]
3254  << endl;
3255  }
3256  Info<< endl;
3257 
3258 
3259  // Reapply with new counts (in nPosOrientation). This will cause
3260  // zones with a negative count to be flipped.
3261  consistentOrientation
3262  (
3263  isMasterFace,
3264  patch,
3265  nMasterFacesPerEdge,
3266  faceToConnectedZone,
3267  nPosOrientation,
3268 
3269  meshFlipMap
3270  );
3271  }
3272  }
3273 
3274 
3275  // Put the faces into the correct zone
3276  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3277 
3278  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
3279  {
3280  label faceZoneI = faceToZone[faceI];
3281 
3282  if (faceZoneI != -1)
3283  {
3284  // Orient face zone to have slave cells in max cell zone.
3285  // Note: logic to use flipMap should be consistent with logic
3286  // to pick up the freeStandingBaffleFaces!
3287 
3288  label ownZone = cellToZone[faceOwner[faceI]];
3289  label neiZone = cellToZone[faceNeighbour[faceI]];
3290 
3291  bool flip;
3292 
3293  if (ownZone == neiZone)
3294  {
3295  // free-standing face. Use geometrically derived orientation
3296  flip = meshFlipMap[faceI];
3297  }
3298  else
3299  {
3300  flip =
3301  (
3302  ownZone == -1
3303  || (neiZone != -1 && ownZone > neiZone)
3304  );
3305  }
3306 
3307  meshMod.setAction
3308  (
3310  (
3311  mesh_.faces()[faceI], // modified face
3312  faceI, // label of face
3313  faceOwner[faceI], // owner
3314  faceNeighbour[faceI], // neighbour
3315  false, // face flip
3316  -1, // patch for face
3317  false, // remove from zone
3318  faceZoneI, // zone for face
3319  flip // face flip in zone
3320  )
3321  );
3322  }
3323  }
3324 
3325 
3326  // Set owner as no-flip
3327  forAll(patches, patchI)
3328  {
3329  const polyPatch& pp = patches[patchI];
3330 
3331  label faceI = pp.start();
3332 
3333  forAll(pp, i)
3334  {
3335  label faceZoneI = faceToZone[faceI];
3336 
3337  if (faceZoneI != -1)
3338  {
3339  label ownZone = cellToZone[faceOwner[faceI]];
3340  label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
3341 
3342  bool flip;
3343 
3344  if (ownZone == neiZone)
3345  {
3346  // free-standing face. Use geometrically derived orientation
3347  flip = meshFlipMap[faceI];
3348  }
3349  else
3350  {
3351  flip =
3352  (
3353  ownZone == -1
3354  || (neiZone != -1 && ownZone > neiZone)
3355  );
3356  }
3357 
3358  meshMod.setAction
3359  (
3361  (
3362  mesh_.faces()[faceI], // modified face
3363  faceI, // label of face
3364  faceOwner[faceI], // owner
3365  -1, // neighbour
3366  false, // face flip
3367  patchI, // patch for face
3368  false, // remove from zone
3369  faceZoneI, // zone for face
3370  flip // face flip in zone
3371  )
3372  );
3373  }
3374  faceI++;
3375  }
3376  }
3377 
3378 
3379  // Put the cells into the correct zone
3380  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
3381 
3382  forAll(cellToZone, cellI)
3383  {
3384  label zoneI = cellToZone[cellI];
3385 
3386  if (zoneI >= 0)
3387  {
3388  meshMod.setAction
3389  (
3391  (
3392  cellI,
3393  false, // removeFromZone
3394  zoneI
3395  )
3396  );
3397  }
3398  }
3399 
3400  // Change the mesh (no inflation, parallel sync)
3401  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
3402 
3403  // Update fields
3404  mesh_.updateMesh(map);
3405 
3406  // Move mesh if in inflation mode
3407  if (map().hasMotionPoints())
3408  {
3409  mesh_.movePoints(map().preMotionPoints());
3410  }
3411  else
3412  {
3413  // Delete mesh volumes.
3414  mesh_.clearOut();
3415  }
3416 
3417  // Reset the instance for if in overwrite mode
3418  mesh_.setInstance(timeName());
3419 
3420  // Print some stats (note: zones are synchronised)
3421  if (mesh_.cellZones().size() > 0)
3422  {
3423  Info<< "CellZones:" << endl;
3424  forAll(mesh_.cellZones(), zoneI)
3425  {
3426  const cellZone& cz = mesh_.cellZones()[zoneI];
3427  Info<< " " << cz.name()
3428  << "\tsize:" << returnReduce(cz.size(), sumOp<label>())
3429  << endl;
3430  }
3431  Info<< endl;
3432  }
3433  if (mesh_.faceZones().size() > 0)
3434  {
3435  Info<< "FaceZones:" << endl;
3436  forAll(mesh_.faceZones(), zoneI)
3437  {
3438  const faceZone& fz = mesh_.faceZones()[zoneI];
3439  Info<< " " << fz.name()
3440  << "\tsize:" << returnReduce(fz.size(), sumOp<label>())
3441  << endl;
3442  }
3443  Info<< endl;
3444  }
3445 
3446  // None of the faces has changed, only the zones. Still...
3447  updateMesh(map, labelList());
3448 
3449  return map;
3450 }
3451 
3452 
3453 // ************************************************************************* //
const vectorField & faceAreas() const
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTime.C:74
labelList getExposedFaces(const labelList &cellsToRemove) const
Get labels of exposed faces.
Definition: removeCells.C:73
const cellZoneMesh & cellZones() const
Return cell zone mesh.
Definition: polyMesh.H:469
A List with indirect addressing.
Definition: IndirectList.H:102
void calc(const argList &args, const Time &runTime, const fvMesh &mesh)
static labelList getInsidePointNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of surfaces with a cellZone that have &#39;insidePoint&#39;.
Transport of region for use in PatchEdgeFaceWave.
OFstream which keeps track of vertices.
Definition: OBJstream.H:53
label countHits() const
Count number of intersections (local)
void findNearestIntersection(const labelList &surfacesToTest, const pointField &start, const pointField &end, labelList &surface1, List< pointIndexHit > &hit1, labelList &region1, labelList &surface2, List< pointIndexHit > &hit2, labelList &region2) const
Find intersection nearest to the endpoints. surface1,2 are.
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:249
label setAction(const topoAction &action)
For compatibility with polyTopoChange: set topological action.
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:463
autoPtr< mapPolyMesh > createBaffles(const labelList &ownPatch, const labelList &neiPatch)
Create baffle for every internal face where ownPatch != -1.
static PackedBoolList getMasterFaces(const polyMesh &)
Get per face whether it is uncoupled or a master of a.
Definition: syncTools.C:156
bool write() const
Write mesh and all data.
label size() const
Return the number of elements in the PtrList.
Definition: PtrListI.H:32
label nFaces() const
const searchableSurfaces & geometry() const
static label findRegion(const polyMesh &, const labelList &cellRegion, const vector &perturbVec, const point &p)
Find region point is in. Uses optional perturbation to re-test.
A list of face labels.
Definition: faceSet.H:48
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
void baffleAndSplitMesh(const bool handleSnapProblems, const snapParameters &snapParams, const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField &perpendicularAngle, const bool mergeFreeStanding, const scalar freeStandingAngle, const dictionary &motionDict, Time &runTime, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const point &keepPoint)
Split off unreachable areas of mesh.
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
void setInstance(const fileName &)
Set the instance for mesh files.
Definition: polyMeshIO.C:32
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
fileName path() const
Return path.
Definition: Time.H:281
const word & name() const
Return name.
labelList meshEdges(const edgeList &allEdges, const labelListList &cellEdges, const labelList &faceCells) const
Return labels of patch edges in the global edge list using.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:61
void checkZoneFaces() const
Debug helper: check faceZones are not on processor patches.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:390
label whichZone(const label objectIndex) const
Given a global object index, return the zone it is in.
Definition: ZoneMesh.C:232
dimensioned< scalar > mag(const dimensioned< Type > &)
An STL-conforming const_iterator.
Definition: HashTable.H:470
A subset of mesh cells.
Definition: cellZone.H:61
const edgeList & edges() const
Return mesh edges. Uses calcEdges.
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
const labelListList & pointEdges() const
const labelList & tetBasePtIs() const
Return the tetBasePtIs.
Definition: polyMesh.C:864
void flip()
Reverse orientation.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: fvMesh.C:728
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
label nRegions() const
Return total number of regions.
Definition: regionSplit.H:218
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
autoPtr< mapPolyMesh > splitMesh(const label nBufferLayers, const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const point &keepPoint)
Split off (with optional buffer layers) unreachable areas.
scalar f1
Definition: createFields.H:28
const Type & second() const
Return second.
Definition: Pair.H:99
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
A bit-packed bool list.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void updateMesh(const mapPolyMesh &, const labelList &changedFaces)
Update for external change to mesh. changedFaces are in new mesh.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
Given list of cells to remove insert all the topology changes.
Definition: removeCells.H:59
void setSize(const label)
Dummy setSize function.
Definition: FixedListI.H:177
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static labelList getNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of named surfaces (surfaces with faceZoneName)
label nCells() const
const cellList & cells() const
autoPtr< mapPolyMesh > splitMeshRegions(const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const point &keepPoint)
Split mesh. Keep part containing point.
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
label nEdges() const
const vectorField & cellCentres() const
messageStream Info
void printMeshInfo(const bool, const string &) const
Print some mesh stats.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
bool updateEdge(const polyMesh &mesh, const indirectPrimitivePatch &patch, const label edgeI, const label faceI, const patchFaceOrientation &faceInfo, const scalar tol, TrackingData &td)
Influence of face on edge.
Simple container to keep together snap specific information.
static labelList getClosedNamedSurfaces(const PtrList< surfaceZonesInfo > &surfList, const searchableSurfaces &allGeometry, const labelList &surfaces)
Get indices of surfaces with a cellZone that are closed and.
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:380
T & first()
Return the first element of the list.
Definition: UListI.H:117
const labelListList & faceEdges() const
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
virtual void updateMesh(const mapPolyMesh &mpm)
Update mesh corresponding to the given map.
Definition: fvMesh.C:804
Class containing data for face removal.
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:317
patches[0]
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
const wordList & names() const
Names of surfaces.
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
label n
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
void findInside(const labelList &surfacesToTest, const pointField &pt, labelList &insideSurfaces) const
Detect if a point is &#39;inside&#39; (closed) surfaces.
static labelList addFaceZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:68
static void mapCombineScatter(const List< commsStruct > &comms, Container &Values, const int tag, const label comm)
Scatter data. Reverse of combineGather.
void setRefinement(const localPointRegion &regionSide, polyTopoChange &)
Play commands into polyTopoChange to duplicate points. Gets.
const Field< PointType > & localPoints() const
Return pointField of points in patch.
static const char nl
Definition: Ostream.H:260
Wave propagation of information along patch. Every iteration information goes through one layer of fa...
const double e
Elementary charge.
Definition: doubleFloat.H:78
#define WarningInFunction
Report a warning using Foam::Warning.
const Field< PointType > & faceCentres() const
Return face centres for patch.
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
virtual void sync(const polyMesh &mesh)
Sync faceSet across coupled patches.
Definition: faceSet.C:113
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1035
Foam::polyBoundaryMesh.
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#define forAll(list, i)
Definition: UList.H:421
virtual Ostream & write(const char)
Write character.
Definition: OBJstream.C:82
autoPtr< mapPolyMesh > createZoneBaffles(const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, List< labelPair > &)
Create baffles for faces straddling zoned surfaces. Return.
virtual Ostream & write(const token &)=0
Write next token to stream.
dimensionedScalar cos(const dimensionedScalar &ds)
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:429
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
Unit conversion functions.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
Class describing modification of a face.
static labelList addCellZonesToMesh(const PtrList< surfaceZonesInfo > &surfList, const labelList &namedSurfaces, polyMesh &mesh)
label size() const
Return number of elements in table.
List< Key > toc() const
Return the table of contents.
Definition: HashTable.C:201
static void mapCombineGather(const List< commsStruct > &comms, Container &Values, const CombineOp &cop, const int tag, const label comm)
const labelListList & pointFaces() const
errorManip< error > abort(error &err)
Definition: errorManip.H:131
const word & name() const
Return name.
Definition: IOobject.H:260
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1073
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:404
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
const word & name() const
Return name.
Definition: zone.H:150
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:88
label nTotalPoints() const
Return total number of points in decomposed mesh. Not.
Transport of orientation for use in PatchEdgeFaceWave.
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
A list of faces which address into the list of points.
static writeType writeLevel()
Get/set write level.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
label nInternalFaces() const
labelList intersectedFaces() const
Get faces with intersection.
const Map< label > & meshPointMap() const
Per point that is to be duplicated the local index.
error FatalError
const PtrList< surfaceZonesInfo > & surfZones() const
Duplicate points.
const boundBox & bounds() const
Return mesh bounding box.
Definition: polyMesh.H:427
ZoneMesh< faceZone, polyMesh > faceZoneMesh
A ZoneMesh with the type faceZone.
List< label > labelList
A List of labels.
Definition: labelList.H:56
void clearOut()
Clear all geometry and addressing.
Definition: fvMesh.C:231
autoPtr< mapPolyMesh > dupNonManifoldPoints()
Find boundary points that connect to more than one cell.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
const labelListList & faceEdges() const
Return face-edge addressing.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
const Type & first() const
Return first.
Definition: Pair.H:87
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
void reverse(UList< T > &, const label n)
Definition: UListI.H:322
Direct mesh changes based on v1.3 polyTopoChange syntax.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:106
const List< Face > & localFaces() const
Return patch faces addressing into local point list.
Class describing modification of a cell.
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:97
word timeName() const
Replacement for Time::timeName() : return oldInstance (if.
void checkData()
Debugging: check that all faces still obey start()>end()
Takes mesh with &#39;baffles&#39; (= boundary faces sharing points). Determines for selected points on bounda...
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1060
iterator find(const label &)
Find and return an iterator set at the hashedEntry.
bool updateEdge(const polyMesh &mesh, const indirectPrimitivePatch &patch, const label edgeI, const label faceI, const patchEdgeFaceRegion &faceInfo, const scalar tol, TrackingData &td)
Influence of face on edge.
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){pointMap[i]=i;}for(label i=0;i< nPoints;i++){if(f[i] > 0.0){hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei){if(edges[ei].mag(points)< SMALL){label start=pointMap[edges[ei].start()];while(start!=pointMap[start]){start=pointMap[start];}label end=pointMap[edges[ei].end()];while(end!=pointMap[end]){end=pointMap[end];}label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;}}cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){cellShape &cs=cellShapes[celli];forAll(cs, i){cs[i]=pointMap[cs[i]];}cs.collapse();}label bcIDs[11]={-1, 0, 2, 4,-1, 5,-1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&symmetryPolyPatch::typeName,&wedgePolyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&symmetryPolyPatch::typeName,&oldCyclicPolyPatch::typeName};enum patchTypeNames{PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={"piston","valve","liner","cylinderHead","axis","wedge","inflow","outflow","presin","presout","symmetryPlane","cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
const labelList & surfaces() const
bool mkDir(const fileName &, mode_t=0777)
Make a directory and return an error if it could not be created.
Definition: POSIX.C:420
autoPtr< mapPolyMesh > mergeBaffles(const List< labelPair > &)
Merge baffles. Gets pairs of faces.
label nEdges() const
Return number of edges in patch.
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1079
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:209
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:117
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition: regionSplit.H:116
label nPoints() const
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1200
autoPtr< mapPolyMesh > zonify(const point &keepPoint, const bool allowFreeStandingZoneFaces)
Put faces/cells into zones according to surface specification.
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
const vectorField & faceCentres() const
static const label labelMax
Definition: label.H:62
bool insert(const Key &key)
Insert a new entry.
Definition: HashSet.H:116