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);
136 for (label celli = 0; celli <
mesh.
nCells(); ++celli)
155 void Foam::snappyVoxelMeshDriver::markSurfaceRefinement
163 const refinementSurfaces&
s = meshRefiner_.surfaces();
166 label geomi =
s.surfaces()[surfi];
167 const searchableSurface& geom =
s.geometry()[geomi];
169 if (isA<triSurface>(geom))
171 const triSurface& ts = refCast<const triSurface>(geom);
174 for (
const labelledTri& tri : ts)
176 label regioni = tri.region();
177 label globalRegioni =
s.regionOffset()[surfi]+regioni;
178 const boundBox triBb(tri.box(
points));
181 label level =
s.minLevel()[globalRegioni];
207 void Foam::snappyVoxelMeshDriver::findVoxels
214 voxels.setSize(locationsOutsideMesh.size());
216 forAll(locationsOutsideMesh, loci)
218 const point& pt = locationsOutsideMesh[loci];
221 if (voxeli == -1 || voxelLevel[voxeli] ==
labelMax)
224 << pt <<
" is outside mesh with bounding box " 229 voxels[loci] = voxeli;
235 void Foam::snappyVoxelMeshDriver::floodFill
237 const label startVoxeli,
238 const label newLevel,
242 DynamicList<labelVector> front;
245 DynamicList<labelVector> newFront;
249 for (
const auto& voxel : front)
252 if (voxelLevel[voxeli] == -1)
254 voxelLevel[voxeli] = 0;
265 if (newFront.empty())
269 front.transfer(newFront);
274 void Foam::snappyVoxelMeshDriver::max
284 for (label
k = 0;
k < n_[2];
k++)
286 const label start1 = voxeli;
287 for (label j = 0; j < n_[1]; j++)
289 const label start0 = voxeli;
290 for (label i = 0; i < n_[0]; i++)
299 voxeli = start0 + off[1];
301 voxeli = start1 + off[2];
313 for (
const auto level : voxelLevel)
325 for (label
k = 0;
k < n_[2];
k++)
327 const label start1 = voxeli;
328 for (label j = 0; j < n_[1]; j++)
330 const label start0 = voxeli;
331 for (label i = 0; i < n_[0]; i++)
333 label level = voxelLevel[voxeli];
335 if (level != -1 && level !=
labelMax)
341 voxeli = start0 + off[1];
343 voxeli = start1 + off[2];
352 Foam::snappyVoxelMeshDriver::snappyVoxelMeshDriver
359 meshRefiner_(meshRefiner),
360 globalToMasterPatch_(globalToMasterPatch),
361 globalToSlavePatch_(globalToSlavePatch),
362 bb_(meshRefiner_.
mesh().bounds())
385 <<
"Cell size estimate :" <<
nl 387 <<
setw(2) << label(0) <<
setw(oldWidth)
388 <<
" : " << level0Len <<
nl 390 <<
setw(2) << maxLevel <<
setw(oldWidth)
391 <<
" : " << level0Len/
pow(2.0, maxLevel) <<
nl 399 round(meshSpan.x()/level0Len),
400 round(meshSpan.y()/level0Len),
401 round(meshSpan.z()/level0Len)
403 label nTot = n_.x()*n_.y()*n_.z();
404 while (nTot < 1000000)
407 nTot = n_.x()*n_.y()*n_.z();
410 Info<<
"Voxellating initial mesh : " << n_ <<
nl <<
endl;
412 tmp<pointField> tcc(voxelCentres());
415 Info<<
"Voxel refinement :" <<
nl 416 <<
" Initial : (" << nTot <<
')' <<
endl;
419 isInside(cc, isVoxelInMesh);
424 globalRegion_.
setSize(nTot, -1);
427 forAll(isVoxelInMesh, voxeli)
429 if (!isVoxelInMesh[voxeli])
432 globalRegion_[voxeli] = -1;
438 Info<<
" After removing outside cells : " << count(voxelLevel_)
449 const refinementParameters& refineParams
452 const scalar level0Len = meshRefiner_.meshCutter().level0EdgeLength();
454 tmp<pointField> tcc(voxelCentres());
458 isInside(cc, isVoxelInMesh);
463 markSurfaceRefinement(voxelLevel_, globalRegion_);
467 Info<<
" After surface refinement : " <<
count(voxelLevel_)
473 const pointField& outsidePoints = refineParams.locationsOutsideMesh();
481 labelList outsideOldLevel(outsideMeshVoxels.size(), -1);
482 forAll(outsideMeshVoxels, loci)
484 label voxeli = outsideMeshVoxels[loci];
487 outsideOldLevel[loci] = voxelLevel_[outsideMeshVoxels[loci]];
488 if (outsideOldLevel[loci] >= 0)
491 << outsidePoints[loci]
492 <<
" is inside mesh or close to surface" <<
endl;
503 refineParams.locationsInMesh(),
507 forAll(insideMeshVoxels, loci)
509 label voxeli = insideMeshVoxels[loci];
512 if (voxelLevel_[voxeli] != -1)
515 << refineParams.locationsInMesh()[loci]
516 <<
" is marked as a surface voxel " << voxeli
517 <<
" with cell level " << voxelLevel_[voxeli] <<
endl;
522 floodFill(voxeli, 0, voxelLevel_);
529 Info<<
" After keeping inside voxels : " <<
count(voxelLevel_)
536 forAll(outsideMeshVoxels, loci)
538 label voxeli = outsideMeshVoxels[loci];
539 if (voxeli >= 0 && voxelLevel_[voxeli] != outsideOldLevel[loci])
542 << outsidePoints[loci]
543 <<
" is reachable from an inside location" <<
nl 544 <<
"Either your locations are too close to the" 545 <<
" geometry or there might be a leak in the" 546 <<
" geometry" <<
endl;
554 meshRefiner_.shells().findHigherLevel(cc, voxelLevel_, maxLevel);
557 max(maxLevel, voxelLevel_);
564 Info<<
" After shell refinement : " << levelCounts <<
endl;
568 const vector meshSpan(bb_.span());
576 forAll(levelCounts, leveli)
578 const scalar
s = level0Len/
pow(2.0, leveli);
579 const scalar nCellsPerVoxel
585 cellCount += levelCounts[leveli]*nCellsPerVoxel;
587 Info<<
"Estimated cell count : " << cellCount <<
endl;
static void listCombineGather(UList< T > &values, const CombineOp &cop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Combines List elements.
static labelVector offset(const labelVector &nDivs)
Change in combined voxel index for change in components.
const refinementSurfaces & surfaces() const
Reference to surface search engines.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
constexpr char nl
The newline '\n' character (0x0a)
const shellSurfaces & shells() const
Reference to refinement shells (regions)
Ostream & endl(Ostream &os)
Add newline and flush stream.
label maxLevel() const
Highest shell level.
label k
Boltzmann constant.
scalar level0EdgeLength() const
Typical edge length between unrefined points.
void doRefine(const refinementParameters &refineParams)
#define forAll(list, i)
Loop across all elements in list.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
static label index(const labelVector &nDivs, const labelVector &voxel)
Find cells. Returns number of cells found.
vectorField pointField
pointField is a vectorField.
void setSize(const label n)
Alias for resize()
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
OSstream Sout
OSstream wrapped stdout (std::cout)
const globalMeshData & globalData() const
Return parallel info (demand-driven)
static labelVector index3(const labelVector &nDivs, const label voxeli)
Combined voxel index to individual indices.
const refinementFeatures & features() const
Reference to feature edge mesh.
const hexRef8 & meshCutter() const
Reference to meshcutting engine.
Istream and Ostream manipulators taking arguments.
static point centre(const boundBox &bb, const labelVector &nDivs, const labelVector &voxel)
Voxel index to voxel centre.
defineTypeNameAndDebug(combustionModel, 0)
label nTotalCells() const noexcept
Total global number of mesh cells.
Helper class which maintains intersections of (changing) mesh with (static) surfaces.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
label nCells() const noexcept
Number of mesh cells.
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
static void fill(Container &elems, const boundBox &bb, const labelVector &nDivs, const boundBox &subBb, const Type val)
Fill voxels indicated by bounding box.
vector span() const
The bounding box span (from minimum to maximum)
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Vector< label > labelVector
Vector of labels.
Omanip< int > setw(const int i)
virtual int width() const override
Get width of output field.
List< label > labelList
A List of labels.
A class for managing temporary objects.
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))
List< bool > boolList
A List of bools.
boundBox cellBb(const label celli) const
The bounding box for given cell index.
const labelListList & levels() const
Per featureEdgeMesh the list of level.
const labelList & maxLevel() const
From global region number to refinement level.