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