50void Foam::snappyVoxelMeshDriver::addNeighbours
52 const labelList& cellLevel,
53 const labelVector& voxel,
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];
112void 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]];
167void 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];
219void 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;
247void 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);
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];
371 meshRefiner_(meshRefiner),
372 globalToMasterPatch_(globalToMasterPatch),
373 globalToSlavePatch_(globalToSlavePatch),
374 bb_(meshRefiner_.
mesh().bounds())
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
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);
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_)
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_);
572 const labelList levelCounts(count(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;
Istream and Ostream manipulators taking arguments.
void setSize(const label n)
Alias for resize()
virtual int width() const
Get width of output field.
virtual char fill() const =0
Get padding character.
scalar centre() const
Mid-point location, zero for an empty list.
static void listCombineGather(const List< commsStruct > &comms, List< T > &values, const CombineOp &cop, const int tag, const label comm)
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
void size(const label n)
Older name for setAddressableSize.
const Cmpt & z() const
Access to the vector z component.
const Cmpt & y() const
Access to the vector y component.
const Cmpt & x() const
Access to the vector x component.
static const boundBox invertedBox
A large inverted boundBox: min/max == +/- ROOTVGREAT.
vector span() const
The bounding box span (from minimum to maximum)
label nTotalCells() const noexcept
Return total number of cells in decomposed mesh.
scalar level0EdgeLength() const
Typical edge length between unrefined points.
const vector & offset() const noexcept
Offset vector (from patch faces to destination mesh objects)
Helper class which maintains intersections of (changing) mesh with (static) surfaces.
const hexRef8 & meshCutter() const
Reference to meshcutting engine.
const refinementSurfaces & surfaces() const
Reference to surface search engines.
const shellSurfaces & shells() const
Reference to refinement shells (regions)
const refinementFeatures & features() const
Reference to feature edge mesh.
static const complex max
complex (VGREAT,VGREAT)
virtual const faceList & faces() const
Return raw faces.
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
const globalMeshData & globalData() const
Return parallel info.
virtual const pointField & points() const
Return raw points.
label nCells() const noexcept
Number of mesh cells.
const cellList & cells() const
const labelListList & levels() const
Per featureEdgeMesh the list of level.
Simple container to keep together refinement specific information.
const pointField & locationsOutsideMesh() const
Optional points which are checked to be outside the mesh.
const pointField & locationsInMesh() const
Areas to keep.
const labelList & maxLevel() const
From global region number to refinement level.
label maxLevel() const
Highest shell level.
Equivalent of snappyRefineDriver but operating on a voxel mesh.
void doRefine(const refinementParameters &refineParams)
splitCell * master() const
A class for managing temporary objects.
static labelVector index3(const labelVector &nDivs, const label voxeli)
Combined voxel index to individual indices.
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
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))
#define WarningInFunction
Report a warning using Foam::Warning.
unsigned int count(const UList< bool > &bools, const bool val=true)
Count number of 'true' entries.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
List< label > labelList
A List of labels.
List< cell > cellList
A List of cells.
Vector< label > labelVector
Vector of labels.
messageStream Info
Information stream (stdout output on master, null elsewhere)
vectorField pointField
pointField is a vectorField.
vector point
Point is a vector.
OSstream Sout
OSstream wrapped stdout (std::cout)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Omanip< int > setw(const int i)
Ostream & endl(Ostream &os)
Add newline and flush stream.
List< bool > boolList
A List of bools.
List< face > faceList
A List of faces.
constexpr char nl
The newline '\n' character (0x0a)
#define forAll(list, i)
Loop across all elements in list.