zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: snappyHexMesh: add -dry-run argument. See #972.

Browse Source
This commit is contained in:
mattijs
2018-11-26 12:23:01 +00:00
parent f8985566e9
commit 0e5325bff9
48 changed files with 2253 additions and 328 deletions
Download Patch File Download Diff File Expand all files Collapse all files

403
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
View File

@ -62,7 +62,7 @@ Description
#include "fvMeshTools.H"
#include "profiling.H"
#include "processorMeshes.H"
#include "vtkSetWriter.H"
#include "snappyVoxelMeshDriver.H"
using namespace Foam;
@ -325,7 +325,8 @@ autoPtr<refinementSurfaces> createRefinementSurfaces
maxLevel,
gapLevel,
scalarField(nRegions, -GREAT), //perpendicularAngle,
patchInfo
patchInfo,
false //dryRun
)
);
@ -537,8 +538,61 @@ void extractSurface
}
label checkAlignment(const polyMesh& mesh, const scalar tol, Ostream& os)
{
// Check all edges aligned with one of the coordinate axes
const faceList& faces = mesh.faces();
const pointField& points = mesh.points();
label nUnaligned = 0;
forAll(faces, facei)
{
const face& f = faces[facei];
forAll(f, fp)
{
label fp1 = f.fcIndex(fp);
const linePointRef e(edge(f[fp], f[fp1]).line(points));
const vector v(e.vec());
const scalar magV(mag(v));
if (magV > ROOTVSMALL)
{
for
(
direction dir = 0;
dir < pTraits<vector>::nComponents;
++dir
)
{
const scalar s(mag(v[dir]));
if (s > magV*tol && s < magV*(1-tol))
{
++nUnaligned;
break;
}
}
}
}
}
reduce(nUnaligned, sumOp<label>());
if (nUnaligned)
{
os << "Initial mesh has " << nUnaligned
<< " edges unaligned with any of the coordinate axes" << nl << endl;
}
return nUnaligned;
}
// Check writing tolerance before doing any serious work
scalar getMergeDistance(const polyMesh& mesh, const scalar mergeTol)
scalar getMergeDistance
(
const polyMesh& mesh,
const scalar mergeTol,
const bool dryRun
)
{
const boundBox& meshBb = mesh.bounds();
scalar mergeDist = mergeTol * meshBb.mag();
@ -550,7 +604,7 @@ scalar getMergeDistance(const polyMesh& mesh, const scalar mergeTol)
<< endl;
// check writing tolerance
if (mesh.time().writeFormat() == IOstream::ASCII)
if (mesh.time().writeFormat() == IOstream::ASCII && !dryRun)
{
const scalar writeTol = std::pow
(
@ -690,6 +744,11 @@ int main(int argc, char *argv[])
"checkGeometry",
"Check all surface geometry for quality"
);
argList::addBoolOption
(
"dry-run",
"Check case set-up only using a single time step"
);
argList::addOption
(
"surfaceSimplify",
@ -718,39 +777,17 @@ int main(int argc, char *argv[])
const bool overwrite = args.found("overwrite");
const bool checkGeometry = args.found("checkGeometry");
const bool surfaceSimplify = args.found("surfaceSimplify");
const bool dryRun = args.optionFound("dry-run");
autoPtr<fvMesh> meshPtr;
if (dryRun)
{
word regionName = fvMesh::defaultRegion;
if (args.readIfPresent("region", regionName))
{
Info<< "Create mesh " << regionName << " for time = "
<< runTime.timeName() << nl << endl;
}
else
{
Info<< "Create mesh for time = "
<< runTime.timeName() << nl << endl;
}
meshPtr.reset
(
new fvMesh
(
IOobject
(
regionName,
runTime.timeName(),
runTime,
IOobject::MUST_READ
)
)
);
Info<< "Operating in dry-run mode to detect set-up errors"
<< nl << endl;
}
fvMesh& mesh = meshPtr();
#include "createNamedMesh.H"
Info<< "Read mesh in = "
<< runTime.cpuTimeIncrement() << " s" << endl;
@ -758,6 +795,13 @@ int main(int argc, char *argv[])
mesh.boundaryMesh().checkParallelSync(true);
meshRefinement::checkCoupledFaceZones(mesh);
if (dryRun)
{
// Check if mesh aligned with cartesian axes
checkAlignment(mesh, 1e-6, Pout); //FatalIOError);
}
// Read meshing dictionary
const word dictName("snappyHexMeshDict");
@ -766,25 +810,36 @@ int main(int argc, char *argv[])
// all surface geometry
const dictionary& geometryDict = meshDict.subDict("geometry");
const dictionary& geometryDict =
meshRefinement::subDict(meshDict, "geometry", dryRun);
// refinement parameters
const dictionary& refineDict = meshDict.subDict("castellatedMeshControls");
const dictionary& refineDict =
meshRefinement::subDict(meshDict, "castellatedMeshControls", dryRun);
// mesh motion and mesh quality parameters
const dictionary& motionDict = meshDict.subDict("meshQualityControls");
const dictionary& motionDict =
meshRefinement::subDict(meshDict, "meshQualityControls", dryRun);
// snap-to-surface parameters
const dictionary& snapDict = meshDict.subDict("snapControls");
const dictionary& snapDict =
meshRefinement::subDict(meshDict, "snapControls", dryRun);
// layer addition parameters
const dictionary& layerDict = meshDict.subDict("addLayersControls");
const dictionary& layerDict =
meshRefinement::subDict(meshDict, "addLayersControls", dryRun);
// absolute merge distance
const scalar mergeDist = getMergeDistance
(
mesh,
meshDict.get<scalar>("mergeTolerance")
meshRefinement::get<scalar>
(
meshDict,
"mergeTolerance",
dryRun
),
dryRun
);
const bool keepPatches(meshDict.lookupOrDefault("keepPatches", false));
@ -975,7 +1030,7 @@ int main(int argc, char *argv[])
const scalar defaultCellSize =
motionDict.get<scalar>("defaultCellSize");
const scalar initialCellSize = ::pow(meshPtr().V()[0], 1.0/3.0);
const scalar initialCellSize = ::pow(mesh.V()[0], 1.0/3.0);
//Info<< "Wanted cell size = " << defaultCellSize << endl;
//Info<< "Current cell size = " << initialCellSize << endl;
@ -1001,8 +1056,14 @@ int main(int argc, char *argv[])
new refinementSurfaces
(
allGeometry,
refineDict.subDict("refinementSurfaces"),
refineDict.lookupOrDefault("gapLevelIncrement", 0)
meshRefinement::subDict
(
refineDict,
"refinementSurfaces",
dryRun
),
refineDict.lookupOrDefault("gapLevelIncrement", 0),
dryRun
)
);
@ -1017,6 +1078,8 @@ int main(int argc, char *argv[])
if (checkGeometry)
{
// Check geometry amongst itself (e.g. intersection, size differences)
// Extract patchInfo
List<wordList> patchTypes(allGeometry.size());
@ -1035,7 +1098,14 @@ int main(int argc, char *argv[])
if (patchInfo.set(globalRegioni))
{
patchTypes[geomi][regioni] =
patchInfo[globalRegioni].get<word>("type");
meshRefinement::get<word>
(
patchInfo[globalRegioni],
"type",
dryRun,
keyType::REGEX,
word::null
);
}
else
{
@ -1058,10 +1128,49 @@ int main(int argc, char *argv[])
true
);
return 0;
if (!dryRun)
{
return 0;
}
}
if (dryRun)
{
// Check geometry to mesh bounding box
Info<< "Checking for geometry size relative to mesh." << endl;
const boundBox& meshBb = mesh.bounds();
forAll(allGeometry, geomi)
{
const searchableSurface& s = allGeometry[geomi];
const boundBox& bb = s.bounds();
scalar ratio = bb.mag() / meshBb.mag();
if (ratio > 100.0 || ratio < 1.0/100.0)
{
Warning
<< " " << allGeometry.names()[geomi]
<< " bounds differ from mesh"
<< " by more than a factor 100:" << nl
<< " bounding box : " << bb << nl
<< " mesh bounding box : " << meshBb
<< endl;
}
if (!meshBb.contains(bb))
{
Warning
<< " " << allGeometry.names()[geomi]
<< " bounds not fully contained in mesh" << nl
<< " bounding box : " << bb << nl
<< " mesh bounding box : " << meshBb
<< endl;
}
}
Info<< endl;
}
// Read refinement shells
// ~~~~~~~~~~~~~~~~~~~~~~
@ -1070,7 +1179,8 @@ int main(int argc, char *argv[])
shellSurfaces shells
(
allGeometry,
refineDict.subDict("refinementRegions")
meshRefinement::subDict(refineDict, "refinementRegions", dryRun),
dryRun
);
Info<< "Read refinement shells in = "
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl;
@ -1093,7 +1203,7 @@ int main(int argc, char *argv[])
Info<< "Reading limit shells." << endl;
}
shellSurfaces limitShells(allGeometry, limitDict);
shellSurfaces limitShells(allGeometry, limitDict, dryRun);
if (!limitDict.empty())
{
@ -1101,6 +1211,29 @@ int main(int argc, char *argv[])
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl;
}
if (dryRun)
{
// Check for use of all geometry
const wordList& allGeomNames = allGeometry.names();
labelHashSet unusedGeometries(identity(allGeomNames.size()));
unusedGeometries.erase(surfaces.surfaces());
unusedGeometries.erase(shells.shells());
unusedGeometries.erase(limitShells.shells());
if (unusedGeometries.size())
{
IOWarningInFunction(geometryDict)
<< "The following geometry entries are not used:" << nl;
for (const label geomi : unusedGeometries)
{
Info<< " " << allGeomNames[geomi] << nl;
}
Info<< endl;
}
}
// Read feature meshes
@ -1110,12 +1243,33 @@ int main(int argc, char *argv[])
refinementFeatures features
(
mesh,
refineDict.lookup("features")
meshRefinement::lookup(refineDict, "features", dryRun),
dryRun
);
Info<< "Read features in = "
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl;
if (dryRun)
{
// Check geometry to mesh bounding box
Info<< "Checking for line geometry size relative to surface geometry."
<< endl;
OStringStream os;
bool hasErrors = features.checkSizes
(
100.0, //const scalar maxRatio,
mesh.bounds(),
true, //const bool report,
os //FatalIOError
);
if (hasErrors)
{
Warning<< os.str() << endl;
}
}
// Refinement engine
// ~~~~~~~~~~~~~~~~~
@ -1134,10 +1288,17 @@ int main(int argc, char *argv[])
surfaces, // for surface intersection refinement
features, // for feature edges/point based refinement
shells, // for volume (inside/outside) refinement
limitShells // limit of volume refinement
limitShells, // limit of volume refinement
labelList(), // initial faces to test
dryRun
);
Info<< "Calculated surface intersections in = "
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl;
if (!dryRun)
{
meshRefiner.updateIntersections(identity(mesh.nFaces()));
Info<< "Calculated surface intersections in = "
<< mesh.time().cpuTimeIncrement() << " s" << nl << endl;
}
// Some stats
meshRefiner.printMeshInfo(debugLevel, "Initial mesh");
@ -1151,10 +1312,10 @@ int main(int argc, char *argv[])
// Refinement parameters
const refinementParameters refineParams(refineDict);
const refinementParameters refineParams(refineDict, dryRun);
// Snap parameters
const snapParameters snapParams(snapDict);
const snapParameters snapParams(snapDict, dryRun);
@ -1464,9 +1625,46 @@ int main(int argc, char *argv[])
// Now do the real work -refinement -snapping -layers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const bool wantRefine(meshDict.get<bool>("castellatedMesh"));
const bool wantSnap(meshDict.get<bool>("snap"));
const bool wantLayers(meshDict.get<bool>("addLayers"));
const bool wantRefine
(
meshRefinement::get<bool>(meshDict, "castellatedMesh", dryRun)
);
const bool wantSnap
(
meshRefinement::get<bool>(meshDict, "snap", dryRun)
);
const bool wantLayers
(
meshRefinement::get<bool>(meshDict, "addLayers", dryRun)
);
if (dryRun)
{
string errorMsg(FatalError.message());
string IOerrorMsg(FatalIOError.message());
if (errorMsg.size() || IOerrorMsg.size())
{
//errorMsg = "[dryRun] " + errorMsg;
//errorMsg.replaceAll("\n", "\n[dryRun] ");
//IOerrorMsg = "[dryRun] " + IOerrorMsg;
//IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
Warning
<< nl
<< "Missing/incorrect required dictionary entries:" << nl
<< nl
<< IOerrorMsg.c_str() << nl
<< errorMsg.c_str() << nl << nl
<< "Exiting dry-run" << nl << endl;
FatalError.clear();
FatalIOError.clear();
return 0;
}
}
const bool mergePatchFaces
(
@ -1492,7 +1690,8 @@ int main(int argc, char *argv[])
distributor,
globalToMasterPatch,
globalToSlavePatch,
setFormatter
setFormatter,
dryRun
);
@ -1518,13 +1717,16 @@ int main(int argc, char *argv[])
fvMeshTools::removeEmptyPatches(mesh, true);
}
writeMesh
(
"Refined mesh",
meshRefiner,
debugLevel,
meshRefinement::writeLevel()
);
if (!dryRun)
{
writeMesh
(
"Refined mesh",
meshRefiner,
debugLevel,
meshRefinement::writeLevel()
);
}
Info<< "Mesh refined in = "
<< timer.cpuTimeIncrement() << " s." << endl;
@ -1540,7 +1742,8 @@ int main(int argc, char *argv[])
(
meshRefiner,
globalToMasterPatch,
globalToSlavePatch
globalToSlavePatch,
dryRun
);
if (!overwrite && !debugLevel)
@ -1568,13 +1771,16 @@ int main(int argc, char *argv[])
fvMeshTools::removeEmptyPatches(mesh, true);
}
writeMesh
(
"Snapped mesh",
meshRefiner,
debugLevel,
meshRefinement::writeLevel()
);
if (!dryRun)
{
writeMesh
(
"Snapped mesh",
meshRefiner,
debugLevel,
meshRefinement::writeLevel()
);
}
Info<< "Mesh snapped in = "
<< timer.cpuTimeIncrement() << " s." << endl;
@ -1587,13 +1793,19 @@ int main(int argc, char *argv[])
cpuTime timer;
// Layer addition parameters
const layerParameters layerParams(layerDict, mesh.boundaryMesh());
const layerParameters layerParams
(
layerDict,
mesh.boundaryMesh(),
dryRun
);
snappyLayerDriver layerDriver
(
meshRefiner,
globalToMasterPatch,
globalToSlavePatch
globalToSlavePatch,
dryRun
);
// Use the maxLocalCells from the refinement parameters
@ -1626,13 +1838,16 @@ int main(int argc, char *argv[])
fvMeshTools::removeEmptyPatches(mesh, true);
}
writeMesh
(
"Layer mesh",
meshRefiner,
debugLevel,
meshRefinement::writeLevel()
);
if (!dryRun)
{
writeMesh
(
"Layer mesh",
meshRefiner,
debugLevel,
meshRefinement::writeLevel()
);
}
Info<< "Layers added in = "
<< timer.cpuTimeIncrement() << " s." << endl;
@ -1647,7 +1862,7 @@ int main(int argc, char *argv[])
// Check final mesh
Info<< "Checking final mesh ..." << endl;
faceSet wrongFaces(mesh, "wrongFaces", mesh.nFaces()/100);
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces, dryRun);
const label nErrors = returnReduce
(
wrongFaces.size(),
@ -1736,6 +1951,34 @@ int main(int argc, char *argv[])
Info<< "Finished meshing in = "
<< runTime.elapsedCpuTime() << " s." << endl;
if (dryRun)
{
string errorMsg(FatalError.message());
string IOerrorMsg(FatalIOError.message());
if (errorMsg.size() || IOerrorMsg.size())
{
//errorMsg = "[dryRun] " + errorMsg;
//errorMsg.replaceAll("\n", "\n[dryRun] ");
//IOerrorMsg = "[dryRun] " + IOerrorMsg;
//IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
Perr<< nl
<< "Missing/incorrect required dictionary entries:" << nl
<< nl
<< IOerrorMsg.c_str() << nl
<< errorMsg.c_str() << nl << nl
<< "Exiting dry-run" << nl << endl;
FatalError.clear();
FatalIOError.clear();
return 0;
}
}
Info<< "End\n" << endl;
return 0;

6
src/OpenFOAM/db/error/error.C
View File

@ -203,6 +203,12 @@ Foam::string Foam::error::message() const
}
void Foam::error::clear() const
{
return messageStreamPtr_->reset();
}
void Foam::error::exit(const int errNo)
{
if (!throwExceptions_ && JobInfo::constructed)

3
src/OpenFOAM/db/error/error.H
View File

@ -114,6 +114,9 @@ public:
string message() const;
//- Clear any messages
void clear() const;
inline const string& functionName() const
{
return functionName_;

12
src/dynamicMesh/motionSmoother/motionSmoother.C
View File

@ -41,7 +41,8 @@ Foam::motionSmoother::motionSmoother
pointMesh& pMesh,
indirectPrimitivePatch& pp,
const labelList& adaptPatchIDs,
const dictionary& paramDict
const dictionary& paramDict,
const bool dryRun
)
:
motionSmootherData(pMesh),
@ -54,7 +55,8 @@ Foam::motionSmoother::motionSmoother
motionSmootherData::scale_,
motionSmootherData::oldPoints_,
adaptPatchIDs,
paramDict
paramDict,
dryRun
)
{}
@ -65,7 +67,8 @@ Foam::motionSmoother::motionSmoother
indirectPrimitivePatch& pp,
const labelList& adaptPatchIDs,
const pointVectorField& displacement,
const dictionary& paramDict
const dictionary& paramDict,
const bool dryRun
)
:
motionSmootherData(displacement),
@ -78,7 +81,8 @@ Foam::motionSmoother::motionSmoother
motionSmootherData::scale_,
motionSmootherData::oldPoints_,
adaptPatchIDs,
paramDict
paramDict,
dryRun
)
{}

6
src/dynamicMesh/motionSmoother/motionSmoother.H
View File

@ -107,7 +107,8 @@ public:
pointMesh& pMesh,
indirectPrimitivePatch& pp, //!< 'outside' points
const labelList& adaptPatchIDs, //!< patches forming 'outside'
const dictionary& paramDict
const dictionary& paramDict,
const bool dryRun = false
);
//- Construct from mesh, patches to work on and smoothing parameters and
@ -118,7 +119,8 @@ public:
indirectPrimitivePatch& pp, //!< 'outside' points
const labelList& adaptPatchIDs, //!< patches forming 'outside'
const pointVectorField& displacement,
const dictionary& paramDict
const dictionary& paramDict,
const bool dryRun = false
);
};

26
src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
View File

@ -328,7 +328,8 @@ Foam::motionSmootherAlgo::motionSmootherAlgo
pointScalarField& scale,
pointField& oldPoints,
const labelList& adaptPatchIDs,
const dictionary& paramDict
const dictionary& paramDict,
const bool dryRun
)
:
mesh_(mesh),
@ -339,6 +340,7 @@ Foam::motionSmootherAlgo::motionSmootherAlgo
oldPoints_(oldPoints),
adaptPatchIDs_(adaptPatchIDs),
paramDict_(paramDict),
dryRun_(dryRun),
isInternalPoint_(mesh_.nPoints(), true)
{
updateMesh();
@ -856,9 +858,14 @@ bool Foam::motionSmootherAlgo::scaleMesh
}
}
const scalar errorReduction = paramDict.get<scalar>("errorReduction");
const label nSmoothScale = paramDict.get<label>("nSmoothScale");
const scalar errorReduction = get<scalar>
(
paramDict, "errorReduction", dryRun_, keyType::REGEX_RECURSIVE
);
const label nSmoothScale = get<label>
(
paramDict, "nSmoothScale", dryRun_, keyType::REGEX_RECURSIVE
);
// Note: displacement_ should already be synced already from setDisplacement
// but just to make sure.
@ -885,7 +892,16 @@ bool Foam::motionSmootherAlgo::scaleMesh
// Check. Returns parallel number of incorrect faces.
faceSet wrongFaces(mesh_, "wrongFaces", mesh_.nFaces()/100+100);
checkMesh(false, mesh_, meshQualityDict, checkFaces, baffles, wrongFaces);
checkMesh
(
false,
mesh_,
meshQualityDict,
checkFaces,
baffles,
wrongFaces,
dryRun_
);
if (returnReduce(wrongFaces.size(), sumOp<label>()) <= nAllowableErrors)
{

45
src/dynamicMesh/motionSmoother/motionSmootherAlgo.H
View File

@ -156,9 +156,12 @@ class motionSmootherAlgo
// displacement on.
const labelList adaptPatchIDs_;
// Smoothing and checking parameters
//- Smoothing and checking parameters
dictionary paramDict_;
//- In test/dry-run mode?
const bool dryRun_;
//- Is mesh point on boundary or not
bitSet isInternalPoint_;
@ -303,7 +306,8 @@ public:
pointScalarField& scale,
pointField& oldPoints,
const labelList& adaptPatchIDs, // patches forming 'outside'
const dictionary& paramDict
const dictionary& paramDict,
const bool dryRun = false
);
@ -442,7 +446,8 @@ public:
const bool report,
const polyMesh& mesh,
const dictionary& dict,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun = false
);
//- Check (subset of mesh) with mesh settings in dict.
@ -454,7 +459,8 @@ public:
const polyMesh& mesh,
const dictionary& dict,
const labelList& checkFaces,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun = false
);
//- Check (subset of mesh including baffles) with mesh settings
@ -467,7 +473,8 @@ public:
const dictionary& dict,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun = false
);
//- Check part of mesh with mesh settings in dict.
@ -480,7 +487,8 @@ public:
const polyMeshGeometry&,
const pointField&,
const labelList& checkFaces,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun = false
);
//- Check part of mesh including baffles with mesh settings in dict.
@ -494,7 +502,8 @@ public:
const pointField&,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun = false
);
// Helper functions to manipulate displacement vector.
@ -508,6 +517,28 @@ public:
const scalarField& edgeWeight,
GeometricField<Type, pointPatchField, pointMesh>& newFld
) const;
//- Wrapper around dictionary::get which does not exit
template<class Type>
static Type get
(
const dictionary& dict,
const word& keyword,
const bool noExit,
enum keyType::option matchOpt,
const Type& defaultValue = Zero
);
//
// //- Wrapper around dictionary::get which does not exit
// template<class Type>
// static Type get
// (
// const dictionary& dict,
// const word& keyword,
// const bool noExit,
// bool recursive,
// bool patternMatch
// );
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

145
src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
View File

@ -35,7 +35,8 @@ bool Foam::motionSmootherAlgo::checkMesh
const polyMesh& mesh,
const dictionary& dict,
const labelList& checkFaces,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun
)
{
List<labelPair> emptyBaffles;
@ -46,7 +47,8 @@ bool Foam::motionSmootherAlgo::checkMesh
dict,
checkFaces,
emptyBaffles,
wrongFaces
wrongFaces,
dryRun
);
}
@ -57,54 +59,60 @@ bool Foam::motionSmootherAlgo::checkMesh
const dictionary& dict,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun
)
{
const scalar maxNonOrtho
(
dict.get<scalar>("maxNonOrtho", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "maxNonOrtho", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minVol
(
dict.get<scalar>("minVol", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minVol", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minTetQuality
(
dict.get<scalar>("minTetQuality", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minTetQuality", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar maxConcave
(
dict.get<scalar>("maxConcave", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "maxConcave", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minArea
(
dict.get<scalar>("minArea", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minArea", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar maxIntSkew
(
dict.get<scalar>("maxInternalSkewness", keyType::REGEX_RECURSIVE)
get<scalar>
(
dict, "maxInternalSkewness", dryRun, keyType::REGEX_RECURSIVE
)
);
const scalar maxBounSkew
(
dict.get<scalar>("maxBoundarySkewness", keyType::REGEX_RECURSIVE)
get<scalar>
(
dict, "maxBoundarySkewness", dryRun, keyType::REGEX_RECURSIVE
)
);
const scalar minWeight
(
dict.get<scalar>("minFaceWeight", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minFaceWeight", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minVolRatio
(
dict.get<scalar>("minVolRatio", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minVolRatio", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minTwist
(
dict.get<scalar>("minTwist", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minTwist", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minTriangleTwist
(
dict.get<scalar>("minTriangleTwist", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minTriangleTwist", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minFaceFlatness
(
dict.lookupOrDefault<scalar>
@ -114,8 +122,37 @@ bool Foam::motionSmootherAlgo::checkMesh
);
const scalar minDet
(
dict.get<scalar>("minDeterminant", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minDeterminant", dryRun, keyType::REGEX_RECURSIVE)
);
if (dryRun)
{
string errorMsg(FatalError.message());
string IOerrorMsg(FatalIOError.message());
if (errorMsg.size() || IOerrorMsg.size())
{
//errorMsg = "[dryRun] " + errorMsg;
//errorMsg.replaceAll("\n", "\n[dryRun] ");
//IOerrorMsg = "[dryRun] " + IOerrorMsg;
//IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
IOWarningInFunction(dict)
<< nl
<< "Missing/incorrect required dictionary entries:" << nl
<< nl
<< IOerrorMsg.c_str() << nl
<< errorMsg.c_str() << nl << nl
<< "Exiting dry-run" << nl << endl;
FatalError.clear();
FatalIOError.clear();
}
return false;
}
label nWrongFaces = 0;
Info<< "Checking faces in error :" << endl;
@ -422,7 +459,8 @@ bool Foam::motionSmootherAlgo::checkMesh
const bool report,
const polyMesh& mesh,
const dictionary& dict,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun
)
{
return checkMesh
@ -431,7 +469,8 @@ bool Foam::motionSmootherAlgo::checkMesh
mesh,
dict,
identity(mesh.nFaces()),
wrongFaces
wrongFaces,
dryRun
);
}
@ -442,7 +481,8 @@ bool Foam::motionSmootherAlgo::checkMesh
const polyMeshGeometry& meshGeom,
const pointField& points,
const labelList& checkFaces,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun
)
{
List<labelPair> emptyBaffles;
@ -455,7 +495,8 @@ bool Foam::motionSmootherAlgo::checkMesh
points,
checkFaces,
emptyBaffles,
wrongFaces
wrongFaces,
dryRun
);
}
@ -468,64 +509,92 @@ bool Foam::motionSmootherAlgo::checkMesh
const pointField& points,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet& wrongFaces
labelHashSet& wrongFaces,
const bool dryRun
)
{
const scalar maxNonOrtho
(
dict.get<scalar>("maxNonOrtho", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "maxNonOrtho", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minVol
(
dict.get<scalar>("minVol", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minVol", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minTetQuality
(
dict.get<scalar>("minTetQuality", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minTetQuality", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar maxConcave
(
dict.get<scalar>("maxConcave", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "maxConcave", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minArea
(
dict.get<scalar>("minArea", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minArea", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar maxIntSkew
(
dict.get<scalar>("maxInternalSkewness", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "maxInternalSkewness", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar maxBounSkew
(
dict.get<scalar>("maxBoundarySkewness", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "maxBoundarySkewness", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minWeight
(
dict.get<scalar>("minFaceWeight", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minFaceWeight", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minVolRatio
(
dict.get<scalar>("minVolRatio", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minVolRatio", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minTwist
(
dict.get<scalar>("minTwist", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minTwist", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minTriangleTwist
(
dict.get<scalar>("minTriangleTwist", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minTriangleTwist", dryRun, keyType::REGEX_RECURSIVE)
);
const scalar minFaceFlatness
scalar minFaceFlatness = -1.0;
dict.readIfPresent
(
dict.lookupOrDefault<scalar>
(
"minFaceFlatness", -1, keyType::REGEX_RECURSIVE
)
"minFaceFlatness",
minFaceFlatness,
keyType::REGEX_RECURSIVE
);
const scalar minDet
(
dict.get<scalar>("minDeterminant", keyType::REGEX_RECURSIVE)
get<scalar>(dict, "minDeterminant", dryRun, keyType::REGEX_RECURSIVE)
);
if (dryRun)
{
string errorMsg(FatalError.message());
string IOerrorMsg(FatalIOError.message());
if (errorMsg.size() || IOerrorMsg.size())
{
errorMsg = "[dryRun] " + errorMsg;
errorMsg.replaceAll("\n", "\n[dryRun] ");
IOerrorMsg = "[dryRun] " + IOerrorMsg;
IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
Perr<< nl
<< "Missing/incorrect required dictionary entries:" << nl
<< nl
<< IOerrorMsg.c_str() << nl
<< errorMsg.c_str() << nl << nl
<< "Exiting dry-run" << nl << endl;
FatalError.clear();
FatalIOError.clear();
}
return false;
}
label nWrongFaces = 0;
Info<< "Checking faces in error :" << endl;

31
src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C
View File

@ -287,4 +287,35 @@ void Foam::motionSmootherAlgo::testSyncField
}
template<class Type>
Type Foam::motionSmootherAlgo::get
(
const dictionary& dict,
const word& keyword,
const bool noExit,
enum keyType::option matchOpt,
const Type& defaultValue
)
{
Type val(defaultValue);
if
(
!dict.readEntry
(
keyword,
val,
matchOpt,
!noExit
)
)
{
FatalIOError
<< "Entry '" << keyword << "' not found in dictionary "
<< dict.name() << endl;
}
return val;
}
// ************************************************************************* //

1
src/mesh/snappyHexMesh/Make/files
View File

@ -2,6 +2,7 @@ snappyHexMeshDriver/snappyLayerDriver.C
snappyHexMeshDriver/snappySnapDriver.C
snappyHexMeshDriver/snappySnapDriverFeature.C
snappyHexMeshDriver/snappyRefineDriver.C
snappyHexMeshDriver/snappyVoxelMeshDriver.C
snappyHexMeshDriver/layerParameters/layerParameters.C
snappyHexMeshDriver/refinementParameters/refinementParameters.C

3
src/mesh/snappyHexMesh/Make/options
View File

@ -1,10 +1,12 @@
EXE_INC = \
/* -g -DFULLDEBUG -O0 */ \
-I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/overset/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
@ -14,4 +16,5 @@ LIB_LIBS = \
-llagrangian \
-lmeshTools \
-lfvMotionSolvers \
-loverset \
-ldistributed

17
src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C
View File

@ -109,10 +109,11 @@ Foam::displacementMotionSolverMeshMover::displacementMotionSolverMeshMover
(
const dictionary& dict,
const List<labelPair>& baffles,
pointVectorField& pointDisplacement
pointVectorField& pointDisplacement,
const bool dryRun
)
:
externalDisplacementMeshMover(dict, baffles, pointDisplacement),
externalDisplacementMeshMover(dict, baffles, pointDisplacement, dryRun),
solverPtr_
(
@ -178,7 +179,8 @@ Foam::displacementMotionSolverMeshMover::displacementMotionSolverMeshMover
scale_,
oldPoints_,
adaptPatchIDs_,
dict
dict,
dryRun
),
fieldSmoother_(mesh())
@ -205,12 +207,15 @@ bool Foam::displacementMotionSolverMeshMover::move
// Note that this has to update the pointDisplacement boundary conditions
// as well, not just the internal field.
{
const label nSmoothPatchThickness
const label nSmoothPatchThickness = meshRefinement::get<label>
(
moveDict.get<label>("nSmoothThickness")
moveDict, "nSmoothThickness", dryRun_, keyType::REGEX
);
const word minThicknessName(moveDict.get<word>("minThicknessName"));
const word minThicknessName = meshRefinement::get<word>
(
moveDict, "minThicknessName", dryRun_, keyType::REGEX, word::null
);
scalarField zeroMinThickness;

3
src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H
View File

@ -112,7 +112,8 @@ public:
(
const dictionary& dict,
const List<labelPair>& baffles,
pointVectorField& pointDisplacemen
pointVectorField& pointDisplacement,
const bool dryRun
);

11
src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C
View File

@ -117,11 +117,13 @@ Foam::externalDisplacementMeshMover::externalDisplacementMeshMover
(
const dictionary& dict,
const List<labelPair>& baffles,
pointVectorField& pointDisplacement
pointVectorField& pointDisplacement,
const bool dryRun
)
:
baffles_(baffles),
pointDisplacement_(pointDisplacement)
pointDisplacement_(pointDisplacement),
dryRun_(dryRun)
{}
@ -133,7 +135,8 @@ Foam::externalDisplacementMeshMover::New
const word& type,
const dictionary& dict,
const List<labelPair>& baffles,
pointVectorField& pointDisplacement
pointVectorField& pointDisplacement,
const bool dryRun
)
{
Info<< "Selecting externalDisplacementMeshMover " << type << endl;
@ -152,7 +155,7 @@ Foam::externalDisplacementMeshMover::New
return autoPtr<externalDisplacementMeshMover>
(
cstrIter()(dict, baffles, pointDisplacement)
cstrIter()(dict, baffles, pointDisplacement, dryRun)
);
}

14
src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H
View File

@ -66,6 +66,9 @@ protected:
//- Reference to point motion field
pointVectorField& pointDisplacement_;
//- In dry-run mode?
const bool dryRun_;
// Protected Member functions
@ -110,9 +113,10 @@ public:
(
const dictionary& dict,
const List<labelPair>& baffles,
pointVectorField& pointDisplacement
pointVectorField& pointDisplacement,
const bool dryRun
),
(dict, baffles, pointDisplacement)
(dict, baffles, pointDisplacement, dryRun)
);
@ -124,7 +128,8 @@ public:
(
const dictionary& dict,
const List<labelPair>& baffles,
pointVectorField& pointDisplacement
pointVectorField& pointDisplacement,
const bool dryRun
);
@ -136,7 +141,8 @@ public:
const word& type,
const dictionary& dict,
const List<labelPair>& baffles,
pointVectorField& pointDisplacement
pointVectorField& pointDisplacement,
const bool dryRun = false
);

79
src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C
View File

@ -138,38 +138,55 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
//- Smooth surface normals
const label nSmoothSurfaceNormals
(
coeffDict.get<label>("nSmoothSurfaceNormals")
meshRefinement::get<label>
(
coeffDict,
"nSmoothSurfaceNormals",
dryRun_
)
);
// Note: parameter name changed
// "minMedianAxisAngle" -> "minMedialAxisAngle" (DEC-2013)
// but not previously reported.
scalar minMedialAxisAngleCos = Foam::cos
scalar minMedialAxisAngle(Zero);
if
(
degToRad
!coeffDict.readCompat
(
coeffDict.getCompat<scalar>
(
"minMedialAxisAngle", {{ "minMedianAxisAngle", 1712 }}
)
"minMedialAxisAngle",
{{ "minMedianAxisAngle", 1712 }},
minMedialAxisAngle,
keyType::REGEX,
!dryRun_
)
);
)
{
FatalIOError
<< "Entry '" << "minMedialAxisAngle"
<< "' not found in dictionary " << coeffDict.name() << endl;
}
const scalar minMedialAxisAngleCos(Foam::cos(degToRad(minMedialAxisAngle)));
//- Feature angle when to stop adding layers
const scalar featureAngle = coeffDict.get<scalar>("featureAngle");
const scalar featureAngle
(
meshRefinement::get<scalar>(coeffDict, "featureAngle", dryRun_)
);
//- When to slip along wall
const scalar slipFeatureAngle =
coeffDict.lookupOrDefault<scalar>
(
"slipFeatureAngle",
0.5*featureAngle
);
(
coeffDict.found("slipFeatureAngle")
? coeffDict.get<scalar>("slipFeatureAngle")
: 0.5*featureAngle
);
//- Smooth internal normals
const label nSmoothNormals
(
coeffDict.get<label>("nSmoothNormals")
meshRefinement::get<label>(coeffDict, "nSmoothNormals", dryRun_)
);
//- Number of edges walking out
@ -1222,10 +1239,11 @@ Foam::medialAxisMeshMover::medialAxisMeshMover
(
const dictionary& dict,
const List<labelPair>& baffles,
pointVectorField& pointDisplacement
pointVectorField& pointDisplacement,
const bool dryRun
)
:
externalDisplacementMeshMover(dict, baffles, pointDisplacement),
externalDisplacementMeshMover(dict, baffles, pointDisplacement, dryRun),
adaptPatchIDs_(getFixedValueBCs(pointDisplacement)),
adaptPatchPtr_(getPatch(mesh(), adaptPatchIDs_)),
scale_
@ -1251,7 +1269,8 @@ Foam::medialAxisMeshMover::medialAxisMeshMover
scale_,
oldPoints_,
adaptPatchIDs_,
dict
dict,
dryRun
),
fieldSmoother_(mesh()),
dispVec_
@ -1341,23 +1360,15 @@ void Foam::medialAxisMeshMover::calculateDisplacement
// ~~~~~~~~~~~~~
//- (lambda-mu) smoothing of internal displacement
const label nSmoothDisplacement = coeffDict.lookupOrDefault
(
"nSmoothDisplacement",
0
);
const label nSmoothDisplacement =
coeffDict.lookupOrDefault("nSmoothDisplacement", 0);
//- Layer thickness too big
const scalar maxThicknessToMedialRatio
(
coeffDict.get<scalar>("maxThicknessToMedialRatio")
);
const scalar maxThicknessToMedialRatio =
coeffDict.get<scalar>("maxThicknessToMedialRatio");
//- Feature angle when to stop adding layers
const scalar featureAngle
(
coeffDict.get<scalar>("featureAngle")
);
const scalar featureAngle = coeffDict.get<scalar>("featureAngle");
//- Stop layer growth where mesh wraps around sharp edge
scalar layerTerminationAngle = coeffDict.lookupOrDefault<scalar>
@ -1368,10 +1379,8 @@ void Foam::medialAxisMeshMover::calculateDisplacement
scalar minCosLayerTermination = Foam::cos(degToRad(layerTerminationAngle));
//- Smoothing wanted patch thickness
const label nSmoothPatchThickness
(
coeffDict.get<label>("nSmoothThickness")
);
const label nSmoothPatchThickness =
coeffDict.get<label>("nSmoothThickness");
//- Number of edges walking out
const label nMedialAxisIter = coeffDict.lookupOrDefault<label>

3
src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H
View File

@ -207,7 +207,8 @@ public:
(
const dictionary& dict,
const List<labelPair>& baffles,
pointVectorField& pointDisplacement
pointVectorField& pointDisplacement,
const bool dryRun
);

142
src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
View File

@ -487,7 +487,11 @@ void Foam::meshRefinement::checkData()
neiBoundaryFc
);
}
// Check meshRefinement
const labelList& surfIndex = surfaceIndex();
{
// Get boundary face centre and level. Coupled aware.
labelList neiLevel(nBnd);
@ -549,14 +553,14 @@ void Foam::meshRefinement::checkData()
// Check
forAll(surfaceHit, facei)
{
if (surfaceIndex_[facei] != surfaceHit[facei])
if (surfIndex[facei] != surfaceHit[facei])
{
if (mesh_.isInternalFace(facei))
{
WarningInFunction
<< "Internal face:" << facei
<< " fc:" << mesh_.faceCentres()[facei]
<< " cached surfaceIndex_:" << surfaceIndex_[facei]
<< " cached surfaceIndex_:" << surfIndex[facei]
<< " current:" << surfaceHit[facei]
<< " ownCc:"
<< mesh_.cellCentres()[mesh_.faceOwner()[facei]]
@ -566,14 +570,14 @@ void Foam::meshRefinement::checkData()
}
else if
(
surfaceIndex_[facei]
surfIndex[facei]
!= neiHit[facei-mesh_.nInternalFaces()]
)
{
WarningInFunction
<< "Boundary face:" << facei
<< " fc:" << mesh_.faceCentres()[facei]
<< " cached surfaceIndex_:" << surfaceIndex_[facei]
<< " cached surfaceIndex_:" << surfIndex[facei]
<< " current:" << surfaceHit[facei]
<< " ownCc:"
<< mesh_.cellCentres()[mesh_.faceOwner()[facei]]
@ -996,7 +1000,8 @@ Foam::label Foam::meshRefinement::splitFacesUndo
false, // report
mesh_,
motionDict,
errorFaces
errorFaces,
dryRun_
);
if (!hasErrors)
{
@ -1208,7 +1213,9 @@ Foam::meshRefinement::meshRefinement
const refinementSurfaces& surfaces,
const refinementFeatures& features,
const shellSurfaces& shells,
const shellSurfaces& limitShells
const shellSurfaces& limitShells,
const labelUList& checkFaces,
const bool dryRun
)
:
mesh_(mesh),
@ -1219,6 +1226,7 @@ Foam::meshRefinement::meshRefinement
features_(features),
shells_(shells),
limitShells_(limitShells),
dryRun_(dryRun),
meshCutter_
(
mesh,
@ -1241,12 +1249,35 @@ Foam::meshRefinement::meshRefinement
userFaceData_(0)
{
// recalculate intersections for all faces
updateIntersections(identity(mesh_.nFaces()));
updateIntersections(checkFaces);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::labelList& Foam::meshRefinement::surfaceIndex() const
{
if (surfaceIndex_.size() != mesh_.nFaces())
{
const_cast<meshRefinement&>(*this).updateIntersections
(
identity(mesh_.nFaces())
);
}
return surfaceIndex_;
}
Foam::labelList& Foam::meshRefinement::surfaceIndex()
{
if (surfaceIndex_.size() != mesh_.nFaces())
{
updateIntersections(identity(mesh_.nFaces()));
}
return surfaceIndex_;
}
Foam::label Foam::meshRefinement::countHits() const
{
// Stats on edges to test. Count proc faces only once.
@ -1254,9 +1285,11 @@ Foam::label Foam::meshRefinement::countHits() const
label nHits = 0;
forAll(surfaceIndex_, facei)
const labelList& surfIndex = surfaceIndex();
forAll(surfIndex, facei)
{
if (surfaceIndex_[facei] >= 0 && isMasterFace.test(facei))
if (surfIndex[facei] >= 0 && isMasterFace.test(facei))
{
++nHits;
}
@ -1520,9 +1553,11 @@ Foam::labelList Foam::meshRefinement::intersectedFaces() const
{
label nBoundaryFaces = 0;
forAll(surfaceIndex_, facei)
const labelList& surfIndex = surfaceIndex();
forAll(surfIndex, facei)
{
if (surfaceIndex_[facei] != -1)
if (surfIndex[facei] != -1)
{
nBoundaryFaces++;
}
@ -1531,9 +1566,9 @@ Foam::labelList Foam::meshRefinement::intersectedFaces() const
labelList surfaceFaces(nBoundaryFaces);
nBoundaryFaces = 0;
forAll(surfaceIndex_, facei)
forAll(surfIndex, facei)
{
if (surfaceIndex_[facei] != -1)
if (surfIndex[facei] != -1)
{
surfaceFaces[nBoundaryFaces++] = facei;
}
@ -1550,9 +1585,11 @@ Foam::labelList Foam::meshRefinement::intersectedPoints() const
bitSet isBoundaryPoint(mesh_.nPoints());
label nBoundaryPoints = 0;
forAll(surfaceIndex_, facei)
const labelList& surfIndex = surfaceIndex();
forAll(surfIndex, facei)
{
if (surfaceIndex_[facei] != -1)
if (surfIndex[facei] != -1)
{
const face& f = faces[facei];
@ -3094,6 +3131,9 @@ void Foam::meshRefinement::write
if (writeFlags && !(writeFlags & NOWRITEREFINEMENT))
{
meshCutter_.write();
// Force calculation before writing
(void)surfaceIndex();
surfaceIndex_.write();
}
@ -3157,4 +3197,76 @@ void Foam::meshRefinement::writeLevel(const writeType flags)
//}
const Foam::dictionary& Foam::meshRefinement::subDict
(
const dictionary& dict,
const word& keyword,
const bool noExit
)
{
if (noExit)
{
// Find non-recursive with patterns
const dictionary::const_searcher finder
(
dict.csearch
(
keyword,
keyType::REGEX
)
);
if (!finder.found())
{
FatalIOErrorInFunction(dict)
<< "Entry '" << keyword << "' not found in dictionary "
<< dict.name();
return dictionary::null;
}
else
{
return finder.dict();
}
}
else
{
return dict.subDict(keyword);
}
}
Foam::ITstream& Foam::meshRefinement::lookup
(
const dictionary& dict,
const word& keyword,
const bool noExit
)
{
if (noExit)
{
const dictionary::const_searcher finder
(
dict.csearch(keyword, keyType::REGEX)
);
if (!finder.found())
{
FatalIOErrorInFunction(dict)
<< "Entry '" << keyword << "' not found in dictionary "
<< dict.name();
// Fake entry
return dict.first()->stream();
}
else
{
return finder.ref().stream();
}
}
else
{
return dict.lookup(keyword);
}
}
// ************************************************************************* //

44
src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
View File

@ -165,6 +165,9 @@ private:
//- All limit-refinement interaction
const shellSurfaces& limitShells_;
//- Are we operating in test mode?
const bool dryRun_;
//- Refinement engine
hexRef8 meshCutter_;
@ -819,7 +822,9 @@ public:
const refinementSurfaces&,
const refinementFeatures&,
const shellSurfaces&, // omnidirectional refinement
const shellSurfaces& // limit refinement
const shellSurfaces&, // limit refinement
const labelUList& checkFaces, // initial faces to check
const bool dryRun
);
@ -885,15 +890,9 @@ public:
}
//- Per start-end edge the index of the surface hit
const labelList& surfaceIndex() const
{
return surfaceIndex_;
}
const labelList& surfaceIndex() const;
labelList& surfaceIndex()
{
return surfaceIndex_;
}
labelList& surfaceIndex();
//- For faces originating from processor faces store the original
// patch
@ -1531,6 +1530,33 @@ public:
const EnumContainer& namedEnum,
const wordList& words
);
//- Wrapper around dictionary::get which does not exit
template<class Type>
static Type get
(
const dictionary& dict,
const word& keyword,
const bool noExit,
enum keyType::option matchOpt = keyType::REGEX,
const Type& defaultValue = Zero
);
//- Wrapper around dictionary::subDict which does not exit
static const dictionary& subDict
(
const dictionary& dict,
const word& keyword,
const bool noExit
);
//- Wrapper around dictionary::lookup which does not exit
static ITstream& lookup
(
const dictionary& dict,
const word& keyword,
const bool noExit
);
};

2
src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
View File

@ -1378,7 +1378,7 @@ Foam::label Foam::meshRefinement::markInternalGapRefinement
<< mesh_.globalData().nTotalFaces() << endl;
transportData::trackData td(surfaceIndex_);
transportData::trackData td(surfaceIndex());
FaceCellWave<transportData, transportData::trackData> deltaCalc
(

6
src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
View File

@ -414,7 +414,8 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
false, // report
mesh_,
motionDict,
errorFaces
errorFaces,
dryRun_
);
//if (checkEdgeConnectivity)
@ -898,7 +899,8 @@ Foam::label Foam::meshRefinement::mergeEdgesUndo
false, // report
mesh_,
motionDict,
errorFaces
errorFaces,
dryRun_
);
//if (checkEdgeConnectivity)
//{

17
src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
View File

@ -633,8 +633,8 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
if (checkCollapse)
{
motionDict.readEntry("minArea", minArea);
motionDict.readEntry("maxNonOrtho", maxNonOrtho);
minArea = get<scalar>(motionDict, "minArea", dryRun_);
maxNonOrtho = get<scalar>(motionDict, "maxNonOrtho", dryRun_);
Info<< "markFacesOnProblemCells :"
<< " Deleting all-anchor surface cells only if"
@ -1123,7 +1123,14 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
// Check initial mesh
Info<< "Checking initial mesh ..." << endl;
labelHashSet wrongFaces(mesh_.nFaces()/100);
motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
motionSmoother::checkMesh
(
false,
mesh_,
motionDict,
wrongFaces,
dryRun_
);
const label nInitErrors = returnReduce
(
wrongFaces.size(),
@ -1235,7 +1242,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
// nWrongFaces = nNewWrongFaces;
//}
scalar minArea(motionDict.get<scalar>("minArea"));
scalar minArea(get<scalar>(motionDict, "minArea", dryRun_));
if (minArea > -SMALL)
{
polyMeshGeometry::checkFaceArea
@ -1262,7 +1269,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
nWrongFaces = nNewWrongFaces;
}
scalar minDet(motionDict.get<scalar>("minDeterminant"));
scalar minDet(get<scalar>(motionDict, "minDeterminant", dryRun_));
if (minDet > -1)
{
polyMeshGeometry::checkCellDeterminant

6
src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
View File

@ -913,9 +913,11 @@ Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
label nTest = 0;
forAll(surfaceIndex_, faceI)
const labelList& surfIndex = surfaceIndex();
forAll(surfIndex, faceI)
{
if (surfaceIndex_[faceI] != -1)
if (surfIndex[faceI] != -1)
{
label own = mesh_.faceOwner()[faceI];

31
src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C
View File

@ -325,4 +325,35 @@ void Foam::meshRefinement::weightedSum
}
template<class Type>
Type Foam::meshRefinement::get
(
const dictionary& dict,
const word& keyword,
const bool noExit,
enum keyType::option matchOpt,
const Type& defaultValue
)
{
Type val(defaultValue);
if
(
!dict.readEntry
(
keyword,
val,
matchOpt,
!noExit
)
)
{
FatalIOError
<< "Entry '" << keyword << "' not found in dictionary "
<< dict.name() << endl;
}
return val;
}
// ************************************************************************* //

103
src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C
View File

@ -28,6 +28,7 @@ License
#include "Tuple2.H"
#include "DynamicField.H"
#include "featureEdgeMesh.H"
#include "meshRefinement.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -41,7 +42,17 @@ void Foam::refinementFeatures::read
{
const dictionary& dict = featDicts[featI];
fileName featFileName(dict.lookup("file"));
fileName featFileName
(
meshRefinement::get<fileName>
(
dict,
"file",
dryRun_,
keyType::REGEX,
fileName::null
)
);
// Try reading extendedEdgeMesh first
@ -233,7 +244,19 @@ void Foam::refinementFeatures::read
else
{
// Look up 'level' for single level
levels_[featI] = labelList(1, dict.get<label>("level"));
levels_[featI] =
labelList
(
1,
meshRefinement::get<label>
(
dict,
"level",
dryRun_,
keyType::REGEX,
0
)
);
distances_[featI] = scalarField(1, Zero);
}
@ -439,14 +462,16 @@ Foam::refinementFeatures::regionEdgeTrees() const
Foam::refinementFeatures::refinementFeatures
(
const objectRegistry& io,
const PtrList<dictionary>& featDicts
const PtrList<dictionary>& featDicts,
const bool dryRun
)
:
PtrList<extendedFeatureEdgeMesh>(featDicts.size()),
distances_(featDicts.size()),
levels_(featDicts.size()),
edgeTrees_(featDicts.size()),
pointTrees_(featDicts.size())
pointTrees_(featDicts.size()),
dryRun_(dryRun)
{
// Read features
read(io, featDicts);
@ -531,6 +556,76 @@ Foam::refinementFeatures::refinementFeatures
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::refinementFeatures::checkSizes
(
const scalar maxRatio,
const boundBox& meshBb,
const bool report,
Ostream& os
) const
{
if (report)
{
os<< "Checking for size." << endl;
}
bool hasError = false;
forAll(*this, i)
{
const extendedFeatureEdgeMesh& em = operator[](i);
const boundBox bb(em.points(), true);
for (label j = i+1; j < size(); j++)
{
const extendedFeatureEdgeMesh& em2 = operator[](j);
const boundBox bb2(em2.points(), true);
scalar ratio = bb.mag()/bb2.mag();
if (ratio > maxRatio || ratio < 1.0/maxRatio)
{
hasError = true;
if (report)
{
os << " " << em.name()
<< " bounds differ from " << em2.name()
<< " by more than a factor 100:" << nl
<< " bounding box : " << bb << nl
<< " bounding box : " << bb2
<< endl;
}
}
}
}
forAll(*this, i)
{
const extendedFeatureEdgeMesh& em = operator[](i);
const boundBox bb(em.points(), true);
if (!meshBb.contains(bb))
{
if (report)
{
os << " " << em.name()
<< " bounds not fully contained in mesh"<< nl
<< " bounding box : " << bb << nl
<< " mesh bounding box : " << meshBb
<< endl;
}
}
}
if (report)
{
os<< endl;
}
return returnReduce(hasError, orOp<bool>());
}
void Foam::refinementFeatures::findNearestEdge
(
const pointField& samples,

16
src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H
View File

@ -67,6 +67,9 @@ class refinementFeatures
//- Features points
PtrList<indexedOctree<treeDataPoint>> pointTrees_;
//- Dry-run mode?
const bool dryRun_;
//- Region edge trees (demand driven)
mutable autoPtr<PtrList<indexedOctree<treeDataEdge>>>
regionEdgeTreesPtr_;
@ -112,7 +115,8 @@ public:
refinementFeatures
(
const objectRegistry& io,
const PtrList<dictionary>& featDicts
const PtrList<dictionary>& featDicts,
const bool dryRun = false
);
@ -189,6 +193,16 @@ public:
labelList& maxLevel
) const;
// Other
//- Check sizes - return true if error and stream to os
bool checkSizes
(
const scalar maxRatio,
const boundBox& meshBb,
const bool report,
Ostream& os
) const;
};

40
src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
View File

@ -32,6 +32,8 @@ License
#include "searchableSurfacesQueries.H"
#include "UPtrList.H"
#include "volumeType.H"
// For dictionary::get wrapper
#include "meshRefinement.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -124,14 +126,16 @@ Foam::refinementSurfaces::refinementSurfaces
(
const searchableSurfaces& allGeometry,
const dictionary& surfacesDict,
const label gapLevelIncrement
const label gapLevelIncrement,
const bool dryRun
)
:
allGeometry_(allGeometry),
surfaces_(surfacesDict.size()),
names_(surfacesDict.size()),
surfZones_(surfacesDict.size()),
regionOffset_(surfacesDict.size())
regionOffset_(surfacesDict.size()),
dryRun_(dryRun)
{
// Wildcard specification : loop over all surface, all regions
// and try to find a match.
@ -195,7 +199,18 @@ Foam::refinementSurfaces::refinementSurfaces
names_[surfI] = geomName;
surfaces_[surfI] = geomI;
const labelPair refLevel(dict.lookup("level"));
const labelPair refLevel
(
meshRefinement::get<labelPair>
(
dict,
"level",
dryRun_,
keyType::REGEX,
labelPair(0, 0)
)
);
globalMinLevel[surfI] = refLevel[0];
globalMaxLevel[surfI] = refLevel[1];
globalLevelIncr[surfI] = dict.lookupOrDefault
@ -281,7 +296,18 @@ Foam::refinementSurfaces::refinementSurfaces
regionNames[regionI]
);
const labelPair refLevel(regionDict.lookup("level"));
const labelPair refLevel
(
meshRefinement::get<labelPair>
(
regionDict,
"level",
dryRun_,
keyType::REGEX,
labelPair(0, 0)
)
);
regionMinLevel[surfI].insert(regionI, refLevel[0]);
regionMaxLevel[surfI].insert(regionI, refLevel[1]);
@ -476,7 +502,8 @@ Foam::refinementSurfaces::refinementSurfaces
const labelList& maxLevel,
const labelList& gapLevel,
const scalarField& perpendicularAngle,
PtrList<dictionary>& patchInfo
PtrList<dictionary>& patchInfo,
const bool dryRun
)
:
allGeometry_(allGeometry),
@ -488,7 +515,8 @@ Foam::refinementSurfaces::refinementSurfaces
maxLevel_(maxLevel),
gapLevel_(gapLevel),
perpendicularAngle_(perpendicularAngle),
patchInfo_(patchInfo.size())
patchInfo_(patchInfo.size()),
dryRun_(dryRun)
{
forAll(patchInfo_, pI)
{

9
src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
View File

@ -98,6 +98,9 @@ class refinementSurfaces
//- From global region number to patchType
PtrList<dictionary> patchInfo_;
//- Are we operating in test mode?
const bool dryRun_;
// Private Member Functions
@ -127,7 +130,8 @@ public:
(
const searchableSurfaces& allGeometry,
const dictionary&,
const label gapLevelIncrement
const label gapLevelIncrement,
const bool dryRun
);
//- Construct from components
@ -142,7 +146,8 @@ public:
const labelList& maxLevel,
const labelList& gapLevel,
const scalarField& perpendicularAngle,
PtrList<dictionary>& patchInfo
PtrList<dictionary>& patchInfo,
const bool dryRun
);

6
src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
View File

@ -569,10 +569,12 @@ void Foam::shellSurfaces::findLevel
Foam::shellSurfaces::shellSurfaces
(
const searchableSurfaces& allGeometry,
const dictionary& shellsDict
const dictionary& shellsDict,
const bool dryRun
)
:
allGeometry_(allGeometry)
allGeometry_(allGeometry),
dryRun_(dryRun)
{
// Wilcard specification : loop over all surfaces and try to find a match.

6
src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H
View File

@ -74,6 +74,9 @@ private:
//- Reference to all geometry.
const searchableSurfaces& allGeometry_;
//- Dry-run mode?
const bool dryRun_;
//- Indices of surfaces that are shells
labelList shells_;
@ -173,7 +176,8 @@ public:
shellSurfaces
(
const searchableSurfaces& allGeometry,
const dictionary& shellsDict
const dictionary& shellsDict,
const bool dryRun = false
);

27
src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C
View File

@ -96,12 +96,13 @@ Foam::scalar Foam::layerParameters::layerExpansionRatio
Foam::layerParameters::layerParameters
(
const dictionary& dict,
const polyBoundaryMesh& boundaryMesh
const polyBoundaryMesh& boundaryMesh,
const bool dryRun
)
:
dict_(dict),
numLayers_(boundaryMesh.size(), -1),
relativeSizes_(dict.lookup("relativeSizes")),
relativeSizes_(meshRefinement::get<bool>(dict, "relativeSizes", dryRun)),
layerSpec_(ILLEGAL),
firstLayerThickness_(boundaryMesh.size(), -123),
finalLayerThickness_(boundaryMesh.size(), -123),
@ -110,9 +111,9 @@ Foam::layerParameters::layerParameters
minThickness_
(
boundaryMesh.size(),
dict.get<scalar>("minThickness")
meshRefinement::get<scalar>(dict, "minThickness", dryRun)
),
featureAngle_(dict.get<scalar>("featureAngle")),
featureAngle_(meshRefinement::get<scalar>(dict, "featureAngle", dryRun)),
mergePatchFacesAngle_
(
dict.lookupOrDefault<scalar>
@ -125,16 +126,16 @@ Foam::layerParameters::layerParameters
(
dict.lookupOrDefault("concaveAngle", defaultConcaveAngle)
),
nGrow_(dict.get<label>("nGrow")),
nGrow_(meshRefinement::get<label>(dict, "nGrow", dryRun)),
maxFaceThicknessRatio_
(
dict.get<scalar>("maxFaceThicknessRatio")
meshRefinement::get<scalar>(dict, "maxFaceThicknessRatio", dryRun)
),
nBufferCellsNoExtrude_
(
dict.get<label>("nBufferCellsNoExtrude")
meshRefinement::get<label>(dict, "nBufferCellsNoExtrude", dryRun)
),
nLayerIter_(dict.get<label>("nLayerIter")),
nLayerIter_(meshRefinement::get<label>(dict, "nLayerIter", dryRun)),
nRelaxedIter_(labelMax),
additionalReporting_(dict.lookupOrDefault("additionalReporting", false)),
meshShrinker_
@ -144,7 +145,8 @@ Foam::layerParameters::layerParameters
"meshShrinker",
medialAxisMeshMover::typeName
)
)
),
dryRun_(dryRun)
{
// Detect layer specification mode
@ -255,7 +257,12 @@ Foam::layerParameters::layerParameters
}
const dictionary& layersDict = dict.subDict("layers");
const dictionary& layersDict = meshRefinement::subDict
(
dict,
"layers",
dryRun
);
for (const entry& dEntry : layersDict)
{

10
src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H
View File

@ -132,6 +132,9 @@ private:
word meshShrinker_;
//- In dry-run mode?
const bool dryRun_;
// Private Member Functions
@ -155,7 +158,12 @@ public:
// Constructors
//- Construct from dictionary
layerParameters(const dictionary& dict, const polyBoundaryMesh&);
layerParameters
(
const dictionary& dict,
const polyBoundaryMesh&,
const bool dryRun = false
);
// Member Functions

49
src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C
View File

@ -29,14 +29,28 @@ License
#include "globalIndex.H"
#include "Tuple2.H"
#include "wallPolyPatch.H"
#include "meshRefinement.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::refinementParameters::refinementParameters(const dictionary& dict)
Foam::refinementParameters::refinementParameters
(
const dictionary& dict,
const bool dryRun
)
:
maxGlobalCells_(dict.get<label>("maxGlobalCells")),
maxLocalCells_(dict.get<label>("maxLocalCells")),
minRefineCells_(dict.get<label>("minRefinementCells")),
maxGlobalCells_
(
meshRefinement::get<label>(dict, "maxGlobalCells", dryRun)
),
maxLocalCells_
(
meshRefinement::get<label>(dict, "maxLocalCells", dryRun)
),
minRefineCells_
(
meshRefinement::get<label>(dict, "minRefinementCells", dryRun)
),
planarAngle_
(
dict.lookupOrDefault
@ -45,7 +59,10 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
dict.get<scalar>("resolveFeatureAngle")
)
),
nBufferLayers_(dict.get<label>("nCellsBetweenLevels")),
nBufferLayers_
(
meshRefinement::get<label>(dict, "nCellsBetweenLevels", dryRun)
),
locationsOutsideMesh_
(
dict.lookupOrDefault
@ -55,7 +72,10 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
)
),
faceZoneControls_(dict.subOrEmptyDict("faceZoneControls")),
allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")),
allowFreeStandingZoneFaces_
(
meshRefinement::get<bool>(dict, "allowFreeStandingZoneFaces", dryRun)
),
useTopologicalSnapDetection_
(
dict.lookupOrDefault("useTopologicalSnapDetection", true)
@ -69,9 +89,11 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
(
dict.lookupOrDefault<Switch>("interfaceRefine", true)
),
nErodeCellZone_(dict.lookupOrDefault<label>("nCellZoneErodeIter", 0))
nErodeCellZone_(dict.lookupOrDefault<label>("nCellZoneErodeIter", 0)),
dryRun_(dryRun)
{
point locationInMesh;
List<Tuple2<point, word>> pointsToZone;
if (dict.readIfPresent("locationInMesh", locationInMesh))
{
locationsInMesh_.append(locationInMesh);
@ -84,7 +106,7 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
<< exit(FatalIOError);
}
}
else
else if (dict.readIfPresent("locationsInMesh", pointsToZone))
{
List<Tuple2<point, word>> pointsToZone(dict.lookup("locationsInMesh"));
label nZones = locationsInMesh_.size();
@ -102,9 +124,18 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
nZones++;
}
}
else
{
IOWarningInFunction(dict)
<< "No 'locationInMesh' or 'locationsInMesh' provided"
<< endl;
}
const scalar featAngle(dict.get<scalar>("resolveFeatureAngle"));
const scalar featAngle
(
meshRefinement::get<scalar>(dict, "resolveFeatureAngle", dryRun)
);
if (featAngle < 0 || featAngle > 180)
{

18
src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H
View File

@ -70,7 +70,7 @@ class refinementParameters
scalar curvature_;
//- Planarity criterion
scalar planarAngle_;
const scalar planarAngle_;
//- Number of layers between different refinement levels
const label nBufferLayers_;
@ -93,20 +93,22 @@ class refinementParameters
//- FaceZone faces allowed which have owner and neighbour in same
// cellZone?
Switch allowFreeStandingZoneFaces_;
const Switch allowFreeStandingZoneFaces_;
//- Use old topology based problem-cell removal (cells with 8 points
// on surface)
Switch useTopologicalSnapDetection_;
const Switch useTopologicalSnapDetection_;
//- Allowed load unbalance
scalar maxLoadUnbalance_;
const scalar maxLoadUnbalance_;
Switch handleSnapProblems_;
const Switch handleSnapProblems_;
Switch interfaceRefine_;
const Switch interfaceRefine_;
label nErodeCellZone_;
const label nErodeCellZone_;
const bool dryRun_;
// Private Member Functions
@ -123,7 +125,7 @@ public:
// Constructors
//- Construct from dictionary - new syntax
refinementParameters(const dictionary& dict);
refinementParameters(const dictionary& dict, const bool dryRun = false);
// Member Functions

26
src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.C
View File

@ -24,17 +24,30 @@ License
\*---------------------------------------------------------------------------*/
#include "snapParameters.H"
#include "meshRefinement.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from dictionary
Foam::snapParameters::snapParameters(const dictionary& dict)
Foam::snapParameters::snapParameters(const dictionary& dict, const bool dryRun)
:
nSmoothPatch_(dict.get<label>("nSmoothPatch")),
nSmoothPatch_
(
meshRefinement::get<label>(dict, "nSmoothPatch", dryRun)
),
nSmoothInternal_(dict.lookupOrDefault("nSmoothInternal", 0)),
snapTol_(dict.get<scalar>("tolerance")),
nSmoothDispl_(dict.get<label>("nSolveIter")),
nSnap_(dict.get<label>("nRelaxIter")),
snapTol_
(
meshRefinement::get<scalar>(dict, "tolerance", dryRun)
),
nSmoothDispl_
(
meshRefinement::get<label>(dict, "nSolveIter", dryRun)
),
nSnap_
(
meshRefinement::get<label>(dict, "nRelaxIter", dryRun)
),
nFeatureSnap_(dict.lookupOrDefault("nFeatureSnapIter", -1)),
explicitFeatureSnap_(dict.lookupOrDefault("explicitFeatureSnap", true)),
implicitFeatureSnap_(dict.lookupOrDefault("implicitFeatureSnap", false)),
@ -60,7 +73,8 @@ Foam::snapParameters::snapParameters(const dictionary& dict)
dict.lookupOrDefault("nFaceSplitInterval", labelMin)
),
concaveAngle_(dict.lookupOrDefault("concaveAngle", 45)),
minAreaRatio_(dict.lookupOrDefault("minAreaRatio", 0.3))
minAreaRatio_(dict.lookupOrDefault("minAreaRatio", 0.3)),
dryRun_(dryRun)
{}

11
src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.H
View File

@ -86,13 +86,16 @@ class snapParameters
//- How often needs face splitting be run
label nFaceSplitInterval_;
const label nFaceSplitInterval_;
//- When is angle too concave too split
scalar concaveAngle_;
const scalar concaveAngle_;
//- When is face-split not sufficiently diagonal
scalar minAreaRatio_;
const scalar minAreaRatio_;
//- in dry-run mode
const bool dryRun_;
// Private Member Functions
@ -109,7 +112,7 @@ public:
// Constructors
//- Construct from dictionary
snapParameters(const dictionary& dict);
snapParameters(const dictionary& dict, const bool dryRun = false);
// Member Functions

45
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
View File

@ -2648,7 +2648,8 @@ Foam::label Foam::snappyLayerDriver::checkAndUnmark
meshQualityDict,
identity(newMesh.nFaces()),
baffles,
wrongFaces
wrongFaces,
false // dryRun_
);
Info<< "Detected " << returnReduce(wrongFaces.size(), sumOp<label>())
<< " illegal faces"
@ -3215,12 +3216,14 @@ Foam::snappyLayerDriver::snappyLayerDriver
(
meshRefinement& meshRefiner,
const labelList& globalToMasterPatch,
const labelList& globalToSlavePatch
const labelList& globalToSlavePatch,
const bool dryRun
)
:
meshRefiner_(meshRefiner),
globalToMasterPatch_(globalToMasterPatch),
globalToSlavePatch_(globalToSlavePatch)
globalToSlavePatch_(globalToSlavePatch),
dryRun_(dryRun)
{}
@ -3894,6 +3897,40 @@ void Foam::snappyLayerDriver::addLayers
internalBaffles,
displacement
);
if (dryRun_)
{
string errorMsg(FatalError.message());
string IOerrorMsg(FatalIOError.message());
if (errorMsg.size() || IOerrorMsg.size())
{
//errorMsg = "[dryRun] " + errorMsg;
//errorMsg.replaceAll("\n", "\n[dryRun] ");
//IOerrorMsg = "[dryRun] " + IOerrorMsg;
//IOerrorMsg.replaceAll("\n", "\n[dryRun] ");
IOWarningInFunction(combinedDict)
<< nl
<< "Missing/incorrect required dictionary entries:"
<< nl << nl
<< IOerrorMsg.c_str() << nl
<< errorMsg.c_str() << nl << nl
<< "Exiting dry-run" << nl << endl;
if (Pstream::parRun())
{
Perr<< "\nFOAM parallel run exiting\n" << endl;
Pstream::exit(0);
}
else
{
Perr<< "\nFOAM exiting\n" << endl;
::exit(0);
}
}
}
}
@ -4691,7 +4728,7 @@ void Foam::snappyLayerDriver::doLayers
// Check initial mesh
Info<< "Checking initial mesh ..." << endl;
labelHashSet wrongFaces(mesh.nFaces()/100);
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces, dryRun_);
const label nInitErrors = returnReduce
(
wrongFaces.size(),

6
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H
View File

@ -106,6 +106,9 @@ private:
//- From surface region to patch
const labelList globalToSlavePatch_;
//- Are we operating in test mode?
const bool dryRun_;
// Private Member Functions
@ -618,7 +621,8 @@ public:
(
meshRefinement& meshRefiner,
const labelList& globalToMasterPatch,
const labelList& globalToSlavePatch
const labelList& globalToSlavePatch,
const bool dryRun = false
);

95
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
View File

@ -43,6 +43,7 @@ License
#include "searchableSurfaces.H"
#include "fvMeshSubset.H"
#include "interpolationTable.H"
#include "snappyVoxelMeshDriver.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -61,7 +62,8 @@ Foam::snappyRefineDriver::snappyRefineDriver
fvMeshDistribute& distributor,
const labelUList& globalToMasterPatch,
const labelUList& globalToSlavePatch,
const writer<scalar>& setFormatter
const writer<scalar>& setFormatter,
const bool dryRun
)
:
meshRefiner_(meshRefiner),
@ -69,7 +71,8 @@ Foam::snappyRefineDriver::snappyRefineDriver
distributor_(distributor),
globalToMasterPatch_(globalToMasterPatch),
globalToSlavePatch_(globalToSlavePatch),
setFormatter_(setFormatter)
setFormatter_(setFormatter),
dryRun_(dryRun)
{}
@ -82,6 +85,11 @@ Foam::label Foam::snappyRefineDriver::featureEdgeRefine
const label minRefine
)
{
if (dryRun_)
{
return 0;
}
if (refineParams.minRefineCells() == -1)
{
// Special setting to be able to restart shm on meshes with inconsistent
@ -199,6 +207,11 @@ Foam::label Foam::snappyRefineDriver::smallFeatureRefine
const label maxIter
)
{
if (dryRun_)
{
return 0;
}
if (refineParams.minRefineCells() == -1)
{
// Special setting to be able to restart shm on meshes with inconsistent
@ -209,7 +222,6 @@ Foam::label Foam::snappyRefineDriver::smallFeatureRefine
addProfiling(feature, "snappyHexMesh::refine::smallFeature");
const fvMesh& mesh = meshRefiner_.mesh();
label iter = 0;
// See if any surface has an extendedGapLevel
@ -329,6 +341,11 @@ Foam::label Foam::snappyRefineDriver::surfaceOnlyRefine
const label maxIter
)
{
if (dryRun_)
{
return 0;
}
if (refineParams.minRefineCells() == -1)
{
// Special setting to be able to restart shm on meshes with inconsistent
@ -460,6 +477,11 @@ Foam::label Foam::snappyRefineDriver::gapOnlyRefine
const label maxIter
)
{
if (dryRun_)
{
return 0;
}
if (refineParams.minRefineCells() == -1)
{
// Special setting to be able to restart shm on meshes with inconsistent
@ -706,6 +728,11 @@ Foam::label Foam::snappyRefineDriver::bigGapOnlyRefine
return 0;
}
if (dryRun_)
{
return 0;
}
const fvMesh& mesh = meshRefiner_.mesh();
label iter = 0;
@ -855,6 +882,11 @@ Foam::label Foam::snappyRefineDriver::danglingCellRefine
return 0;
}
if (dryRun_)
{
return 0;
}
addProfiling(dangling, "snappyHexMesh::refine::danglingCell");
const fvMesh& mesh = meshRefiner_.mesh();
@ -1008,6 +1040,11 @@ Foam::label Foam::snappyRefineDriver::refinementInterfaceRefine
return 0;
}
if (dryRun_)
{
return 0;
}
addProfiling(interface, "snappyHexMesh::refine::transition");
const fvMesh& mesh = meshRefiner_.mesh();
@ -1355,6 +1392,11 @@ void Foam::snappyRefineDriver::removeInsideCells
return;
}
if (dryRun_)
{
return;
}
Info<< nl
<< "Removing mesh beyond surface intersections" << nl
<< "------------------------------------------" << nl
@ -1420,6 +1462,11 @@ Foam::label Foam::snappyRefineDriver::shellRefine
const label maxIter
)
{
if (dryRun_)
{
return 0;
}
if (refineParams.minRefineCells() == -1)
{
// Special setting to be able to restart shm on meshes with inconsistent
@ -1601,6 +1648,11 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine
const label maxIter
)
{
if (dryRun_)
{
return 0;
}
addProfiling(shell, "snappyHexMesh::refine::directionalShell");
const fvMesh& mesh = meshRefiner_.mesh();
const shellSurfaces& shells = meshRefiner_.shells();
@ -2323,6 +2375,11 @@ void Foam::snappyRefineDriver::baffleAndSplitMesh
const dictionary& motionDict
)
{
if (dryRun_)
{
return;
}
addProfiling(split, "snappyHexMesh::refine::splitting");
Info<< nl
<< "Splitting mesh at surface intersections" << nl
@ -2454,6 +2511,11 @@ void Foam::snappyRefineDriver::splitAndMergeBaffles
const dictionary& motionDict
)
{
if (dryRun_)
{
return;
}
Info<< nl
<< "Handling cells with snap problems" << nl
<< "---------------------------------" << nl
@ -2656,6 +2718,11 @@ void Foam::snappyRefineDriver::mergePatchFaces
const dictionary& motionDict
)
{
if (dryRun_)
{
return;
}
addProfiling(merge, "snappyHexMesh::refine::merge");
Info<< nl
<< "Merge refined boundary faces" << nl
@ -2719,9 +2786,29 @@ void Foam::snappyRefineDriver::doRefine
const fvMesh& mesh = meshRefiner_.mesh();
// Check that all the keep points are inside the mesh.
refineParams.findCells(true, mesh, refineParams.locationsInMesh());
// Check that all the keep points are inside the mesh.
if (dryRun_)
{
refineParams.findCells(true, mesh, refineParams.locationsOutsideMesh());
}
// Estimate cell sizes
if (dryRun_)
{
snappyVoxelMeshDriver voxelDriver
(
meshRefiner_,
globalToMasterPatch_,
globalToSlavePatch_
);
voxelDriver.doRefine(refineParams);
}
// Refine around feature edges
featureEdgeRefine
(
@ -2892,7 +2979,7 @@ void Foam::snappyRefineDriver::doRefine
}
if (Pstream::parRun())
if (!dryRun_ && Pstream::parRun())
{
Info<< nl
<< "Doing final balancing" << nl

23
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H
View File

@ -39,6 +39,8 @@ SourceFiles
#include "scalarField.H"
#include "Tuple2.H"
#include "writer.H"
#include "DynamicList.H"
#include "labelVector.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -80,9 +82,27 @@ class snappyRefineDriver
//- How to write lines. Used e.g. when writing leak-paths
const writer<scalar>& setFormatter_;
//- Are we operating in test mode?
const bool dryRun_;
// Private Member Functions
void addNeighbours
(
const labelVector& n,
const labelList& cellLevel,
const labelVector& voxel,
const label voxeli,
DynamicList<labelVector>& front
) const;
//- Rough estimate of cell size and cell count
void estimateCellSizeAndCount
(
const refinementParameters& refineParams
) const;
//- Refine all cells pierced by explicit feature edges
label featureEdgeRefine
(
@ -225,7 +245,8 @@ public:
fvMeshDistribute& distributor,
const labelUList& globalToMasterPatch,
const labelUList& globalToSlavePatch,
const writer<scalar>& setFormatter
const writer<scalar>& setFormatter,
const bool dryRun = false
);

20
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
View File

@ -757,12 +757,14 @@ Foam::snappySnapDriver::snappySnapDriver
(
meshRefinement& meshRefiner,
const labelList& globalToMasterPatch,
const labelList& globalToSlavePatch
const labelList& globalToSlavePatch,
const bool dryRun
)
:
meshRefiner_(meshRefiner),
globalToMasterPatch_(globalToMasterPatch),
globalToSlavePatch_(globalToSlavePatch)
globalToSlavePatch_(globalToSlavePatch),
dryRun_(dryRun)
{}
@ -2606,7 +2608,11 @@ void Foam::snappySnapDriver::doSnap
if (snapParams.nFeatureSnap() > 0)
{
doFeatures = true;
nFeatIter = snapParams.nFeatureSnap();
if (!dryRun_)
{
nFeatIter = snapParams.nFeatureSnap();
}
Info<< "Snapping to features in " << nFeatIter
<< " iterations ..." << endl;
@ -2652,7 +2658,8 @@ void Foam::snappySnapDriver::doSnap
pointMesh::New(mesh),
adaptPatchIDs
),
motionDict
motionDict,
dryRun_
)
);
@ -2660,7 +2667,7 @@ void Foam::snappySnapDriver::doSnap
// Check initial mesh
Info<< "Checking initial mesh ..." << endl;
labelHashSet wrongFaces(mesh.nFaces()/100);
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces);
motionSmoother::checkMesh(false, mesh, motionDict, wrongFaces, dryRun_);
const label nInitErrors = returnReduce
(
wrongFaces.size(),
@ -2934,7 +2941,8 @@ void Foam::snappySnapDriver::doSnap
pointMesh::New(mesh),
adaptPatchIDs
),
motionDict
motionDict,
dryRun_
)
);

6
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H
View File

@ -67,6 +67,9 @@ class snappySnapDriver
//- From global surface region to slave side patch
const labelList globalToSlavePatch_;
//- Are we operating in test mode?
const bool dryRun_;
// Private Member Functions
@ -670,7 +673,8 @@ public:
(
meshRefinement& meshRefiner,
const labelList& globalToMasterPatch,
const labelList& globalToSlavePatch
const labelList& globalToSlavePatch,
const bool dryRun = false
);

594
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.C Normal file
View File

@ -0,0 +1,594 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "snappyVoxelMeshDriver.H"
#include "meshRefinement.H"
#include "fvMesh.H"
#include "Time.H"
#include "refinementParameters.H"
#include "refinementSurfaces.H"
#include "refinementFeatures.H"
#include "shellSurfaces.H"
#include "searchableSurfaces.H"
#include "voxelMeshSearch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(snappyVoxelMeshDriver, 0);
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::snappyVoxelMeshDriver::addNeighbours
(
const labelList& cellLevel,
const labelVector& voxel,
const label voxeli,
DynamicList<labelVector>& front
) const
{
const labelVector off(voxelMeshSearch::offset(n_));
for (direction dir = 0; dir < 3; ++dir)
{
if (voxel[dir] > 0)
{
labelVector left(voxel);
left[dir] -= 1;
if (cellLevel[voxeli-off[dir]] == -1)
{
front.append(left);
}
}
if (voxel[dir] < n_[dir]-1)
{
labelVector right(voxel);
right[dir] += 1;
if (cellLevel[voxeli+off[dir]] == -1)
{
front.append(right);
}
}
}
}
// Insert cell level for all volume refinement
Foam::tmp<Foam::pointField> Foam::snappyVoxelMeshDriver::voxelCentres() const
{
tmp<pointField> tcc(tmp<pointField>::New(n_.x()*n_.y()*n_.z()));
pointField& cc = tcc.ref();
const labelVector off(voxelMeshSearch::offset(n_));
label voxeli = voxelMeshSearch::index(n_, labelVector(0, 0, 0));
for (label k = 0; k < n_[2]; k++)
{
const label start1 = voxeli;
for (label j = 0; j < n_[1]; j++)
{
const label start0 = voxeli;
for (label i = 0; i < n_[0]; i++)
{
const labelVector voxel(i, j, k);
cc[voxeli] = voxelMeshSearch::centre(bb_, n_, voxel);
voxeli += off[0];
}
voxeli = start0 + off[1];
}
voxeli = start1 + off[2];
}
return tcc;
}
void Foam::snappyVoxelMeshDriver::isInside
(
const pointField& cc,
boolList& isVoxelInMesh
) const
{
const fvMesh& mesh = meshRefiner_.mesh();
isVoxelInMesh.setSize(cc.size());
if (isVoxelInMesh.size() < mesh.globalData().nTotalCells())
{
forAll(cc, voxeli)
{
const label celli = mesh.findCell
(
cc[voxeli],
polyMesh::FACE_PLANES
);
isVoxelInMesh[voxeli] = (celli != -1);
}
Pstream::listCombineGather(isVoxelInMesh, orEqOp<bool>());
}
else
{
const cellList& cells = mesh.cells();
const faceList& faces = mesh.faces();
const pointField& points = mesh.points();
for (label celli = 0; celli < mesh.nCells(); celli++)
{
const cell& cFaces = cells[celli];
boundBox cellBb(boundBox::invertedBox);
forAll(cFaces, cFacei)
{
const face& f = faces[cFaces[cFacei]];
forAll(f, fp)
{
cellBb.add(points[f[fp]]);
}
}
voxelMeshSearch::fill
(
isVoxelInMesh,
bb_,
n_,
cellBb,
1,
orEqOp<bool>()
);
}
Pstream::listCombineGather(isVoxelInMesh, orEqOp<bool>());
}
}
void Foam::snappyVoxelMeshDriver::markSurfaceRefinement
(
labelList& voxelLevel,
labelList& globalRegion
) const
{
// Insert cell level for all refinementSurfaces
const refinementSurfaces& s = meshRefiner_.surfaces();
forAll(s.surfaces(), surfi)
{
label geomi = s.surfaces()[surfi];
const searchableSurface& geom = s.geometry()[geomi];
//Pout<< "Geometry:" << s.names()[surfi] << endl;
if (isA<triSurface>(geom))
{
const triSurface& ts = refCast<const triSurface>(geom);
const pointField& points = ts.points();
forAll(ts, trii)
{
label regioni = ts[trii].region();
label globalRegioni = s.regionOffset()[surfi]+regioni;
const boundBox triBb(points, ts[trii], false);
// Fill cellLevel
label level = s.minLevel()[globalRegioni];
voxelMeshSearch::fill
(
voxelLevel,
bb_,
n_,
triBb,
level,
maxEqOp<label>()
);
voxelMeshSearch::fill
(
globalRegion,
bb_,
n_,
triBb,
globalRegioni,
maxEqOp<label>()
);
}
}
// else: maybe do intersection tests?
}
}
void Foam::snappyVoxelMeshDriver::findVoxels
(
const labelList& voxelLevel,
const pointField& locationsOutsideMesh,
labelList& voxels
) const
{
voxels.setSize(locationsOutsideMesh.size());
voxels = -1;
forAll(locationsOutsideMesh, loci)
{
const point& pt = locationsOutsideMesh[loci];
label voxeli = voxelMeshSearch::index(bb_, n_, pt, false);
if (voxeli == -1 || voxelLevel[voxeli] == labelMax)
{
WarningInFunction << "Location outside mesh "
<< pt << " is outside mesh with bounding box "
<< bb_ << endl;
}
else
{
voxels[loci] = voxeli;
}
}
}
void Foam::snappyVoxelMeshDriver::floodFill
(
const label startVoxeli,
const label newLevel,
labelList& voxelLevel
) const
{
DynamicList<labelVector> front;
front.append(voxelMeshSearch::index3(n_, startVoxeli));
DynamicList<labelVector> newFront;
while (true)
{
newFront.clear();
for (const auto& voxel : front)
{
label voxeli = voxelMeshSearch::index(n_, voxel);
if (voxelLevel[voxeli] == -1)
{
voxelLevel[voxeli] = 0;
addNeighbours
(
voxelLevel,
voxel,
voxeli,
newFront
);
}
}
if (newFront.empty())
{
break;
}
front.transfer(newFront);
}
}
void Foam::snappyVoxelMeshDriver::max
(
const labelList& maxLevel,
labelList& voxelLevel
) const
{
// Mark voxels with level
const labelVector off(voxelMeshSearch::offset(n_));
label voxeli = voxelMeshSearch::index(n_, labelVector(0, 0, 0));
for (label k = 0; k < n_[2]; k++)
{
const label start1 = voxeli;
for (label j = 0; j < n_[1]; j++)
{
const label start0 = voxeli;
for (label i = 0; i < n_[0]; i++)
{
voxelLevel[voxeli] = Foam::max
(
voxelLevel[voxeli],
maxLevel[voxeli]
);
voxeli += off[0];
}
voxeli = start0 + off[1];
}
voxeli = start1 + off[2];
}
}
Foam::labelList Foam::snappyVoxelMeshDriver::count
(
const labelList& voxelLevel
) const
{
label maxLevel = 0;
for (const auto level : voxelLevel)
{
if (level != labelMax)
{
maxLevel = Foam::max(maxLevel, level);
}
}
labelList count(maxLevel+1, 0);
const labelVector off(voxelMeshSearch::offset(n_));
label voxeli = voxelMeshSearch::index(n_, labelVector(0, 0, 0));
for (label k = 0; k < n_[2]; k++)
{
const label start1 = voxeli;
for (label j = 0; j < n_[1]; j++)
{
const label start0 = voxeli;
for (label i = 0; i < n_[0]; i++)
{
label level = voxelLevel[voxeli];
if (level != -1 && level != labelMax)
{
++count[level];
}
voxeli += off[0];
}
voxeli = start0 + off[1];
}
voxeli = start1 + off[2];
}
return count;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::snappyVoxelMeshDriver::snappyVoxelMeshDriver
(
meshRefinement& meshRefiner,
const labelUList& globalToMasterPatch,
const labelUList& globalToSlavePatch
)
:
meshRefiner_(meshRefiner),
globalToMasterPatch_(globalToMasterPatch),
globalToSlavePatch_(globalToSlavePatch),
bb_(meshRefiner_.mesh().bounds())
{
label maxLevel = labelMin;
// Feature refinement
const labelListList& featLevels = meshRefiner_.features().levels();
forAll(featLevels, feati)
{
maxLevel = Foam::max(maxLevel, Foam::max(featLevels[feati]));
}
// Surface refinement
const labelList& surfaceLevels = meshRefiner_.surfaces().maxLevel();
maxLevel = Foam::max(maxLevel, Foam::max(surfaceLevels));
// Shell refinement
maxLevel = Foam::max(maxLevel, meshRefiner_.shells().maxLevel());
const scalar level0Len = meshRefiner_.meshCutter().level0EdgeLength();
Info<< "Cell size estimate :" << nl
<< " Maximum : " << level0Len << nl
<< " Minimum : " << level0Len/pow(2.0, maxLevel) << nl
<< endl;
// Define voxel mesh with similar dimensions as mesh
const vector meshSpan(bb_.span());
n_ = labelVector
(
round(meshSpan.x()/level0Len),
round(meshSpan.y()/level0Len),
round(meshSpan.z()/level0Len)
);
label nTot = n_.x()*n_.y()*n_.z();
while (nTot < 1000000) //1048576)
{
n_ *= 2;
nTot = n_.x()*n_.y()*n_.z();
}
Info<< "Voxellating initial mesh : " << n_ << nl << endl;
tmp<pointField> tcc(voxelCentres());
const pointField& cc = tcc();
Info<< "Voxel refinement :" << nl
<< " Initial : (" << nTot << ')' << endl;
boolList isVoxelInMesh;
isInside(cc, isVoxelInMesh);
if (Pstream::master())
{
voxelLevel_.setSize(nTot, -1);
globalRegion_.setSize(nTot, -1);
// Remove cells outside initial mesh
forAll(isVoxelInMesh, voxeli)
{
if (!isVoxelInMesh[voxeli])
{
voxelLevel_[voxeli] = labelMax;
globalRegion_[voxeli] = -1;
}
}
//if (debug)
{
Info<< " After removing outside cells : " << count(voxelLevel_)
<< endl;
}
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::snappyVoxelMeshDriver::doRefine
(
const refinementParameters& refineParams
)
{
const scalar level0Len = meshRefiner_.meshCutter().level0EdgeLength();
tmp<pointField> tcc(voxelCentres());
const pointField& cc = tcc();
boolList isVoxelInMesh;
isInside(cc, isVoxelInMesh);
if (Pstream::master())
{
// Mark voxels containing (parts of) triangles
markSurfaceRefinement(voxelLevel_, globalRegion_);
//if (debug)
{
Info<< " After surface refinement : " << count(voxelLevel_)
<< endl;
}
// Find outside locations (and their current refinement level)
const pointField& outsidePoints = refineParams.locationsOutsideMesh();
labelList outsideMeshVoxels;
findVoxels
(
voxelLevel_,
outsidePoints,
outsideMeshVoxels
);
labelList outsideOldLevel(outsideMeshVoxels.size(), -1);
forAll(outsideMeshVoxels, loci)
{
label voxeli = outsideMeshVoxels[loci];
if (voxeli >= 0)
{
outsideOldLevel[loci] = voxelLevel_[outsideMeshVoxels[loci]];
if (outsideOldLevel[loci] >= 0)
{
WarningInFunction << "Location outside mesh "
<< outsidePoints[loci]
<< " is inside mesh or close to surface" << endl;
}
}
}
// Find inside locations
labelList insideMeshVoxels;
findVoxels
(
voxelLevel_,
refineParams.locationsInMesh(),
insideMeshVoxels
);
forAll(insideMeshVoxels, loci)
{
label voxeli = insideMeshVoxels[loci];
if (voxeli != -1)
{
if (voxelLevel_[voxeli] != -1)
{
WarningInFunction << "Location inside mesh "
<< refineParams.locationsInMesh()[loci]
<< " is marked as a surface voxel " << voxeli
<< " with cell level " << voxelLevel_[voxeli] << endl;
}
else
{
// Flood-fill out from voxel
floodFill(voxeli, 0, voxelLevel_);
}
}
}
//if (debug)
{
Info<< " After keeping inside voxels : " << count(voxelLevel_)
<< endl;
}
// Re-check the outside locations to see if they have been bled into
{
forAll(outsideMeshVoxels, loci)
{
label voxeli = outsideMeshVoxels[loci];
if (voxeli >= 0 && voxelLevel_[voxeli] != outsideOldLevel[loci])
{
WarningInFunction << "Location outside mesh "
<< outsidePoints[loci]
<< " is reachable from an inside location" << nl
<< "Either your locations are too close to the"
<< " geometry or there might be a leak in the"
<< " geometry" << endl;
}
}
}
// Shell refinement : find ccs inside higher shells
labelList maxLevel;
meshRefiner_.shells().findHigherLevel(cc, voxelLevel_, maxLevel);
// Assign max of maxLevel and voxelLevel
max(maxLevel, voxelLevel_);
// Determine number of levels
const labelList levelCounts(count(voxelLevel_));
//if (debug)
{
Info<< " After shell refinement : " << levelCounts << endl;
}
const vector meshSpan(bb_.span());
const vector voxel0Size
(
meshSpan[0]/n_[0],
meshSpan[1]/n_[1],
meshSpan[2]/n_[2]
);
label cellCount = 0;
forAll(levelCounts, leveli)
{
const scalar s = level0Len/pow(2.0, leveli);
const scalar nCellsPerVoxel
(
voxel0Size[0]/s
*voxel0Size[1]/s
*voxel0Size[2]/s
);
cellCount += levelCounts[leveli]*nCellsPerVoxel;
}
Info<< "Estimated cell count : " << cellCount << endl;
}
}
// ************************************************************************* //

172
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.H Normal file
View File

@ -0,0 +1,172 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2018 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::snappyVoxelMeshDriver
Description
Equivalent of snappyRefineDriver but operating on a voxel mesh.
Used to estimate cell size count from refinement (currently).
SourceFiles
snappyVoxelMeshDriver.C
\*---------------------------------------------------------------------------*/
#ifndef snappyVoxelMeshDriver_H
#define snappyVoxelMeshDriver_H
#include "labelList.H"
#include "DynamicList.H"
#include "labelVector.H"
#include "boundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declarations
class refinementParameters;
class layerParameters;
class meshRefinement;
/*---------------------------------------------------------------------------*\
Class snappyVoxelMeshDriver Declaration
\*---------------------------------------------------------------------------*/
class snappyVoxelMeshDriver
{
// Private data
//- Mesh+surface
meshRefinement& meshRefiner_;
//- From surface region to patch
const labelList globalToMasterPatch_;
//- From surface region to patch
const labelList globalToSlavePatch_;
//- Bounding box
boundBox bb_;
//- Number of voxels
labelVector n_;
//- Per voxel the refinement level
labelList voxelLevel_;
//- Per insersected voxel the originating surface
labelList globalRegion_;
// Private Member Functions
void addNeighbours
(
const labelList& cellLevel,
const labelVector& voxel,
const label voxeli,
DynamicList<labelVector>& front
) const;
tmp<pointField> voxelCentres() const;
void isInside
(
const pointField& voxelCentres,
boolList& isVoxelInMesh
) const;
void markSurfaceRefinement
(
labelList& voxelLevel,
labelList& globalRegion
) const;
void findVoxels
(
const labelList& voxelLevel,
const pointField& locationsOutsideMesh,
labelList& voxels
) const;
void floodFill
(
const label voxeli,
const label newLevel,
labelList& voxelLevel
) const;
void max
(
const labelList& maxLevel,
labelList& voxelLevel
) const;
labelList count(const labelList& voxelLevel) const;
//- No copy construct
snappyVoxelMeshDriver(const snappyVoxelMeshDriver&) = delete;
//- No copy assignment
void operator=(const snappyVoxelMeshDriver&) = delete;
public:
//- Runtime type information
ClassName("snappyVoxelMeshDriver");
// Constructors
//- Construct from components
snappyVoxelMeshDriver
(
meshRefinement& meshRefiner,
const labelUList& globalToMasterPatch,
const labelUList& globalToSlavePatch
);
// Member Functions
void doRefine(const refinementParameters& refineParams);
//void doLayers(const layerParameters& layerParams);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

72
src/meshTools/searchableSurfaces/searchableSurfaces/searchableSurfaces.C
View File

@ -48,18 +48,51 @@ bool Foam::searchableSurfaces::connected
const pointIndexHit& hit
)
{
const triFace& localFace = s.localFaces()[hit.index()];
const edge& e = s.edges()[edgeI];
const labelList& mp = s.meshPoints();
const edge meshE(mp[e[0]], mp[e[1]]);
forAll(localFace, i)
const triFace& f = s[hit.index()];
forAll(f, i)
{
if (e.otherVertex(localFace[i]) != -1)
if (meshE.otherVertex(f[i]) != -1)
{
return true;
}
}
return false;
// Account for triangle intersection routine going wrong for
// lines in same plane as triangle. Tbd!
vector eVec(meshE.vec(s.points()));
scalar magEVec(mag(eVec));
if (magEVec > ROOTVSMALL)
{
vector n(f.areaNormal(s.points()));
scalar magArea(mag(n));
if (magArea > ROOTVSMALL)
{
n /= magArea;
if (mag(n&(eVec/magEVec)) < SMALL)
{
// Bad intersection
return true;
}
else
{
return false;
}
}
else
{
return true;
}
}
else
{
return true;
}
}
@ -619,20 +652,20 @@ bool Foam::searchableSurfaces::checkIntersection
{
List<pointIndexHit> hits;
if (i == j)
{
// Slightly shorten the edges to prevent finding lots of
// intersections. The fast triangle intersection routine
// has problems with rays passing through a point of the
// triangle.
vectorField d
(
max(tolerance, 10*s0.tolerance())
*(end-start)
);
start += d;
end -= d;
}
//if (i == j)
//{
// // Slightly shorten the edges to prevent finding lots of
// // intersections. The fast triangle intersection routine
// // has problems with rays passing through a point of the
// // triangle. Now done in 'connected' routine. Tbd.
// vectorField d
// (
// max(tolerance, 10*s0.tolerance())
// *(end-start)
// );
// start += d;
// end -= d;
//}
operator[](j).findLineAny(start, end, hits);
@ -654,6 +687,9 @@ bool Foam::searchableSurfaces::checkIntersection
}
}
// tdb. What about distributedTriSurfaceMesh?
//reduce(nHits, sumOp<label>());
if (nHits > 0)
{
if (report)

39
src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C
View File

@ -42,6 +42,43 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::labelVector Foam::voxelMeshSearch::offset
(
const labelVector& nDivs
)
{
return labelVector(1, nDivs.x(), nDivs.x()*nDivs.y());
}
Foam::label Foam::voxelMeshSearch::index
(
const labelVector& nDivs,
const labelVector& voxel
)
{
return voxel.x()+voxel.y()*nDivs.x()+voxel.z()*nDivs.x()*nDivs.y();
}
Foam::labelVector Foam::voxelMeshSearch::index3
(
const labelVector& nDivs,
label voxeli
)
{
const label nxy = nDivs.x()*nDivs.y();
labelVector voxel;
voxel.z() = voxeli/nxy;
voxeli = voxeli % nxy;
voxel.y() = voxeli/nDivs.x();
voxel.x() = voxeli%nDivs.x();
return voxel;
}
Foam::labelVector Foam::voxelMeshSearch::index3
(
const boundBox& bb,
@ -98,7 +135,7 @@ Foam::label Foam::voxelMeshSearch::index
return -1;
}
return v[0] + g[0]*v[1] + g[1]*g.y()*v[2];
return index(g, v);
}

48
src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H
View File

@ -122,6 +122,26 @@ public:
//Voxel helper functions
//- Voxel indices to combined voxel index
static label index
(
const labelVector& nDivs,
const labelVector& voxel
);
//- Change in combined voxel index for change in components
static labelVector offset
(
const labelVector& nDivs
);
//- Combined voxel index to individual indices
static labelVector index3
(
const labelVector& nDivs,
const label voxeli
);
//- Coordinate to voxel indices
static labelVector index3
(
@ -130,6 +150,14 @@ public:
const point& p
);
//- Voxel index to voxel centre
static point centre
(
const boundBox& bb,
const labelVector& nDivs,
const labelVector& voxel
);
//- Coordinate to combined voxel index. If clip makes sure
// components are all inside. If not clip returns -1 if outside bb.
static label index
@ -140,14 +168,6 @@ public:
const bool clip
);
//- Voxel index to voxel centre
static point centre
(
const boundBox& bb,
const labelVector& nDivs,
const labelVector& voxel
);
//- Fill voxels indicated by bounding box
template<class Container, class Type>
static void fill
@ -159,6 +179,18 @@ public:
const Type val
);
//- Fill voxels indicated by bounding box
template<class Container, class Type, class CombineOp>
static void fill
(
Container& elems,
const boundBox& bb,
const labelVector& nDivs,
const boundBox& subBb,
const Type val,
const CombineOp& cop = eqOp<Type>()
);
//- Check if any voxel inside bounding box is set to val or
// not set to val (isNot = true)
template<class Container, class Type>

84
src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearchTemplates.C
View File

@ -52,15 +52,67 @@ void Foam::voxelMeshSearch::fill
minIds[cmpt] = max(minIds[cmpt], 0);
}
for (label i = minIds[0]; i <= maxIds[0]; i++)
const labelVector off(offset(nDivs));
label voxeli = index(nDivs, minIds);
for (label k = minIds[2]; k <= maxIds[2]; k++)
{
const label start1 = voxeli;
for (label j = minIds[1]; j <= maxIds[1]; j++)
{
for (label k = minIds[2]; k <= maxIds[2]; k++)
const label start0 = voxeli;
for (label i = minIds[0]; i <= maxIds[0]; i++)
{
elems[i+j*nDivs.x()+k*nDivs.x()*nDivs.y()] = val;
elems[voxeli] = val;
voxeli += off[0];
}
voxeli = start0 + off[1];
}
voxeli = start1 + off[2];
}
}
template<class Container, class Type, class CombineOp>
void Foam::voxelMeshSearch::fill
(
Container& elems,
const boundBox& bb,
const labelVector& nDivs,
const boundBox& subBb,
const Type val,
const CombineOp& cop
)
{
labelVector minIds(index3(bb, nDivs, subBb.min()));
labelVector maxIds(index3(bb, nDivs, subBb.max()));
for (direction cmpt = 0; cmpt < 3; cmpt++)
{
if (maxIds[cmpt] < 0 || minIds[cmpt] >= nDivs[cmpt])
{
return;
}
// Clip
maxIds[cmpt] = min(maxIds[cmpt], nDivs[cmpt]-1);
minIds[cmpt] = max(minIds[cmpt], 0);
}
const labelVector off(offset(nDivs));
label voxeli = index(nDivs, minIds);
for (label k = minIds[2]; k <= maxIds[2]; k++)
{
const label start1 = voxeli;
for (label j = minIds[1]; j <= maxIds[1]; j++)
{
const label start0 = voxeli;
for (label i = minIds[0]; i <= maxIds[0]; i++)
{
cop(elems[voxeli], val);
voxeli += off[0];
}
voxeli = start0 + off[1];
}
voxeli = start1 + off[2];
}
}
@ -100,19 +152,26 @@ bool Foam::voxelMeshSearch::overlaps
}
for (label i = minIds[0]; i <= maxIds[0]; i++)
const labelVector off(offset(nDivs));
label voxeli = index(nDivs, minIds);
for (label k = minIds[2]; k <= maxIds[2]; k++)
{
const label start1 = voxeli;
for (label j = minIds[1]; j <= maxIds[1]; j++)
{
for (label k = minIds[2]; k <= maxIds[2]; k++)
const label start0 = voxeli;
for (label i = minIds[0]; i <= maxIds[0]; i++)
{
const Type elemVal = elems[i+j*nDivs.x()+k*nDivs.x()*nDivs.y()];
const Type elemVal = elems[voxeli];
if (isNot != (elemVal == val))
{
return true;
}
voxeli += off[0];
}
voxeli = start0 + off[1];
}
voxeli = start1 + off[2];
}
return false;
}
@ -136,19 +195,26 @@ void Foam::voxelMeshSearch::write
<< exit(FatalError);
}
for (label i = 0; i < nDivs[0]; i++)
const labelVector off(offset(nDivs));
label voxeli = index(nDivs, labelVector(0, 0, 0));
for (label k = 0; k < nDivs[2]; k++)
{
const label start1 = voxeli;
for (label j = 0; j < nDivs[1]; j++)
{
for (label k = 0; k < nDivs[2]; k++)
const label start0 = voxeli;
for (label i = 0; i < nDivs[0]; i++)
{
const Type elemVal = elems[i+j*nDivs.x()+k*nDivs.x()*nDivs.y()];
const Type& elemVal = elems[voxeli];
if (isNot != (elemVal == val))
{
os.write(centre(bb, nDivs, labelVector(i, j, k)));
}
voxeli += off[0];
}
voxeli = start0 + off[1];
}
voxeli = start1 + off[2];
}
}