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

Copyright © 2011-2018  OpenFOAM | OPENFOAM® is a registered trademark of OpenCFD Ltd.

Generated by   doxygen 1.8.17