meshRefinementBlock.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2018-2024 OpenCFD Ltd.
9 -------------------------------------------------------------------------------
10 License
11  This file is part of OpenFOAM.
12 
13  OpenFOAM is free software: you can redistribute it and/or modify it
14  under the terms of the GNU General Public License as published by
15  the Free Software Foundation, either version 3 of the License, or
16  (at your option) any later version.
17 
18  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21  for more details.
22 
23  You should have received a copy of the GNU General Public License
24  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25 
26 \*---------------------------------------------------------------------------*/
27 
28 #include "meshRefinement.H"
29 #include "fvMesh.H"
30 #include "Time.H"
31 #include "refinementSurfaces.H"
32 #include "removeCells.H"
33 #include "unitConversion.H"
34 #include "bitSet.H"
35 #include "volFields.H"
36 
37 // Leak path
38 #include "shortestPathSet.H"
39 #include "meshSearch.H"
40 #include "topoDistanceData.H"
41 #include "FaceCellWave.H"
42 #include "removeCells.H"
43 #include "regionSplit.H"
44 
45 #include "volFields.H"
46 #include "wallPoints.H"
47 #include "searchableSurfaces.H"
48 #include "distributedTriSurfaceMesh.H"
49 
50 #include "holeToFace.H"
51 #include "refinementParameters.H"
52 #include "indirectPrimitivePatch.H"
53 #include "OBJstream.H"
54 #include "PatchTools.H"
55 
56 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
57 
58 //Foam::label Foam::meshRefinement::markFakeGapRefinement
59 //(
60 // const scalar planarCos,
61 //
62 // const label nAllowRefine,
63 // const labelList& neiLevel,
64 // const pointField& neiCc,
65 //
66 // labelList& refineCell,
67 // label& nRefine
68 //) const
69 //{
70 // label oldNRefine = nRefine;
71 //
72 //
73 // // Collect candidate faces (i.e. intersecting any surface and
74 // // owner/neighbour not yet refined.
75 // const labelList testFaces(getRefineCandidateFaces(refineCell));
76 //
77 // // Collect segments
78 // pointField start(testFaces.size());
79 // pointField end(testFaces.size());
80 // labelList minLevel(testFaces.size());
81 //
82 // calcCellCellRays
83 // (
84 // neiCc,
85 // neiLevel,
86 // testFaces,
87 // start,
88 // end,
89 // minLevel
90 // );
91 //
92 //
93 // // Re-use the gap shooting methods. This needs:
94 // // - shell gapLevel : faked. Set to 0,labelMax
95 // // - surface gapLevel : faked by overwriting
96 //
97 //
98 // List<FixedList<label, 3>>& surfGapLevel = const_cast
99 // <
100 // List<FixedList<label, 3>>&
101 // >(surfaces_.extendedGapLevel());
102 //
103 // List<volumeType>& surfGapMode = const_cast
104 // <
105 // List<volumeType>&
106 // >(surfaces_.extendedGapMode());
107 //
108 // const List<FixedList<label, 3>> surfOldLevel(surfGapLevel);
109 // const List<volumeType> surfOldMode(surfGapMode);
110 //
111 // // Set the extended gap levels
112 // forAll(surfaces_.gapLevel(), regioni)
113 // {
114 // surfGapLevel[regioni] = FixedList<label, 3>
115 // ({
116 // 3,
117 // -1,
118 // surfaces_.gapLevel()[regioni]+1
119 // });
120 // }
121 // surfGapMode = volumeType::MIXED;
122 //
123 //Pout<< "gapLevel was:" << surfOldLevel << endl;
124 //Pout<< "gapLevel now:" << surfGapLevel << endl;
125 //Pout<< "gapMode was:" << surfOldMode << endl;
126 //Pout<< "gapMode now:" << surfGapMode << endl;
127 //Pout<< "nRefine was:" << oldNRefine << endl;
128 //
129 //
130 //
131 // List<List<FixedList<label, 3>>>& shellGapLevel = const_cast
132 // <
133 // List<List<FixedList<label, 3>>>&
134 // >(shells_.extendedGapLevel());
135 //
136 // List<List<volumeType>>& shellGapMode = const_cast
137 // <
138 // List<List<volumeType>>&
139 // >(shells_.extendedGapMode());
140 //
141 // const List<List<FixedList<label, 3>>> shellOldLevel(shellGapLevel);
142 // const List<List<volumeType>> shellOldMode(shellGapMode);
143 //
144 // // Set the extended gap levels
145 // forAll(shellGapLevel, shelli)
146 // {
147 // shellGapLevel[shelli] = FixedList<label, 3>({3, -1, labelMax});
148 // shellGapMode[shelli] = volumeType::MIXED;
149 // }
150 //Pout<< "shellLevel was:" << shellOldLevel << endl;
151 //Pout<< "shellLevel now:" << shellGapLevel << endl;
152 //
153 // const label nAdditionalRefined = markSurfaceGapRefinement
154 // (
155 // planarCos,
156 //
157 // nAllowRefine,
158 // neiLevel,
159 // neiCc,
160 //
161 // refineCell,
162 // nRefine
163 // );
164 //
165 //Pout<< "nRefine now:" << nRefine << endl;
166 //
167 // // Restore
168 // surfGapLevel = surfOldLevel;
169 // surfGapMode = surfOldMode;
170 // shellGapLevel = shellOldLevel;
171 // shellGapMode = shellOldMode;
172 //
173 // return nAdditionalRefined;
174 //}
175 
176 
177 void Foam::meshRefinement::markOutsideFaces
178 (
179  const labelList& cellLevel,
180  const labelList& neiLevel,
181  const labelList& refineCell,
182  bitSet& isOutsideFace
183 ) const
184 {
185  // Get faces:
186  // - on outside of cell set
187  // - inbetween same cell level (i.e. quads)
188 
189  isOutsideFace.setSize(mesh_.nFaces());
190  isOutsideFace = Zero;
191 
192  forAll(mesh_.faceNeighbour(), facei)
193  {
194  label own = mesh_.faceOwner()[facei];
195  label nei = mesh_.faceNeighbour()[facei];
196  if
197  (
198  (cellLevel[own] == cellLevel[nei])
199  && (
200  (refineCell[own] != -1)
201  != (refineCell[nei] != -1)
202  )
203  )
204  {
205  isOutsideFace.set(facei);
206  }
207  }
208  {
209 
210  const label nBnd = mesh_.nBoundaryFaces();
211 
212  labelList neiRefineCell(nBnd);
213  syncTools::swapBoundaryCellList(mesh_, refineCell, neiRefineCell);
214  for (label bFacei = 0; bFacei < nBnd; ++bFacei)
215  {
216  label facei = mesh_.nInternalFaces()+bFacei;
217  label own = mesh_.faceOwner()[facei];
218 
219  if
220  (
221  (cellLevel[own] == neiLevel[bFacei])
222  && (
223  (refineCell[own] != -1)
224  != (neiRefineCell[bFacei] != -1)
225  )
226  )
227  {
228  isOutsideFace.set(facei);
229  }
230  }
231  }
232 }
233 
234 
235 Foam::label Foam::meshRefinement::countFaceDirs
236 (
237  const bitSet& isOutsideFace,
238  const label celli
239 ) const
240 {
241  const cell& cFaces = mesh_.cells()[celli];
242  const vectorField& faceAreas = mesh_.faceAreas();
243 
244  Vector<bool> haveDirs(vector::uniform(false));
245 
246  forAll(cFaces, cFacei)
247  {
248  const label facei = cFaces[cFacei];
249 
250  if (isOutsideFace[facei])
251  {
252  const vector& n = faceAreas[facei];
253  scalar magSqrN = magSqr(n);
254 
255  if (magSqrN > ROOTVSMALL)
256  {
257  for
258  (
259  direction dir = 0;
260  dir < pTraits<vector>::nComponents;
261  dir++
262  )
263  {
264  if (Foam::sqr(n[dir]) > 0.99*magSqrN)
265  {
266  haveDirs[dir] = true;
267  break;
268  }
269  }
270  }
271  }
272  }
273 
274  label nDirs = 0;
275  forAll(haveDirs, dir)
276  {
277  if (haveDirs[dir])
278  {
279  nDirs++;
280  }
281  }
282  return nDirs;
283 }
284 
285 
286 void Foam::meshRefinement::growSet
287 (
288  const labelList& neiLevel,
289  const bitSet& isOutsideFace,
290  labelList& refineCell,
291  label& nRefine
292 ) const
293 {
294  // Get cells with three or more outside faces
295  const cellList& cells = mesh_.cells();
296  forAll(cells, celli)
297  {
298  if (refineCell[celli] == -1)
299  {
300  if (countFaceDirs(isOutsideFace, celli) == 3)
301  {
302  // Mark cell with any value
303  refineCell[celli] = 0;
304  nRefine++;
305  }
306  }
307  }
308 }
309 
310 
311 //void Foam::meshRefinement::markMultiRegionCell
312 //(
313 // const label celli,
314 // const FixedList<label, 3>& surface,
315 //
316 // Map<FixedList<label, 3>>& cellToRegions,
317 // bitSet& isMultiRegion
318 //) const
319 //{
320 // if (!isMultiRegion[celli])
321 // {
322 // Map<FixedList<label, 3>>::iterator fnd = cellToRegions.find(celli);
323 // if (!fnd.good())
324 // {
325 // cellToRegions.insert(celli, surface);
326 // }
327 // else if (fnd() != surface)
328 // {
329 // // Found multiple intersections on cell
330 // isMultiRegion.set(celli);
331 // }
332 // }
333 //}
334 
335 
336 //void Foam::meshRefinement::detectMultiRegionCells
337 //(
338 // const labelListList& faceZones,
339 // const labelList& testFaces,
340 //
341 // const labelList& surface1,
342 // const List<pointIndexHit>& hit1,
343 // const labelList& region1,
344 //
345 // const labelList& surface2,
346 // const List<pointIndexHit>& hit2,
347 // const labelList& region2,
348 //
349 // bitSet& isMultiRegion
350 //) const
351 //{
352 // isMultiRegion.clear();
353 // isMultiRegion.setSize(mesh_.nCells());
354 //
355 // Map<FixedList<label, 3>> cellToRegions(testFaces.size());
356 //
357 // forAll(testFaces, i)
358 // {
359 // const pointIndexHit& info1 = hit1[i];
360 // if (info1.hit())
361 // {
362 // const label facei = testFaces[i];
363 // const labelList& fz1 = faceZones[surface1[i]];
364 // const FixedList<label, 3> surfaceInfo1
365 // ({
366 // surface1[i],
367 // region1[i],
368 // (fz1.size() ? fz1[info1.index()] : region1[i])
369 // });
370 //
371 // markMultiRegionCell
372 // (
373 // mesh_.faceOwner()[facei],
374 // surfaceInfo1,
375 // cellToRegions,
376 // isMultiRegion
377 // );
378 // if (mesh_.isInternalFace(facei))
379 // {
380 // markMultiRegionCell
381 // (
382 // mesh_.faceNeighbour()[facei],
383 // surfaceInfo1,
384 // cellToRegions,
385 // isMultiRegion
386 // );
387 // }
388 //
389 // const pointIndexHit& info2 = hit2[i];
390 //
391 // if (info2.hit() && info1 != info2)
392 // {
393 // const labelList& fz2 = faceZones[surface2[i]];
394 // const FixedList<label, 3> surfaceInfo2
395 // ({
396 // surface2[i],
397 // region2[i],
398 // (fz2.size() ? fz2[info2.index()] : region2[i])
399 // });
400 //
401 // markMultiRegionCell
402 // (
403 // mesh_.faceOwner()[facei],
404 // surfaceInfo2,
405 // cellToRegions,
406 // isMultiRegion
407 // );
408 // if (mesh_.isInternalFace(facei))
409 // {
410 // markMultiRegionCell
411 // (
412 // mesh_.faceNeighbour()[facei],
413 // surfaceInfo2,
414 // cellToRegions,
415 // isMultiRegion
416 // );
417 // }
418 // }
419 // }
420 // }
421 //
422 //
423 // if (debug&meshRefinement::MESH)
424 // {
425 // volScalarField multiCell
426 // (
427 // IOobject
428 // (
429 // "multiCell",
430 // mesh_.time().timeName(),
431 // mesh_,
432 // IOobject::NO_READ,
433 // IOobject::NO_WRITE,
434 // IOobject::NO_REGISTER
435 // ),
436 // mesh_,
437 // dimensionedScalar
438 // (
439 // "zero",
440 // dimensionSet(0, 1, 0, 0, 0),
441 // 0.0
442 // )
443 // );
444 // forAll(isMultiRegion, celli)
445 // {
446 // if (isMultiRegion[celli])
447 // {
448 // multiCell[celli] = 1.0;
449 // }
450 // }
451 //
452 // Info<< "Writing all multi cells to " << multiCell.name() << endl;
453 // multiCell.write();
454 // }
455 //}
456 
457 
458 void Foam::meshRefinement::markProximityCandidateFaces
459 (
460  const labelList& blockedSurfaces,
461  const List<scalarList>& regionToBlockSize,
462  const labelList& neiLevel,
463  const pointField& neiCc,
464 
465  labelList& testFaces
466 ) const
467 {
468  labelListList faceZones(surfaces_.surfaces().size());
469  {
470  // If triSurface do additional zoning based on connectivity
471  for (const label surfi : blockedSurfaces)
472  {
473  const label geomi = surfaces_.surfaces()[surfi];
474  const searchableSurface& s = surfaces_.geometry()[geomi];
475  if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
476  {
477  const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
478  const labelListList& edFaces = surf.edgeFaces();
479  boolList isOpenEdge(edFaces.size(), false);
480  forAll(edFaces, i)
481  {
482  if (edFaces[i].size() == 1)
483  {
484  isOpenEdge[i] = true;
485  }
486  }
487 
488  labelList faceZone;
489  const label nZones = surf.markZones(isOpenEdge, faceZone);
490  if (nZones > 1)
491  {
492  faceZones[surfi].transfer(faceZone);
493  }
494  }
495  }
496  }
497 
498 
499  // Clever limiting of the number of iterations (= max cells in the channel)
500  // since it has too many problematic issues, e.g. with volume refinement
501  // and the real check uses the proper distance anyway just disable.
502  const label nIters = mesh_.globalData().nTotalCells();
503 
504 
505  // Collect candidate faces (i.e. intersecting any surface and
506  // owner/neighbour not yet refined)
507  //const labelList testFaces(getRefineCandidateFaces(refineCell));
508 
509  // Collect segments
510  pointField start(testFaces.size());
511  pointField end(testFaces.size());
512  labelList minLevel(testFaces.size());
513 
514  calcCellCellRays
515  (
516  neiCc,
517  neiLevel,
518  testFaces,
519  start,
520  end,
521  minLevel
522  );
523  // TBD. correct nIters for actual cellLevel (since e.g. volume refinement
524  // might add to cell level). Unfortunately we only have minLevel,
525  // not maxLevel. Problem: what if volume refinement only in middle of
526  // channel? Say channel 1m wide with a 0.1m sphere of refinement
527  // Workaround: have dummy surface with e.g. maxLevel 100 to
528  // force nIters to be high enough.
529 
530 
531  // Test for all intersections (with surfaces of higher gap level than
532  // minLevel) and cache per cell the max surface level and the local normal
533  // on that surface.
534 
535  labelList surface1;
536  List<pointIndexHit> hit1;
537  labelList region1;
538  vectorField normal1;
539 
540  labelList surface2;
541  List<pointIndexHit> hit2;
542  labelList region2;
543  vectorField normal2;
544 
545  surfaces_.findNearestIntersection
546  (
547  blockedSurfaces,
548  start,
549  end,
550 
551  surface1,
552  hit1,
553  region1, // local region
554  normal1,
555 
556  surface2,
557  hit2,
558  region2, // local region
559  normal2
560  );
561 
562 
563  // Detect cells that are using multiple surface regions
564  //bitSet isMultiRegion;
565  //detectMultiRegionCells
566  //(
567  // faceZones,
568  // testFaces,
569  //
570  // surface1,
571  // hit1,
572  // region1,
573  //
574  // surface2,
575  // hit2,
576  // region2,
577  //
578  // isMultiRegion
579  //);
580 
581 
582  label n = 0;
583  forAll(testFaces, i)
584  {
585  if (hit1[i].hit())
586  {
587  n++;
588  }
589  }
590 
591  List<wallPoints> faceDist(n);
592  labelList changedFaces(n);
593  n = 0;
594 
595  DynamicList<point> originLocation(2);
596  DynamicList<scalar> originDistSqr(2);
597  DynamicList<FixedList<label, 3>> originSurface(2);
598  //DynamicList<point> originNormal(2);
599 
600 
601  //- To avoid walking through surfaces we mark all faces that have been
602  // intersected. We can either mark only those faces intersecting
603  // blockedSurfaces (i.e. with a 'blockLevel') or mark faces intersecting
604  // any (refinement) surface (this includes e.g. faceZones). This is
605  // much easier since that information is already cached
606  // (meshRefinement::intersectedFaces())
607 
608  //bitSet isBlockedFace(mesh_.nFaces());
609  forAll(testFaces, i)
610  {
611  if (hit1[i].hit())
612  {
613  const label facei = testFaces[i];
614  //isBlockedFace.set(facei);
615  const point& fc = mesh_.faceCentres()[facei];
616  const labelList& fz1 = faceZones[surface1[i]];
617 
618  originLocation.clear();
619  originDistSqr.clear();
620  //originNormal.clear();
621  originSurface.clear();
622 
623  originLocation.append(hit1[i].point());
624  originDistSqr.append(hit1[i].point().distSqr(fc));
625  //originNormal.append(normal1[i]);
626  originSurface.append
627  (
628  FixedList<label, 3>
629  ({
630  surface1[i],
631  region1[i],
632  (fz1.size() ? fz1[hit1[i].index()] : region1[i])
633  })
634  );
635 
636  if (hit2[i].hit() && hit1[i] != hit2[i])
637  {
638  const labelList& fz2 = faceZones[surface2[i]];
639  originLocation.append(hit2[i].point());
640  originDistSqr.append(hit2[i].point().distSqr(fc));
641  //originNormal.append(normal2[i]);
642  originSurface.append
643  (
644  FixedList<label, 3>
645  ({
646  surface2[i],
647  region2[i],
648  (fz2.size() ? fz2[hit2[i].index()] : region2[i])
649  })
650  );
651  }
652 
653  // Collect all seed data. Currently walking does not look at
654  // surface direction - if so pass in surface normal as well
655  faceDist[n] = wallPoints
656  (
657  originLocation, // origin
658  originDistSqr, // distance to origin
659  originSurface // surface+region+zone
660  //originNormal // normal at origin
661  );
662  changedFaces[n] = facei;
663  n++;
664  }
665  }
666 
667 
668  // Clear intersection info
669  surface1.clear();
670  hit1.clear();
671  region1.clear();
672  normal1.clear();
673  surface2.clear();
674  hit2.clear();
675  region2.clear();
676  normal2.clear();
677 
678 
679  List<wallPoints> allFaceInfo(mesh_.nFaces());
680  List<wallPoints> allCellInfo(mesh_.nCells());
681 
682  // Any refinement surface (even a faceZone) should stop the gap walking.
683  // This is exactly the information which is cached in the surfaceIndex_
684  // field.
685  const bitSet isBlockedFace(intersectedFaces());
686 
687  wallPoints::trackData td(isBlockedFace, regionToBlockSize);
688  FaceCellWave<wallPoints, wallPoints::trackData> wallDistCalc
689  (
690  mesh_,
691  changedFaces,
692  faceDist,
693  allFaceInfo,
694  allCellInfo,
695  0, // max iterations
696  td
697  );
698  wallDistCalc.iterate(nIters);
699 
700 
701  // Extract faces of visited cells
702 
703  bitSet isProximityFace(mesh_.nFaces(), false);
704 
705  // Make sure to include initial intersected ones
706  isProximityFace.set(testFaces);
707 
708  forAll(allCellInfo, celli)
709  {
710  if (allCellInfo[celli].valid(wallDistCalc.data()))
711  {
712  isProximityFace.set(mesh_.cells()[celli]);
713  }
714  }
715 
716  syncTools::syncFaceList
717  (
718  mesh_,
719  isProximityFace,
720  orEqOp<unsigned int>(),
721  0u
722  );
723 
724  testFaces = isProximityFace.toc();
725 }
726 
727 
728 Foam::label Foam::meshRefinement::markFakeGapRefinement
729 (
730  const scalar planarCos,
731  const labelList& blockedSurfaces,
732  const label nAllowRefine,
733  const labelList& neiLevel,
734  const pointField& neiCc,
735 
736  labelList& refineCell,
737  label& nRefine
738 ) const
739 {
740  // Re-work the blockLevel of the blockedSurfaces into a length scale
741  // and a number of cells to cross
742  List<scalarList> regionToBlockSize(surfaces_.surfaces().size());
743  scalar maxBlockSize(-1);
744  for (const label surfi : blockedSurfaces)
745  {
746  const label geomi = surfaces_.surfaces()[surfi];
747  const searchableSurface& s = surfaces_.geometry()[geomi];
748  const label nRegions = s.regions().size();
749 
750  regionToBlockSize[surfi].setSize(nRegions, -1);
751  for (label regioni = 0; regioni < nRegions; regioni++)
752  {
753  const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
754  const label bLevel = surfaces_.blockLevel()[globalRegioni];
755 
756  // If valid blockLevel specified
757  if (bLevel != labelMax)
758  {
759  regionToBlockSize[surfi][regioni] =
760  meshCutter_.level0EdgeLength()/pow(2.0, bLevel);
761  }
762 
763  maxBlockSize = max(maxBlockSize, regionToBlockSize[surfi][regioni]);
764 
765  //const label mLevel = surfaces_.maxLevel()[globalRegioni];
769  //if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
770  //{
771  // const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
772  //}
773  }
774  }
775 
776  // Collect candidate faces (i.e. intersecting any surface and
777  // owner/neighbour not yet refined) and their cells
778  labelList cellMap;
779  {
780  labelList testFaces(getRefineCandidateFaces(refineCell));
781 
782  // Extend the set by faceCellFace wave upto given 3*blockLevel size
783  markProximityCandidateFaces
784  (
785  blockedSurfaces,
786  regionToBlockSize,
787  neiLevel,
788  neiCc,
789  testFaces
790  );
791 
792 
793  bitSet isTestCell(mesh_.nCells());
794  for (const label facei : testFaces)
795  {
796  const label own = mesh_.faceOwner()[facei];
797  if (refineCell[own] == -1)
798  {
799  isTestCell.set(own);
800  }
801  if (mesh_.isInternalFace(facei))
802  {
803  const label nei = mesh_.faceNeighbour()[facei];
804  if (refineCell[nei] == -1)
805  {
806  isTestCell.set(nei);
807  }
808  }
809  }
810  cellMap = isTestCell.sortedToc();
811  }
812 
813 
814  // Shoot rays in all 6 directions
815  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
816  //
817  // Almost exact copy of meshRefinement::markInternalGapRefinement.
818  // Difference is only in the length scale (now based on blockLevel)
819 
820 
821  // Find per cell centre the nearest point and normal on the surfaces
822  List<pointIndexHit> nearInfo;
823  vectorField nearNormal;
824  labelList nearSurface;
825  labelList nearRegion;
826  {
827  surfaces_.findNearestRegion
828  (
829  blockedSurfaces,
830  pointField(mesh_.cellCentres(), cellMap),
831  scalarField(cellMap.size(), maxBlockSize), // use uniform for now
832  nearSurface,
833  nearInfo,
834  nearRegion,
835  nearNormal
836  );
837  }
838 
839 
840  DynamicList<label> map(nearInfo.size());
841  DynamicField<point> rayStart(nearInfo.size());
842  DynamicField<point> rayEnd(nearInfo.size());
843  DynamicField<scalar> gapSize(nearInfo.size());
844 
845  DynamicField<point> rayStart2(nearInfo.size());
846  DynamicField<point> rayEnd2(nearInfo.size());
847  DynamicField<scalar> gapSize2(nearInfo.size());
848 
849  label nTestCells = 0;
850 
851  forAll(nearInfo, i)
852  {
853  if (nearInfo[i].hit())
854  {
855  const label globalRegioni = surfaces_.globalRegion
856  (
857  nearSurface[i],
858  nearRegion[i]
859  );
860 
861  // Combine info from shell and surface
862  const label bLevel = surfaces_.blockLevel()[globalRegioni];
863  const FixedList<label, 3> gapInfo
864  ({
865  100,
866  bLevel,
867  100
868  });
869  const volumeType gapMode(volumeType::MIXED);
870 
871  // Store wanted number of cells in gap
872  const label celli = cellMap[i];
873  const label cLevel = meshCutter_.cellLevel()[celli];
874 
875  // Construct one or more rays to test for oppositeness
876  const label nRays = generateRays
877  (
878  false, // ignore actual surface normal
879  nearInfo[i].hitPoint(),
880  nearNormal[i],
881  gapInfo,
882  gapMode,
883 
884  mesh_.cellCentres()[celli],
885  cLevel,
886 
887  rayStart,
888  rayEnd,
889  gapSize,
890 
891  rayStart2,
892  rayEnd2,
893  gapSize2
894  );
895  if (nRays > 0)
896  {
897  nTestCells++;
898  for (label j = 0; j < nRays; j++)
899  {
900  map.append(i);
901  }
902  }
903  }
904  }
905 
906  Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
907  << " cells for testing out of "
908  << mesh_.globalData().nTotalCells() << endl;
909  map.shrink();
910  rayStart.shrink();
911  rayEnd.shrink();
912  gapSize.shrink();
913 
914  rayStart2.shrink();
915  rayEnd2.shrink();
916  gapSize2.shrink();
917 
918  cellMap = labelUIndList(cellMap, map)();
919  nearNormal = UIndirectList<vector>(nearNormal, map)();
920  nearInfo.clear();
921  nearSurface.clear();
922  nearRegion.clear();
923 
924 
925  // Do intersections in pairs
926  labelList surf1;
927  labelList region1;
928  List<pointIndexHit> hit1;
929  vectorField normal1;
930  surfaces_.findNearestIntersection
931  (
932  rayStart,
933  rayEnd,
934  surf1,
935  region1,
936  hit1,
937  normal1
938  );
939 
940  labelList surf2;
941  labelList region2;
942  List<pointIndexHit> hit2;
943  vectorField normal2;
944  surfaces_.findNearestIntersection
945  (
946  rayStart2,
947  rayEnd2,
948  surf2,
949  region2,
950  hit2,
951  normal2
952  );
953 
954  // Extract cell based gap size
955  const label oldNRefine = nRefine;
956  forAll(surf1, i)
957  {
958  if (surf1[i] != -1 && surf2[i] != -1)
959  {
960  // Found intersections with surface. Check for
961  // - small gap
962  // - coplanar normals
963 
964  const label celli = cellMap[i];
965  if (celli != -1 && refineCell[celli] == -1)
966  {
967  const scalar d2 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());
968 
969  const scalar maxGapSize
970  (
971  max
972  (
973  regionToBlockSize[surf1[i]][region1[i]],
974  regionToBlockSize[surf2[i]][region2[i]]
975  )
976  );
977 
978  if
979  (
980  (mag(normal1[i]&normal2[i]) > planarCos)
981  && (d2 < Foam::sqr(maxGapSize))
982  )
983  {
984  if
985  (
986  !markForRefine
987  (
988  0, // mark level
989  nAllowRefine,
990  refineCell[celli],
991  nRefine
992  )
993  )
994  {
995  if (debug)
996  {
997  Pout<< "Stopped refining since reaching my cell"
998  << " limit of " << mesh_.nCells()+7*nRefine
999  << endl;
1000  }
1001  break;
1002  }
1003  }
1004  }
1005  }
1006  }
1007 
1008  if
1009  (
1010  returnReduce(nRefine, sumOp<label>())
1011  > returnReduce(nAllowRefine, sumOp<label>())
1012  )
1013  {
1014  Info<< "Reached refinement limit." << endl;
1015  }
1016 
1017  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1018 }
1019 
1020 
1021 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeGapCells
1022 (
1023  const scalar planarAngle,
1024  const labelList& minSurfaceLevel,
1025  const labelList& globalToMasterPatch,
1026  const label growIter
1027 )
1028 {
1029  // Swap neighbouring cell centres and cell level
1030  labelList neiLevel(mesh_.nBoundaryFaces());
1031  pointField neiCc(mesh_.nBoundaryFaces());
1032  calcNeighbourData(neiLevel, neiCc);
1033 
1034  labelList refineCell(mesh_.nCells(), -1);
1035  label nRefine = 0;
1036  //markProximityRefinement
1037  //(
1038  // Foam::cos(degToRad(planarAngle)),
1039  //
1040  // minSurfaceLevel, // surface min level
1041  // labelList(minSurfaceLevel.size(), labelMax), // surfaceGapLevel
1042  //
1043  // labelMax/Pstream::nProcs(), //nAllowRefine,
1044  // neiLevel,
1045  // neiCc,
1046  //
1047  // refineCell,
1048  // nRefine
1049  //);
1050 
1051 
1052  // Determine minimum blockLevel per surface
1053  Map<label> surfToBlockLevel;
1054 
1055  forAll(surfaces_.surfaces(), surfi)
1056  {
1057  const label geomi = surfaces_.surfaces()[surfi];
1058  const searchableSurface& s = surfaces_.geometry()[geomi];
1059  const label nRegions = s.regions().size();
1060 
1061  label minBlockLevel = labelMax;
1062  for (label regioni = 0; regioni < nRegions; regioni++)
1063  {
1064  const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
1065  minBlockLevel = min
1066  (
1067  minBlockLevel,
1068  surfaces_.blockLevel()[globalRegioni]
1069  );
1070  }
1071 
1072  if (minBlockLevel < labelMax)
1073  {
1074  surfToBlockLevel.insert(surfi, minBlockLevel);
1075  }
1076  }
1077 
1078 
1079  //markProximityRefinementWave
1080  //(
1081  // Foam::cos(degToRad(planarAngle)),
1082  // surfToBlockLevel.sortedToc(),
1083  //
1084  // labelMax/Pstream::nProcs(), //nAllowRefine,
1085  // neiLevel,
1086  // neiCc,
1087  //
1088  // refineCell,
1089  // nRefine
1090  //);
1091 
1092  markFakeGapRefinement
1093  (
1094  Foam::cos(degToRad(planarAngle)),
1095  surfToBlockLevel.sortedToc(),
1096 
1097  labelMax/Pstream::nProcs(), //nAllowRefine,
1098  neiLevel,
1099  neiCc,
1100 
1101  refineCell,
1102  nRefine
1103  );
1104 
1105 
1106  Info<< "Marked for blocking due to close opposite surfaces : "
1107  << returnReduce(nRefine, sumOp<label>()) << " cells." << endl;
1108 
1109  // Remove outliers, i.e. cells with all points exposed
1110  if (growIter)
1111  {
1112  labelList oldRefineCell(refineCell);
1113 
1114  // Pass1: extend the set to fill in gaps
1115  bitSet isOutsideFace;
1116  for (label iter = 0; iter < growIter; iter++)
1117  {
1118  // Get outside faces
1119  markOutsideFaces
1120  (
1121  meshCutter_.cellLevel(),
1122  neiLevel,
1123  refineCell,
1124  isOutsideFace
1125  );
1126  // Extend with cells with three outside faces
1127  growSet(neiLevel, isOutsideFace, refineCell, nRefine);
1128  }
1129 
1130 
1131  // Pass2: erode back to original set if pass1 didn't help
1132  for (label iter = 0; iter < growIter; iter++)
1133  {
1134  // Get outside faces. Ignore cell level.
1135  markOutsideFaces
1136  (
1137  labelList(mesh_.nCells(), 0),
1138  labelList(neiLevel.size(), 0),
1139  refineCell,
1140  isOutsideFace
1141  );
1142 
1143  // Unmark cells with three or more outside faces
1144  for (label celli = 0; celli < mesh_.nCells(); celli++)
1145  {
1146  if (refineCell[celli] != -1 && oldRefineCell[celli] == -1)
1147  {
1148  if (countFaceDirs(isOutsideFace, celli) >= 3)
1149  {
1150  refineCell[celli] = -1;
1151  --nRefine;
1152  }
1153  }
1154  }
1155  }
1156 
1157  Info<< "Marked for blocking after filtering : "
1158  << returnReduce(nRefine, sumOp<label>()) << " cells." << endl;
1159  }
1160 
1161 
1162  // Determine patch for every mesh face
1163  const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
1164  labelList unnamedSurfaces(surfaceZonesInfo::getUnnamedSurfaces(surfZones));
1165  const label defaultRegion(surfaces_.globalRegion(unnamedSurfaces[0], 0));
1166 
1167  const labelList nearestRegion
1168  (
1169  nearestIntersection
1170  (
1171  unnamedSurfaces,
1172  defaultRegion
1173  )
1174  );
1175 
1176  // Pack
1177  labelList cellsToRemove(nRefine);
1178  nRefine = 0;
1179 
1180  forAll(refineCell, cellI)
1181  {
1182  if (refineCell[cellI] != -1)
1183  {
1184  cellsToRemove[nRefine++] = cellI;
1185  }
1186  }
1187 
1188  // Remove cells
1189  removeCells cellRemover(mesh_);
1190  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
1191 
1192  labelList exposedPatches(exposedFaces.size());
1193  forAll(exposedFaces, i)
1194  {
1195  label facei = exposedFaces[i];
1196  exposedPatches[i] = globalToMasterPatch[nearestRegion[facei]];
1197  }
1198 
1199  return doRemoveCells
1200  (
1201  cellsToRemove,
1202  exposedFaces,
1203  exposedPatches,
1204  cellRemover
1205  );
1206 }
1207 
1208 
1209 void Foam::meshRefinement::selectIntersectedFaces
1210 (
1211  const labelList& selectedSurfaces,
1212  boolList& isBlockedFace
1213 ) const
1214 {
1215  // Like meshRefinement::selectSeparatedCoupledFaces. tbd: convert to bitSet
1216 
1217  // Check that no connection between inside and outside points
1218  isBlockedFace.setSize(mesh_.nFaces(), false);
1219 
1220  // Block off separated couples.
1221  selectSeparatedCoupledFaces(isBlockedFace);
1222 
1223  // Block off intersections with selected surfaces
1224 
1225  // Mark per face (for efficiency)
1226  boolList isSelectedSurf(surfaces_.surfaces().size(), false);
1227  UIndirectList<bool>(isSelectedSurf, selectedSurfaces) = true;
1228 
1229  forAll(surfaceIndex_, facei)
1230  {
1231  const label surfi = surfaceIndex_[facei];
1232  if (surfi != -1 && isSelectedSurf[surfi])
1233  {
1234  isBlockedFace[facei] = true;
1235  }
1236  }
1237 }
1238 
1239 
1240 //Foam::labelList Foam::meshRefinement::detectLeakCells
1241 //(
1242 // const boolList& isBlockedFace,
1243 // const labelList& leakFaces,
1244 // const labelList& seedCells
1245 //) const
1246 //{
1247 // int dummyTrackData = 0;
1248 // List<topoDistanceData<label>> allFaceInfo(mesh_.nFaces());
1249 // List<topoDistanceData<label>> allCellInfo(mesh_.nCells());
1250 //
1251 // // Block faces
1252 // forAll(isBlockedFace, facei)
1253 // {
1254 // if (isBlockedFace[facei])
1255 // {
1256 // allFaceInfo[facei] = topoDistanceData<label>(labelMax, 123);
1257 // }
1258 // }
1259 // for (const label facei : leakFaces)
1260 // {
1261 // allFaceInfo[facei] = topoDistanceData<label>(labelMax, 123);
1262 // }
1263 //
1264 // // Walk from inside cell
1265 // DynamicList<topoDistanceData<label>> faceDist;
1266 // DynamicList<label> cFaces1;
1267 // for (const label celli : seedCells)
1268 // {
1269 // if (celli != -1)
1270 // {
1271 // const labelList& cFaces = mesh_.cells()[celli];
1272 // faceDist.reserve(cFaces.size());
1273 // cFaces1.reserve(cFaces.size());
1274 //
1275 // for (const label facei : cFaces)
1276 // {
1277 // if (!allFaceInfo[facei].valid(dummyTrackData))
1278 // {
1279 // cFaces1.append(facei);
1280 // faceDist.append(topoDistanceData<label>(0, 123));
1281 // }
1282 // }
1283 // }
1284 // }
1285 //
1286 // // Walk through face-cell wave till all cells are reached
1287 // FaceCellWave<topoDistanceData<label>> wallDistCalc
1288 // (
1289 // mesh_,
1290 // cFaces1,
1291 // faceDist,
1292 // allFaceInfo,
1293 // allCellInfo,
1294 // mesh_.globalData().nTotalCells()+1 // max iterations
1295 // );
1296 //
1297 // label nRemove = 0;
1298 // for (const label facei : leakFaces)
1299 // {
1300 // const label own = mesh_.faceOwner()[facei];
1301 // if (!allCellInfo[own].valid(dummyTrackData))
1302 // {
1303 // nRemove++;
1304 // }
1305 // if (mesh_.isInternalFace(facei))
1306 // {
1307 // const label nei = mesh_.faceNeighbour()[facei];
1308 // if (!allCellInfo[nei].valid(dummyTrackData))
1309 // {
1310 // nRemove++;
1311 // }
1312 // }
1313 // }
1314 //
1315 // labelList cellsToRemove(nRemove);
1316 // nRemove = 0;
1317 // for (const label facei : leakFaces)
1318 // {
1319 // const label own = mesh_.faceOwner()[facei];
1320 // if (!allCellInfo[own].valid(dummyTrackData))
1321 // {
1322 // cellsToRemove[nRemove++] = own;
1323 // }
1324 // if (mesh_.isInternalFace(facei))
1325 // {
1326 // const label nei = mesh_.faceNeighbour()[facei];
1327 // if (!allCellInfo[nei].valid(dummyTrackData))
1328 // {
1329 // cellsToRemove[nRemove++] = nei;
1330 // }
1331 // }
1332 // }
1333 //
1334 // if (debug)
1335 // {
1336 // volScalarField fld
1337 // (
1338 // IOobject
1339 // (
1340 // "cellsToKeep",
1341 // mesh_.time().timeName(),
1342 // mesh_,
1343 // IOobject::NO_READ,
1344 // IOobject::NO_WRITE
1345 // ),
1346 // mesh_,
1347 // dimensionedScalar(dimless, Zero)
1348 // );
1349 // forAll(allCellInfo, celli)
1350 // {
1351 // if (allCellInfo[celli].valid(dummyTrackData))
1352 // {
1353 // fld[celli] = allCellInfo[celli].distance();
1354 // }
1355 // }
1356 // forAll(fld.boundaryField(), patchi)
1357 // {
1358 // const polyPatch& pp = mesh_.boundaryMesh()[patchi];
1359 // SubList<topoDistanceData<label>> p(pp.patchSlice(allFaceInfo));
1360 // scalarField pfld
1361 // (
1362 // fld.boundaryField()[patchi].size(),
1363 // Zero
1364 // );
1365 // forAll(pfld, i)
1366 // {
1367 // if (p[i].valid(dummyTrackData))
1368 // {
1369 // pfld[i] = 1.0*p[i].distance();
1370 // }
1371 // }
1372 // fld.boundaryFieldRef()[patchi] == pfld;
1373 // }
1374 // //Note: do not swap cell values so do not do
1375 // //fld.correctBoundaryConditions();
1376 // Pout<< "Writing distance field for initial cells "
1377 // << seedCells << " to " << fld.objectPath() << endl;
1378 // fld.write();
1379 // }
1380 //
1381 // return cellsToRemove;
1382 //}
1383 //
1384 //
1385 //Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeLeakCells
1386 //(
1387 // const labelList& globalToMasterPatch,
1388 // const labelList& globalToSlavePatch,
1389 // const pointField& locationsInMesh,
1390 // const wordList& zonesInMesh,
1391 // const pointField& locationsOutsideMesh,
1392 // const labelList& selectedSurfaces
1393 //)
1394 //{
1395 // boolList isBlockedFace;
1396 // selectIntersectedFaces(selectedSurfaces, isBlockedFace);
1397 //
1398 // // Determine cell regions
1399 // const regionSplit cellRegion(mesh_, isBlockedFace);
1400 //
1401 // // Detect locationsInMesh regions
1402 // labelList insideCells(locationsInMesh.size(), -1);
1403 // labelList insideRegions(locationsInMesh.size(), -1);
1404 // forAll(locationsInMesh, i)
1405 // {
1406 // insideCells[i] = findCell
1407 // (
1408 // mesh_,
1409 // mergeDistance_*vector::one, //perturbVec,
1410 // locationsInMesh[i]
1411 // );
1412 // if (insideCells[i] != -1)
1413 // {
1414 // insideRegions[i] = cellRegion[insideCells[i]];
1415 // }
1416 // reduce(insideRegions[i], maxOp<label>());
1417 //
1418 // if (insideRegions[i] == -1)
1419 // {
1420 // // See if we can perturb a bit
1421 // insideCells[i] = findCell
1422 // (
1423 // mesh_,
1424 // mergeDistance_*vector::one, //perturbVec,
1425 // locationsInMesh[i]+mergeDistance_*vector::one
1426 // );
1427 // if (insideCells[i] != -1)
1428 // {
1429 // insideRegions[i] = cellRegion[insideCells[i]];
1430 // }
1431 // reduce(insideRegions[i], maxOp<label>());
1432 //
1433 // if (insideRegions[i] == -1)
1434 // {
1435 // FatalErrorInFunction
1436 // << "Cannot find locationInMesh " << locationsInMesh[i]
1437 // << " on any processor" << exit(FatalError);
1438 // }
1439 // }
1440 // }
1441 //
1442 //
1443 // // Check that all the locations outside the
1444 // // mesh do not conflict with those inside
1445 //
1446 // bool haveLeak = false;
1447 // forAll(locationsOutsideMesh, i)
1448 // {
1449 // // Find the region containing the point
1450 // label regioni = findRegion
1451 // (
1452 // mesh_,
1453 // cellRegion,
1454 // mergeDistance_*vector::one, //perturbVec,
1455 // locationsOutsideMesh[i]
1456 // );
1457 //
1458 // if (regioni != -1)
1459 // {
1460 // // Check for locationsOutsideMesh overlapping with inside ones
1461 // if (insideRegions.find(regioni) != -1)
1462 // {
1463 // haveLeak = true;
1464 // WarningInFunction
1465 // << "Outside location " << locationsOutsideMesh[i]
1466 // << " in region " << regioni
1467 // << " is connected to one of the inside points "
1468 // << locationsInMesh << endl;
1469 // }
1470 // }
1471 // }
1472 //
1473 //
1474 // autoPtr<mapPolyMesh> mapPtr;
1475 // if (returnReduceOr(haveLeak))
1476 // {
1477 // // Use shortestPathSet to provide a minimum set of faces needed
1478 // // to close hole. Tbd: maybe directly use wave?
1479 // meshSearch searchEngine(mesh_);
1480 // shortestPathSet leakPath
1481 // (
1482 // "leakPath",
1483 // mesh_,
1484 // searchEngine,
1485 // coordSet::coordFormatNames[coordSet::coordFormat::DISTANCE],
1486 // true,
1487 // 50, // tbd. Number of iterations. This is the maximum
1488 // // number of faces in the leak hole
1489 //
1490 // //pbm.groupPatchIDs()["wall"], // patch to grow from
1491 // meshedPatches(), // patch to grow from
1492 //
1493 // locationsInMesh,
1494 // locationsOutsideMesh,
1495 // isBlockedFace
1496 // );
1497 //
1498 //
1499 // // Use leak path to find minimum set of cells to delete
1500 // const labelList cellsToRemove
1501 // (
1502 // detectLeakCells
1503 // (
1504 // isBlockedFace,
1505 // leakPath.leakFaces(),
1506 // insideCells
1507 // )
1508 // );
1509 //
1510 // // Re-do intersections to find nearest unnamed surface
1511 // const label defaultRegion
1512 // (
1513 // surfaces().globalRegion
1514 // (
1515 // selectedSurfaces[0],
1516 // 0
1517 // )
1518 // );
1519 //
1520 // const labelList nearestRegion
1521 // (
1522 // nearestIntersection
1523 // (
1524 // selectedSurfaces,
1525 // defaultRegion
1526 // )
1527 // );
1528 //
1529 //
1530 // // Remove cells
1531 // removeCells cellRemover(mesh_);
1532 // const labelList exposedFaces
1533 // (
1534 // cellRemover.getExposedFaces(cellsToRemove)
1535 // );
1536 //
1537 // labelList exposedPatches(exposedFaces.size());
1538 // forAll(exposedFaces, i)
1539 // {
1540 // label facei = exposedFaces[i];
1541 // exposedPatches[i] = globalToMasterPatch[nearestRegion[facei]];
1542 // }
1543 //
1544 // mapPtr = doRemoveCells
1545 // (
1546 // cellsToRemove,
1547 // exposedFaces,
1548 // exposedPatches,
1549 // cellRemover
1550 // );
1551 //
1552 //
1553 // // Put the exposed faces into a special faceZone
1554 // {
1555 // // Add to "frozenFaces" zone
1556 // faceZoneMesh& faceZones = mesh_.faceZones();
1557 //
1558 // // Get current frozen faces (if any)
1559 // bitSet isFrozenFace(mesh_.nFaces());
1560 // label zonei = faceZones.findZoneID("frozenFaces");
1561 // if (zonei != -1)
1562 // {
1563 // const bitSet oldSet(mesh_.nFaces(), faceZones[zonei]);
1564 // isFrozenFace.set(oldSet);
1565 // }
1566 //
1567 // // Add newly exposed faces (if not yet in any faceZone!)
1568 // const labelList exposed
1569 // (
1570 // renumber
1571 // (
1572 // mapPtr().reverseFaceMap(),
1573 // exposedFaces
1574 // )
1575 // );
1576 // bitSet isZonedFace(mesh_.nFaces(), faceZones.zoneMap().toc());
1577 // for (const label facei : exposed)
1578 // {
1579 // if (!isZonedFace[facei])
1580 // {
1581 // isFrozenFace.set(facei);
1582 // }
1583 // }
1584 //
1585 // syncTools::syncFaceList
1586 // (
1587 // mesh_,
1588 // isFrozenFace,
1589 // orEqOp<unsigned int>(),
1590 // 0u
1591 // );
1592 //
1593 // // Add faceZone if non existing
1594 // faceZones.clearAddressing();
1595 // if (zonei == -1)
1596 // {
1597 // zonei = faceZones.size();
1598 //
1599 // faceZones.emplace_back
1600 // (
1601 // "frozenFaces", // name
1602 // zonei, // index
1603 // faceZones
1604 // );
1605 // }
1606 //
1607 // // Update faceZone with new contents
1608 // labelList frozenFaces(isFrozenFace.toc());
1609 // boolList frozenFlip(frozenFaces.size(), false);
1610 //
1611 // faceZones[zonei].resetAddressing
1612 // (
1613 // std::move(frozenFaces),
1614 // std::move(frozenFlip)
1615 // );
1616 // }
1617 //
1618 //
1619 // //// Put the exposed points into a special pointZone
1620 // //if (false)
1621 // //{
1622 // // const labelList meshFaceIDs
1623 // // (
1624 // // renumber
1625 // // (
1626 // // mapPtr().reverseFaceMap(),
1627 // // exposedFaces
1628 // // )
1629 // // );
1630 // // const uindirectPrimitivePatch pp
1631 // // (
1632 // // UIndirectList<face>(mesh_.faces(), meshFaceIDs),
1633 // // mesh_.points()
1634 // // );
1635 // //
1636 // // // Count number of faces per edge
1637 // // const labelListList& edgeFaces = pp.edgeFaces();
1638 // // labelList nEdgeFaces(edgeFaces.size());
1639 // // forAll(edgeFaces, edgei)
1640 // // {
1641 // // nEdgeFaces[edgei] = edgeFaces[edgei].size();
1642 // // }
1643 // //
1644 // // // Sync across processor patches
1645 // // if (Pstream::parRun())
1646 // // {
1647 // // const globalMeshData& globalData = mesh_.globalData();
1648 // // const mapDistribute& map = globalData.globalEdgeSlavesMap();
1649 // // const indirectPrimitivePatch& cpp =
1650 // // globalData.coupledPatch();
1651 // //
1652 // // // Match pp edges to coupled edges
1653 // // labelList patchEdges;
1654 // // labelList coupledEdges;
1655 // // PackedBoolList sameEdgeOrientation;
1656 // // PatchTools::matchEdges
1657 // // (
1658 // // pp,
1659 // // cpp,
1660 // // patchEdges,
1661 // // coupledEdges,
1662 // // sameEdgeOrientation
1663 // // );
1664 // //
1665 // //
1666 // // // Convert patch-edge data into cpp-edge data
1667 // // labelList coupledNEdgeFaces(map.constructSize(), Zero);
1668 // // UIndirectList<label>(coupledNEdgeFaces, coupledEdges) =
1669 // // UIndirectList<label>(nEdgeFaces, patchEdges);
1670 // //
1671 // // // Synchronise
1672 // // globalData.syncData
1673 // // (
1674 // // coupledNEdgeFaces,
1675 // // globalData.globalEdgeSlaves(),
1676 // // globalData.globalEdgeTransformedSlaves(),
1677 // // map,
1678 // // plusEqOp<label>()
1679 // // );
1680 // //
1681 // // // Convert back from cpp-edge to patch-edge
1682 // // UIndirectList<label>(nEdgeFaces, patchEdges) =
1683 // // UIndirectList<label>(coupledNEdgeFaces, coupledEdges);
1684 // // }
1685 // //
1686 // // // Freeze all internal points
1687 // // bitSet isFrozenPoint(mesh_.nPoints());
1688 // // forAll(nEdgeFaces, edgei)
1689 // // {
1690 // // if (nEdgeFaces[edgei] != 1)
1691 // // {
1692 // // const edge& e = pp.edges()[edgei];
1693 // // isFrozenPoint.set(pp.meshPoints()[e[0]]);
1694 // // isFrozenPoint.set(pp.meshPoints()[e[1]]);
1695 // // }
1696 // // }
1697 // // // Add to "frozenPoints" zone
1698 // // pointZoneMesh& pointZones = mesh_.pointZones();
1699 // // pointZones.clearAddressing();
1700 // // label zonei = pointZones.findZoneID("frozenPoints");
1701 // // if (zonei != -1)
1702 // // {
1703 // // const bitSet oldSet(mesh_.nPoints(), pointZones[zonei]);
1704 // // // Add to isFrozenPoint
1705 // // isFrozenPoint.set(oldSet);
1706 // // }
1707 // //
1708 // // syncTools::syncPointList
1709 // // (
1710 // // mesh_,
1711 // // isFrozenPoint,
1712 // // orEqOp<unsigned int>(),
1713 // // 0u
1714 // // );
1715 // //
1716 // // if (zonei == -1)
1717 // // {
1718 // // zonei = pointZones.size();
1719 // //
1720 // // pointZones.emplace_back
1721 // // (
1722 // // "frozenPoints", // name
1723 // // isFrozenPoint.sortedToc(), // addressing
1724 // // zonei, // index
1725 // // pointZones // pointZoneMesh
1726 // // );
1727 // // }
1728 // //}
1729 //
1730 //
1731 // if (debug&meshRefinement::MESH)
1732 // {
1733 // const_cast<Time&>(mesh_.time())++;
1734 //
1735 // Pout<< "Writing current mesh to time "
1736 // << timeName() << endl;
1737 // write
1738 // (
1739 // meshRefinement::debugType(debug),
1740 // meshRefinement::writeType
1741 // (
1742 // meshRefinement::writeLevel()
1743 // | meshRefinement::WRITEMESH
1744 // ),
1745 // mesh_.time().path()/timeName()
1746 // );
1747 // Pout<< "Dumped mesh in = "
1748 // << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
1749 // }
1750 // }
1751 // return mapPtr;
1752 //}
1753 
1754 
1755 //Foam::autoPtr<Foam::mapDistribute> Foam::meshRefinement::nearestFace
1756 //(
1757 // const globalIndex& globalSeedFaces,
1758 // const labelList& seedFaces,
1759 // const labelList& closureFaces
1760 //) const
1761 //{
1762 // // Takes a set of faces for which we have information (seedFaces,
1763 // // globalSeedFaces - these are e.g. intersected faces) and walks out
1764 // // across nonSeedFace. Returns for
1765 // // every nonSeedFace the nearest seed face (in global indexing).
1766 // // Used e.g. in hole closing. Assumes that seedFaces, closureFaces
1767 // // are a small subset of the master
1768 // // so are not large - it uses edge addressing on the closureFaces
1769 //
1770 //
1771 // if (seedFaces.size() != globalSeedFaces.localSize())
1772 // {
1773 // FatalErrorInFunction << "problem : seedFaces:" << seedFaces.size()
1774 // << " globalSeedFaces:" << globalSeedFaces.localSize()
1775 // << exit(FatalError);
1776 // }
1777 //
1778 // //// Mark mesh points that are used by any closureFaces. This is for
1779 // //// initial filtering
1780 // //bitSet isNonSeedPoint(mesh.nPoints());
1781 // //for (const label facei : closureFaces)
1782 // //{
1783 // // isNonSeedPoint.set(mesh_.faces()[facei]);
1784 // //}
1785 // //syncTools::syncPointList
1786 // //(
1787 // // mesh_,
1788 // // isNonSeedPoint,
1789 // // orEqOp<unsigned int>(),
1790 // // 0u
1791 // //);
1792 //
1793 // // Make patch
1794 // const uindirectPrimitivePatch pp
1795 // (
1796 // IndirectList<face>(mesh_.faces(), closureFaces),
1797 // mesh_.points()
1798 // );
1799 // const edgeList& edges = pp.edges();
1800 // const labelList& mp = pp.meshPoints();
1801 // const label nBndEdges = pp.nEdges() - pp.nInternalEdges();
1802 //
1803 // // For all faces in seedFaces mark the edge with a face. No special
1804 // // handling for multiple faces sharing the edge - first one wins
1805 // EdgeMap<label> edgeMap(pp.nEdges());
1806 // for (const label facei : seedFaces)
1807 // {
1808 // const label globalFacei = globalSeedFaces.toGlobal(facei);
1809 // const face& f = mesh_.faces()[facei];
1810 // forAll(f, fp)
1811 // {
1812 // label nextFp = f.fcIndex(fp);
1813 // edgeMap.insert(edge(f[fp], f[nextFp]), globalFacei);
1814 // }
1815 // }
1816 // syncTools::syncEdgeMap(mesh_, edgeMap, maxEqOp<label>());
1817 //
1818 //
1819 //
1820 // // Seed
1821 // DynamicList<label> initialEdges(2*nBndEdges);
1822 // DynamicList<edgeTopoDistanceData<label, uindirectPrimitivePatch>>
1823 // initialEdgesInfo(2*nBndEdges);
1824 // forAll(edges, edgei)
1825 // {
1826 // const edge& e = edges[edgei];
1827 // const edge meshE = edge(mp[e[0]], mp[e[1]]);
1828 //
1829 // const auto iter = edgeMap.cfind(meshE);
1830 // if (iter.good())
1831 // {
1832 // initialEdges.append(edgei);
1833 // initialEdgesInfo.append
1834 // (
1835 // edgeTopoDistanceData<label, uindirectPrimitivePatch>
1836 // (
1837 // 0, // distance
1838 // iter.val() // globalFacei
1839 // )
1840 // );
1841 // }
1842 // }
1843 //
1844 // // Data on all edges and faces
1845 // List<edgeTopoDistanceData<label, uindirectPrimitivePatch>> allEdgeInfo
1846 // (
1847 // pp.nEdges()
1848 // );
1849 // List<edgeTopoDistanceData<label, uindirectPrimitivePatch>> allFaceInfo
1850 // (
1851 // pp.size()
1852 // );
1853 //
1854 // // Walk
1855 // PatchEdgeFaceWave
1856 // <
1857 // uindirectPrimitivePatch,
1858 // edgeTopoDistanceData<label, uindirectPrimitivePatch>
1859 // > calc
1860 // (
1861 // mesh_,
1862 // pp,
1863 // initialEdges,
1864 // initialEdgesInfo,
1865 // allEdgeInfo,
1866 // allFaceInfo,
1867 // returnReduce(pp.nEdges(), sumOp<label>())
1868 // );
1869 //
1870 //
1871 // // Per non-seed face the seed face
1872 // labelList closureToSeed(pp.size(), -1);
1873 // forAll(allFaceInfo, facei)
1874 // {
1875 // if (allFaceInfo[facei].valid(calc.data()))
1876 // {
1877 // closureToSeed[facei] = allFaceInfo[facei].data();
1878 // }
1879 // }
1880 //
1881 // List<Map<label>> compactMap;
1882 // return autoPtr<mapDistribute>::New
1883 // (
1884 // globalSeedFaces,
1885 // closureToSeed,
1886 // compactMap
1887 // );
1888 //}
1889 
1890 
1891 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::blockLeakFaces
1892 (
1893  const labelList& globalToMasterPatch,
1894  const labelList& globalToSlavePatch,
1895  const pointField& locationsInMesh,
1896  const wordList& zonesInMesh,
1897  const pointField& locationsOutsideMesh,
1898  const labelList& selectedSurfaces
1899 )
1900 {
1901  // Problem: this is based on cached intersection information. This might
1902  // not include the surface we actually want to use. In which case the
1903  // surface would not be seen as intersected!
1904  boolList isBlockedFace;
1905  selectIntersectedFaces(selectedSurfaces, isBlockedFace);
1906 
1907  // Determine cell regions
1908  const regionSplit cellRegion(mesh_, isBlockedFace);
1909 
1910  // Detect locationsInMesh regions
1911  labelList insideCells(locationsInMesh.size(), -1);
1912  labelList insideRegions(locationsInMesh.size(), -1);
1913  forAll(locationsInMesh, i)
1914  {
1915  insideCells[i] = findCell
1916  (
1917  mesh_,
1918  mergeDistance_*vector::one, //perturbVec,
1919  locationsInMesh[i]
1920  );
1921  if (insideCells[i] != -1)
1922  {
1923  insideRegions[i] = cellRegion[insideCells[i]];
1924  }
1925  reduce(insideRegions[i], maxOp<label>());
1926 
1927  if (insideRegions[i] == -1)
1928  {
1929  // See if we can perturb a bit
1930  insideCells[i] = findCell
1931  (
1932  mesh_,
1933  mergeDistance_*vector::one, //perturbVec,
1934  locationsInMesh[i]+mergeDistance_*vector::one
1935  );
1936  if (insideCells[i] != -1)
1937  {
1938  insideRegions[i] = cellRegion[insideCells[i]];
1939  }
1940  reduce(insideRegions[i], maxOp<label>());
1941 
1942  if (insideRegions[i] == -1)
1943  {
1944  FatalErrorInFunction
1945  << "Cannot find locationInMesh " << locationsInMesh[i]
1946  << " on any processor" << exit(FatalError);
1947  }
1948  }
1949  }
1950 
1951 
1952  // Check that all the locations outside the
1953  // mesh do not conflict with those inside
1954 
1955  bool haveLeak = false;
1956  forAll(locationsOutsideMesh, i)
1957  {
1958  // Find the region containing the point
1959  label regioni = findRegion
1960  (
1961  mesh_,
1962  cellRegion,
1963  mergeDistance_*vector::one, //perturbVec,
1964  locationsOutsideMesh[i]
1965  );
1966 
1967  if (regioni != -1)
1968  {
1969  // Check for locationsOutsideMesh overlapping with inside ones
1970  if (insideRegions.find(regioni) != -1)
1971  {
1972  haveLeak = true;
1973  WarningInFunction
1974  << "Outside location " << locationsOutsideMesh[i]
1975  << " in region " << regioni
1976  << " is connected to one of the inside points "
1977  << locationsInMesh << endl;
1978  }
1979  }
1980  }
1981 
1982 
1983  autoPtr<mapPolyMesh> mapPtr;
1984  if (returnReduceOr(haveLeak))
1985  {
1986  // Use holeToFace to provide a minimum set of faces needed
1987  // to close hole.
1988 
1989  const List<pointField> allLocations
1990  (
1991  refinementParameters::zonePoints
1992  (
1993  locationsInMesh,
1994  zonesInMesh,
1995  locationsOutsideMesh
1996  )
1997  );
1998 
1999  const labelList blockingFaces(findIndices(isBlockedFace, true));
2000 
2001  labelList closureFaces;
2002  labelList closureToBlocked;
2003  autoPtr<mapDistribute> closureMapPtr;
2004  {
2005  const globalIndex globalBlockingFaces(blockingFaces.size());
2006 
2007  closureMapPtr = holeToFace::calcClosure
2008  (
2009  mesh_,
2010  allLocations,
2011  blockingFaces,
2012  globalBlockingFaces,
2013  true, // allow erosion
2014 
2015  closureFaces,
2016  closureToBlocked
2017  );
2018 
2019  if (debug)
2020  {
2021  Pout<< "meshRefinement::blockLeakFaces :"
2022  << " found closure faces:" << closureFaces.size()
2023  << " map:" << bool(closureMapPtr) << endl;
2024  }
2025 
2026  if (!closureMapPtr)
2027  {
2028  FatalErrorInFunction
2029  << "have leak but did not find any closure faces"
2030  << exit(FatalError);
2031  }
2032  }
2033 
2034  // Baffle faces
2035  labelList ownPatch(mesh_.nFaces(), -1);
2036  labelList neiPatch(mesh_.nFaces(), -1);
2037 
2038  // Keep (global) boundary faces in their patch
2039  {
2040  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2041  for (label patchi = 0; patchi < patches.nNonProcessor(); ++patchi)
2042  {
2043  const polyPatch& pp = patches[patchi];
2044 
2045  forAll(pp, i)
2046  {
2047  ownPatch[pp.start()+i] = patchi;
2048  neiPatch[pp.start()+i] = patchi;
2049  }
2050  }
2051  }
2052 
2053  const faceZoneMesh& fzs = mesh_.faceZones();
2054 
2055  // Mark zone per face
2056  labelList faceToZone(mesh_.nFaces(), -1);
2057  boolList faceToFlip(mesh_.nFaces(), false);
2058  forAll(fzs, zonei)
2059  {
2060  const labelList& addressing = fzs[zonei];
2061  const boolList& flipMap = fzs[zonei].flipMap();
2062 
2063  forAll(addressing, i)
2064  {
2065  faceToZone[addressing[i]] = zonei;
2066  faceToFlip[addressing[i]] = flipMap[i];
2067  }
2068  }
2069 
2070 
2071  // Fetch patch and zone info from blockingFaces
2072  labelList packedOwnPatch(labelUIndList(ownPatch, blockingFaces));
2073  closureMapPtr->distribute(packedOwnPatch);
2074  labelList packedNeiPatch(labelUIndList(neiPatch, blockingFaces));
2075  closureMapPtr->distribute(packedNeiPatch);
2076  labelList packedZone(labelUIndList(faceToZone, blockingFaces));
2077  closureMapPtr->distribute(packedZone);
2078  boolList packedFlip(UIndirectList<bool>(faceToFlip, blockingFaces));
2079  closureMapPtr->distribute(packedFlip);
2080  forAll(closureFaces, i)
2081  {
2082  const label facei = closureFaces[i];
2083  const label sloti = closureToBlocked[i];
2084 
2085  if (sloti != -1)
2086  {
2087  // TBD. how to orient own/nei patch?
2088  ownPatch[facei] = packedOwnPatch[sloti];
2089  neiPatch[facei] = packedNeiPatch[sloti];
2090  faceToZone[facei] = packedZone[sloti];
2091  faceToFlip[facei] = packedFlip[sloti];
2092  }
2093  }
2094 
2095 
2096  // Add faces to faceZone. For now do this outside of createBaffles
2097  // below
2098  {
2099  List<DynamicList<label>> zoneToFaces(fzs.size());
2100  List<DynamicList<bool>> zoneToFlip(fzs.size());
2101 
2102  // Add any to-be-patched face
2103  forAll(faceToZone, facei)
2104  {
2105  const label zonei = faceToZone[facei];
2106  if (zonei != -1)
2107  {
2108  zoneToFaces[zonei].append(facei);
2109  zoneToFlip[zonei].append(faceToFlip[facei]);
2110  }
2111  }
2112 
2113  forAll(zoneToFaces, zonei)
2114  {
2115  surfaceZonesInfo::addFaceZone
2116  (
2117  fzs[zonei].name(),
2118  zoneToFaces[zonei],
2119  zoneToFlip[zonei],
2120  mesh_
2121  );
2122  }
2123  }
2124 
2125  // Put the points of closureFaces into a special pointZone
2126  {
2127  const uindirectPrimitivePatch pp
2128  (
2129  UIndirectList<face>(mesh_.faces(), closureFaces),
2130  mesh_.points()
2131  );
2132 
2133  // Count number of faces per edge
2134  const labelList nEdgeFaces(countEdgeFaces(pp));
2135 
2136  // Freeze all internal points
2137  bitSet isFrozenPoint(mesh_.nPoints());
2138  forAll(nEdgeFaces, edgei)
2139  {
2140  if (nEdgeFaces[edgei] != 1)
2141  {
2142  const edge& e = pp.edges()[edgei];
2143  isFrozenPoint.set(pp.meshPoints()[e[0]]);
2144  isFrozenPoint.set(pp.meshPoints()[e[1]]);
2145  }
2146  }
2147 
2148  // Lookup/add pointZone and include its points
2149  pointZoneMesh& pointZones =
2150  const_cast<pointZoneMesh&>(mesh_.pointZones());
2151  const label zonei = addPointZone("frozenPoints");
2152  const bitSet oldSet(mesh_.nPoints(), pointZones[zonei]);
2153  isFrozenPoint.set(oldSet);
2154 
2155  syncTools::syncPointList
2156  (
2157  mesh_,
2158  isFrozenPoint,
2159  orEqOp<unsigned int>(),
2160  0u
2161  );
2162 
2163  // Override addressing
2164  pointZones.clearAddressing();
2165  pointZones[zonei] = isFrozenPoint.sortedToc();
2166  }
2167 
2168 
2169 
2170  // Create baffles for faces
2171  mapPtr = createBaffles(ownPatch, neiPatch);
2172 
2174  //{
2175  // // Add newly exposed faces (if not yet in any faceZone!)
2176  // const labelList exposed
2177  // (
2178  // renumber
2179  // (
2180  // mapPtr().reverseFaceMap(),
2181  // blockingFaces
2182  // )
2183  // );
2184  //
2185  // surfaceZonesInfo::addFaceZone
2186  // (
2187  // "frozenFaces",
2188  // exposed,
2189  // boolList(exposed.size(), false),
2190  // mesh_
2191  // );
2192  //}
2193 
2194 
2195  if (debug&meshRefinement::MESH)
2196  {
2197  const_cast<Time&>(mesh_.time())++;
2198 
2199  Pout<< "Writing current mesh to time "
2200  << timeName() << endl;
2201  write
2202  (
2203  meshRefinement::debugType(debug),
2204  meshRefinement::writeType
2205  (
2206  meshRefinement::writeLevel()
2207  | meshRefinement::WRITEMESH
2208  ),
2209  mesh_.time().path()/timeName()
2210  );
2211  Pout<< "Dumped mesh in = "
2212  << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
2213  }
2214  }
2215  return mapPtr;
2216 }
2217 
2218 
2219 
2220 /*
2221 Foam::label Foam::meshRefinement::markProximityRefinementWave
2222 (
2223  const scalar planarCos,
2224  const labelList& blockedSurfaces,
2225  const label nAllowRefine,
2226  const labelList& neiLevel,
2227  const pointField& neiCc,
2228 
2229  labelList& refineCell,
2230  label& nRefine
2231 ) const
2232 {
2233  labelListList faceZones(surfaces_.surfaces().size());
2234  {
2235  // If triSurface do additional zoning based on connectivity
2236  for (const label surfi : blockedSurfaces)
2237  {
2238  const label geomi = surfaces_.surfaces()[surfi];
2239  const searchableSurface& s = surfaces_.geometry()[geomi];
2240  if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
2241  {
2242  const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
2243  const labelListList& edFaces = surf.edgeFaces();
2244  boolList isOpenEdge(edFaces.size(), false);
2245  forAll(edFaces, i)
2246  {
2247  if (edFaces[i].size() == 1)
2248  {
2249  isOpenEdge[i] = true;
2250  }
2251  }
2252 
2253  labelList faceZone;
2254  const label nZones = surf.markZones(isOpenEdge, faceZone);
2255  if (nZones > 1)
2256  {
2257  faceZones[surfi].transfer(faceZone);
2258  }
2259  }
2260  }
2261  }
2262 
2263 
2264  // Re-work the blockLevel of the blockedSurfaces into a length scale
2265  // and a number of cells to cross
2266  List<scalarList> regionToBlockSize(surfaces_.surfaces().size());
2267  for (const label surfi : blockedSurfaces)
2268  {
2269  const label geomi = surfaces_.surfaces()[surfi];
2270  const searchableSurface& s = surfaces_.geometry()[geomi];
2271  const label nRegions = s.regions().size();
2272  regionToBlockSize[surfi].setSize(nRegions);
2273  for (label regioni = 0; regioni < nRegions; regioni++)
2274  {
2275  const label globalRegioni = surfaces_.globalRegion(surfi, regioni);
2276  const label bLevel = surfaces_.blockLevel()[globalRegioni];
2277  regionToBlockSize[surfi][regioni] =
2278  meshCutter_.level0EdgeLength()/pow(2.0, bLevel);
2279 
2280  //const label mLevel = surfaces_.maxLevel()[globalRegioni];
2284  //if (isA<triSurfaceMesh>(s) && !isA<distributedTriSurfaceMesh>(s))
2285  //{
2286  // const triSurfaceMesh& surf = refCast<const triSurfaceMesh>(s);
2287  //}
2288 
2289  //nIters = max(nIters, (2<<(mLevel-bLevel)));
2290  }
2291  }
2292 
2293  // Clever limiting of the number of iterations (= max cells in the channel)
2294  // since it has too many problematic issues, e.g. with volume refinement
2295  // and the real check uses the proper distance anyway just disable.
2296  const label nIters = mesh_.globalData().nTotalCells();
2297 
2298 
2299  // Collect candidate faces (i.e. intersecting any surface and
2300  // owner/neighbour not yet refined)
2301  const labelList testFaces(getRefineCandidateFaces(refineCell));
2302 
2303  // Collect segments
2304  pointField start(testFaces.size());
2305  pointField end(testFaces.size());
2306  labelList minLevel(testFaces.size());
2307 
2308  calcCellCellRays
2309  (
2310  neiCc,
2311  neiLevel,
2312  testFaces,
2313  start,
2314  end,
2315  minLevel
2316  );
2317  // TBD. correct nIters for actual cellLevel (since e.g. volume refinement
2318  // might add to cell level). Unfortunately we only have minLevel,
2319  // not maxLevel. Problem: what if volume refinement only in middle of
2320  // channel? Say channel 1m wide with a 0.1m sphere of refinement
2321  // Workaround: have dummy surface with e.g. maxLevel 100 to
2322  // force nIters to be high enough.
2323 
2324 
2325  // Test for all intersections (with surfaces of higher gap level than
2326  // minLevel) and cache per cell the max surface level and the local normal
2327  // on that surface.
2328 
2329  labelList surface1;
2330  List<pointIndexHit> hit1;
2331  labelList region1;
2332  vectorField normal1;
2333 
2334  labelList surface2;
2335  List<pointIndexHit> hit2;
2336  labelList region2;
2337  vectorField normal2;
2338 
2339  surfaces_.findNearestIntersection
2340  (
2341  blockedSurfaces,
2342  start,
2343  end,
2344 
2345  surface1,
2346  hit1,
2347  region1, // local region
2348  normal1,
2349 
2350  surface2,
2351  hit2,
2352  region2, // local region
2353  normal2
2354  );
2355 
2356 
2357  // Detect cells that are using multiple surface regions
2358  //bitSet isMultiRegion;
2359  //detectMultiRegionCells
2360  //(
2361  // faceZones,
2362  // testFaces,
2363  //
2364  // surface1,
2365  // hit1,
2366  // region1,
2367  //
2368  // surface2,
2369  // hit2,
2370  // region2,
2371  //
2372  // isMultiRegion
2373  //);
2374 
2375 
2376  label n = 0;
2377  forAll(testFaces, i)
2378  {
2379  if (hit1[i].hit())
2380  {
2381  n++;
2382  }
2383  }
2384 
2385  List<wallPoints> faceDist(n);
2386  labelList changedFaces(n);
2387  n = 0;
2388 
2389  DynamicList<point> originLocation(2);
2390  DynamicList<scalar> originDistSqr(2);
2391  DynamicList<FixedList<label, 3>> originSurface(2);
2392  //DynamicList<point> originNormal(2);
2393 
2394 
2395  //- To avoid walking through surfaces we mark all faces that have been
2396  // intersected. We can either mark only those faces intersecting
2397  // blockedSurfaces (i.e. with a 'blockLevel') or mark faces intersecting
2398  // any (refinement) surface (this includes e.g. faceZones). This is
2399  // much easier since that information is already cached
2400  // (meshRefinement::intersectedFaces())
2401 
2402  //bitSet isBlockedFace(mesh_.nFaces());
2403  forAll(testFaces, i)
2404  {
2405  if (hit1[i].hit())
2406  {
2407  const label facei = testFaces[i];
2408  //isBlockedFace.set(facei);
2409  const point& fc = mesh_.faceCentres()[facei];
2410  const labelList& fz1 = faceZones[surface1[i]];
2411 
2412  originLocation.clear();
2413  originDistSqr.clear();
2414  //originNormal.clear();
2415  originSurface.clear();
2416 
2417  originLocation.append(hit1[i].point());
2418  originDistSqr.append(hit1[i].point().distSqr(fc));
2419  //originNormal.append(normal1[i]);
2420  originSurface.append
2421  (
2422  FixedList<label, 3>
2423  ({
2424  surface1[i],
2425  region1[i],
2426  (fz1.size() ? fz1[hit1[i].index()] : region1[i])
2427  })
2428  );
2429 
2430  if (hit2[i].hit() && hit1[i] != hit2[i])
2431  {
2432  const labelList& fz2 = faceZones[surface2[i]];
2433  originLocation.append(hit2[i].point());
2434  originDistSqr.append(hit2[i].point().distSqr(fc));
2435  //originNormal.append(normal2[i]);
2436  originSurface.append
2437  (
2438  FixedList<label, 3>
2439  ({
2440  surface2[i],
2441  region2[i],
2442  (fz2.size() ? fz2[hit2[i].index()] : region2[i])
2443  })
2444  );
2445  }
2446 
2447  // Collect all seed data. Currently walking does not look at
2448  // surface direction - if so pass in surface normal as well
2449  faceDist[n] = wallPoints
2450  (
2451  originLocation, // origin
2452  originDistSqr, // distance to origin
2453  originSurface // surface+region+zone
2454  //originNormal // normal at origin
2455  );
2456  changedFaces[n] = facei;
2457  n++;
2458  }
2459  }
2460 
2461 
2462  // Clear intersection info
2463  surface1.clear();
2464  hit1.clear();
2465  region1.clear();
2466  normal1.clear();
2467  surface2.clear();
2468  hit2.clear();
2469  region2.clear();
2470  normal2.clear();
2471 
2472 
2473  List<wallPoints> allFaceInfo(mesh_.nFaces());
2474  List<wallPoints> allCellInfo(mesh_.nCells());
2475 
2476  // Any refinement surface (even a faceZone) should stop the gap walking.
2477  // This is exactly the information which is cached in the surfaceIndex_
2478  // field.
2479  const bitSet isBlockedFace(intersectedFaces());
2480 
2481  wallPoints::trackData td(isBlockedFace, regionToBlockSize);
2482  FaceCellWave<wallPoints, wallPoints::trackData> wallDistCalc
2483  (
2484  mesh_,
2485  changedFaces,
2486  faceDist,
2487  allFaceInfo,
2488  allCellInfo,
2489  0, // max iterations
2490  td
2491  );
2492  wallDistCalc.iterate(nIters);
2493 
2494 
2495  if (debug&meshRefinement::MESH)
2496  {
2497  // Dump current nearest opposite surfaces
2498  volScalarField distance
2499  (
2500  IOobject
2501  (
2502  "gapSize",
2503  mesh_.time().timeName(),
2504  mesh_,
2505  IOobject::NO_READ,
2506  IOobject::NO_WRITE,
2507  IOobject::NO_REGISTER
2508  ),
2509  mesh_,
2510  dimensionedScalar
2511  (
2512  "zero",
2513  dimLength, //dimensionSet(0, 1, 0, 0, 0),
2514  -1
2515  )
2516  );
2517 
2518  forAll(allCellInfo, celli)
2519  {
2520  if (allCellInfo[celli].valid(wallDistCalc.data()))
2521  {
2522  const point& cc = mesh_.cellCentres()[celli];
2523  // Nearest surface points
2524  const List<point>& origin = allCellInfo[celli].origin();
2525 
2526  // Find 'opposite' pair with minimum distance
2527  for (label i = 0; i < origin.size(); i++)
2528  {
2529  for (label j = i + 1; j < origin.size(); j++)
2530  {
2531  if (((cc-origin[i]) & (cc-origin[j])) < 0)
2532  {
2533  const scalar d(mag(origin[i]-origin[j]));
2534  if (distance[celli] < 0)
2535  {
2536  distance[celli] = d;
2537  }
2538  else
2539  {
2540  distance[celli] = min(distance[celli], d);
2541  }
2542  }
2543  }
2544  }
2545  }
2546  }
2547  distance.correctBoundaryConditions();
2548 
2549  Info<< "Writing measured gap distance to "
2550  << distance.name() << endl;
2551  distance.write();
2552  }
2553 
2554 
2555 
2556  // Detect tight gaps:
2557  // - cell is inbetween the two surfaces
2558  // - two surfaces are planarish
2559  // - two surfaces are not too far apart
2560  // (number of walking iterations is a too-coarse measure)
2561 
2562  scalarField smallGapDistance(mesh_.nCells(), 0.0);
2563  label nMulti = 0;
2564  label nSmallGap = 0;
2565 
2566  //OBJstream str(mesh_.time().timePath()/"multiRegion.obj");
2567 
2568 
2569  forAll(allCellInfo, celli)
2570  {
2571  if (allCellInfo[celli].valid(wallDistCalc.data()))
2572  {
2573  const point& cc = mesh_.cellCentres()[celli];
2574 
2575  const List<point>& origin = allCellInfo[celli].origin();
2576  const List<FixedList<label, 3>>& surface =
2577  allCellInfo[celli].surface();
2578 
2579  // Find pair with minimum distance
2580  for (label i = 0; i < origin.size(); i++)
2581  {
2582  for (label j = i + 1; j < origin.size(); j++)
2583  {
2584  //if (isMultiRegion[celli])
2585  //{
2586  // // The intersection locations are too inaccurate
2587  // // (since not proper nearest, just a cell-cell ray
2588  // // intersection) so include always
2589  //
2590  // smallGapDistance[celli] =
2591  // max(smallGapDistance[celli], maxDist);
2592  //
2593  //
2594  // str.writeLine(cc, origin[i]);
2595  // str.writeLine(cc, origin[j]);
2596  //
2597  // nMulti++;
2598  //}
2599  //else
2600  if (((cc-origin[i]) & (cc-origin[j])) < 0)
2601  {
2602  const label surfi = surface[i][0];
2603  const label regioni = surface[i][1];
2604 
2605  const label surfj = surface[j][0];
2606  const label regionj = surface[j][1];
2607 
2608  const scalar maxSize = max
2609  (
2610  regionToBlockSize[surfi][regioni],
2611  regionToBlockSize[surfj][regionj]
2612  );
2613 
2614  if
2615  (
2616  magSqr(origin[i]-origin[j])
2617  < Foam::sqr(2*maxSize)
2618  )
2619  {
2620  const scalar maxDist
2621  (
2622  max
2623  (
2624  mag(cc-origin[i]),
2625  mag(cc-origin[j])
2626  )
2627  );
2628 
2629  smallGapDistance[celli] =
2630  max(smallGapDistance[celli], maxDist);
2631  nSmallGap++;
2632  }
2633  }
2634  }
2635  }
2636  }
2637  }
2638 
2639 
2640  if (debug)
2641  {
2642  Info<< "Marked for blocking due to intersecting multiple surfaces : "
2643  << returnReduce(nMulti, sumOp<label>()) << " cells." << endl;
2644  Info<< "Marked for blocking due to close opposite surfaces : "
2645  << returnReduce(nSmallGap, sumOp<label>()) << " cells." << endl;
2646  }
2647 
2648  if (debug&meshRefinement::MESH)
2649  {
2650  volScalarField distance
2651  (
2652  IOobject
2653  (
2654  "smallGapDistance",
2655  mesh_.time().timeName(),
2656  mesh_,
2657  IOobject::NO_READ,
2658  IOobject::NO_WRITE,
2659  IOobject::NO_REGISTER
2660  ),
2661  mesh_,
2662  dimensionedScalar(dimLength, Zero)
2663  );
2664  distance.field() = smallGapDistance;
2665  distance.correctBoundaryConditions();
2666 
2667  Info<< "Writing all small-gap cells to "
2668  << distance.name() << endl;
2669  distance.write();
2670  }
2671 
2672 
2673  // Mark refinement
2674  const label oldNRefine = nRefine;
2675  forAll(smallGapDistance, celli)
2676  {
2677  if (smallGapDistance[celli] > SMALL)
2678  {
2679  if
2680  (
2681  !markForRefine
2682  (
2683  0, // mark level
2684  nAllowRefine,
2685  refineCell[celli],
2686  nRefine
2687  )
2688  )
2689  {
2690  if (debug)
2691  {
2692  Pout<< "Stopped refining since reaching my cell"
2693  << " limit of " << mesh_.nCells()+7*nRefine
2694  << endl;
2695  }
2696  break;
2697  }
2698  }
2699  }
2700 
2701  if
2702  (
2703  returnReduce(nRefine, sumOp<label>())
2704  > returnReduce(nAllowRefine, sumOp<label>())
2705  )
2706  {
2707  Info<< "Reached refinement limit." << endl;
2708  }
2709 
2710  return returnReduce(nRefine-oldNRefine, sumOp<label>());
2711 }
2712 */
2713 
2714 // ************************************************************************* //
Foam::regionSplit
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition: regionSplit.H:136
Foam::syncTools::syncFaceList
static void syncFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop, const bool parRun=UPstream::parRun())
Synchronize values on all mesh faces.
Definition: syncTools.H:465
Foam::meshRefinement::markOutsideFaces
void markOutsideFaces(const labelList &cellLevel, const labelList &neiLevel, const labelList &refineCell, bitSet &isOutsideFace) const
Mark faces on interface between set and rest.
Definition: meshRefinementBlock.C:171
Foam::UList::size
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:116
forAll
forAll(testFaces, i)
Definition: meshRefinementBlock.C:2402
Foam::ZoneMesh::clearAddressing
void clearAddressing()
Clear addressing.
Definition: ZoneMesh.C:905
Foam::direction
uint8_t direction
Definition: direction.H:46
Foam::bitSet::set
void set(const bitSet &bitset)
Set specified bits from another bitset.
Definition: bitSetI.H:498
allCellInfo
List< wallPoints > allCellInfo(mesh_.nCells())
Foam::cellList
List< cell > cellList
List of cell.
Definition: cellListFwd.H:39
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:125
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::List::transfer
void transfer(List< T > &list)
Transfer the contents of the argument List into this list and annul the argument list.
Definition: List.C:326
Foam::FatalError
error FatalError
Error stream (stdout output on all processes), with additional &#39;FOAM FATAL ERROR&#39; header text and sta...
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:608
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition: List.H:521
Foam::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1122
isBlockedFace
const bitSet isBlockedFace(intersectedFaces())
Foam::mapDistribute::distribute
void distribute(List< T > &fld, const bool dummyTransform=true, const int tag=UPstream::msgType()) const
Distribute List data using default commsType, default flip/negate operator.
Definition: mapDistributeTemplates.C:172
Foam::syncTools::swapBoundaryCellList
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData, const bool parRun=UPstream::parRun())
Extract and swap to obtain neighbour cell values for all boundary faces.
Definition: syncToolsTemplates.C:1404
td
wallPoints::trackData td(isBlockedFace, regionToBlockSize)
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:40
Foam::List::set
std::enable_if< std::is_same< bool, TypeT >::value, bool >::type set(const label i, bool val=true)
A bitSet::set() method for a list of bool.
Definition: List.H:493
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:44
unitConversion.H
Unit conversion functions.
Foam::nl
constexpr char nl
The newline &#39;\n&#39; character (0x0a)
Definition: Ostream.H:50
Foam::meshRefinement::removeGapCells
autoPtr< mapPolyMesh > removeGapCells(const scalar planarAngle, const labelList &minSurfaceLevel, const labelList &globalToMasterPatch, const label growIter)
Detect gapRefinement cells and remove them.
Definition: meshRefinementBlock.C:1018
Foam::PackedList::setSize
void setSize(const label n, unsigned int val=0u)
Alias for resize()
Definition: PackedList.H:807
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: IndirectList.H:65
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:531
Foam::meshRefinement::writeLevel
static writeType writeLevel()
Get/set write level.
Definition: meshRefinement.C:3775
Foam::findIndices
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
Foam::meshRefinement::countFaceDirs
label countFaceDirs(const bitSet &isOutsideFace, const label celli) const
Count number of faces on cell that are in set.
Definition: meshRefinementBlock.C:229
Foam::expressions::Detail::nComponents
::Foam::direction nComponents(const expressions::valueTypeCode) noexcept
The number of components associated with given valueTypeCode.
Definition: exprTraits.C:40
Foam::searchableSurface
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
Definition: searchableSurface.H:65
Foam::primitiveMesh::nFaces
label nFaces() const noexcept
Number of mesh faces.
Definition: primitiveMeshI.H:83
Foam::returnReduce
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
Definition: PstreamReduceOps.H:444
Foam::labelListList
List< labelList > labelListList
List of labelList.
Definition: labelList.H:38
Foam::Time
Class to control time during OpenFOAM simulations that is also the top-level objectRegistry.
Definition: Time.H:69
Foam::PrimitivePatch::meshPoints
const labelList & meshPoints() const
Return labelList of mesh points in patch.
Definition: PrimitivePatch.C:323
oldNRefine
const label oldNRefine
Definition: meshRefinementBlock.C:2673
Foam::HashTable::insert
bool insert(const Key &key, const T &obj)
Copy insert a new entry, not overwriting existing entries.
Definition: HashTableI.H:152
timeName
word timeName
Definition: getTimeIndex.H:3
Foam::vtk::write
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
Definition: foamVtkOutputTemplates.C:29
Foam::PrimitivePatch
A list of faces which address into the list of points.
Definition: PrimitivePatch.H:70
Foam::globalIndex
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:61
Foam::UPstream::nProcs
static label nProcs(const label communicator=worldComm)
Number of ranks in parallel run (for given communicator). It is 1 for serial run. ...
Definition: UPstream.H:1077
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:38
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:258
Foam::name
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
Definition: exprTraits.C:127
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
Foam::edge
An edge is a list of two vertex labels. This can correspond to a directed graph edge or an edge on a ...
Definition: edge.H:59
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:130
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:46
Foam::UPtrList::size
label size() const noexcept
The number of entries in the list.
Definition: UPtrListI.H:106
Foam::meshRefinement::blockLeakFaces
autoPtr< mapPolyMesh > blockLeakFaces(const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const pointField &locationsInMesh, const wordList &zonesInMesh, const pointField &locationsOutsideMesh, const labelList &selectedSurfaces)
Baffle faces to break any leak from inside to outside.
Definition: meshRefinementBlock.C:1888
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1116
Foam::meshRefinement::growSet
void growSet(const labelList &neiLevel, const bitSet &isOutsideFace, labelList &refineCell, label &nRefine) const
Add one layer of cells to set.
Definition: meshRefinementBlock.C:280
Foam::PrimitivePatch::edges
const edgeList & edges() const
Return list of edges, address into LOCAL point list.
Definition: PrimitivePatch.C:176
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const noexcept
Number of internal faces.
Definition: primitiveMeshI.H:71
Foam::vector
Vector< scalar > vector
Definition: vector.H:57
Foam::refineCell
Container with cells to refine. Refinement given as single direction.
Definition: refineCell.H:52
wallDistCalc
FaceCellWave< wallPoints, wallPoints::trackData > wallDistCalc(mesh_, changedFaces, faceDist, allFaceInfo, allCellInfo, 0, td)
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:26
Foam::surfaceZonesInfo::getUnnamedSurfaces
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
Definition: surfaceZonesInfo.C:257
Foam::syncTools::syncPointList
static void syncPointList(const polyMesh &mesh, List< T > &pointValues, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
Definition: syncToolsTemplates.C:735
Foam::polyBoundaryMesh
A polyBoundaryMesh is a polyPatch list with additional search methods and registered IO...
Definition: polyBoundaryMesh.H:60
Foam::volumeType::MIXED
A location that is partly inside and outside.
Definition: volumeType.H:67
Foam::ensightOutput::debug
int debug
Static debugging option.
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:201
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:68
Foam::point
vector point
Point is a vector.
Definition: point.H:37
Foam::bitSet
A bitSet stores bits (elements with only two states) in packed internal format and supports a variety...
Definition: bitSet.H:59
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:450
Foam::HashTable::sortedToc
List< Key > sortedToc() const
The table of contents (the keys) in sorted order.
Definition: HashTable.C:163
Foam::UIndirectList
A List with indirect addressing. Like IndirectList but does not store addressing. ...
Definition: faMatrix.H:52
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:54
Foam::meshRefinement::selectIntersectedFaces
void selectIntersectedFaces(const labelList &surfaces, boolList &isBlockedFace) const
Faces currently on boundary or intersected by surface.
Definition: meshRefinementBlock.C:1206
Foam::refinementParameters::zonePoints
static List< pointField > zonePoints(const pointField &locationsInMesh, const wordList &zonesInMesh, const pointField &locationsOutsideMesh)
Helper: per zone (entry in zonesInMesh) the locations with.
Definition: refinementParameters.C:306
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
Foam::labelMax
constexpr label labelMax
Definition: label.H:55
Foam::meshRefinement::writeType
writeType
Enumeration for what to write. Used as a bit-pattern.
Definition: meshRefinement.H:114
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:48
Foam::holeToFace::calcClosure
static autoPtr< mapDistribute > calcClosure(const polyMesh &mesh, const List< pointField > &zonePoints, const labelList &blockedFaces, const globalIndex &globalBlockedFaces, const bool erode, labelList &closureFaces, labelList &closureToBlocked)
Optional direct use to generate the set of faces and the method to.
Definition: holeToFace.C:1190
Foam::meshRefinement::debugType
debugType
Enumeration for what to debug. Used as a bit-pattern.
Definition: meshRefinement.H:92
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:62
Foam::polyPatch
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:69
Foam::polyBoundaryMesh::nNonProcessor
label nNonProcessor() const
The number of patches before the first processor patch.
Definition: polyBoundaryMesh.C:608
s
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
Definition: gmvOutputSpray.H:25
Foam::primitiveMesh::nBoundaryFaces
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
Definition: primitiveMeshI.H:77
Foam::returnReduceOr
bool returnReduceOr(const bool value, const label comm=UPstream::worldComm)
Perform logical (or) MPI Allreduce on a copy. Uses UPstream::reduceOr.
Definition: PstreamReduceOps.H:480
Foam::degToRad
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Definition: unitConversion.H:44
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:60
allFaceInfo
List< wallPoints > allFaceInfo(mesh_.nFaces())
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::reduce
void reduce(T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Reduce inplace (cf. MPI Allreduce) using linear/tree communication schedule.
Definition: PstreamReduceOps.H:49
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::surfaceZonesInfo::addFaceZone
static label addFaceZone(const word &name, const labelList &addressing, const boolList &flipMap, polyMesh &mesh)
Definition: surfaceZonesInfo.C:521
pp
uindirectPrimitivePatch pp(UIndirectList< face >(mesh.faces(), faceLabels), mesh.points())
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127