Loading...
Searching...
No Matches
meshRefinementGapRefine.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2015 OpenFOAM Foundation
9 Copyright (C) 2015-2022 OpenCFD Ltd.
10-------------------------------------------------------------------------------
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 hit1[i].hitPoint().distSqr(hit2[i].hitPoint())
492 < Foam::sqr(gapSize[i])
493 )
494 )
495 {
496 if
497 (
498 !markForRefine
499 (
500 surf1[i],
501 nAllowRefine,
502 refineCell[cellI],
503 nRefine
504 )
505 )
506 {
507 break;
508 }
509 }
510 }
511 }
512
513 if
514 (
515 returnReduce(nRefine, sumOp<label>())
516 > returnReduce(nAllowRefine, sumOp<label>())
517 )
518 {
519 Info<< "Reached refinement limit." << endl;
520 }
521 }
522
523 return returnReduce(nRefine-oldNRefine, sumOp<label>());
524}
525
526
527//Foam::meshRefinement::findNearestOppositeOp::findNearestOppositeOp
528//(
529// const indexedOctree<treeDataTriSurface>& tree,
530// const point& oppositePoint,
531// const vector& oppositeNormal,
532// const scalar minCos
533//)
534//:
535// tree_(tree),
536// oppositePoint_(oppositePoint),
537// oppositeNormal_(oppositeNormal),
538// minCos_(minCos)
539//{}
540//
541//
542//void Foam::meshRefinement::findNearestOppositeOp::operator()
543//(
544// const labelUList& indices,
545// const point& sample,
546// scalar& nearestDistSqr,
547// label& minIndex,
548// point& nearestPoint
549//) const
550//{
551// const treeDataTriSurface& shape = tree_.shapes();
552// const triSurface& patch = shape.patch();
553// const pointField& points = patch.points();
554//
555// forAll(indices, i)
556// {
557// const label index = indices[i];
558// const labelledTri& f = patch[index];
559//
560// pointHit nearHit = f.nearestPoint(sample, points);
561// scalar distSqr = sqr(nearHit.distance());
562//
563// if (distSqr < nearestDistSqr)
564// {
565// // Nearer. Check if
566// // - a bit way from other hit
567// // - in correct search cone
568// vector d(nearHit.point()-oppositePoint_);
569// scalar normalDist(d&oppositeNormal_);
570//
571// if (normalDist > Foam::sqr(SMALL) && normalDist/mag(d) > minCos_)
572// {
573// nearestDistSqr = distSqr;
574// minIndex = index;
575// nearestPoint = nearHit.point();
576// }
577// }
578// }
579//}
580//
581//
582//void Foam::meshRefinement::searchCone
583//(
584// const label surfI,
585// labelList& nearMap, // cells
586// scalarField& nearGap, // gap size
587// List<pointIndexHit>& nearInfo, // nearest point on surface
588// List<pointIndexHit>& oppositeInfo // detected point on gap (or miss)
589//) const
590//{
591// const labelList& cellLevel = meshCutter_.cellLevel();
592// const pointField& cellCentres = mesh_.cellCentres();
593// const scalar edge0Len = meshCutter_.level0EdgeLength();
594//
595// const labelList& surfaceIndices = surfaces_.surfaces();
596// const List<FixedList<label, 3>>& extendedGapLevel =
597// surfaces_.extendedGapLevel();
598// const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
599//
600//
601// label geomI = surfaceIndices[surfI];
602// const searchableSurface& geom = surfaces_.geometry()[geomI];
603//
604// const triSurfaceMesh& s = refCast<const triSurfaceMesh>(geom);
605// const indexedOctree<treeDataTriSurface>& tree = s.tree();
606//
607//
608// const scalar searchCos = Foam::cos(degToRad(30.0));
609//
610// // Normals for ray shooting and inside/outside detection
611// vectorField nearNormal;
612// geom.getNormal(nearInfo, nearNormal);
613// // Regions
614// labelList nearRegion;
615// geom.getRegion(nearInfo, nearRegion);
616//
617//
618// // Now loop over all near points and search in the half cone
619// labelList map(nearInfo.size());
620// label compactI = 0;
621//
622// oppositeInfo.setSize(nearInfo.size());
623//
624// forAll(nearInfo, i)
625// {
626// label globalRegionI =
627// surfaces_.globalRegion(surfI, nearRegion[i]);
628//
629// // Get updated gap information now we have the region
630// label nGapCells = extendedGapLevel[globalRegionI][0];
631// label minLevel = extendedGapLevel[globalRegionI][1];
632// label maxLevel = extendedGapLevel[globalRegionI][2];
633// volumeType mode = extendedGapMode[globalRegionI];
634//
635// label cellI = nearMap[i];
636// label cLevel = cellLevel[cellI];
637//
638// if (cLevel >= minLevel && cLevel < maxLevel)
639// {
640// scalar cellSize = edge0Len/pow(2.0, cLevel);
641//
642// // Update gap size
643// nearGap[i] = nGapCells*cellSize;
644//
645// const point& nearPt = nearInfo[i].hitPoint();
646// vector v(cellCentres[cellI]-nearPt);
647// scalar magV = mag(v);
648//
649// // Like with ray shooting we want to
650// // - find triangles up to nearGap away on the wanted side of the
651// // surface
652// // - find triangles up to 0.5*cellSize away on the unwanted side
653// // of the surface. This is for cells straddling the surface
654// // where
655// // the cell centre might be on the wrong side of the surface
656//
657// // Tbd: check that cell centre is inbetween the gap hits
658// // (only if the cell is far enough away)
659//
660// scalar posNormalSize = 0.0;
661// scalar negNormalSize = 0.0;
662//
663// if (mode == volumeType::OUTSIDE)
664// {
665// posNormalSize = nearGap[i];
666// if (magV < 0.5*cellSize)
667// {
668// negNormalSize = 0.5*cellSize;
669// }
670// }
671// else if (mode == volumeType::INSIDE)
672// {
673// if (magV < 0.5*cellSize)
674// {
675// posNormalSize = 0.5*cellSize;
676// }
677// negNormalSize = nearGap[i];
678// }
679// else
680// {
681// posNormalSize = nearGap[i];
682// negNormalSize = nearGap[i];
683// }
684//
685// // Test with positive normal
686// oppositeInfo[compactI] = tree.findNearest
687// (
688// nearPt,
689// sqr(posNormalSize),
690// findNearestOppositeOp
691// (
692// tree,
693// nearPt,
694// nearNormal[i],
695// searchCos
696// )
697// );
698//
699// if (oppositeInfo[compactI].hit())
700// {
701// map[compactI++] = i;
702// }
703// else
704// {
705// // Test with negative normal
706// oppositeInfo[compactI] = tree.findNearest
707// (
708// nearPt,
709// sqr(negNormalSize),
710// findNearestOppositeOp
711// (
712// tree,
713// nearPt,
714// -nearNormal[i],
715// searchCos
716// )
717// );
718//
719// if (oppositeInfo[compactI].hit())
720// {
721// map[compactI++] = i;
722// }
723// }
724// }
725// }
726//
727// Info<< "Selected " << returnReduce(compactI, sumOp<label>())
728// << " hits on the correct side out of "
729// << returnReduce(map.size(), sumOp<label>()) << endl;
730// map.setSize(compactI);
731// oppositeInfo.setSize(compactI);
732//
733// nearMap = labelUIndList(nearMap, map)();
734// nearGap = UIndirectList<scalar>(nearGap, map)();
735// nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
736// nearNormal = UIndirectList<vector>(nearNormal, map)();
737//
738// // Exclude hits which aren't opposite enough. E.g. you might find
739// // a point on a perpendicular wall - but this does not constitute a gap.
740// vectorField oppositeNormal;
741// geom.getNormal(oppositeInfo, oppositeNormal);
742//
743// compactI = 0;
744// forAll(oppositeInfo, i)
745// {
746// if ((nearNormal[i] & oppositeNormal[i]) < -0.707)
747// {
748// map[compactI++] = i;
749// }
750// }
751//
752// Info<< "Selected " << returnReduce(compactI, sumOp<label>())
753// << " hits opposite the nearest out of "
754// << returnReduce(map.size(), sumOp<label>()) << endl;
755// map.setSize(compactI);
756//
757// nearMap = labelUIndList(nearMap, map)();
758// nearGap = UIndirectList<scalar>(nearGap, map)();
759// nearInfo = UIndirectList<pointIndexHit>(nearInfo, map)();
760// oppositeInfo = UIndirectList<pointIndexHit>(oppositeInfo, map)();
761//}
762
763
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 =
1406 hit1[i].hitPoint().distSqr(hit2[i].hitPoint());
1407
1408 if
1409 (
1410 cellI != -1
1411 && (mag(normal1[i]&normal2[i]) > planarCos)
1412 && (d2 < Foam::sqr(gapSize[i]))
1413 )
1414 {
1415 detectedGapSize[cellI] = min
1416 (
1417 detectedGapSize[cellI],
1418 Foam::sqrt(d2)
1419 );
1420 }
1421 }
1422 }
1423
1424 // Spread it
1425 if (spreadGapSize)
1426 {
1427 scalarField boundaryGapSize;
1428 syncTools::swapBoundaryCellList
1429 (
1430 mesh_,
1431 detectedGapSize,
1432 boundaryGapSize
1433 );
1434
1435 // Field on cells and faces
1436 List<transportData> cellData(mesh_.nCells());
1437 List<transportData> faceData(mesh_.nFaces());
1438
1439 // Start of walk
1440 const pointField& faceCentres = mesh_.faceCentres();
1441
1442 DynamicList<label> frontFaces(mesh_.nFaces());
1443 DynamicList<transportData> frontData(mesh_.nFaces());
1444 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1445 {
1446 label own = mesh_.faceOwner()[faceI];
1447 label nei = mesh_.faceNeighbour()[faceI];
1448
1449 scalar minSize = min
1450 (
1451 detectedGapSize[own],
1452 detectedGapSize[nei]
1453 );
1454
1455 if (minSize < GREAT)
1456 {
1457 frontFaces.append(faceI);
1458 frontData.append
1459 (
1460 transportData
1461 (
1462 faceCentres[faceI],
1463 minSize,
1464 0.0
1465 )
1466 );
1467 }
1468 }
1469 for
1470 (
1471 label faceI = mesh_.nInternalFaces();
1472 faceI < mesh_.nFaces();
1473 faceI++
1474 )
1475 {
1476 label own = mesh_.faceOwner()[faceI];
1477
1478 scalar minSize = min
1479 (
1480 detectedGapSize[own],
1481 boundaryGapSize[faceI-mesh_.nInternalFaces()]
1482 );
1483
1484 if (minSize < GREAT)
1485 {
1486 frontFaces.append(faceI);
1487 frontData.append
1488 (
1489 transportData
1490 (
1491 faceCentres[faceI],
1492 minSize,
1493 0.0
1494 )
1495 );
1496 }
1497 }
1498
1499 Info<< "Selected "
1500 << returnReduce(frontFaces.size(), sumOp<label>())
1501 << " faces for spreading gap size out of "
1502 << mesh_.globalData().nTotalFaces() << endl;
1503
1504
1505 transportData::trackData td(surfaceIndex());
1506
1507 FaceCellWave<transportData, transportData::trackData> deltaCalc
1508 (
1509 mesh_,
1510 frontFaces,
1511 frontData,
1512 faceData,
1513 cellData,
1514 mesh_.globalData().nTotalCells()+1,
1515 td
1516 );
1517
1518
1519 forAll(cellMap, i)
1520 {
1521 label cellI = cellMap[i];
1522 if
1523 (
1524 cellI != -1
1525 && cellData[cellI].valid(deltaCalc.data())
1526 && numGapCells[cellI] != -1
1527 )
1528 {
1529 // Update transported gap size
1530 detectedGapSize[cellI] = min
1531 (
1532 detectedGapSize[cellI],
1533 cellData[cellI].data()
1534 );
1535 }
1536 }
1537 }
1538
1539
1540 // Use it
1541 forAll(cellMap, i)
1542 {
1543 label cellI = cellMap[i];
1544
1545 if (cellI != -1 && numGapCells[cellI] != -1)
1546 {
1547 // Needed gap size
1548 label cLevel = cellLevel[cellI];
1549 scalar cellSize =
1550 meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
1551 scalar neededGapSize = numGapCells[cellI]*cellSize;
1552
1553 if (neededGapSize > detectedGapSize[cellI])
1554 {
1555 if
1556 (
1557 !markForRefine
1558 (
1559 123,
1560 nAllowRefine,
1561 refineCell[cellI],
1562 nRefine
1563 )
1564 )
1565 {
1566 break;
1567 }
1568 }
1569 }
1570 }
1571
1572
1573 if
1574 (
1575 returnReduce(nRefine, sumOp<label>())
1576 > returnReduce(nAllowRefine, sumOp<label>())
1577 )
1578 {
1579 Info<< "Reached refinement limit." << endl;
1580 }
1581 }
1582
1583 return returnReduce(nRefine-oldNRefine, sumOp<label>());
1584}
1585
1586
1587Foam::label Foam::meshRefinement::markSmallFeatureRefinement
1588(
1589 const scalar planarCos,
1590 const label nAllowRefine,
1591 const labelList& neiLevel,
1592 const pointField& neiCc,
1593
1594 labelList& refineCell,
1595 label& nRefine
1596) const
1597{
1598 const labelList& cellLevel = meshCutter_.cellLevel();
1599 const labelList& surfaceIndices = surfaces_.surfaces();
1600 const List<FixedList<label, 3>>& extendedGapLevel =
1601 surfaces_.extendedGapLevel();
1602 const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();
1603 const boolList& extendedGapSelf = surfaces_.gapSelf();
1604
1605 label oldNRefine = nRefine;
1606
1607 // Check that we're using any gap refinement
1608 labelList shellMaxLevel(shells_.maxGapLevel());
1609
1610 if (max(shellMaxLevel) == 0)
1611 {
1612 return 0;
1613 }
1614
1615 //- Force calculation of tetBasePt
1616 (void)mesh_.tetBasePtIs();
1617 (void)mesh_.cellTree();
1618
1619
1620 forAll(surfaceIndices, surfI)
1621 {
1622 label geomI = surfaceIndices[surfI];
1623 const searchableSurface& geom = surfaces_.geometry()[geomI];
1624
1625
1626 // Get the element index in a roundabout way. Problem is e.g.
1627 // distributed surface where local indices differ from global
1628 // ones (needed for getRegion call)
1629
1630 pointField ctrs;
1631 labelList region;
1632 vectorField normal;
1633 {
1634 // Representative local coordinates and bounding sphere
1635 scalarField radiusSqr;
1636 geom.boundingSpheres(ctrs, radiusSqr);
1637
1638 List<pointIndexHit> info;
1639 geom.findNearest(ctrs, radiusSqr, info);
1640
1641 forAll(info, i)
1642 {
1643 if (!info[i].hit())
1644 {
1645 FatalErrorInFunction
1646 << "fc:" << ctrs[i]
1647 << " radius:" << radiusSqr[i]
1648 << exit(FatalError);
1649 }
1650 }
1651
1652 geom.getRegion(info, region);
1653 geom.getNormal(info, normal);
1654 }
1655
1656 // Do test to see whether triangles are inside/outside shell with
1657 // applicable specification (minLevel <= celllevel < maxLevel)
1658 List<FixedList<label, 3>> shellGapInfo;
1659 List<volumeType> shellGapMode;
1660 labelList gapShell;
1661 shells_.findHigherGapLevel
1662 (
1663 ctrs,
1664 labelList(ctrs.size(), Zero),
1665
1666 gapShell,
1667 shellGapInfo,
1668 shellGapMode
1669 );
1670
1671
1672 DynamicList<label> map(ctrs.size());
1673 DynamicList<label> cellMap(ctrs.size());
1674
1675 DynamicField<point> rayStart(ctrs.size());
1676 DynamicField<point> rayEnd(ctrs.size());
1677 DynamicField<scalar> gapSize(ctrs.size());
1678
1679 label nTestCells = 0;
1680
1681 forAll(ctrs, i)
1682 {
1683 if (shellGapInfo[i][2] > 0)
1684 {
1685 label globalRegionI = surfaces_.globalRegion(surfI, region[i]);
1686
1687 // Combine info from shell and surface
1688 FixedList<label, 3> gapInfo;
1689 volumeType gapMode;
1690 mergeGapInfo
1691 (
1692 shellGapInfo[i],
1693 shellGapMode[i],
1694
1695 extendedGapLevel[globalRegionI],
1696 extendedGapMode[globalRegionI],
1697
1698 gapInfo,
1699 gapMode
1700 );
1701
1702 //- Option 1: use octree nearest searching inside polyMesh
1703 //label cellI = mesh_.findCell(pt, polyMesh::CELL_TETS);
1704
1705 //- Option 2: use octree 'inside' searching inside polyMesh. Is
1706 // much faster.
1707 label cellI = -1;
1708 const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
1709 if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
1710 {
1711 cellI = tree.findInside(ctrs[i]);
1712 }
1713
1714 if (cellI != -1 && refineCell[cellI] == -1)
1715 {
1716 // Construct one or two rays to test for oppositeness
1717 // Note that we always want to use the surface normal
1718 // and not the vector from cell centre to surface point
1719
1720 label nRays = generateRays
1721 (
1722 ctrs[i],
1723 normal[i],
1724 gapInfo,
1725 gapMode,
1726
1727 cellLevel[cellI],
1728
1729 rayStart,
1730 rayEnd
1731 );
1732
1733 if (nRays > 0)
1734 {
1735 nTestCells++;
1736 for (label j = 0; j < nRays; j++)
1737 {
1738 cellMap.append(cellI);
1739 map.append(i);
1740 }
1741 }
1742 }
1743 }
1744 }
1745
1746 Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
1747 << " cells containing triangle centres out of "
1748 << mesh_.globalData().nTotalCells() << endl;
1749 map.shrink();
1750 cellMap.shrink();
1751 rayStart.shrink();
1752 rayEnd.shrink();
1753
1754 ctrs.clear();
1755 region.clear();
1756 shellGapInfo.clear();
1757 shellGapMode.clear();
1758 normal = UIndirectList<vector>(normal, map)();
1759
1760 // Do intersections.
1761 labelList surfaceHit;
1762 vectorField surfaceNormal;
1763 surfaces_.findNearestIntersection
1764 (
1765 rayStart,
1766 rayEnd,
1767 surfaceHit,
1768 surfaceNormal
1769 );
1770
1771
1772 label nOldRefine = 0;
1773
1774
1775 forAll(surfaceHit, i)
1776 {
1777 // Combine selfProx of shell and surfaces. Ignore regions for
1778 // now
1779 const label shelli = gapShell[map[i]];
1780 bool selfProx = true;
1781 if (shelli != -1)
1782 {
1783 selfProx = shells_.gapSelf()[shelli][0];
1784 }
1785 if (surfI != -1 && selfProx)
1786 {
1787 const label globalRegioni = surfaces_.globalRegion(surfI, 0);
1788 selfProx = extendedGapSelf[globalRegioni];
1789 }
1790
1791 if
1792 (
1793 surfaceHit[i] != -1
1794 && (surfaceHit[i] != surfI || selfProx)
1795 )
1796 {
1797 // Found intersection with surface. Check coplanar normals.
1798 label cellI = cellMap[i];
1799
1800 if (mag(normal[i]&surfaceNormal[i]) > planarCos)
1801 {
1802 if
1803 (
1804 !markForRefine
1805 (
1806 surfaceHit[i],
1807 nAllowRefine,
1808 refineCell[cellI],
1809 nRefine
1810 )
1811 )
1812 {
1813 break;
1814 }
1815 }
1816 }
1817 }
1818
1819 Info<< "For surface " << geom.name() << " found "
1820 << returnReduce(nRefine-nOldRefine, sumOp<label>())
1821 << " cells in small gaps" << endl;
1822
1823 if
1824 (
1825 returnReduce(nRefine, sumOp<label>())
1826 > returnReduce(nAllowRefine, sumOp<label>())
1827 )
1828 {
1829 Info<< "Reached refinement limit." << endl;
1830 }
1831 }
1832
1833 return returnReduce(nRefine-oldNRefine, sumOp<label>());
1834}
1835
1836
1838Foam::label Foam::meshRefinement::markSurfaceFieldRefinement
1839(
1840 const label nAllowRefine,
1841 const labelList& neiLevel,
1842 const pointField& neiCc,
1843
1844 labelList& refineCell,
1845 label& nRefine
1846) const
1847{
1848 const labelList& cellLevel = meshCutter_.cellLevel();
1849 const labelList& surfaceIndices = surfaces_.surfaces();
1850
1851 label oldNRefine = nRefine;
1852
1853 //- Force calculation of tetBasePt
1854 (void)mesh_.tetBasePtIs();
1855 (void)mesh_.cellTree();
1856 const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
1857
1858
1859 forAll(surfaceIndices, surfI)
1860 {
1861 label geomI = surfaceIndices[surfI];
1862 const searchableSurface& geom = surfaces_.geometry()[geomI];
1863
1864 // Get the element index in a roundabout way. Problem is e.g.
1865 // distributed surface where local indices differ from global
1866 // ones (needed for getRegion call)
1867
1868 pointField ctrs;
1869 labelList region;
1870 labelList minLevelField;
1871 {
1872 // Representative local coordinates and bounding sphere
1873 scalarField radiusSqr;
1874 geom.boundingSpheres(ctrs, radiusSqr);
1875
1876 List<pointIndexHit> info;
1877 geom.findNearest(ctrs, radiusSqr, info);
1878
1879 forAll(info, i)
1880 {
1881 if (!info[i].hit())
1882 {
1883 FatalErrorInFunction
1884 << "fc:" << ctrs[i]
1885 << " radius:" << radiusSqr[i]
1886 << exit(FatalError);
1887 }
1888 }
1889
1890 geom.getRegion(info, region);
1891 geom.getField(info, minLevelField);
1892 }
1893
1894 if (minLevelField.size() != geom.size())
1895 {
1896 Pout<< "** no minLevelField" << endl;
1897 continue;
1898 }
1899
1900
1901 label nOldRefine = 0;
1902
1903 forAll(ctrs, i)
1904 {
1905 label cellI = -1;
1906 if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
1907 {
1908 cellI = tree.findInside(ctrs[i]);
1909 }
1910
1911 if
1912 (
1913 cellI != -1
1914 && refineCell[cellI] == -1
1915 && minLevelField[i] > cellLevel[cellI]
1916 )
1917 {
1918 if
1919 (
1920 !markForRefine
1921 (
1922 surfI,
1923 nAllowRefine,
1924 refineCell[cellI],
1925 nRefine
1926 )
1927 )
1928 {
1929 break;
1930 }
1931 }
1932 }
1933
1934 Info<< "For surface " << geom.name() << " found "
1935 << returnReduce(nRefine-nOldRefine, sumOp<label>())
1936 << " cells containing cached refinement field" << endl;
1937
1938 if
1939 (
1940 returnReduce(nRefine, sumOp<label>())
1941 > returnReduce(nAllowRefine, sumOp<label>())
1942 )
1943 {
1944 Info<< "Reached refinement limit." << endl;
1945 }
1946 }
1947
1948 return returnReduce(nRefine-oldNRefine, sumOp<label>());
1949}
1951
1952
1953// ************************************************************************* //
n
label n
Definition TABSMDCalcMethod2.H:31
Foam::DynamicField
Dynamically sized Field. Similar to DynamicList, but inheriting from a Field instead of a List.
Definition DynamicField.H:62
Foam::DynamicList
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition DynamicList.H:68
Foam::FaceCellWave
Wave propagation of information through grid. Every iteration information goes through one layer of c...
Definition FaceCellWave.H:209
Foam::FixedList
A 1D vector of objects of type <T> with a fixed length <N>.
Definition FixedList.H:73
Foam::List
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition List.H:72
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition ListI.H:133
Foam::UIndirectList
A List with indirect addressing. Like IndirectList but does not store addressing.
Definition IndirectList.H:76
Foam::indexedOctree
Non-pointer based hierarchical recursive searching.
Definition indexedOctree.H:298
Foam::refineCell
Container with cells to refine. Refinement given as single direction.
Definition refineCell.H:53
Foam::searchableSurface
Base class of (analytical or triangulated) surface. Encapsulates all the search routines....
Definition searchableSurface.H:68
Foam::syncTools::swapBoundaryCellList
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData, const bool parRun=UPstream::parRun())
Extract and swap to obtain neighbour cell values for all boundary faces.
Definition syncToolsTemplates.C:1425
Foam::transportData::trackData
Class used to pass additional data in.
Definition transportData.H:60
Foam::transportData
Holds information (coordinate and distance). Walks out 0.5*distance.
Definition transportData.H:53
Foam::volumeType
An enumeration wrapper for classification of a location as being inside/outside of a volume.
Definition volumeType.H:56
Foam::volumeType::OUTSIDE
@ OUTSIDE
A location outside the volume.
Definition volumeType.H:66
Foam::volumeType::MIXED
@ MIXED
A location that is partly inside and outside.
Definition volumeType.H:67
Foam::volumeType::INSIDE
@ INSIDE
A location inside the volume.
Definition volumeType.H:65
FatalErrorInFunction
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition error.H:600
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))
td
wallPoints::trackData td(isBlockedFace, regionToBlockSize)
oldNRefine
const label oldNRefine
Definition meshRefinementBlock.C:2673
Foam::SVG::end
const char * end
Definition SVGTools.H:223
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition hashSets.C:40
Foam::labelList
List< label > labelList
A List of labels.
Definition List.H:62
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition dimensionedSymmTensor.C:44
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition IndirectList.H:65
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 POSIX.C:775
Foam::scalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Definition primitiveFieldsFwd.H:46
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition dimensionedScalar.C:68
Foam::returnReduce
T returnReduce(const T &value, BinaryOp bop, const int tag=UPstream::msgType(), const int communicator=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
Definition PstreamReduceOps.H:333
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition Ostream.H:519
Foam::sqrt
dimensionedScalar sqrt(const dimensionedScalar &ds)
Definition dimensionedScalar.C:137
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:26
Foam::identity
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i), works like std::iota() but returning a...
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition primitiveFieldsFwd.H:48
Foam::point
vector point
Point is a vector.
Definition point.H:37
Foam::boolList
List< bool > boolList
A List of bools.
Definition List.H:60
Foam::Zero
static constexpr const zero Zero
Global zero (0).
Definition zero.H:127
Foam::FatalError
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition pointFieldFwd.H:38
Foam::vector
Vector< scalar > vector
Definition vector.H:57
Foam::exit
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition errorManip.H:125
tree
Tree tree(triangles.begin(), triangles.end())
e
volScalarField & e
Definition createFields.H:11
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition stdFoam.H:299