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