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. Ignore regions for
457  // now
458  const label cellI = cellMap[i];
459  const label shelli =
460  (
461  cellToCompact[cellI] != -1
462  ? gapShell[cellToCompact[cellI]]
463  : -1
464  );
465 
466  bool selfProx = true;
467  if (shelli != -1)
468  {
469  selfProx = shells_.gapSelf()[shelli][0];
470  }
471  if (surf1[i] != -1 && selfProx)
472  {
473  const label globalRegioni = surfaces_.globalRegion(surf1[i], 0);
474  selfProx = extendedGapSelf[globalRegioni];
475  }
476 
477  if
478  (
479  surf1[i] != -1
480  && surf2[i] != -1
481  && (surf2[i] != surf1[i] || selfProx)
482  )
483  {
484  // Found intersection with surface. Check opposite normal.
485  label cellI = cellMap[i];
486 
487  if
488  (
489  cellI != -1
490  && (mag(normal1[i]&normal2[i]) > planarCos)
491  && (
492  magSqr(hit1[i].hitPoint()-hit2[i].hitPoint())
493  < Foam::sqr(gapSize[i])
494  )
495  )
496  {
497  if
498  (
499  !markForRefine
500  (
501  surf1[i],
502  nAllowRefine,
503  refineCell[cellI],
504  nRefine
505  )
506  )
507  {
508  break;
509  }
510  }
511  }
512  }
513 
514  if
515  (
516  returnReduce(nRefine, sumOp<label>())
517  > returnReduce(nAllowRefine, sumOp<label>())
518  )
519  {
520  Info<< "Reached refinement limit." << endl;
521  }
522  }
523 
524  return returnReduce(nRefine-oldNRefine, sumOp<label>());
525 }
526 
527 
528 //Foam::meshRefinement::findNearestOppositeOp::findNearestOppositeOp
529 //(
530 // const indexedOctree<treeDataTriSurface>& tree,
531 // const point& oppositePoint,
532 // const vector& oppositeNormal,
533 // const scalar minCos
534 //)
535 //:
536 // tree_(tree),
537 // oppositePoint_(oppositePoint),
538 // oppositeNormal_(oppositeNormal),
539 // minCos_(minCos)
540 //{}
541 //
542 //
543 //void Foam::meshRefinement::findNearestOppositeOp::operator()
544 //(
545 // const labelUList& indices,
546 // const point& sample,
547 // scalar& nearestDistSqr,
548 // label& minIndex,
549 // point& nearestPoint
550 //) const
551 //{
552 // const treeDataTriSurface& shape = tree_.shapes();
553 // const triSurface& patch = shape.patch();
554 // const pointField& points = patch.points();
555 //
556 // forAll(indices, i)
557 // {
558 // const label index = indices[i];
559 // const labelledTri& f = patch[index];
560 //
561 // pointHit nearHit = f.nearestPoint(sample, points);
562 // scalar distSqr = sqr(nearHit.distance());
563 //
564 // if (distSqr < nearestDistSqr)
565 // {
566 // // Nearer. Check if
567 // // - a bit way from other hit
568 // // - in correct search cone
569 // vector d(nearHit.rawPoint()-oppositePoint_);
570 // scalar normalDist(d&oppositeNormal_);
571 //
572 // if (normalDist > Foam::sqr(SMALL) && normalDist/mag(d) > minCos_)
573 // {
574 // nearestDistSqr = distSqr;
575 // minIndex = index;
576 // nearestPoint = nearHit.rawPoint();
577 // }
578 // }
579 // }
580 //}
581 //
582 //
583 //void Foam::meshRefinement::searchCone
584 //(
585 // const label surfI,
586 // labelList& nearMap, // cells
587 // scalarField& nearGap, // gap size
588 // List<pointIndexHit>& nearInfo, // nearest point on surface
589 // List<pointIndexHit>& oppositeInfo // detected point on gap (or miss)
590 //) const
591 //{
592 // const labelList& cellLevel = meshCutter_.cellLevel();
593 // const pointField& cellCentres = mesh_.cellCentres();
594 // const scalar edge0Len = meshCutter_.level0EdgeLength();
595 //
596 // const labelList& surfaceIndices = surfaces_.surfaces();
597 // const List<FixedList<label, 3>>& extendedGapLevel =
598 // surfaces_.extendedGapLevel();
599 // const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
600 //
601 //
602 // label geomI = surfaceIndices[surfI];
603 // const searchableSurface& geom = surfaces_.geometry()[geomI];
604 //
605 // const triSurfaceMesh& s = refCast<const triSurfaceMesh>(geom);
606 // const indexedOctree<treeDataTriSurface>& tree = s.tree();
607 //
608 //
609 // const scalar searchCos = Foam::cos(degToRad(30.0));
610 //
611 // // Normals for ray shooting and inside/outside detection
612 // vectorField nearNormal;
613 // geom.getNormal(nearInfo, nearNormal);
614 // // Regions
615 // labelList nearRegion;
616 // geom.getRegion(nearInfo, nearRegion);
617 //
618 //
619 // // Now loop over all near points and search in the half cone
620 // labelList map(nearInfo.size());
621 // label compactI = 0;
622 //
623 // oppositeInfo.setSize(nearInfo.size());
624 //
625 // forAll(nearInfo, i)
626 // {
627 // label globalRegionI =
628 // surfaces_.globalRegion(surfI, nearRegion[i]);
629 //
630 // // Get updated gap information now we have the region
631 // label nGapCells = extendedGapLevel[globalRegionI][0];
632 // label minLevel = extendedGapLevel[globalRegionI][1];
633 // label maxLevel = extendedGapLevel[globalRegionI][2];
634 // volumeType mode = extendedGapMode[globalRegionI];
635 //
636 // label cellI = nearMap[i];
637 // label cLevel = cellLevel[cellI];
638 //
639 // if (cLevel >= minLevel && cLevel < maxLevel)
640 // {
641 // scalar cellSize = edge0Len/pow(2.0, cLevel);
642 //
643 // // Update gap size
644 // nearGap[i] = nGapCells*cellSize;
645 //
646 // const point& nearPt = nearInfo[i].hitPoint();
647 // vector v(cellCentres[cellI]-nearPt);
648 // scalar magV = mag(v);
649 //
650 // // Like with ray shooting we want to
651 // // - find triangles up to nearGap away on the wanted side of the
652 // // surface
653 // // - find triangles up to 0.5*cellSize away on the unwanted side
654 // // of the surface. This is for cells straddling the surface
655 // // where
656 // // the cell centre might be on the wrong side of the surface
657 //
658 // // Tbd: check that cell centre is inbetween the gap hits
659 // // (only if the cell is far enough away)
660 //
661 // scalar posNormalSize = 0.0;
662 // scalar negNormalSize = 0.0;
663 //
664 // if (mode == volumeType::OUTSIDE)
665 // {
666 // posNormalSize = nearGap[i];
667 // if (magV < 0.5*cellSize)
668 // {
669 // negNormalSize = 0.5*cellSize;
670 // }
671 // }
672 // else if (mode == volumeType::INSIDE)
673 // {
674 // if (magV < 0.5*cellSize)
675 // {
676 // posNormalSize = 0.5*cellSize;
677 // }
678 // negNormalSize = nearGap[i];
679 // }
680 // else
681 // {
682 // posNormalSize = nearGap[i];
683 // negNormalSize = nearGap[i];
684 // }
685 //
686 // // Test with positive normal
687 // oppositeInfo[compactI] = tree.findNearest
688 // (
689 // nearPt,
690 // sqr(posNormalSize),
691 // findNearestOppositeOp
692 // (
693 // tree,
694 // nearPt,
695 // nearNormal[i],
696 // searchCos
697 // )
698 // );
699 //
700 // if (oppositeInfo[compactI].hit())
701 // {
702 // map[compactI++] = i;
703 // }
704 // else
705 // {
706 // // Test with negative normal
707 // oppositeInfo[compactI] = tree.findNearest
708 // (
709 // nearPt,
710 // sqr(negNormalSize),
711 // findNearestOppositeOp
712 // (
713 // tree,
714 // nearPt,
715 // -nearNormal[i],
716 // searchCos
717 // )
718 // );
719 //
720 // if (oppositeInfo[compactI].hit())
721 // {
722 // map[compactI++] = i;
723 // }
724 // }
725 // }
726 // }
727 //
728 // Info<< "Selected " << returnReduce(compactI, sumOp<label>())
729 // << " hits on the correct side out of "
730 // << returnReduce(map.size(), sumOp<label>()) << endl;
731 // map.setSize(compactI);
732 // oppositeInfo.setSize(compactI);
733 //
734 // nearMap = labelUIndList(nearMap, map)();
735 // nearGap = UIndirectList<scalar>(nearGap, map)();
736 // nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
737 // nearNormal = UIndirectList<vector>(nearNormal, map)();
738 //
739 // // Exclude hits which aren't opposite enough. E.g. you might find
740 // // a point on a perpendicular wall - but this does not constitute a gap.
741 // vectorField oppositeNormal;
742 // geom.getNormal(oppositeInfo, oppositeNormal);
743 //
744 // compactI = 0;
745 // forAll(oppositeInfo, i)
746 // {
747 // if ((nearNormal[i] & oppositeNormal[i]) < -0.707)
748 // {
749 // map[compactI++] = i;
750 // }
751 // }
752 //
753 // Info<< "Selected " << returnReduce(compactI, sumOp<label>())
754 // << " hits opposite the nearest out of "
755 // << returnReduce(map.size(), sumOp<label>()) << endl;
756 // map.setSize(compactI);
757 //
758 // nearMap = labelUIndList(nearMap, map)();
759 // nearGap = UIndirectList<scalar>(nearGap, map)();
760 // nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
761 // oppositeInfo = UIndirectList<pointIndexHit>(oppositeInfo, map)();
762 //}
763 
764 
765 Foam::label Foam::meshRefinement::generateRays
766 (
767  const point& nearPoint,
768  const vector& nearNormal,
769  const FixedList<label, 3>& gapInfo,
770  const volumeType& mode,
771 
772  const label cLevel,
773 
774  DynamicField<point>& start,
775  DynamicField<point>& end
776 ) const
777 {
778  label nOldRays = start.size();
779 
780  if (cLevel >= gapInfo[1] && cLevel < gapInfo[2] && gapInfo[0] > 0)
781  {
782  scalar cellSize = meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
783 
784  // Calculate gap size
785  scalar nearGap = gapInfo[0]*cellSize;
786 
787  const vector& n = nearNormal;
788 
789  // Situation 'C' above: cell too close. Use surface
790  // -normal and -point to shoot rays
791 
792  if (mode == volumeType::OUTSIDE)
793  {
794  start.append(nearPoint+1e-6*n);
795  end.append(nearPoint+nearGap*n);
796  }
797  else if (mode == volumeType::INSIDE)
798  {
799  start.append(nearPoint-1e-6*n);
800  end.append(nearPoint-nearGap*n);
801  }
802  else if (mode == volumeType::MIXED)
803  {
804  start.append(nearPoint+1e-6*n);
805  end.append(nearPoint+nearGap*n);
806 
807  start.append(nearPoint-1e-6*n);
808  end.append(nearPoint-nearGap*n);
809  }
810  }
811 
812  return start.size()-nOldRays;
813 }
814 
815 
816 Foam::label Foam::meshRefinement::generateRays
817 (
818  const bool useSurfaceNormal,
819 
820  const point& nearPoint,
821  const vector& nearNormal,
822  const FixedList<label, 3>& gapInfo,
823  const volumeType& mode,
824 
825  const point& cc,
826  const label cLevel,
827 
828  DynamicField<point>& start,
829  DynamicField<point>& end,
830  DynamicField<scalar>& gapSize,
831 
832  DynamicField<point>& start2,
833  DynamicField<point>& end2,
834  DynamicField<scalar>& gapSize2
835 ) const
836 {
837  // We want to handle the following cases:
838  // - surface: small gap (marked with 'surface'). gap might be
839  // on inside or outside of surface.
840  // - A: cell well inside the gap.
841  // - B: cell well outside the gap.
842  // - C: cell straddling the gap. cell centre might be inside
843  // or outside
844  //
845  // +---+
846  // | B |
847  // +---+
848  //
849  // +------+
850  // | |
851  // | C |
852  // --------|------|----surface
853  // +------+
854  //
855  // +---+
856  // | A |
857  // +---+
858  //
859  //
860  // --------------------surface
861  //
862  // So:
863  // - find nearest point on surface
864  // - in situation A,B decide if on wanted side of surface
865  // - detect if locally a gap (and the cell inside the gap) by
866  // shooting a ray from the point on the surface in the direction
867  // of
868  // - A,B: the cell centre
869  // - C: the surface normal and/or negative surface normal
870  // and see we hit anything
871  //
872  // Variations of this scheme:
873  // - always shoot in the direction of the surface normal. This needs
874  // then an additional check to make sure the cell centre is
875  // somewhere inside the gap
876  // - instead of ray shooting use a 'constrained' nearest search
877  // by e.g. looking inside a search cone (implemented in searchCone).
878  // The problem with this constrained nearest is that it still uses
879  // the absolute nearest point on each triangle and only afterwards
880  // checks if it is inside the search cone.
881 
882 
883  // Decide which near points are good:
884  // - with updated minLevel and maxLevel and nearGap make sure
885  // the cell is still a candidate
886  // NOTE: inside the gap the nearest point on the surface will
887  // be HALF the gap size - otherwise we would have found
888  // a point on the opposite side
889  // - if the mode is both sides
890  // - or if the hit is inside the current cell (situation 'C',
891  // magV < 0.5cellSize)
892  // - or otherwise if on the correct side
893 
894  label nOldRays = start.size();
895 
896  if (cLevel >= gapInfo[1] && cLevel < gapInfo[2] && gapInfo[0] > 0)
897  {
898  scalar cellSize = meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
899 
900  // Calculate gap size
901  scalar nearGap = gapInfo[0]*cellSize;
902 
903  // Distance to nearest
904  vector v(cc-nearPoint);
905  scalar magV = mag(v);
906 
907  if (useSurfaceNormal || magV < 0.5*cellSize)
908  {
909  const vector& n = nearNormal;
910 
911  // Situation 'C' above: cell too close. Use surface
912  // -normal and -point to shoot rays
913 
914  if (mode == volumeType::OUTSIDE)
915  {
916  start.append(nearPoint+1e-6*n);
917  end.append(nearPoint+nearGap*n);
918  gapSize.append(nearGap);
919  // Second vector so we get pairs of intersections
920  start2.append(nearPoint+1e-6*n);
921  end2.append(nearPoint-1e-6*n);
922  gapSize2.append(gapSize.last());
923  }
924  else if (mode == volumeType::INSIDE)
925  {
926  start.append(nearPoint-1e-6*n);
927  end.append(nearPoint-nearGap*n);
928  gapSize.append(nearGap);
929  // Second vector so we get pairs of intersections
930  start2.append(nearPoint-1e-6*n);
931  end2.append(nearPoint+1e-6*n);
932  gapSize2.append(gapSize.last());
933  }
934  else if (mode == volumeType::MIXED)
935  {
936  // Do both rays:
937  // Outside
938  {
939  start.append(nearPoint+1e-6*n);
940  end.append(nearPoint+nearGap*n);
941  gapSize.append(nearGap);
942  // Second vector so we get pairs of intersections
943  start2.append(nearPoint+1e-6*n);
944  end2.append(nearPoint-1e-6*n);
945  gapSize2.append(gapSize.last());
946  }
947  // Inside
948  {
949  start.append(nearPoint-1e-6*n);
950  end.append(nearPoint-nearGap*n);
951  gapSize.append(nearGap);
952  // Second vector so we get pairs of intersections
953  start2.append(nearPoint-1e-6*n);
954  end2.append(nearPoint+1e-6*n);
955  gapSize2.append(gapSize.last());
956  }
957  }
958  }
959  else
960  {
961  // Situation 'A' or 'B' above: cell well away. Test if
962  // cell on correct side of surface and shoot ray through
963  // cell centre. Note: no need to shoot ray in other
964  // direction since we're trying to detect cell inside
965  // the gap.
966 
967  scalar s = (v&nearNormal);
968 
969  if
970  (
971  (mode == volumeType::MIXED)
972  || (mode == volumeType::OUTSIDE && s > SMALL)
973  || (mode == volumeType::INSIDE && s < -SMALL)
974  )
975  {
977  //vector n(v/(magV+ROOTVSMALL));
978  //
979  //start.append(cc);
980  //end.append(cc+nearGap*n);
981  //gapSize.append(nearGap);
982  //
983  //start2.append(cc);
984  //end2.append(cc-nearGap*n);
985  //gapSize2.append(nearGap);
986 
987 
990  //start.append(cc);
991  //end.append(cc+nearGap*vector(1, 0, 0));
992  //gapSize.append(nearGap);
993  //
994  //start2.append(cc);
995  //end2.append(cc-nearGap*vector(1, 0, 0));
996  //gapSize2.append(nearGap);
997  //
999  //start.append(cc);
1000  //end.append(cc+nearGap*vector(0, 1, 0));
1001  //gapSize.append(nearGap);
1002  //
1003  //start2.append(cc);
1004  //end2.append(cc-nearGap*vector(0, 1, 0));
1005  //gapSize2.append(nearGap);
1006  //
1008  //start.append(cc);
1009  //end.append(cc+nearGap*vector(0, 0, 1));
1010  //gapSize.append(nearGap);
1011  //
1012  //start2.append(cc);
1013  //end2.append(cc-nearGap*vector(0, 0, 1));
1014  //gapSize2.append(nearGap);
1015 
1016 
1017  // 3 axes aligned with normal
1018 
1019  // Use vector through cell centre
1020  vector n(v/(magV+ROOTVSMALL));
1021 
1022  // Get second vector. Make sure it is sufficiently perpendicular
1023  vector e2(1, 0, 0);
1024  scalar s = (e2 & n);
1025  if (mag(s) < 0.9)
1026  {
1027  e2 -= s*n;
1028  }
1029  else
1030  {
1031  e2 = vector(0, 1, 0);
1032  e2 -= (e2 & n)*n;
1033  }
1034  e2 /= mag(e2);
1035 
1036  // Third vector
1037  vector e3 = n ^ e2;
1038 
1039 
1040  // Rays in first direction
1041  start.append(cc);
1042  end.append(cc+nearGap*n);
1043  gapSize.append(nearGap);
1044 
1045  start2.append(cc);
1046  end2.append(cc-nearGap*n);
1047  gapSize2.append(nearGap);
1048 
1049  // Rays in second direction
1050  start.append(cc);
1051  end.append(cc+nearGap*e2);
1052  gapSize.append(nearGap);
1053 
1054  start2.append(cc);
1055  end2.append(cc-nearGap*e2);
1056  gapSize2.append(nearGap);
1057 
1058  // Rays in third direction
1059  start.append(cc);
1060  end.append(cc+nearGap*e3);
1061  gapSize.append(nearGap);
1062 
1063  start2.append(cc);
1064  end2.append(cc-nearGap*e3);
1065  gapSize2.append(nearGap);
1066  }
1067  }
1068  }
1069 
1070  return start.size()-nOldRays;
1071 }
1072 
1073 
1074 void Foam::meshRefinement::selectGapCandidates
1075 (
1076  const labelList& refineCell,
1077  const label nRefine,
1078 
1079  labelList& cellMap,
1080  labelList& gapShell,
1081  List<FixedList<label, 3>>& shellGapInfo,
1082  List<volumeType>& shellGapMode
1083 ) const
1084 {
1085  const labelList& cellLevel = meshCutter_.cellLevel();
1086  const pointField& cellCentres = mesh_.cellCentres();
1087 
1088  // Collect cells to test
1089  cellMap.setSize(cellLevel.size()-nRefine);
1090  label compactI = 0;
1091 
1092  forAll(cellLevel, cellI)
1093  {
1094  if (refineCell[cellI] == -1)
1095  {
1096  cellMap[compactI++] = cellI;
1097  }
1098  }
1099  Info<< "Selected " << returnReduce(compactI, sumOp<label>())
1100  << " unmarked cells out of "
1101  << mesh_.globalData().nTotalCells() << endl;
1102  cellMap.setSize(compactI);
1103 
1104  // Do test to see whether cells are inside/outside shell with
1105  // applicable specification (minLevel <= celllevel < maxLevel)
1106  shells_.findHigherGapLevel
1107  (
1108  pointField(cellCentres, cellMap),
1109  labelUIndList(cellLevel, cellMap)(),
1110 
1111  gapShell,
1112  shellGapInfo,
1113  shellGapMode
1114  );
1115 
1116  // Compact out hits
1117 
1118  labelList map(shellGapInfo.size());
1119  compactI = 0;
1120  forAll(shellGapInfo, i)
1121  {
1122  if (shellGapInfo[i][2] > 0)
1123  {
1124  map[compactI++] = i;
1125  }
1126  }
1127 
1128  Info<< "Selected " << returnReduce(compactI, sumOp<label>())
1129  << " cells inside gap shells out of "
1130  << mesh_.globalData().nTotalCells() << endl;
1131 
1132  map.setSize(compactI);
1133  cellMap = labelUIndList(cellMap, map)();
1134  gapShell = labelUIndList(gapShell, map)();
1135  shellGapInfo = UIndirectList<FixedList<label, 3>>(shellGapInfo, map)();
1136  shellGapMode = UIndirectList<volumeType>(shellGapMode, map)();
1137 }
1138 
1139 
1140 void Foam::meshRefinement::mergeGapInfo
1141 (
1142  const FixedList<label, 3>& shellGapInfo,
1143  const volumeType shellGapMode,
1144  const FixedList<label, 3>& surfGapInfo,
1145  const volumeType surfGapMode,
1146 
1147  FixedList<label, 3>& gapInfo,
1148  volumeType& gapMode
1149 ) const
1150 {
1151  if (surfGapInfo[0] == 0)
1152  {
1153  gapInfo = shellGapInfo;
1154  gapMode = shellGapMode;
1155  }
1156  else if (shellGapInfo[0] == 0)
1157  {
1158  gapInfo = surfGapInfo;
1159  gapMode = surfGapMode;
1160  }
1161  else
1162  {
1163  // Both specify a level. Does surface level win? Or does information
1164  // need to be merged?
1165 
1166  //gapInfo[0] = max(surfGapInfo[0], shellGapInfo[0]);
1167  //gapInfo[1] = min(surfGapInfo[1], shellGapInfo[1]);
1168  //gapInfo[2] = max(surfGapInfo[2], shellGapInfo[2]);
1169  gapInfo = surfGapInfo;
1170  gapMode = surfGapMode;
1171  }
1172 }
1173 
1174 
1175 Foam::label Foam::meshRefinement::markInternalGapRefinement
1176 (
1177  const scalar planarCos,
1178  const bool spreadGapSize,
1179  const label nAllowRefine,
1180 
1181  labelList& refineCell,
1182  label& nRefine,
1183  labelList& numGapCells,
1184  scalarField& detectedGapSize
1185 ) const
1186 {
1187  detectedGapSize.setSize(mesh_.nCells());
1188  detectedGapSize = GREAT;
1189  numGapCells.setSize(mesh_.nCells());
1190  numGapCells = -1;
1191 
1192  const labelList& cellLevel = meshCutter_.cellLevel();
1193  const pointField& cellCentres = mesh_.cellCentres();
1194  const scalar edge0Len = meshCutter_.level0EdgeLength();
1195 
1196  const List<FixedList<label, 3>>& extendedGapLevel =
1197  surfaces_.extendedGapLevel();
1198  const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
1199  const boolList& extendedGapSelf = surfaces_.gapSelf();
1200 
1201  // Get the gap level for the shells
1202  const labelList maxLevel(shells_.maxGapLevel());
1203 
1204  label oldNRefine = nRefine;
1205 
1206  if (max(maxLevel) > 0)
1207  {
1208  // Collect cells to test
1209  labelList cellMap;
1210  labelList gapShell;
1211  List<FixedList<label, 3>> shellGapInfo;
1212  List<volumeType> shellGapMode;
1213  selectGapCandidates
1214  (
1215  refineCell,
1216  nRefine,
1217 
1218  cellMap,
1219  gapShell,
1220  shellGapInfo,
1221  shellGapMode
1222  );
1223 
1224  // Find nearest point and normal on the surfaces
1225  List<pointIndexHit> nearInfo;
1226  vectorField nearNormal;
1227  labelList nearSurface;
1228  labelList nearRegion;
1229  {
1230  // Now we have both the cell-level and the gap size information. Use
1231  // this to calculate the gap size
1232  scalarField gapSize(cellMap.size());
1233  forAll(cellMap, i)
1234  {
1235  label cellI = cellMap[i];
1236  scalar cellSize = edge0Len/pow(2.0, cellLevel[cellI]);
1237  gapSize[i] = shellGapInfo[i][0]*cellSize;
1238  }
1239 
1240  surfaces_.findNearestRegion
1241  (
1242  identity(surfaces_.surfaces().size()),
1243  pointField(cellCentres, cellMap),
1244  sqr(gapSize),
1245  nearSurface,
1246  nearInfo,
1247  nearRegion,
1248  nearNormal
1249  );
1250  }
1251 
1252 
1253 
1254  DynamicList<label> map(nearInfo.size());
1255  DynamicField<point> rayStart(nearInfo.size());
1256  DynamicField<point> rayEnd(nearInfo.size());
1257  DynamicField<scalar> gapSize(nearInfo.size());
1258 
1259  DynamicField<point> rayStart2(nearInfo.size());
1260  DynamicField<point> rayEnd2(nearInfo.size());
1261  DynamicField<scalar> gapSize2(nearInfo.size());
1262 
1263  label nTestCells = 0;
1264 
1265  forAll(nearInfo, i)
1266  {
1267  if (nearInfo[i].hit())
1268  {
1269  label globalRegionI = surfaces_.globalRegion
1270  (
1271  nearSurface[i],
1272  nearRegion[i]
1273  );
1274 
1275  // Combine info from shell and surface
1276  FixedList<label, 3> gapInfo;
1277  volumeType gapMode;
1278  mergeGapInfo
1279  (
1280  shellGapInfo[i],
1281  shellGapMode[i],
1282 
1283  extendedGapLevel[globalRegionI],
1284  extendedGapMode[globalRegionI],
1285 
1286  gapInfo,
1287  gapMode
1288  );
1289 
1290  // Store wanted number of cells in gap
1291  label cellI = cellMap[i];
1292  label cLevel = cellLevel[cellI];
1293  if (cLevel >= gapInfo[1] && cLevel < gapInfo[2])
1294  {
1295  numGapCells[cellI] = max(numGapCells[cellI], gapInfo[0]);
1296  }
1297 
1298  // Construct one or more rays to test for oppositeness
1299  label nRays = generateRays
1300  (
1301  false,
1302  nearInfo[i].hitPoint(),
1303  nearNormal[i],
1304  gapInfo,
1305  gapMode,
1306 
1307  cellCentres[cellI],
1308  cLevel,
1309 
1310  rayStart,
1311  rayEnd,
1312  gapSize,
1313 
1314  rayStart2,
1315  rayEnd2,
1316  gapSize2
1317  );
1318  if (nRays > 0)
1319  {
1320  nTestCells++;
1321  for (label j = 0; j < nRays; j++)
1322  {
1323  map.append(i);
1324  }
1325  }
1326  }
1327  }
1328 
1329  Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
1330  << " cells for testing out of "
1331  << mesh_.globalData().nTotalCells() << endl;
1332  map.shrink();
1333  rayStart.shrink();
1334  rayEnd.shrink();
1335  gapSize.shrink();
1336 
1337  rayStart2.shrink();
1338  rayEnd2.shrink();
1339  gapSize2.shrink();
1340 
1341  cellMap = labelUIndList(cellMap, map)();
1342  nearNormal = UIndirectList<vector>(nearNormal, map)();
1343  shellGapInfo.clear();
1344  shellGapMode.clear();
1345  nearInfo.clear();
1346  nearSurface.clear();
1347  nearRegion.clear();
1348 
1349 
1350  // Do intersections in pairs
1351  labelList surf1;
1352  List<pointIndexHit> hit1;
1353  vectorField normal1;
1354  surfaces_.findNearestIntersection
1355  (
1356  rayStart,
1357  rayEnd,
1358  surf1,
1359  hit1,
1360  normal1
1361  );
1362 
1363  labelList surf2;
1364  List<pointIndexHit> hit2;
1365  vectorField normal2;
1366  surfaces_.findNearestIntersection
1367  (
1368  rayStart2,
1369  rayEnd2,
1370  surf2,
1371  hit2,
1372  normal2
1373  );
1374 
1375  // Extract cell based gap size
1376  forAll(surf1, i)
1377  {
1378  // Combine selfProx of shell and surfaces. Ignore regions for
1379  // now
1380  const label shelli = gapShell[map[i]];
1381 
1382  bool selfProx = true;
1383  if (shelli != -1)
1384  {
1385  selfProx = shells_.gapSelf()[shelli][0];
1386  }
1387  if (surf1[i] != -1 && selfProx)
1388  {
1389  const label globalRegioni = surfaces_.globalRegion(surf1[i], 0);
1390  selfProx = extendedGapSelf[globalRegioni];
1391  }
1392 
1393  if
1394  (
1395  surf1[i] != -1
1396  && surf2[i] != -1
1397  && (surf2[i] != surf1[i] || selfProx)
1398  )
1399  {
1400  // Found intersections with surface. Check for
1401  // - small gap
1402  // - coplanar normals
1403 
1404  const label cellI = cellMap[i];
1405 
1406  const scalar d2 = magSqr(hit1[i].hitPoint()-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  // Field on cells and faces
1428  List<transportData> cellData(mesh_.nCells());
1429  List<transportData> faceData(mesh_.nFaces());
1430 
1431  // Start of walk
1432  const pointField& faceCentres = mesh_.faceCentres();
1433 
1434  DynamicList<label> frontFaces(mesh_.nFaces());
1435  DynamicList<transportData> frontData(mesh_.nFaces());
1436  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1437  {
1438  label own = mesh_.faceOwner()[faceI];
1439  label nei = mesh_.faceNeighbour()[faceI];
1440 
1441  scalar minSize = min
1442  (
1443  detectedGapSize[own],
1444  detectedGapSize[nei]
1445  );
1446 
1447  if (minSize < GREAT)
1448  {
1449  frontFaces.append(faceI);
1450  frontData.append
1451  (
1452  transportData
1453  (
1454  faceCentres[faceI],
1455  minSize,
1456  0.0
1457  )
1458  );
1459  }
1460  }
1461  for
1462  (
1463  label faceI = mesh_.nInternalFaces();
1464  faceI < mesh_.nFaces();
1465  faceI++
1466  )
1467  {
1468  label own = mesh_.faceOwner()[faceI];
1469 
1470  if (detectedGapSize[own] < GREAT)
1471  {
1472  frontFaces.append(faceI);
1473  frontData.append
1474  (
1475  transportData
1476  (
1477  faceCentres[faceI],
1478  detectedGapSize[own],
1479  0.0
1480  )
1481  );
1482  }
1483  }
1484 
1485  Info<< "Selected "
1486  << returnReduce(frontFaces.size(), sumOp<label>())
1487  << " faces for spreading gap size out of "
1488  << mesh_.globalData().nTotalFaces() << endl;
1489 
1490 
1491  transportData::trackData td(surfaceIndex());
1492 
1493  FaceCellWave<transportData, transportData::trackData> deltaCalc
1494  (
1495  mesh_,
1496  frontFaces,
1497  frontData,
1498  faceData,
1499  cellData,
1500  mesh_.globalData().nTotalCells()+1,
1501  td
1502  );
1503 
1504 
1505  forAll(cellMap, i)
1506  {
1507  label cellI = cellMap[i];
1508  if
1509  (
1510  cellI != -1
1511  && cellData[cellI].valid(deltaCalc.data())
1512  && numGapCells[cellI] != -1
1513  )
1514  {
1515  // Update transported gap size
1516  detectedGapSize[cellI] = min
1517  (
1518  detectedGapSize[cellI],
1519  cellData[cellI].data()
1520  );
1521  }
1522  }
1523  }
1524 
1525 
1526  // Use it
1527  forAll(cellMap, i)
1528  {
1529  label cellI = cellMap[i];
1530 
1531  if (cellI != -1 && numGapCells[cellI] != -1)
1532  {
1533  // Needed gap size
1534  label cLevel = cellLevel[cellI];
1535  scalar cellSize =
1536  meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
1537  scalar neededGapSize = numGapCells[cellI]*cellSize;
1538 
1539  if (neededGapSize > detectedGapSize[cellI])
1540  {
1541  if
1542  (
1543  !markForRefine
1544  (
1545  123,
1546  nAllowRefine,
1547  refineCell[cellI],
1548  nRefine
1549  )
1550  )
1551  {
1552  break;
1553  }
1554  }
1555  }
1556  }
1557 
1558 
1559  if
1560  (
1561  returnReduce(nRefine, sumOp<label>())
1562  > returnReduce(nAllowRefine, sumOp<label>())
1563  )
1564  {
1565  Info<< "Reached refinement limit." << endl;
1566  }
1567  }
1568 
1569  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1570 }
1571 
1572 
1573 Foam::label Foam::meshRefinement::markSmallFeatureRefinement
1574 (
1575  const scalar planarCos,
1576  const label nAllowRefine,
1577  const labelList& neiLevel,
1578  const pointField& neiCc,
1579 
1580  labelList& refineCell,
1581  label& nRefine
1582 ) const
1583 {
1584  const labelList& cellLevel = meshCutter_.cellLevel();
1585  const labelList& surfaceIndices = surfaces_.surfaces();
1586  const List<FixedList<label, 3>>& extendedGapLevel =
1587  surfaces_.extendedGapLevel();
1588  const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
1589  const boolList& extendedGapSelf = surfaces_.gapSelf();
1590 
1591  label oldNRefine = nRefine;
1592 
1593  // Check that we're using any gap refinement
1594  labelList shellMaxLevel(shells_.maxGapLevel());
1595 
1596  if (max(shellMaxLevel) == 0)
1597  {
1598  return 0;
1599  }
1600 
1601  //- Force calculation of tetBasePt
1602  (void)mesh_.tetBasePtIs();
1603  (void)mesh_.cellTree();
1604 
1605 
1606  forAll(surfaceIndices, surfI)
1607  {
1608  label geomI = surfaceIndices[surfI];
1609  const searchableSurface& geom = surfaces_.geometry()[geomI];
1610 
1611 
1612  // Get the element index in a roundabout way. Problem is e.g.
1613  // distributed surface where local indices differ from global
1614  // ones (needed for getRegion call)
1615 
1616  pointField ctrs;
1617  labelList region;
1618  vectorField normal;
1619  {
1620  // Representative local coordinates and bounding sphere
1621  scalarField radiusSqr;
1622  geom.boundingSpheres(ctrs, radiusSqr);
1623 
1624  List<pointIndexHit> info;
1625  geom.findNearest(ctrs, radiusSqr, info);
1626 
1627  forAll(info, i)
1628  {
1629  if (!info[i].hit())
1630  {
1631  FatalErrorInFunction
1632  << "fc:" << ctrs[i]
1633  << " radius:" << radiusSqr[i]
1634  << exit(FatalError);
1635  }
1636  }
1637 
1638  geom.getRegion(info, region);
1639  geom.getNormal(info, normal);
1640  }
1641 
1642  // Do test to see whether triangles are inside/outside shell with
1643  // applicable specification (minLevel <= celllevel < maxLevel)
1644  List<FixedList<label, 3>> shellGapInfo;
1645  List<volumeType> shellGapMode;
1646  labelList gapShell;
1647  shells_.findHigherGapLevel
1648  (
1649  ctrs,
1650  labelList(ctrs.size(), Zero),
1651 
1652  gapShell,
1653  shellGapInfo,
1654  shellGapMode
1655  );
1656 
1657 
1658  DynamicList<label> map(ctrs.size());
1659  DynamicList<label> cellMap(ctrs.size());
1660 
1661  DynamicField<point> rayStart(ctrs.size());
1662  DynamicField<point> rayEnd(ctrs.size());
1663  DynamicField<scalar> gapSize(ctrs.size());
1664 
1665  label nTestCells = 0;
1666 
1667  forAll(ctrs, i)
1668  {
1669  if (shellGapInfo[i][2] > 0)
1670  {
1671  label globalRegionI = surfaces_.globalRegion(surfI, region[i]);
1672 
1673  // Combine info from shell and surface
1674  FixedList<label, 3> gapInfo;
1675  volumeType gapMode;
1676  mergeGapInfo
1677  (
1678  shellGapInfo[i],
1679  shellGapMode[i],
1680 
1681  extendedGapLevel[globalRegionI],
1682  extendedGapMode[globalRegionI],
1683 
1684  gapInfo,
1685  gapMode
1686  );
1687 
1688  //- Option 1: use octree nearest searching inside polyMesh
1689  //label cellI = mesh_.findCell(pt, polyMesh::CELL_TETS);
1690 
1691  //- Option 2: use octree 'inside' searching inside polyMesh. Is
1692  // much faster.
1693  label cellI = -1;
1694  const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
1695  if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
1696  {
1697  cellI = tree.findInside(ctrs[i]);
1698  }
1699 
1700  if (cellI != -1 && refineCell[cellI] == -1)
1701  {
1702  // Construct one or two rays to test for oppositeness
1703  // Note that we always want to use the surface normal
1704  // and not the vector from cell centre to surface point
1705 
1706  label nRays = generateRays
1707  (
1708  ctrs[i],
1709  normal[i],
1710  gapInfo,
1711  gapMode,
1712 
1713  cellLevel[cellI],
1714 
1715  rayStart,
1716  rayEnd
1717  );
1718 
1719  if (nRays > 0)
1720  {
1721  nTestCells++;
1722  for (label j = 0; j < nRays; j++)
1723  {
1724  cellMap.append(cellI);
1725  map.append(i);
1726  }
1727  }
1728  }
1729  }
1730  }
1731 
1732  Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
1733  << " cells containing triangle centres out of "
1734  << mesh_.globalData().nTotalCells() << endl;
1735  map.shrink();
1736  cellMap.shrink();
1737  rayStart.shrink();
1738  rayEnd.shrink();
1739 
1740  ctrs.clear();
1741  region.clear();
1742  shellGapInfo.clear();
1743  shellGapMode.clear();
1744  normal = UIndirectList<vector>(normal, map)();
1745 
1746  // Do intersections.
1747  labelList surfaceHit;
1748  vectorField surfaceNormal;
1749  surfaces_.findNearestIntersection
1750  (
1751  rayStart,
1752  rayEnd,
1753  surfaceHit,
1754  surfaceNormal
1755  );
1756 
1757 
1758  label nOldRefine = 0;
1759 
1760 
1761  forAll(surfaceHit, i)
1762  {
1763  // Combine selfProx of shell and surfaces. Ignore regions for
1764  // now
1765  const label shelli = gapShell[map[i]];
1766  bool selfProx = true;
1767  if (shelli != -1)
1768  {
1769  selfProx = shells_.gapSelf()[shelli][0];
1770  }
1771  if (surfI != -1 && selfProx)
1772  {
1773  const label globalRegioni = surfaces_.globalRegion(surfI, 0);
1774  selfProx = extendedGapSelf[globalRegioni];
1775  }
1776 
1777  if
1778  (
1779  surfaceHit[i] != -1
1780  && (surfaceHit[i] != surfI || selfProx)
1781  )
1782  {
1783  // Found intersection with surface. Check coplanar normals.
1784  label cellI = cellMap[i];
1785 
1786  if (mag(normal[i]&surfaceNormal[i]) > planarCos)
1787  {
1788  if
1789  (
1790  !markForRefine
1791  (
1792  surfaceHit[i],
1793  nAllowRefine,
1794  refineCell[cellI],
1795  nRefine
1796  )
1797  )
1798  {
1799  break;
1800  }
1801  }
1802  }
1803  }
1804 
1805  Info<< "For surface " << geom.name() << " found "
1806  << returnReduce(nRefine-nOldRefine, sumOp<label>())
1807  << " cells in small gaps" << endl;
1808 
1809  if
1810  (
1811  returnReduce(nRefine, sumOp<label>())
1812  > returnReduce(nAllowRefine, sumOp<label>())
1813  )
1814  {
1815  Info<< "Reached refinement limit." << endl;
1816  }
1817  }
1818 
1819  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1820 }
1821 
1822 
1823 // ************************************************************************* //
shellSurfaces.H
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:71
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::primitiveMesh::nInternalFaces
label nInternalFaces() const
Number of internal faces.
Definition: primitiveMeshI.H:78
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:814
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:69
Foam::List::append
void append(const T &val)
Append an element at the end of the list.
Definition: ListI.H:182
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:350
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::nearInfo
Tuple2< scalar, label > nearInfo
Private class for finding nearest.
Definition: sampledMeshedSurface.C:76
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
Number of mesh cells.
Definition: primitiveMeshI.H:96
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
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:1076
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::primitiveMesh::nBoundaryFaces
label nBoundaryFaces() const
Number of boundary faces (== nFaces - nInternalFaces)
Definition: primitiveMeshI.H:84
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:175
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:372
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::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:144
Foam::primitiveMesh::isInternalFace
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
Definition: primitiveMeshI.H:102
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:212
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:1082
Foam::volumeType::OUTSIDE
A location outside the volume.
Definition: volumeType.H:69
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: UIndirectList.H:58