50 void Foam::snappyVoxelMeshDriver::addNeighbours
55 DynamicList<labelVector>& front
66 if (cellLevel[voxeli-off[dir]] == -1)
71 if (voxel[dir] < n_[dir]-1)
75 if (cellLevel[voxeli+off[dir]] == -1)
92 for (label
k = 0;
k < n_[2];
k++)
94 const label start1 = voxeli;
95 for (label j = 0; j < n_[1]; j++)
97 const label start0 = voxeli;
98 for (label i = 0; i < n_[0]; i++)
104 voxeli = start0 + off[1];
106 voxeli = start1 + off[2];
112 void Foam::snappyVoxelMeshDriver::isInside
118 const fvMesh&
mesh = meshRefiner_.mesh();
120 isVoxelInMesh.setSize(cc.size());
130 isVoxelInMesh[voxeli] = (celli != -1);
140 for (label celli = 0; celli <
mesh.
nCells(); celli++)
142 const cell& cFaces =
cells[celli];
146 const face&
f = faces[cFaces[cFacei]];
167 void Foam::snappyVoxelMeshDriver::markSurfaceRefinement
175 const refinementSurfaces&
s = meshRefiner_.surfaces();
178 label geomi =
s.surfaces()[surfi];
179 const searchableSurface& geom =
s.geometry()[geomi];
181 if (isA<triSurface>(geom))
183 const triSurface& ts = refCast<const triSurface>(geom);
188 label regioni = ts[trii].region();
189 label globalRegioni =
s.regionOffset()[surfi]+regioni;
190 const boundBox triBb(
points, ts[trii],
false);
193 label level =
s.minLevel()[globalRegioni];
219 void Foam::snappyVoxelMeshDriver::findVoxels
226 voxels.setSize(locationsOutsideMesh.size());
228 forAll(locationsOutsideMesh, loci)
230 const point& pt = locationsOutsideMesh[loci];
233 if (voxeli == -1 || voxelLevel[voxeli] ==
labelMax)
236 << pt <<
" is outside mesh with bounding box "
241 voxels[loci] = voxeli;
247 void Foam::snappyVoxelMeshDriver::floodFill
249 const label startVoxeli,
250 const label newLevel,
254 DynamicList<labelVector> front;
257 DynamicList<labelVector> newFront;
261 for (
const auto& voxel : front)
264 if (voxelLevel[voxeli] == -1)
266 voxelLevel[voxeli] = 0;
277 if (newFront.empty())
281 front.transfer(newFront);
286 void Foam::snappyVoxelMeshDriver::max
296 for (label
k = 0;
k < n_[2];
k++)
298 const label start1 = voxeli;
299 for (label j = 0; j < n_[1]; j++)
301 const label start0 = voxeli;
302 for (label i = 0; i < n_[0]; i++)
311 voxeli = start0 + off[1];
313 voxeli = start1 + off[2];
325 for (
const auto level : voxelLevel)
337 for (label
k = 0;
k < n_[2];
k++)
339 const label start1 = voxeli;
340 for (label j = 0; j < n_[1]; j++)
342 const label start0 = voxeli;
343 for (label i = 0; i < n_[0]; i++)
345 label level = voxelLevel[voxeli];
347 if (level != -1 && level !=
labelMax)
353 voxeli = start0 + off[1];
355 voxeli = start1 + off[2];
364 Foam::snappyVoxelMeshDriver::snappyVoxelMeshDriver
371 meshRefiner_(meshRefiner),
372 globalToMasterPatch_(globalToMasterPatch),
373 globalToSlavePatch_(globalToSlavePatch),
374 bb_(meshRefiner_.mesh().bounds())
379 const labelListList& featLevels = meshRefiner_.features().levels();
386 const labelList& surfaceLevels = meshRefiner_.surfaces().maxLevel();
390 maxLevel =
Foam::max(maxLevel, meshRefiner_.shells().maxLevel());
392 const scalar level0Len = meshRefiner_.meshCutter().level0EdgeLength();
397 <<
"Cell size estimate :" <<
nl
399 <<
setw(2) << label(0) <<
setw(oldWidth)
400 <<
" : " << level0Len <<
nl
402 <<
setw(2) << maxLevel <<
setw(oldWidth)
403 <<
" : " << level0Len/
pow(2.0, maxLevel) <<
nl
408 const vector meshSpan(bb_.span());
411 round(meshSpan.x()/level0Len),
412 round(meshSpan.y()/level0Len),
413 round(meshSpan.z()/level0Len)
415 label nTot = n_.x()*n_.y()*n_.z();
416 while (nTot < 1000000)
419 nTot = n_.x()*n_.y()*n_.z();
422 Info<<
"Voxellating initial mesh : " << n_ <<
nl <<
endl;
427 Info<<
"Voxel refinement :" <<
nl
428 <<
" Initial : (" << nTot <<
')' <<
endl;
431 isInside(cc, isVoxelInMesh);
435 voxelLevel_.setSize(nTot, -1);
436 globalRegion_.setSize(nTot, -1);
439 forAll(isVoxelInMesh, voxeli)
441 if (!isVoxelInMesh[voxeli])
444 globalRegion_[voxeli] = -1;
450 Info<<
" After removing outside cells : " <<
count(voxelLevel_)
464 const scalar level0Len = meshRefiner_.meshCutter().level0EdgeLength();
470 isInside(cc, isVoxelInMesh);
475 markSurfaceRefinement(voxelLevel_, globalRegion_);
479 Info<<
" After surface refinement : " <<
count(voxelLevel_)
493 labelList outsideOldLevel(outsideMeshVoxels.size(), -1);
494 forAll(outsideMeshVoxels, loci)
496 label voxeli = outsideMeshVoxels[loci];
499 outsideOldLevel[loci] = voxelLevel_[outsideMeshVoxels[loci]];
500 if (outsideOldLevel[loci] >= 0)
503 << outsidePoints[loci]
504 <<
" is inside mesh or close to surface" <<
endl;
519 forAll(insideMeshVoxels, loci)
521 label voxeli = insideMeshVoxels[loci];
524 if (voxelLevel_[voxeli] != -1)
528 <<
" is marked as a surface voxel " << voxeli
529 <<
" with cell level " << voxelLevel_[voxeli] <<
endl;
534 floodFill(voxeli, 0, voxelLevel_);
541 Info<<
" After keeping inside voxels : " <<
count(voxelLevel_)
548 forAll(outsideMeshVoxels, loci)
550 label voxeli = outsideMeshVoxels[loci];
551 if (voxeli >= 0 && voxelLevel_[voxeli] != outsideOldLevel[loci])
554 << outsidePoints[loci]
555 <<
" is reachable from an inside location" <<
nl
556 <<
"Either your locations are too close to the"
557 <<
" geometry or there might be a leak in the"
558 <<
" geometry" <<
endl;
566 meshRefiner_.shells().findHigherLevel(cc, voxelLevel_, maxLevel);
569 max(maxLevel, voxelLevel_);
576 Info<<
" After shell refinement : " << levelCounts <<
endl;
580 const vector meshSpan(bb_.span());
588 forAll(levelCounts, leveli)
590 const scalar
s = level0Len/
pow(2.0, leveli);
591 const scalar nCellsPerVoxel
597 cellCount += levelCounts[leveli]*nCellsPerVoxel;
599 Info<<
"Estimated cell count : " << cellCount <<
endl;