meshRefinementGapRefine.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) 2015 OpenFOAM Foundation
9  Copyright (C) 2015-2022 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "meshRefinement.H"
30 #include "Time.H"
31 #include "refinementSurfaces.H"
32 #include "refinementFeatures.H"
33 #include "shellSurfaces.H"
34 #include "triSurfaceMesh.H"
35 #include "treeDataCell.H"
36 #include "searchableSurfaces.H"
37 #include "DynamicField.H"
38 #include "transportData.H"
39 #include "FaceCellWave.H"
40 #include "volFields.H"
41 #include "zeroGradientFvPatchFields.H"
42 
43 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
44 
45 Foam::label Foam::meshRefinement::markSurfaceGapRefinement
46 (
47  const scalar planarCos,
48 
49  const label nAllowRefine,
50  const labelList& neiLevel,
51  const pointField& neiCc,
52 
53  labelList& refineCell,
54  label& nRefine
55 ) const
56 {
57  const labelList& cellLevel = meshCutter_.cellLevel();
58  const pointField& cellCentres = mesh_.cellCentres();
59 
60  // Get the gap level for the shells
61  const labelList maxLevel(shells_.maxGapLevel());
62 
63  label oldNRefine = nRefine;
64 
65  if (max(maxLevel) > 0)
66  {
67  // Use cached surfaceIndex_ to detect if any intersection. If so
68  // re-intersect to determine level wanted.
69 
70  // Collect candidate faces
71  // ~~~~~~~~~~~~~~~~~~~~~~~
72 
73  labelList testFaces(getRefineCandidateFaces(refineCell));
74 
75  // Collect segments
76  // ~~~~~~~~~~~~~~~~
77 
78  pointField start(testFaces.size());
79  pointField end(testFaces.size());
80  {
81  labelList minLevel(testFaces.size());
82  calcCellCellRays
83  (
84  neiCc,
85  neiLevel,
86  testFaces,
87  start,
88  end,
89  minLevel
90  );
91  }
92 
93 
94  // Collect cells to test for inside/outside in shell
95  labelList cellToCompact(mesh_.nCells(), -1);
96  labelList bFaceToCompact(mesh_.nBoundaryFaces(), -1);
97  labelList gapShell;
98  List<FixedList<label, 3>> shellGapInfo;
99  List<volumeType> shellGapMode;
100  {
101  DynamicField<point> compactToCc(mesh_.nCells()/10);
102  DynamicList<label> compactToLevel(compactToCc.capacity());
103  forAll(testFaces, i)
104  {
105  label faceI = testFaces[i];
106  label own = mesh_.faceOwner()[faceI];
107  if (cellToCompact[own] == -1)
108  {
109  cellToCompact[own] = compactToCc.size();
110  compactToCc.append(cellCentres[own]);
111  compactToLevel.append(cellLevel[own]);
112  }
113  if (mesh_.isInternalFace(faceI))
114  {
115  label nei = mesh_.faceNeighbour()[faceI];
116  if (cellToCompact[nei] == -1)
117  {
118  cellToCompact[nei] = compactToCc.size();
119  compactToCc.append(cellCentres[nei]);
120  compactToLevel.append(cellLevel[nei]);
121  }
122  }
123  else
124  {
125  label bFaceI = faceI - mesh_.nInternalFaces();
126  if (bFaceToCompact[bFaceI] == -1)
127  {
128  bFaceToCompact[bFaceI] = compactToCc.size();
129  compactToCc.append(neiCc[bFaceI]);
130  compactToLevel.append(neiLevel[bFaceI]);
131  }
132  }
133  }
134 
135  shells_.findHigherGapLevel
136  (
137  compactToCc,
138  compactToLevel,
139 
140  gapShell,
141  shellGapInfo,
142  shellGapMode
143  );
144  }
145 
146 
147  //const fileName dir(mesh_.time().path()/timeName());
148  //if (debug)
149  //{
150  // mkDir(dir);
151  // OBJstream insideStr(dir/"insideShell.obj");
152  // OBJstream outsideStr(dir/"outsideShell.obj");
153  // Pout<< "Writing points to:" << nl
154  // << " inside : " << insideStr.name() << nl
155  // << " outside: " << outsideStr.name() << nl
156  // << endl;
157  //
158  // forAll(cellToCompact, celli)
159  // {
160  // const label compacti = cellToCompact[celli];
161  //
162  // if (compacti != -1)
163  // {
164  // if (gapShell[compacti] != -1)
165  // {
166  // insideStr.write(mesh_.cellCentres()[celli]);
167  // }
168  // else
169  // {
170  // outsideStr.write(mesh_.cellCentres()[celli]);
171  // }
172  // }
173  // }
174  // forAll(bFaceToCompact, bFacei)
175  // {
176  // const label compacti = bFaceToCompact[bFacei];
177  // if (compacti != -1)
178  // {
179  // if (gapShell[compacti] != -1)
180  // {
181  // insideStr.write(neiCc[bFacei]);
182  // }
183  // else
184  // {
185  // outsideStr.write(neiCc[bFacei]);
186  // }
187  // }
188  // }
189  //}
190 
191 
192  const List<FixedList<label, 3>>& extendedGapLevel =
193  surfaces_.extendedGapLevel();
194  const List<volumeType>& extendedGapMode =
195  surfaces_.extendedGapMode();
196  const boolList& extendedGapSelf = surfaces_.gapSelf();
197 
198  labelList ccSurface1;
199  List<pointIndexHit> ccHit1;
200  labelList ccRegion1;
201  vectorField ccNormal1;
202  {
203  labelList ccSurface2;
204  List<pointIndexHit> ccHit2;
205  labelList ccRegion2;
206  vectorField ccNormal2;
207 
208  surfaces_.findNearestIntersection
209  (
210  identity(surfaces_.surfaces().size()),
211  start,
212  end,
213 
214  ccSurface1,
215  ccHit1,
216  ccRegion1,
217  ccNormal1,
218 
219  ccSurface2,
220  ccHit2,
221  ccRegion2,
222  ccNormal2
223  );
224  }
225 
226  start.clear();
227  end.clear();
228 
229  DynamicField<point> rayStart(2*ccSurface1.size());
230  DynamicField<point> rayEnd(2*ccSurface1.size());
231  DynamicField<scalar> gapSize(2*ccSurface1.size());
232 
233  DynamicField<point> rayStart2(2*ccSurface1.size());
234  DynamicField<point> rayEnd2(2*ccSurface1.size());
235  DynamicField<scalar> gapSize2(2*ccSurface1.size());
236 
237  DynamicList<label> cellMap(2*ccSurface1.size());
238  DynamicList<label> compactMap(2*ccSurface1.size());
239 
240  forAll(ccSurface1, i)
241  {
242  label surfI = ccSurface1[i];
243 
244  if (surfI != -1)
245  {
246  label globalRegionI =
247  surfaces_.globalRegion(surfI, ccRegion1[i]);
248 
249  label faceI = testFaces[i];
250  const point& surfPt = ccHit1[i].hitPoint();
251 
252  label own = mesh_.faceOwner()[faceI];
253  if
254  (
255  cellToCompact[own] != -1
256  && shellGapInfo[cellToCompact[own]][2] > 0
257  )
258  {
259  // Combine info from shell and surface
260  label compactI = cellToCompact[own];
261  FixedList<label, 3> gapInfo;
262  volumeType gapMode;
263  mergeGapInfo
264  (
265  shellGapInfo[compactI],
266  shellGapMode[compactI],
267  extendedGapLevel[globalRegionI],
268  extendedGapMode[globalRegionI],
269 
270  gapInfo,
271  gapMode
272  );
273 
274  const point& cc = cellCentres[own];
275  label nRays = generateRays
276  (
277  false,
278  surfPt,
279  ccNormal1[i],
280  gapInfo,
281  gapMode,
282  surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
283  cellLevel[own],
284 
285  rayStart,
286  rayEnd,
287  gapSize,
288 
289  rayStart2,
290  rayEnd2,
291  gapSize2
292  );
293  for (label j = 0; j < nRays; j++)
294  {
295  cellMap.append(own);
296  compactMap.append(i);
297  }
298  }
299  if (mesh_.isInternalFace(faceI))
300  {
301  label nei = mesh_.faceNeighbour()[faceI];
302  if
303  (
304  cellToCompact[nei] != -1
305  && shellGapInfo[cellToCompact[nei]][2] > 0
306  )
307  {
308  // Combine info from shell and surface
309  label compactI = cellToCompact[nei];
310  FixedList<label, 3> gapInfo;
311  volumeType gapMode;
312  mergeGapInfo
313  (
314  shellGapInfo[compactI],
315  shellGapMode[compactI],
316  extendedGapLevel[globalRegionI],
317  extendedGapMode[globalRegionI],
318 
319  gapInfo,
320  gapMode
321  );
322 
323  const point& cc = cellCentres[nei];
324  label nRays = generateRays
325  (
326  false,
327  surfPt,
328  ccNormal1[i],
329  gapInfo,
330  gapMode,
331  surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
332  cellLevel[nei],
333 
334  rayStart,
335  rayEnd,
336  gapSize,
337 
338  rayStart2,
339  rayEnd2,
340  gapSize2
341  );
342  for (label j = 0; j < nRays; j++)
343  {
344  cellMap.append(nei);
345  compactMap.append(i);
346  }
347  }
348  }
349  else
350  {
351  // Note: on coupled face. What cell are we going to
352  // refine? We've got the neighbouring cell centre
353  // and level but we cannot mark it for refinement on
354  // this side...
355  label bFaceI = faceI - mesh_.nInternalFaces();
356 
357  if
358  (
359  bFaceToCompact[bFaceI] != -1
360  && shellGapInfo[bFaceToCompact[bFaceI]][2] > 0
361  )
362  {
363  // Combine info from shell and surface
364  label compactI = bFaceToCompact[bFaceI];
365  FixedList<label, 3> gapInfo;
366  volumeType gapMode;
367  mergeGapInfo
368  (
369  shellGapInfo[compactI],
370  shellGapMode[compactI],
371  extendedGapLevel[globalRegionI],
372  extendedGapMode[globalRegionI],
373 
374  gapInfo,
375  gapMode
376  );
377 
378  const point& cc = neiCc[bFaceI];
379  label nRays = generateRays
380  (
381  false,
382  surfPt,
383  ccNormal1[i],
384  gapInfo,
385  gapMode,
386  surfPt+((cc-surfPt)&ccNormal1[i])*ccNormal1[i],
387  neiLevel[bFaceI],
388 
389  rayStart,
390  rayEnd,
391  gapSize,
392 
393  rayStart2,
394  rayEnd2,
395  gapSize2
396  );
397  for (label j = 0; j < nRays; j++)
398  {
399  cellMap.append(-1); // See above.
400  compactMap.append(i);
401  }
402  }
403  }
404  }
405  }
406 
407  Info<< "Shooting " << returnReduce(rayStart.size(), sumOp<label>())
408  << " rays from " << returnReduce(testFaces.size(), sumOp<label>())
409  << " intersected faces" << endl;
410 
411  rayStart.shrink();
412  rayEnd.shrink();
413  gapSize.shrink();
414 
415  rayStart2.shrink();
416  rayEnd2.shrink();
417  gapSize2.shrink();
418 
419  cellMap.shrink();
420  compactMap.shrink();
421 
422  testFaces.clear();
423  ccSurface1.clear();
424  ccHit1.clear();
425  ccRegion1.clear();
426  ccNormal1 = UIndirectList<vector>(ccNormal1, compactMap)();
427 
428 
429  // Do intersections in pairs
430  labelList surf1;
431  List<pointIndexHit> hit1;
432  vectorField normal1;
433  surfaces_.findNearestIntersection
434  (
435  rayStart,
436  rayEnd,
437  surf1,
438  hit1,
439  normal1
440  );
441 
442  labelList surf2;
443  List<pointIndexHit> hit2;
444  vectorField normal2;
445  surfaces_.findNearestIntersection
446  (
447  rayStart2,
448  rayEnd2,
449  surf2,
450  hit2,
451  normal2
452  );
453 
454  forAll(surf1, i)
455  {
456  // Combine selfProx of shell and surfaces.
457  // Ignore regions for now
458  const label cellI = cellMap[i];
459 
460  const label shelli =
461  (
462  (cellI != -1 && cellToCompact[cellI] != -1)
463  ? gapShell[cellToCompact[cellI]]
464  : -1
465  );
466 
467  bool selfProx = true;
468  if (shelli != -1)
469  {
470  selfProx = shells_.gapSelf()[shelli][0];
471  }
472  if (surf1[i] != -1 && selfProx)
473  {
474  const label globalRegioni = surfaces_.globalRegion(surf1[i], 0);
475  selfProx = extendedGapSelf[globalRegioni];
476  }
477 
478  if
479  (
480  surf1[i] != -1
481  && surf2[i] != -1
482  && (surf2[i] != surf1[i] || selfProx)
483  )
484  {
485  // Found intersection with surface. Check opposite normal.
486  if
487  (
488  cellI != -1
489  && (mag(normal1[i]&normal2[i]) > planarCos)
490  && (
491  hit1[i].hitPoint().distSqr(hit2[i].hitPoint())
492  < Foam::sqr(gapSize[i])
493  )
494  )
495  {
496  if
497  (
498  !markForRefine
499  (
500  surf1[i],
501  nAllowRefine,
502  refineCell[cellI],
503  nRefine
504  )
505  )
506  {
507  break;
508  }
509  }
510  }
511  }
512 
513  if
514  (
515  returnReduce(nRefine, sumOp<label>())
516  > returnReduce(nAllowRefine, sumOp<label>())
517  )
518  {
519  Info<< "Reached refinement limit." << endl;
520  }
521  }
522 
523  return returnReduce(nRefine-oldNRefine, sumOp<label>());
524 }
525 
526 
527 //Foam::meshRefinement::findNearestOppositeOp::findNearestOppositeOp
528 //(
529 // const indexedOctree<treeDataTriSurface>& tree,
530 // const point& oppositePoint,
531 // const vector& oppositeNormal,
532 // const scalar minCos
533 //)
534 //:
535 // tree_(tree),
536 // oppositePoint_(oppositePoint),
537 // oppositeNormal_(oppositeNormal),
538 // minCos_(minCos)
539 //{}
540 //
541 //
542 //void Foam::meshRefinement::findNearestOppositeOp::operator()
543 //(
544 // const labelUList& indices,
545 // const point& sample,
546 // scalar& nearestDistSqr,
547 // label& minIndex,
548 // point& nearestPoint
549 //) const
550 //{
551 // const treeDataTriSurface& shape = tree_.shapes();
552 // const triSurface& patch = shape.patch();
553 // const pointField& points = patch.points();
554 //
555 // forAll(indices, i)
556 // {
557 // const label index = indices[i];
558 // const labelledTri& f = patch[index];
559 //
560 // pointHit nearHit = f.nearestPoint(sample, points);
561 // scalar distSqr = sqr(nearHit.distance());
562 //
563 // if (distSqr < nearestDistSqr)
564 // {
565 // // Nearer. Check if
566 // // - a bit way from other hit
567 // // - in correct search cone
568 // vector d(nearHit.point()-oppositePoint_);
569 // scalar normalDist(d&oppositeNormal_);
570 //
571 // if (normalDist > Foam::sqr(SMALL) && normalDist/mag(d) > minCos_)
572 // {
573 // nearestDistSqr = distSqr;
574 // minIndex = index;
575 // nearestPoint = nearHit.point();
576 // }
577 // }
578 // }
579 //}
580 //
581 //
582 //void Foam::meshRefinement::searchCone
583 //(
584 // const label surfI,
585 // labelList& nearMap, // cells
586 // scalarField& nearGap, // gap size
587 // List<pointIndexHit>& nearInfo, // nearest point on surface
588 // List<pointIndexHit>& oppositeInfo // detected point on gap (or miss)
589 //) const
590 //{
591 // const labelList& cellLevel = meshCutter_.cellLevel();
592 // const pointField& cellCentres = mesh_.cellCentres();
593 // const scalar edge0Len = meshCutter_.level0EdgeLength();
594 //
595 // const labelList& surfaceIndices = surfaces_.surfaces();
596 // const List<FixedList<label, 3>>& extendedGapLevel =
597 // surfaces_.extendedGapLevel();
598 // const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
599 //
600 //
601 // label geomI = surfaceIndices[surfI];
602 // const searchableSurface& geom = surfaces_.geometry()[geomI];
603 //
604 // const triSurfaceMesh& s = refCast<const triSurfaceMesh>(geom);
605 // const indexedOctree<treeDataTriSurface>& tree = s.tree();
606 //
607 //
608 // const scalar searchCos = Foam::cos(degToRad(30.0));
609 //
610 // // Normals for ray shooting and inside/outside detection
611 // vectorField nearNormal;
612 // geom.getNormal(nearInfo, nearNormal);
613 // // Regions
614 // labelList nearRegion;
615 // geom.getRegion(nearInfo, nearRegion);
616 //
617 //
618 // // Now loop over all near points and search in the half cone
619 // labelList map(nearInfo.size());
620 // label compactI = 0;
621 //
622 // oppositeInfo.setSize(nearInfo.size());
623 //
624 // forAll(nearInfo, i)
625 // {
626 // label globalRegionI =
627 // surfaces_.globalRegion(surfI, nearRegion[i]);
628 //
629 // // Get updated gap information now we have the region
630 // label nGapCells = extendedGapLevel[globalRegionI][0];
631 // label minLevel = extendedGapLevel[globalRegionI][1];
632 // label maxLevel = extendedGapLevel[globalRegionI][2];
633 // volumeType mode = extendedGapMode[globalRegionI];
634 //
635 // label cellI = nearMap[i];
636 // label cLevel = cellLevel[cellI];
637 //
638 // if (cLevel >= minLevel && cLevel < maxLevel)
639 // {
640 // scalar cellSize = edge0Len/pow(2.0, cLevel);
641 //
642 // // Update gap size
643 // nearGap[i] = nGapCells*cellSize;
644 //
645 // const point& nearPt = nearInfo[i].hitPoint();
646 // vector v(cellCentres[cellI]-nearPt);
647 // scalar magV = mag(v);
648 //
649 // // Like with ray shooting we want to
650 // // - find triangles up to nearGap away on the wanted side of the
651 // // surface
652 // // - find triangles up to 0.5*cellSize away on the unwanted side
653 // // of the surface. This is for cells straddling the surface
654 // // where
655 // // the cell centre might be on the wrong side of the surface
656 //
657 // // Tbd: check that cell centre is inbetween the gap hits
658 // // (only if the cell is far enough away)
659 //
660 // scalar posNormalSize = 0.0;
661 // scalar negNormalSize = 0.0;
662 //
663 // if (mode == volumeType::OUTSIDE)
664 // {
665 // posNormalSize = nearGap[i];
666 // if (magV < 0.5*cellSize)
667 // {
668 // negNormalSize = 0.5*cellSize;
669 // }
670 // }
671 // else if (mode == volumeType::INSIDE)
672 // {
673 // if (magV < 0.5*cellSize)
674 // {
675 // posNormalSize = 0.5*cellSize;
676 // }
677 // negNormalSize = nearGap[i];
678 // }
679 // else
680 // {
681 // posNormalSize = nearGap[i];
682 // negNormalSize = nearGap[i];
683 // }
684 //
685 // // Test with positive normal
686 // oppositeInfo[compactI] = tree.findNearest
687 // (
688 // nearPt,
689 // sqr(posNormalSize),
690 // findNearestOppositeOp
691 // (
692 // tree,
693 // nearPt,
694 // nearNormal[i],
695 // searchCos
696 // )
697 // );
698 //
699 // if (oppositeInfo[compactI].hit())
700 // {
701 // map[compactI++] = i;
702 // }
703 // else
704 // {
705 // // Test with negative normal
706 // oppositeInfo[compactI] = tree.findNearest
707 // (
708 // nearPt,
709 // sqr(negNormalSize),
710 // findNearestOppositeOp
711 // (
712 // tree,
713 // nearPt,
714 // -nearNormal[i],
715 // searchCos
716 // )
717 // );
718 //
719 // if (oppositeInfo[compactI].hit())
720 // {
721 // map[compactI++] = i;
722 // }
723 // }
724 // }
725 // }
726 //
727 // Info<< "Selected " << returnReduce(compactI, sumOp<label>())
728 // << " hits on the correct side out of "
729 // << returnReduce(map.size(), sumOp<label>()) << endl;
730 // map.setSize(compactI);
731 // oppositeInfo.setSize(compactI);
732 //
733 // nearMap = labelUIndList(nearMap, map)();
734 // nearGap = UIndirectList<scalar>(nearGap, map)();
735 // nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
736 // nearNormal = UIndirectList<vector>(nearNormal, map)();
737 //
738 // // Exclude hits which aren't opposite enough. E.g. you might find
739 // // a point on a perpendicular wall - but this does not constitute a gap.
740 // vectorField oppositeNormal;
741 // geom.getNormal(oppositeInfo, oppositeNormal);
742 //
743 // compactI = 0;
744 // forAll(oppositeInfo, i)
745 // {
746 // if ((nearNormal[i] & oppositeNormal[i]) < -0.707)
747 // {
748 // map[compactI++] = i;
749 // }
750 // }
751 //
752 // Info<< "Selected " << returnReduce(compactI, sumOp<label>())
753 // << " hits opposite the nearest out of "
754 // << returnReduce(map.size(), sumOp<label>()) << endl;
755 // map.setSize(compactI);
756 //
757 // nearMap = labelUIndList(nearMap, map)();
758 // nearGap = UIndirectList<scalar>(nearGap, map)();
759 // nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
760 // oppositeInfo = UIndirectList<pointIndexHit>(oppositeInfo, map)();
761 //}
762 
763 
764 Foam::label Foam::meshRefinement::generateRays
765 (
766  const point& nearPoint,
767  const vector& nearNormal,
768  const FixedList<label, 3>& gapInfo,
769  const volumeType& mode,
770 
771  const label cLevel,
772 
773  DynamicField<point>& start,
774  DynamicField<point>& end
775 ) const
776 {
777  label nOldRays = start.size();
778 
779  if (cLevel >= gapInfo[1] && cLevel < gapInfo[2] && gapInfo[0] > 0)
780  {
781  scalar cellSize = meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
782 
783  // Calculate gap size
784  scalar nearGap = gapInfo[0]*cellSize;
785 
786  const vector& n = nearNormal;
787 
788  // Situation 'C' above: cell too close. Use surface
789  // -normal and -point to shoot rays
790 
791  if (mode == volumeType::OUTSIDE)
792  {
793  start.append(nearPoint+1e-6*n);
794  end.append(nearPoint+nearGap*n);
795  }
796  else if (mode == volumeType::INSIDE)
797  {
798  start.append(nearPoint-1e-6*n);
799  end.append(nearPoint-nearGap*n);
800  }
801  else if (mode == volumeType::MIXED)
802  {
803  start.append(nearPoint+1e-6*n);
804  end.append(nearPoint+nearGap*n);
805 
806  start.append(nearPoint-1e-6*n);
807  end.append(nearPoint-nearGap*n);
808  }
809  }
810 
811  return start.size()-nOldRays;
812 }
813 
814 
815 Foam::label Foam::meshRefinement::generateRays
816 (
817  const bool useSurfaceNormal,
818 
819  const point& nearPoint,
820  const vector& nearNormal,
821  const FixedList<label, 3>& gapInfo,
822  const volumeType& mode,
823 
824  const point& cc,
825  const label cLevel,
826 
827  DynamicField<point>& start,
828  DynamicField<point>& end,
829  DynamicField<scalar>& gapSize,
830 
831  DynamicField<point>& start2,
832  DynamicField<point>& end2,
833  DynamicField<scalar>& gapSize2
834 ) const
835 {
836  // We want to handle the following cases:
837  // - surface: small gap (marked with 'surface'). gap might be
838  // on inside or outside of surface.
839  // - A: cell well inside the gap.
840  // - B: cell well outside the gap.
841  // - C: cell straddling the gap. cell centre might be inside
842  // or outside
843  //
844  // +---+
845  // | B |
846  // +---+
847  //
848  // +------+
849  // | |
850  // | C |
851  // --------|------|----surface
852  // +------+
853  //
854  // +---+
855  // | A |
856  // +---+
857  //
858  //
859  // --------------------surface
860  //
861  // So:
862  // - find nearest point on surface
863  // - in situation A,B decide if on wanted side of surface
864  // - detect if locally a gap (and the cell inside the gap) by
865  // shooting a ray from the point on the surface in the direction
866  // of
867  // - A,B: the cell centre
868  // - C: the surface normal and/or negative surface normal
869  // and see we hit anything
870  //
871  // Variations of this scheme:
872  // - always shoot in the direction of the surface normal. This needs
873  // then an additional check to make sure the cell centre is
874  // somewhere inside the gap
875  // - instead of ray shooting use a 'constrained' nearest search
876  // by e.g. looking inside a search cone (implemented in searchCone).
877  // The problem with this constrained nearest is that it still uses
878  // the absolute nearest point on each triangle and only afterwards
879  // checks if it is inside the search cone.
880 
881 
882  // Decide which near points are good:
883  // - with updated minLevel and maxLevel and nearGap make sure
884  // the cell is still a candidate
885  // NOTE: inside the gap the nearest point on the surface will
886  // be HALF the gap size - otherwise we would have found
887  // a point on the opposite side
888  // - if the mode is both sides
889  // - or if the hit is inside the current cell (situation 'C',
890  // magV < 0.5cellSize)
891  // - or otherwise if on the correct side
892 
893  label nOldRays = start.size();
894 
895  if (cLevel >= gapInfo[1] && cLevel < gapInfo[2] && gapInfo[0] > 0)
896  {
897  scalar cellSize = meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
898 
899  // Calculate gap size
900  scalar nearGap = gapInfo[0]*cellSize;
901 
902  // Distance to nearest
903  vector v(cc-nearPoint);
904  scalar magV = mag(v);
905 
906  if (useSurfaceNormal || magV < 0.5*cellSize)
907  {
908  const vector& n = nearNormal;
909 
910  // Situation 'C' above: cell too close. Use surface
911  // -normal and -point to shoot rays
912 
913  if (mode == volumeType::OUTSIDE)
914  {
915  start.append(nearPoint+1e-6*n);
916  end.append(nearPoint+nearGap*n);
917  gapSize.append(nearGap);
918  // Second vector so we get pairs of intersections
919  start2.append(nearPoint+1e-6*n);
920  end2.append(nearPoint-1e-6*n);
921  gapSize2.append(gapSize.last());
922  }
923  else if (mode == volumeType::INSIDE)
924  {
925  start.append(nearPoint-1e-6*n);
926  end.append(nearPoint-nearGap*n);
927  gapSize.append(nearGap);
928  // Second vector so we get pairs of intersections
929  start2.append(nearPoint-1e-6*n);
930  end2.append(nearPoint+1e-6*n);
931  gapSize2.append(gapSize.last());
932  }
933  else if (mode == volumeType::MIXED)
934  {
935  // Do both rays:
936  // Outside
937  {
938  start.append(nearPoint+1e-6*n);
939  end.append(nearPoint+nearGap*n);
940  gapSize.append(nearGap);
941  // Second vector so we get pairs of intersections
942  start2.append(nearPoint+1e-6*n);
943  end2.append(nearPoint-1e-6*n);
944  gapSize2.append(gapSize.last());
945  }
946  // Inside
947  {
948  start.append(nearPoint-1e-6*n);
949  end.append(nearPoint-nearGap*n);
950  gapSize.append(nearGap);
951  // Second vector so we get pairs of intersections
952  start2.append(nearPoint-1e-6*n);
953  end2.append(nearPoint+1e-6*n);
954  gapSize2.append(gapSize.last());
955  }
956  }
957  }
958  else
959  {
960  // Situation 'A' or 'B' above: cell well away. Test if
961  // cell on correct side of surface and shoot ray through
962  // cell centre. Note: no need to shoot ray in other
963  // direction since we're trying to detect cell inside
964  // the gap.
965 
966  scalar s = (v&nearNormal);
967 
968  if
969  (
970  (mode == volumeType::MIXED)
971  || (mode == volumeType::OUTSIDE && s > SMALL)
972  || (mode == volumeType::INSIDE && s < -SMALL)
973  )
974  {
976  //vector n(v/(magV+ROOTVSMALL));
977  //
978  //start.append(cc);
979  //end.append(cc+nearGap*n);
980  //gapSize.append(nearGap);
981  //
982  //start2.append(cc);
983  //end2.append(cc-nearGap*n);
984  //gapSize2.append(nearGap);
985 
986 
989  //start.append(cc);
990  //end.append(cc+nearGap*vector(1, 0, 0));
991  //gapSize.append(nearGap);
992  //
993  //start2.append(cc);
994  //end2.append(cc-nearGap*vector(1, 0, 0));
995  //gapSize2.append(nearGap);
996  //
998  //start.append(cc);
999  //end.append(cc+nearGap*vector(0, 1, 0));
1000  //gapSize.append(nearGap);
1001  //
1002  //start2.append(cc);
1003  //end2.append(cc-nearGap*vector(0, 1, 0));
1004  //gapSize2.append(nearGap);
1005  //
1007  //start.append(cc);
1008  //end.append(cc+nearGap*vector(0, 0, 1));
1009  //gapSize.append(nearGap);
1010  //
1011  //start2.append(cc);
1012  //end2.append(cc-nearGap*vector(0, 0, 1));
1013  //gapSize2.append(nearGap);
1014 
1015 
1016  // 3 axes aligned with normal
1017 
1018  // Use vector through cell centre
1019  vector n(v/(magV+ROOTVSMALL));
1020 
1021  // Get second vector. Make sure it is sufficiently perpendicular
1022  vector e2(1, 0, 0);
1023  scalar s = (e2 & n);
1024  if (mag(s) < 0.9)
1025  {
1026  e2 -= s*n;
1027  }
1028  else
1029  {
1030  e2 = vector(0, 1, 0);
1031  e2 -= (e2 & n)*n;
1032  }
1033  e2 /= mag(e2);
1034 
1035  // Third vector
1036  vector e3 = n ^ e2;
1037 
1038 
1039  // Rays in first direction
1040  start.append(cc);
1041  end.append(cc+nearGap*n);
1042  gapSize.append(nearGap);
1043 
1044  start2.append(cc);
1045  end2.append(cc-nearGap*n);
1046  gapSize2.append(nearGap);
1047 
1048  // Rays in second direction
1049  start.append(cc);
1050  end.append(cc+nearGap*e2);
1051  gapSize.append(nearGap);
1052 
1053  start2.append(cc);
1054  end2.append(cc-nearGap*e2);
1055  gapSize2.append(nearGap);
1056 
1057  // Rays in third direction
1058  start.append(cc);
1059  end.append(cc+nearGap*e3);
1060  gapSize.append(nearGap);
1061 
1062  start2.append(cc);
1063  end2.append(cc-nearGap*e3);
1064  gapSize2.append(nearGap);
1065  }
1066  }
1067  }
1068 
1069  return start.size()-nOldRays;
1070 }
1071 
1072 
1073 void Foam::meshRefinement::selectGapCandidates
1074 (
1075  const labelList& refineCell,
1076  const label nRefine,
1077 
1078  labelList& cellMap,
1079  labelList& gapShell,
1080  List<FixedList<label, 3>>& shellGapInfo,
1081  List<volumeType>& shellGapMode
1082 ) const
1083 {
1084  const labelList& cellLevel = meshCutter_.cellLevel();
1085  const pointField& cellCentres = mesh_.cellCentres();
1086 
1087  // Collect cells to test
1088  cellMap.setSize(cellLevel.size()-nRefine);
1089  label compactI = 0;
1090 
1091  forAll(cellLevel, cellI)
1092  {
1093  if (refineCell[cellI] == -1)
1094  {
1095  cellMap[compactI++] = cellI;
1096  }
1097  }
1098  Info<< "Selected " << returnReduce(compactI, sumOp<label>())
1099  << " unmarked cells out of "
1100  << mesh_.globalData().nTotalCells() << endl;
1101  cellMap.setSize(compactI);
1102 
1103  // Do test to see whether cells are inside/outside shell with
1104  // applicable specification (minLevel <= celllevel < maxLevel)
1105  shells_.findHigherGapLevel
1106  (
1107  pointField(cellCentres, cellMap),
1108  labelUIndList(cellLevel, cellMap)(),
1109 
1110  gapShell,
1111  shellGapInfo,
1112  shellGapMode
1113  );
1114 
1115  // Compact out hits
1116 
1117  labelList map(shellGapInfo.size());
1118  compactI = 0;
1119  forAll(shellGapInfo, i)
1120  {
1121  if (shellGapInfo[i][2] > 0)
1122  {
1123  map[compactI++] = i;
1124  }
1125  }
1126 
1127  Info<< "Selected " << returnReduce(compactI, sumOp<label>())
1128  << " cells inside gap shells out of "
1129  << mesh_.globalData().nTotalCells() << endl;
1130 
1131  map.setSize(compactI);
1132  cellMap = labelUIndList(cellMap, map)();
1133  gapShell = labelUIndList(gapShell, map)();
1134  shellGapInfo = UIndirectList<FixedList<label, 3>>(shellGapInfo, map)();
1135  shellGapMode = UIndirectList<volumeType>(shellGapMode, map)();
1136 }
1137 
1138 
1139 void Foam::meshRefinement::mergeGapInfo
1140 (
1141  const FixedList<label, 3>& shellGapInfo,
1142  const volumeType shellGapMode,
1143  const FixedList<label, 3>& surfGapInfo,
1144  const volumeType surfGapMode,
1145 
1146  FixedList<label, 3>& gapInfo,
1147  volumeType& gapMode
1148 ) const
1149 {
1150  if (surfGapInfo[0] == 0)
1151  {
1152  gapInfo = shellGapInfo;
1153  gapMode = shellGapMode;
1154  }
1155  else if (shellGapInfo[0] == 0)
1156  {
1157  gapInfo = surfGapInfo;
1158  gapMode = surfGapMode;
1159  }
1160  else
1161  {
1162  // Both specify a level. Does surface level win? Or does information
1163  // need to be merged?
1164 
1165  //gapInfo[0] = max(surfGapInfo[0], shellGapInfo[0]);
1166  //gapInfo[1] = min(surfGapInfo[1], shellGapInfo[1]);
1167  //gapInfo[2] = max(surfGapInfo[2], shellGapInfo[2]);
1168  gapInfo = surfGapInfo;
1169  gapMode = surfGapMode;
1170  }
1171 }
1172 
1173 
1174 Foam::label Foam::meshRefinement::markInternalGapRefinement
1175 (
1176  const scalar planarCos,
1177  const bool spreadGapSize,
1178  const label nAllowRefine,
1179 
1180  labelList& refineCell,
1181  label& nRefine,
1182  labelList& numGapCells,
1183  scalarField& detectedGapSize
1184 ) const
1185 {
1186  detectedGapSize.setSize(mesh_.nCells());
1187  detectedGapSize = GREAT;
1188  numGapCells.setSize(mesh_.nCells());
1189  numGapCells = -1;
1190 
1191  const labelList& cellLevel = meshCutter_.cellLevel();
1192  const pointField& cellCentres = mesh_.cellCentres();
1193  const scalar edge0Len = meshCutter_.level0EdgeLength();
1194 
1195  const List<FixedList<label, 3>>& extendedGapLevel =
1196  surfaces_.extendedGapLevel();
1197  const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
1198  const boolList& extendedGapSelf = surfaces_.gapSelf();
1199 
1200  // Get the gap level for the shells
1201  const labelList maxLevel(shells_.maxGapLevel());
1202 
1203  label oldNRefine = nRefine;
1204 
1205  if (max(maxLevel) > 0)
1206  {
1207  // Collect cells to test
1208  labelList cellMap;
1209  labelList gapShell;
1210  List<FixedList<label, 3>> shellGapInfo;
1211  List<volumeType> shellGapMode;
1212  selectGapCandidates
1213  (
1214  refineCell,
1215  nRefine,
1216 
1217  cellMap,
1218  gapShell,
1219  shellGapInfo,
1220  shellGapMode
1221  );
1222 
1223  // Find nearest point and normal on the surfaces
1224  List<pointIndexHit> nearInfo;
1225  vectorField nearNormal;
1226  labelList nearSurface;
1227  labelList nearRegion;
1228  {
1229  // Now we have both the cell-level and the gap size information. Use
1230  // this to calculate the gap size
1231  scalarField gapSize(cellMap.size());
1232  forAll(cellMap, i)
1233  {
1234  label cellI = cellMap[i];
1235  scalar cellSize = edge0Len/pow(2.0, cellLevel[cellI]);
1236  gapSize[i] = shellGapInfo[i][0]*cellSize;
1237  }
1238 
1239  surfaces_.findNearestRegion
1240  (
1241  identity(surfaces_.surfaces().size()),
1242  pointField(cellCentres, cellMap),
1243  sqr(gapSize),
1244  nearSurface,
1245  nearInfo,
1246  nearRegion,
1247  nearNormal
1248  );
1249  }
1250 
1251 
1252 
1253  DynamicList<label> map(nearInfo.size());
1254  DynamicField<point> rayStart(nearInfo.size());
1255  DynamicField<point> rayEnd(nearInfo.size());
1256  DynamicField<scalar> gapSize(nearInfo.size());
1257 
1258  DynamicField<point> rayStart2(nearInfo.size());
1259  DynamicField<point> rayEnd2(nearInfo.size());
1260  DynamicField<scalar> gapSize2(nearInfo.size());
1261 
1262  label nTestCells = 0;
1263 
1264  forAll(nearInfo, i)
1265  {
1266  if (nearInfo[i].hit())
1267  {
1268  label globalRegionI = surfaces_.globalRegion
1269  (
1270  nearSurface[i],
1271  nearRegion[i]
1272  );
1273 
1274  // Combine info from shell and surface
1275  FixedList<label, 3> gapInfo;
1276  volumeType gapMode;
1277  mergeGapInfo
1278  (
1279  shellGapInfo[i],
1280  shellGapMode[i],
1281 
1282  extendedGapLevel[globalRegionI],
1283  extendedGapMode[globalRegionI],
1284 
1285  gapInfo,
1286  gapMode
1287  );
1288 
1289  // Store wanted number of cells in gap
1290  label cellI = cellMap[i];
1291  label cLevel = cellLevel[cellI];
1292  if (cLevel >= gapInfo[1] && cLevel < gapInfo[2])
1293  {
1294  numGapCells[cellI] = max(numGapCells[cellI], gapInfo[0]);
1295  }
1296 
1297  // Construct one or more rays to test for oppositeness
1298  label nRays = generateRays
1299  (
1300  false,
1301  nearInfo[i].hitPoint(),
1302  nearNormal[i],
1303  gapInfo,
1304  gapMode,
1305 
1306  cellCentres[cellI],
1307  cLevel,
1308 
1309  rayStart,
1310  rayEnd,
1311  gapSize,
1312 
1313  rayStart2,
1314  rayEnd2,
1315  gapSize2
1316  );
1317  if (nRays > 0)
1318  {
1319  nTestCells++;
1320  for (label j = 0; j < nRays; j++)
1321  {
1322  map.append(i);
1323  }
1324  }
1325  }
1326  }
1327 
1328  Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
1329  << " cells for testing out of "
1330  << mesh_.globalData().nTotalCells() << endl;
1331  map.shrink();
1332  rayStart.shrink();
1333  rayEnd.shrink();
1334  gapSize.shrink();
1335 
1336  rayStart2.shrink();
1337  rayEnd2.shrink();
1338  gapSize2.shrink();
1339 
1340  cellMap = labelUIndList(cellMap, map)();
1341  nearNormal = UIndirectList<vector>(nearNormal, map)();
1342  shellGapInfo.clear();
1343  shellGapMode.clear();
1344  nearInfo.clear();
1345  nearSurface.clear();
1346  nearRegion.clear();
1347 
1348 
1349  // Do intersections in pairs
1350  labelList surf1;
1351  List<pointIndexHit> hit1;
1352  vectorField normal1;
1353  surfaces_.findNearestIntersection
1354  (
1355  rayStart,
1356  rayEnd,
1357  surf1,
1358  hit1,
1359  normal1
1360  );
1361 
1362  labelList surf2;
1363  List<pointIndexHit> hit2;
1364  vectorField normal2;
1365  surfaces_.findNearestIntersection
1366  (
1367  rayStart2,
1368  rayEnd2,
1369  surf2,
1370  hit2,
1371  normal2
1372  );
1373 
1374  // Extract cell based gap size
1375  forAll(surf1, i)
1376  {
1377  // Combine selfProx of shell and surfaces. Ignore regions for
1378  // now
1379  const label shelli = gapShell[map[i]];
1380 
1381  bool selfProx = true;
1382  if (shelli != -1)
1383  {
1384  selfProx = shells_.gapSelf()[shelli][0];
1385  }
1386  if (surf1[i] != -1 && selfProx)
1387  {
1388  const label globalRegioni = surfaces_.globalRegion(surf1[i], 0);
1389  selfProx = extendedGapSelf[globalRegioni];
1390  }
1391 
1392  if
1393  (
1394  surf1[i] != -1
1395  && surf2[i] != -1
1396  && (surf2[i] != surf1[i] || selfProx)
1397  )
1398  {
1399  // Found intersections with surface. Check for
1400  // - small gap
1401  // - coplanar normals
1402 
1403  const label cellI = cellMap[i];
1404 
1405  const scalar d2 =
1406  hit1[i].hitPoint().distSqr(hit2[i].hitPoint());
1407 
1408  if
1409  (
1410  cellI != -1
1411  && (mag(normal1[i]&normal2[i]) > planarCos)
1412  && (d2 < Foam::sqr(gapSize[i]))
1413  )
1414  {
1415  detectedGapSize[cellI] = min
1416  (
1417  detectedGapSize[cellI],
1418  Foam::sqrt(d2)
1419  );
1420  }
1421  }
1422  }
1423 
1424  // Spread it
1425  if (spreadGapSize)
1426  {
1427  scalarField boundaryGapSize;
1428  syncTools::swapBoundaryCellList
1429  (
1430  mesh_,
1431  detectedGapSize,
1432  boundaryGapSize
1433  );
1434 
1435  // Field on cells and faces
1436  List<transportData> cellData(mesh_.nCells());
1437  List<transportData> faceData(mesh_.nFaces());
1438 
1439  // Start of walk
1440  const pointField& faceCentres = mesh_.faceCentres();
1441 
1442  DynamicList<label> frontFaces(mesh_.nFaces());
1443  DynamicList<transportData> frontData(mesh_.nFaces());
1444  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1445  {
1446  label own = mesh_.faceOwner()[faceI];
1447  label nei = mesh_.faceNeighbour()[faceI];
1448 
1449  scalar minSize = min
1450  (
1451  detectedGapSize[own],
1452  detectedGapSize[nei]
1453  );
1454 
1455  if (minSize < GREAT)
1456  {
1457  frontFaces.append(faceI);
1458  frontData.append
1459  (
1460  transportData
1461  (
1462  faceCentres[faceI],
1463  minSize,
1464  0.0
1465  )
1466  );
1467  }
1468  }
1469  for
1470  (
1471  label faceI = mesh_.nInternalFaces();
1472  faceI < mesh_.nFaces();
1473  faceI++
1474  )
1475  {
1476  label own = mesh_.faceOwner()[faceI];
1477 
1478  scalar minSize = min
1479  (
1480  detectedGapSize[own],
1481  boundaryGapSize[faceI-mesh_.nInternalFaces()]
1482  );
1483 
1484  if (minSize < GREAT)
1485  {
1486  frontFaces.append(faceI);
1487  frontData.append
1488  (
1489  transportData
1490  (
1491  faceCentres[faceI],
1492  minSize,
1493  0.0
1494  )
1495  );
1496  }
1497  }
1498 
1499  Info<< "Selected "
1500  << returnReduce(frontFaces.size(), sumOp<label>())
1501  << " faces for spreading gap size out of "
1502  << mesh_.globalData().nTotalFaces() << endl;
1503 
1504 
1505  transportData::trackData td(surfaceIndex());
1506 
1507  FaceCellWave<transportData, transportData::trackData> deltaCalc
1508  (
1509  mesh_,
1510  frontFaces,
1511  frontData,
1512  faceData,
1513  cellData,
1514  mesh_.globalData().nTotalCells()+1,
1515  td
1516  );
1517 
1518 
1519  forAll(cellMap, i)
1520  {
1521  label cellI = cellMap[i];
1522  if
1523  (
1524  cellI != -1
1525  && cellData[cellI].valid(deltaCalc.data())
1526  && numGapCells[cellI] != -1
1527  )
1528  {
1529  // Update transported gap size
1530  detectedGapSize[cellI] = min
1531  (
1532  detectedGapSize[cellI],
1533  cellData[cellI].data()
1534  );
1535  }
1536  }
1537  }
1538 
1539 
1540  // Use it
1541  forAll(cellMap, i)
1542  {
1543  label cellI = cellMap[i];
1544 
1545  if (cellI != -1 && numGapCells[cellI] != -1)
1546  {
1547  // Needed gap size
1548  label cLevel = cellLevel[cellI];
1549  scalar cellSize =
1550  meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
1551  scalar neededGapSize = numGapCells[cellI]*cellSize;
1552 
1553  if (neededGapSize > detectedGapSize[cellI])
1554  {
1555  if
1556  (
1557  !markForRefine
1558  (
1559  123,
1560  nAllowRefine,
1561  refineCell[cellI],
1562  nRefine
1563  )
1564  )
1565  {
1566  break;
1567  }
1568  }
1569  }
1570  }
1571 
1572 
1573  if
1574  (
1575  returnReduce(nRefine, sumOp<label>())
1576  > returnReduce(nAllowRefine, sumOp<label>())
1577  )
1578  {
1579  Info<< "Reached refinement limit." << endl;
1580  }
1581  }
1582 
1583  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1584 }
1585 
1586 
1587 Foam::label Foam::meshRefinement::markSmallFeatureRefinement
1588 (
1589  const scalar planarCos,
1590  const label nAllowRefine,
1591  const labelList& neiLevel,
1592  const pointField& neiCc,
1593 
1594  labelList& refineCell,
1595  label& nRefine
1596 ) const
1597 {
1598  const labelList& cellLevel = meshCutter_.cellLevel();
1599  const labelList& surfaceIndices = surfaces_.surfaces();
1600  const List<FixedList<label, 3>>& extendedGapLevel =
1601  surfaces_.extendedGapLevel();
1602  const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
1603  const boolList& extendedGapSelf = surfaces_.gapSelf();
1604 
1605  label oldNRefine = nRefine;
1606 
1607  // Check that we're using any gap refinement
1608  labelList shellMaxLevel(shells_.maxGapLevel());
1609 
1610  if (max(shellMaxLevel) == 0)
1611  {
1612  return 0;
1613  }
1614 
1615  //- Force calculation of tetBasePt
1616  (void)mesh_.tetBasePtIs();
1617  (void)mesh_.cellTree();
1618 
1619 
1620  forAll(surfaceIndices, surfI)
1621  {
1622  label geomI = surfaceIndices[surfI];
1623  const searchableSurface& geom = surfaces_.geometry()[geomI];
1624 
1625 
1626  // Get the element index in a roundabout way. Problem is e.g.
1627  // distributed surface where local indices differ from global
1628  // ones (needed for getRegion call)
1629 
1630  pointField ctrs;
1631  labelList region;
1632  vectorField normal;
1633  {
1634  // Representative local coordinates and bounding sphere
1635  scalarField radiusSqr;
1636  geom.boundingSpheres(ctrs, radiusSqr);
1637 
1638  List<pointIndexHit> info;
1639  geom.findNearest(ctrs, radiusSqr, info);
1640 
1641  forAll(info, i)
1642  {
1643  if (!info[i].hit())
1644  {
1645  FatalErrorInFunction
1646  << "fc:" << ctrs[i]
1647  << " radius:" << radiusSqr[i]
1648  << exit(FatalError);
1649  }
1650  }
1651 
1652  geom.getRegion(info, region);
1653  geom.getNormal(info, normal);
1654  }
1655 
1656  // Do test to see whether triangles are inside/outside shell with
1657  // applicable specification (minLevel <= celllevel < maxLevel)
1658  List<FixedList<label, 3>> shellGapInfo;
1659  List<volumeType> shellGapMode;
1660  labelList gapShell;
1661  shells_.findHigherGapLevel
1662  (
1663  ctrs,
1664  labelList(ctrs.size(), Zero),
1665 
1666  gapShell,
1667  shellGapInfo,
1668  shellGapMode
1669  );
1670 
1671 
1672  DynamicList<label> map(ctrs.size());
1673  DynamicList<label> cellMap(ctrs.size());
1674 
1675  DynamicField<point> rayStart(ctrs.size());
1676  DynamicField<point> rayEnd(ctrs.size());
1677  DynamicField<scalar> gapSize(ctrs.size());
1678 
1679  label nTestCells = 0;
1680 
1681  forAll(ctrs, i)
1682  {
1683  if (shellGapInfo[i][2] > 0)
1684  {
1685  label globalRegionI = surfaces_.globalRegion(surfI, region[i]);
1686 
1687  // Combine info from shell and surface
1688  FixedList<label, 3> gapInfo;
1689  volumeType gapMode;
1690  mergeGapInfo
1691  (
1692  shellGapInfo[i],
1693  shellGapMode[i],
1694 
1695  extendedGapLevel[globalRegionI],
1696  extendedGapMode[globalRegionI],
1697 
1698  gapInfo,
1699  gapMode
1700  );
1701 
1702  //- Option 1: use octree nearest searching inside polyMesh
1703  //label cellI = mesh_.findCell(pt, polyMesh::CELL_TETS);
1704 
1705  //- Option 2: use octree 'inside' searching inside polyMesh. Is
1706  // much faster.
1707  label cellI = -1;
1708  const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
1709  if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
1710  {
1711  cellI = tree.findInside(ctrs[i]);
1712  }
1713 
1714  if (cellI != -1 && refineCell[cellI] == -1)
1715  {
1716  // Construct one or two rays to test for oppositeness
1717  // Note that we always want to use the surface normal
1718  // and not the vector from cell centre to surface point
1719 
1720  label nRays = generateRays
1721  (
1722  ctrs[i],
1723  normal[i],
1724  gapInfo,
1725  gapMode,
1726 
1727  cellLevel[cellI],
1728 
1729  rayStart,
1730  rayEnd
1731  );
1732 
1733  if (nRays > 0)
1734  {
1735  nTestCells++;
1736  for (label j = 0; j < nRays; j++)
1737  {
1738  cellMap.append(cellI);
1739  map.append(i);
1740  }
1741  }
1742  }
1743  }
1744  }
1745 
1746  Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
1747  << " cells containing triangle centres out of "
1748  << mesh_.globalData().nTotalCells() << endl;
1749  map.shrink();
1750  cellMap.shrink();
1751  rayStart.shrink();
1752  rayEnd.shrink();
1753 
1754  ctrs.clear();
1755  region.clear();
1756  shellGapInfo.clear();
1757  shellGapMode.clear();
1758  normal = UIndirectList<vector>(normal, map)();
1759 
1760  // Do intersections.
1761  labelList surfaceHit;
1762  vectorField surfaceNormal;
1763  surfaces_.findNearestIntersection
1764  (
1765  rayStart,
1766  rayEnd,
1767  surfaceHit,
1768  surfaceNormal
1769  );
1770 
1771 
1772  label nOldRefine = 0;
1773 
1774 
1775  forAll(surfaceHit, i)
1776  {
1777  // Combine selfProx of shell and surfaces. Ignore regions for
1778  // now
1779  const label shelli = gapShell[map[i]];
1780  bool selfProx = true;
1781  if (shelli != -1)
1782  {
1783  selfProx = shells_.gapSelf()[shelli][0];
1784  }
1785  if (surfI != -1 && selfProx)
1786  {
1787  const label globalRegioni = surfaces_.globalRegion(surfI, 0);
1788  selfProx = extendedGapSelf[globalRegioni];
1789  }
1790 
1791  if
1792  (
1793  surfaceHit[i] != -1
1794  && (surfaceHit[i] != surfI || selfProx)
1795  )
1796  {
1797  // Found intersection with surface. Check coplanar normals.
1798  label cellI = cellMap[i];
1799 
1800  if (mag(normal[i]&surfaceNormal[i]) > planarCos)
1801  {
1802  if
1803  (
1804  !markForRefine
1805  (
1806  surfaceHit[i],
1807  nAllowRefine,
1808  refineCell[cellI],
1809  nRefine
1810  )
1811  )
1812  {
1813  break;
1814  }
1815  }
1816  }
1817  }
1818 
1819  Info<< "For surface " << geom.name() << " found "
1820  << returnReduce(nRefine-nOldRefine, sumOp<label>())
1821  << " cells in small gaps" << endl;
1822 
1823  if
1824  (
1825  returnReduce(nRefine, sumOp<label>())
1826  > returnReduce(nAllowRefine, sumOp<label>())
1827  )
1828  {
1829  Info<< "Reached refinement limit." << endl;
1830  }
1831  }
1832 
1833  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1834 }
1835 
1836 
1838 Foam::label Foam::meshRefinement::markSurfaceFieldRefinement
1839 (
1840  const label nAllowRefine,
1841  const labelList& neiLevel,
1842  const pointField& neiCc,
1843 
1844  labelList& refineCell,
1845  label& nRefine
1846 ) const
1847 {
1848  const labelList& cellLevel = meshCutter_.cellLevel();
1849  const labelList& surfaceIndices = surfaces_.surfaces();
1850 
1851  label oldNRefine = nRefine;
1852 
1853  //- Force calculation of tetBasePt
1854  (void)mesh_.tetBasePtIs();
1855  (void)mesh_.cellTree();
1856  const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
1857 
1858 
1859  forAll(surfaceIndices, surfI)
1860  {
1861  label geomI = surfaceIndices[surfI];
1862  const searchableSurface& geom = surfaces_.geometry()[geomI];
1863 
1864  // Get the element index in a roundabout way. Problem is e.g.
1865  // distributed surface where local indices differ from global
1866  // ones (needed for getRegion call)
1867 
1868  pointField ctrs;
1869  labelList region;
1870  labelList minLevelField;
1871  {
1872  // Representative local coordinates and bounding sphere
1873  scalarField radiusSqr;
1874  geom.boundingSpheres(ctrs, radiusSqr);
1875 
1876  List<pointIndexHit> info;
1877  geom.findNearest(ctrs, radiusSqr, info);
1878 
1879  forAll(info, i)
1880  {
1881  if (!info[i].hit())
1882  {
1883  FatalErrorInFunction
1884  << "fc:" << ctrs[i]
1885  << " radius:" << radiusSqr[i]
1886  << exit(FatalError);
1887  }
1888  }
1889 
1890  geom.getRegion(info, region);
1891  geom.getField(info, minLevelField);
1892  }
1893 
1894  if (minLevelField.size() != geom.size())
1895  {
1896  Pout<< "** no minLevelField" << endl;
1897  continue;
1898  }
1899 
1900 
1901  label nOldRefine = 0;
1902 
1903  forAll(ctrs, i)
1904  {
1905  label cellI = -1;
1906  if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
1907  {
1908  cellI = tree.findInside(ctrs[i]);
1909  }
1910 
1911  if
1912  (
1913  cellI != -1
1914  && refineCell[cellI] == -1
1915  && minLevelField[i] > cellLevel[cellI]
1916  )
1917  {
1918  if
1919  (
1920  !markForRefine
1921  (
1922  surfI,
1923  nAllowRefine,
1924  refineCell[cellI],
1925  nRefine
1926  )
1927  )
1928  {
1929  break;
1930  }
1931  }
1932  }
1933 
1934  Info<< "For surface " << geom.name() << " found "
1935  << returnReduce(nRefine-nOldRefine, sumOp<label>())
1936  << " cells containing cached refinement field" << endl;
1937 
1938  if
1939  (
1940  returnReduce(nRefine, sumOp<label>())
1941  > returnReduce(nAllowRefine, sumOp<label>())
1942  )
1943  {
1944  Info<< "Reached refinement limit." << endl;
1945  }
1946  }
1947 
1948  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1949 }
1951 
1952 
1953 // ************************************************************************* //
Foam::UList::size
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:116
Foam::refinementSurfaces::surfaces
const labelList & surfaces() const
Definition: refinementSurfaces.H:238
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::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::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1122
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::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:44
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: IndirectList.H:65
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:137
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:531
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::primitiveMesh::isInternalFace
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
Definition: primitiveMeshI.H:96
oldNRefine
const label oldNRefine
Definition: meshRefinementBlock.C:2673
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:421
Foam::refinementSurfaces::gapSelf
const boolList & gapSelf() const
From global region number to whether to detect gaps to same.
Definition: refinementSurfaces.H:333
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::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:320
Foam::identity
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
Definition: labelLists.C:44
Foam::shellSurfaces::gapSelf
const boolListList & gapSelf() const
Per shell, per region whether to test for gap with same surface.
Definition: shellSurfaces.H:249
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:46
tree
Tree tree(triangles.begin(), triangles.end())
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1116
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::volumeType::INSIDE
A location inside the volume.
Definition: volumeType.H:65
Foam::refinementSurfaces::extendedGapMode
const List< volumeType > & extendedGapMode() const
From global region number to side of surface to detect.
Definition: refinementSurfaces.H:323
Foam::volumeType::OUTSIDE
A location outside the volume.
Definition: volumeType.H:66
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::primitiveMesh::cellCentres
const vectorField & cellCentres() const
Definition: primitiveMeshCellCentresAndVols.C:78
Foam::shellSurfaces::maxGapLevel
labelList maxGapLevel() const
Highest shell gap level.
Definition: shellSurfaces.C:830
Foam::volumeType::MIXED
A location that is partly inside and outside.
Definition: volumeType.H:67
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::refinementSurfaces::extendedGapLevel
const List< FixedList< label, 3 > > & extendedGapLevel() const
From global region number to specification of gap and its.
Definition: refinementSurfaces.H:315
Foam::point
vector point
Point is a vector.
Definition: point.H:37
Foam::primitiveMesh::nCells
label nCells() const noexcept
Number of mesh cells.
Definition: primitiveMeshI.H:89
Foam::refinementSurfaces::globalRegion
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
Definition: refinementSurfaces.H:376
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:48
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:62
Foam::refinementSurfaces::findNearestIntersection
void findNearestIntersection(const labelList &surfacesToTest, const pointField &start, const pointField &end, labelList &surface1, List< pointIndexHit > &hit1, labelList &region1, labelList &surface2, List< pointIndexHit > &hit2, labelList &region2) const
Find intersection nearest to the endpoints. surface1,2 are.
Definition: refinementSurfaces.C:1658
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::mode
mode_t mode(const fileName &name, const bool followLink=true)
Return the file mode, normally following symbolic links.
Definition: POSIX.C:773
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:60
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::hexRef8::cellLevel
const labelIOList & cellLevel() const
Definition: hexRef8.H:481
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127