50 void Foam::autoDensity::writeOBJ
52 const treeBoundBox& bb,
56 OFstream str(time().
path()/
name +
".obj");
62 for (
const point& pt : bbPoints)
69 str <<
"l " <<
e[0] + 1 <<
' ' <<
e[1] + 1 <<
nl;
74 bool Foam::autoDensity::combinedOverlaps(
const treeBoundBox& box)
const 79 decomposition().overlapsThisProcessor(box)
80 || geometryToConformTo().overlaps(box);
83 return geometryToConformTo().overlaps(box);
87 bool Foam::autoDensity::combinedInside(
const point&
p)
const 92 decomposition().positionOnThisProcessor(
p)
93 && geometryToConformTo().inside(
p);
96 return geometryToConformTo().inside(
p);
108 return geometryToConformTo().wellInside
111 minimumSurfaceDistanceCoeffSqr_*
sqr(sizes)
115 Field<bool> inside(
pts.
size(),
true);
122 geometryToConformTo().wellInside
125 minimumSurfaceDistanceCoeffSqr_*
sqr(sizes)
131 decomposition().positionOnThisProcessor(
pts)
150 inside[i] = (insideA[i] && insideB[i]);
157 bool Foam::autoDensity::combinedWellInside
167 inside = decomposition().positionOnThisProcessor(
p);
174 && geometryToConformTo().wellInside
177 minimumSurfaceDistanceCoeffSqr_*
sqr(size)
184 Foam::label Foam::autoDensity::recurseAndFill
186 DynamicList<Vb::Point>& initialPoints,
187 const treeBoundBox& bb,
196 treeBoundBox subBB = bb.subBbox(i);
198 word newName = recursionName +
"_" +
Foam::name(i);
202 if (combinedOverlaps(subBB))
226 word(newName +
"_overlap")
229 Pout<< newName +
"_overlap " << subBB <<
endl;
232 if (!fillBox(initialPoints, subBB,
true))
249 else if (combinedInside(subBB.centre()))
259 Pout<< newName +
"_inside " << subBB <<
endl;
262 if (!fillBox(initialPoints, subBB,
false))
291 return treeDepth + 1;
295 bool Foam::autoDensity::fillBox
297 DynamicList<Vb::Point>& initialPoints,
298 const treeBoundBox& bb,
302 const conformationSurfaces& geometry = geometryToConformTo();
304 label initialSize = initialPoints.size();
306 scalar maxCellSize = -GREAT;
308 scalar minCellSize = GREAT;
310 scalar maxDensity = 1/
pow3(minCellSize);
312 scalar volumeAdded = 0.0;
318 scalar totalVolume = bb.volume();
320 label trialPoints = 0;
322 bool wellInside =
false;
334 geometry.findSurfaceNearest
346 Pout<<
"box wellInside, no need to sample surface." <<
endl;
353 if (!overlapping && !wellInside)
362 scalarField cornerSizes = cellShapeControls().cellSize(corners);
364 Field<bool> insideCorners = combinedWellInside(corners, cornerSizes);
371 scalar
s = cornerSizes[i];
380 minCellSize =
max(
s, minCellSizeLimit_);
383 if (maxCellSize/minCellSize > maxSizeRatio_)
387 Pout<<
"Abort fill at corner sample stage," 388 <<
" minCellSize " << minCellSize
389 <<
" maxCellSize " << maxCellSize
390 <<
" maxSizeRatio " << maxCellSize/minCellSize
397 if (!insideCorners[i])
404 Pout<<
"Inside box found to have some non-wellInside " 405 <<
"corners, using overlapping fill." 419 label nLine = 6*(surfRes_ - 2);
425 for (label i = 0; i < surfRes_; i++)
429 for (label j = 1; j < surfRes_ - 1 ; j++)
439 delta.x()*(surfRes_ - 1),
453 delta.y()*(surfRes_ - 1),
467 delta.z()*(surfRes_ - 1)
471 lineSizes = cellShapeControls().cellSize(linePoints);
473 Field<bool> insideLines = combinedWellInside
482 scalar
s = lineSizes[i];
491 minCellSize =
max(
s, minCellSizeLimit_);
494 if (maxCellSize/minCellSize > maxSizeRatio_)
498 Pout<<
"Abort fill at surface sample stage, " 499 <<
" minCellSize " << minCellSize
500 <<
" maxCellSize " << maxCellSize
501 <<
" maxSizeRatio " << maxCellSize/minCellSize
516 Pout<<
"Inside box found to have some non-" 517 <<
"wellInside surface points, using " 518 <<
"overlapping fill." 536 volRes_*volRes_*volRes_,
544 for (label i = 0; i < volRes_; i++)
546 for (label j = 0; j < volRes_; j++)
548 for (label
k = 0;
k < volRes_;
k++)
573 scalarField sampleSizes = cellShapeControls().cellSize(samplePoints);
589 scalar
s = sampleSizes[i];
598 minCellSize =
max(
s, minCellSizeLimit_);
601 if (maxCellSize/minCellSize > maxSizeRatio_)
605 Pout<<
"Abort fill at sample stage," 606 <<
" minCellSize " << minCellSize
607 <<
" maxCellSize " << maxCellSize
608 <<
" maxSizeRatio " << maxCellSize/minCellSize
621 Pout<<
"No sample points found inside box" <<
endl;
629 Pout<< scalar(nInside)/scalar(samplePoints.size())
630 <<
" full overlapping box" <<
endl;
633 totalVolume *= scalar(nInside)/scalar(samplePoints.size());
637 Pout<<
"Total volume to fill = " << totalVolume <<
endl;
644 maxDensity = 1/
pow3(
max(minCellSize, SMALL));
652 const point&
p = samplePoints[i];
654 scalar localSize = sampleSizes[i];
656 scalar localDensity = 1/
pow3(localSize);
667 if (localDensity/maxDensity >
rndGen().sample01<scalar>())
669 scalar localVolume = 1/localDensity;
671 if (volumeAdded + localVolume > totalVolume)
677 scalar addProbability =
678 (totalVolume - volumeAdded)/localVolume;
684 Pout<<
"totalVolume " << totalVolume <<
nl 685 <<
"volumeAdded " << volumeAdded <<
nl 686 <<
"localVolume " << localVolume <<
nl 687 <<
"addProbability " << addProbability <<
nl 692 if (addProbability > r)
705 volumeAdded += localVolume;
713 volumeAdded += localVolume;
719 if (volumeAdded < totalVolume)
723 Pout<<
"Adding random points, remaining volume " 724 << totalVolume - volumeAdded
728 maxDensity = 1/
pow3(
max(minCellSize, SMALL));
736 scalar localSize = cellShapeControls().cellSize(
p);
738 bool insidePoint =
false;
747 insidePoint = combinedWellInside(
p, localSize);
752 if (localSize > maxCellSize)
754 maxCellSize = localSize;
757 if (localSize < minCellSize)
759 minCellSize =
max(localSize, minCellSizeLimit_);
761 localSize = minCellSize;
765 maxDensity = 1/
pow3(
max(minCellSize, SMALL));
768 if (maxCellSize/minCellSize > maxSizeRatio_)
772 Pout<<
"Abort fill at random fill stage," 773 <<
" minCellSize " << minCellSize
774 <<
" maxCellSize " << maxCellSize
775 <<
" maxSizeRatio " << maxCellSize/minCellSize
781 initialPoints.resize(initialSize);
786 scalar localDensity = 1/
pow3(
max(localSize, SMALL));
790 if (localDensity/maxDensity >
rndGen().sample01<scalar>())
792 scalar localVolume = 1/localDensity;
794 if (volumeAdded + localVolume > totalVolume)
800 scalar addProbability =
801 (totalVolume - volumeAdded)/localVolume;
807 Pout<<
"totalVolume " << totalVolume <<
nl 808 <<
"volumeAdded " << volumeAdded <<
nl 809 <<
"localVolume " << localVolume <<
nl 810 <<
"addProbability " << addProbability <<
nl 815 if (addProbability > r)
828 volumeAdded += localVolume;
836 volumeAdded += localVolume;
842 globalTrialPoints_ += trialPoints;
847 <<
" locations queried, " << initialPoints.size() - initialSize
848 <<
" points placed, (" 849 << scalar(initialPoints.size() - initialSize)
850 /scalar(
max(trialPoints, 1))
851 <<
" success rate)." <<
nl 852 <<
"minCellSize " << minCellSize
853 <<
", maxCellSize " << maxCellSize
854 <<
", ratio " << maxCellSize/minCellSize
866 const dictionary& initialPointsDict,
869 const conformationSurfaces& geometryToConformTo,
870 const cellShapeControl& cellShapeControls,
871 const autoPtr<backgroundMeshDecomposition>& decomposition
884 globalTrialPoints_(0),
887 detailsDict().getOrDefault<scalar>(
"minCellSizeLimit", 0)
889 minLevels_(detailsDict().
get<label>(
"minLevels")),
890 maxSizeRatio_(detailsDict().
get<scalar>(
"maxSizeRatio")),
891 volRes_(detailsDict().
get<label>(
"sampleResolution")),
894 detailsDict().getOrDefault<label>(
"surfaceSampleResolution", volRes_)
897 if (maxSizeRatio_ <= 1.0)
902 <<
"The maxSizeRatio must be greater than one to be sensible, " 903 <<
"setting to " << maxSizeRatio_
937 Pout<<
" Filling box " << hierBB <<
endl;
940 label treeDepth = recurseAndFill
954 reduce(nInitialPoints, sumOp<label>());
955 reduce(globalTrialPoints_, sumOp<label>());
959 <<
indent << nInitialPoints <<
" points placed" <<
nl 960 <<
indent << globalTrialPoints_ <<
" locations queried" <<
nl 962 << scalar(nInitialPoints)/scalar(
max(globalTrialPoints_, 1))
963 <<
" success rate" <<
nl 966 <<
" levels of recursion (maximum)" List< ReturnType > get(const UPtrList< T > &list, const AccessOp &aop)
List of values generated by applying the access operation to each list item.
virtual const fileName & name() const
Get the name of the output serial stream. (eg, the name of the Fstream file name) ...
void size(const label n)
Older name for setAddressableSize.
Ostream & indent(Ostream &os)
Indent stream.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
constexpr char nl
The newline '\n' character (0x0a)
Ostream & endl(Ostream &os)
Add newline and flush stream.
PointIndexHit< point > pointIndexHit
A PointIndexHit with a 3D point.
static bool & parRun() noexcept
Test if this a parallel run.
label k
Boltzmann constant.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
Macros for easy insertion into run-time selection tables.
#define forAll(list, i)
Loop across all elements in list.
vectorField pointField
pointField is a vectorField.
Type sample01()
Return a sample whose components lie in the range [0,1].
const dimensionedScalar e
Elementary charge.
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
autoDensity(const dictionary &initialPointsDict, const Time &runTime, Random &rndGen, const conformationSurfaces &geometryToConformTo, const cellShapeControl &cellShapeControls, const autoPtr< backgroundMeshDecomposition > &decomposition)
Construct from components.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual List< Vb::Point > initialPoints() const
Return the initial points for the conformalVoronoiMesh.
const backgroundMeshDecomposition & decomposition() const
static const edgeList edges
Edge to point addressing, using octant corner points.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
int debug
Static debugging option.
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
defineTypeNameAndDebug(combustionModel, 0)
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
const conformationSurfaces & geometryToConformTo() const
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
dimensionedScalar pow3(const dimensionedScalar &ds)
vector point
Point is a vector.
#define WarningInFunction
Report a warning using Foam::Warning.
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
Reduce inplace (cf. MPI Allreduce) using specified communication schedule.
messageStream Info
Information stream (stdout output on master, null elsewhere)
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))
void shuffle(UList< T > &list)
Randomise the list order.
const treeBoundBox & procBounds() const
Return the boundBox of this processor.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
treeBoundBox extend(Random &rndGen, const scalar s) const
Return slightly wider bounding box.
static constexpr const zero Zero
Global zero (0)