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