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-2020 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  magSqr(hit1[i].hitPoint()-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.rawPoint()-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.rawPoint();
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 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());
1406 
1407  if
1408  (
1409  cellI != -1
1410  && (mag(normal1[i]&normal2[i]) > planarCos)
1411  && (d2 < Foam::sqr(gapSize[i]))
1412  )
1413  {
1414  detectedGapSize[cellI] = min
1415  (
1416  detectedGapSize[cellI],
1417  Foam::sqrt(d2)
1418  );
1419  }
1420  }
1421  }
1422 
1423  // Spread it
1424  if (spreadGapSize)
1425  {
1426  // Field on cells and faces
1427  List<transportData> cellData(mesh_.nCells());
1428  List<transportData> faceData(mesh_.nFaces());
1429 
1430  // Start of walk
1431  const pointField& faceCentres = mesh_.faceCentres();
1432 
1433  DynamicList<label> frontFaces(mesh_.nFaces());
1434  DynamicList<transportData> frontData(mesh_.nFaces());
1435  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1436  {
1437  label own = mesh_.faceOwner()[faceI];
1438  label nei = mesh_.faceNeighbour()[faceI];
1439 
1440  scalar minSize = min
1441  (
1442  detectedGapSize[own],
1443  detectedGapSize[nei]
1444  );
1445 
1446  if (minSize < GREAT)
1447  {
1448  frontFaces.append(faceI);
1449  frontData.append
1450  (
1451  transportData
1452  (
1453  faceCentres[faceI],
1454  minSize,
1455  0.0
1456  )
1457  );
1458  }
1459  }
1460  for
1461  (
1462  label faceI = mesh_.nInternalFaces();
1463  faceI < mesh_.nFaces();
1464  faceI++
1465  )
1466  {
1467  label own = mesh_.faceOwner()[faceI];
1468 
1469  if (detectedGapSize[own] < GREAT)
1470  {
1471  frontFaces.append(faceI);
1472  frontData.append
1473  (
1474  transportData
1475  (
1476  faceCentres[faceI],
1477  detectedGapSize[own],
1478  0.0
1479  )
1480  );
1481  }
1482  }
1483 
1484  Info<< "Selected "
1485  << returnReduce(frontFaces.size(), sumOp<label>())
1486  << " faces for spreading gap size out of "
1487  << mesh_.globalData().nTotalFaces() << endl;
1488 
1489 
1490  transportData::trackData td(surfaceIndex());
1491 
1492  FaceCellWave<transportData, transportData::trackData> deltaCalc
1493  (
1494  mesh_,
1495  frontFaces,
1496  frontData,
1497  faceData,
1498  cellData,
1499  mesh_.globalData().nTotalCells()+1,
1500  td
1501  );
1502 
1503 
1504  forAll(cellMap, i)
1505  {
1506  label cellI = cellMap[i];
1507  if
1508  (
1509  cellI != -1
1510  && cellData[cellI].valid(deltaCalc.data())
1511  && numGapCells[cellI] != -1
1512  )
1513  {
1514  // Update transported gap size
1515  detectedGapSize[cellI] = min
1516  (
1517  detectedGapSize[cellI],
1518  cellData[cellI].data()
1519  );
1520  }
1521  }
1522  }
1523 
1524 
1525  // Use it
1526  forAll(cellMap, i)
1527  {
1528  label cellI = cellMap[i];
1529 
1530  if (cellI != -1 && numGapCells[cellI] != -1)
1531  {
1532  // Needed gap size
1533  label cLevel = cellLevel[cellI];
1534  scalar cellSize =
1535  meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
1536  scalar neededGapSize = numGapCells[cellI]*cellSize;
1537 
1538  if (neededGapSize > detectedGapSize[cellI])
1539  {
1540  if
1541  (
1542  !markForRefine
1543  (
1544  123,
1545  nAllowRefine,
1546  refineCell[cellI],
1547  nRefine
1548  )
1549  )
1550  {
1551  break;
1552  }
1553  }
1554  }
1555  }
1556 
1557 
1558  if
1559  (
1560  returnReduce(nRefine, sumOp<label>())
1561  > returnReduce(nAllowRefine, sumOp<label>())
1562  )
1563  {
1564  Info<< "Reached refinement limit." << endl;
1565  }
1566  }
1567 
1568  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1569 }
1570 
1571 
1572 Foam::label Foam::meshRefinement::markSmallFeatureRefinement
1573 (
1574  const scalar planarCos,
1575  const label nAllowRefine,
1576  const labelList& neiLevel,
1577  const pointField& neiCc,
1578 
1579  labelList& refineCell,
1580  label& nRefine
1581 ) const
1582 {
1583  const labelList& cellLevel = meshCutter_.cellLevel();
1584  const labelList& surfaceIndices = surfaces_.surfaces();
1585  const List<FixedList<label, 3>>& extendedGapLevel =
1586  surfaces_.extendedGapLevel();
1587  const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
1588  const boolList& extendedGapSelf = surfaces_.gapSelf();
1589 
1590  label oldNRefine = nRefine;
1591 
1592  // Check that we're using any gap refinement
1593  labelList shellMaxLevel(shells_.maxGapLevel());
1594 
1595  if (max(shellMaxLevel) == 0)
1596  {
1597  return 0;
1598  }
1599 
1600  //- Force calculation of tetBasePt
1601  (void)mesh_.tetBasePtIs();
1602  (void)mesh_.cellTree();
1603 
1604 
1605  forAll(surfaceIndices, surfI)
1606  {
1607  label geomI = surfaceIndices[surfI];
1608  const searchableSurface& geom = surfaces_.geometry()[geomI];
1609 
1610 
1611  // Get the element index in a roundabout way. Problem is e.g.
1612  // distributed surface where local indices differ from global
1613  // ones (needed for getRegion call)
1614 
1615  pointField ctrs;
1616  labelList region;
1617  vectorField normal;
1618  {
1619  // Representative local coordinates and bounding sphere
1620  scalarField radiusSqr;
1621  geom.boundingSpheres(ctrs, radiusSqr);
1622 
1623  List<pointIndexHit> info;
1624  geom.findNearest(ctrs, radiusSqr, info);
1625 
1626  forAll(info, i)
1627  {
1628  if (!info[i].hit())
1629  {
1630  FatalErrorInFunction
1631  << "fc:" << ctrs[i]
1632  << " radius:" << radiusSqr[i]
1633  << exit(FatalError);
1634  }
1635  }
1636 
1637  geom.getRegion(info, region);
1638  geom.getNormal(info, normal);
1639  }
1640 
1641  // Do test to see whether triangles are inside/outside shell with
1642  // applicable specification (minLevel <= celllevel < maxLevel)
1643  List<FixedList<label, 3>> shellGapInfo;
1644  List<volumeType> shellGapMode;
1645  labelList gapShell;
1646  shells_.findHigherGapLevel
1647  (
1648  ctrs,
1649  labelList(ctrs.size(), Zero),
1650 
1651  gapShell,
1652  shellGapInfo,
1653  shellGapMode
1654  );
1655 
1656 
1657  DynamicList<label> map(ctrs.size());
1658  DynamicList<label> cellMap(ctrs.size());
1659 
1660  DynamicField<point> rayStart(ctrs.size());
1661  DynamicField<point> rayEnd(ctrs.size());
1662  DynamicField<scalar> gapSize(ctrs.size());
1663 
1664  label nTestCells = 0;
1665 
1666  forAll(ctrs, i)
1667  {
1668  if (shellGapInfo[i][2] > 0)
1669  {
1670  label globalRegionI = surfaces_.globalRegion(surfI, region[i]);
1671 
1672  // Combine info from shell and surface
1673  FixedList<label, 3> gapInfo;
1674  volumeType gapMode;
1675  mergeGapInfo
1676  (
1677  shellGapInfo[i],
1678  shellGapMode[i],
1679 
1680  extendedGapLevel[globalRegionI],
1681  extendedGapMode[globalRegionI],
1682 
1683  gapInfo,
1684  gapMode
1685  );
1686 
1687  //- Option 1: use octree nearest searching inside polyMesh
1688  //label cellI = mesh_.findCell(pt, polyMesh::CELL_TETS);
1689 
1690  //- Option 2: use octree 'inside' searching inside polyMesh. Is
1691  // much faster.
1692  label cellI = -1;
1693  const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
1694  if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
1695  {
1696  cellI = tree.findInside(ctrs[i]);
1697  }
1698 
1699  if (cellI != -1 && refineCell[cellI] == -1)
1700  {
1701  // Construct one or two rays to test for oppositeness
1702  // Note that we always want to use the surface normal
1703  // and not the vector from cell centre to surface point
1704 
1705  label nRays = generateRays
1706  (
1707  ctrs[i],
1708  normal[i],
1709  gapInfo,
1710  gapMode,
1711 
1712  cellLevel[cellI],
1713 
1714  rayStart,
1715  rayEnd
1716  );
1717 
1718  if (nRays > 0)
1719  {
1720  nTestCells++;
1721  for (label j = 0; j < nRays; j++)
1722  {
1723  cellMap.append(cellI);
1724  map.append(i);
1725  }
1726  }
1727  }
1728  }
1729  }
1730 
1731  Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
1732  << " cells containing triangle centres out of "
1733  << mesh_.globalData().nTotalCells() << endl;
1734  map.shrink();
1735  cellMap.shrink();
1736  rayStart.shrink();
1737  rayEnd.shrink();
1738 
1739  ctrs.clear();
1740  region.clear();
1741  shellGapInfo.clear();
1742  shellGapMode.clear();
1743  normal = UIndirectList<vector>(normal, map)();
1744 
1745  // Do intersections.
1746  labelList surfaceHit;
1747  vectorField surfaceNormal;
1748  surfaces_.findNearestIntersection
1749  (
1750  rayStart,
1751  rayEnd,
1752  surfaceHit,
1753  surfaceNormal
1754  );
1755 
1756 
1757  label nOldRefine = 0;
1758 
1759 
1760  forAll(surfaceHit, i)
1761  {
1762  // Combine selfProx of shell and surfaces. Ignore regions for
1763  // now
1764  const label shelli = gapShell[map[i]];
1765  bool selfProx = true;
1766  if (shelli != -1)
1767  {
1768  selfProx = shells_.gapSelf()[shelli][0];
1769  }
1770  if (surfI != -1 && selfProx)
1771  {
1772  const label globalRegioni = surfaces_.globalRegion(surfI, 0);
1773  selfProx = extendedGapSelf[globalRegioni];
1774  }
1775 
1776  if
1777  (
1778  surfaceHit[i] != -1
1779  && (surfaceHit[i] != surfI || selfProx)
1780  )
1781  {
1782  // Found intersection with surface. Check coplanar normals.
1783  label cellI = cellMap[i];
1784 
1785  if (mag(normal[i]&surfaceNormal[i]) > planarCos)
1786  {
1787  if
1788  (
1789  !markForRefine
1790  (
1791  surfaceHit[i],
1792  nAllowRefine,
1793  refineCell[cellI],
1794  nRefine
1795  )
1796  )
1797  {
1798  break;
1799  }
1800  }
1801  }
1802  }
1803 
1804  Info<< "For surface " << geom.name() << " found "
1805  << returnReduce(nRefine-nOldRefine, sumOp<label>())
1806  << " cells in small gaps" << endl;
1807 
1808  if
1809  (
1810  returnReduce(nRefine, sumOp<label>())
1811  > returnReduce(nAllowRefine, sumOp<label>())
1812  )
1813  {
1814  Info<< "Reached refinement limit." << endl;
1815  }
1816  }
1817 
1818  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1819 }
1820 
1821 
1822 // ************************************************************************* //
shellSurfaces.H
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
volFields.H
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition: primitiveFieldsFwd.H:52
Foam::hexRef8::cellLevel
const labelIOList & cellLevel() const
Definition: hexRef8.H:397
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
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::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::refinementSurfaces::surfaces
const labelList & surfaces() const
Definition: refinementSurfaces.H:183
Foam::refinementSurfaces::gapSelf
const boolList & gapSelf() const
From global region number to whether to detect gaps to same.
Definition: refinementSurfaces.H:249
Foam::shellSurfaces::maxGapLevel
labelList maxGapLevel() const
Highest shell gap level.
Definition: shellSurfaces.C:826
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:65
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:175
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
refinementFeatures.H
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::mode
mode_t mode(const fileName &name, const bool followLink=true)
Return the file mode, normally following symbolic links.
Definition: MSwindows.C:564
searchableSurfaces.H
n
label n
Definition: TABSMDCalcMethod2.H:31
Foam::primitiveMesh::nCells
label nCells() const noexcept
Number of mesh cells.
Definition: primitiveMeshI.H:96
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
refinementSurfaces.H
Foam::volumeType::MIXED
A location that is partly inside and outside.
Definition: volumeType.H:70
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1107
Foam::primitiveMesh::nBoundaryFaces
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
Definition: primitiveMeshI.H:84
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::refinementSurfaces::extendedGapLevel
const List< FixedList< label, 3 > > & extendedGapLevel() const
From global region number to specification of gap and its.
Definition: refinementSurfaces.H:236
Foam::FatalError
error FatalError
Foam::refinementSurfaces::globalRegion
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
Definition: refinementSurfaces.H:270
meshRefinement.H
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
treeDataCell.H
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
Time.H
Foam::primitiveMesh::cellCentres
const vectorField & cellCentres() const
Definition: primitiveMeshCellCentresAndVols.C:84
Foam::primitiveMesh::isInternalFace
bool isInternalFace(const label faceIndex) const noexcept
Return true if given face label is internal to the mesh.
Definition: primitiveMeshI.H:103
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
DynamicField.H
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:1192
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
transportData.H
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const noexcept
Number of internal faces.
Definition: primitiveMeshI.H:78
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::refinementSurfaces::extendedGapMode
const List< volumeType > & extendedGapMode() const
From global region number to side of surface to detect.
Definition: refinementSurfaces.H:242
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
FaceCellWave.H
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::identity
labelList identity(const label len, label start=0)
Create identity map of the given length with (map[i] == i)
Definition: labelList.C:38
Foam::volumeType::INSIDE
A location inside the volume.
Definition: volumeType.H:68
Foam::shellSurfaces::gapSelf
const boolListList & gapSelf() const
Per shell, per region whether to test for gap with same surface.
Definition: shellSurfaces.H:204
Foam::point
vector point
Point is a vector.
Definition: point.H:43
triSurfaceMesh.H
zeroGradientFvPatchFields.H
Foam::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1113
Foam::volumeType::OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: UIndirectList.H:58