diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C index dd8bfc5246b18e42d041018b29df9af7a6364d55..f99aaab8577e9f9c9fef2a0cc2029dff92045f8d 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C @@ -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,11 +1128,50 @@ 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; diff --git a/src/OpenFOAM/db/error/error.C b/src/OpenFOAM/db/error/error.C index 4cc8b0a8fb7606822ffc69d3cf8480f10ef38c11..1ff7a969192b52d0eb568a3bd3f3fe574b992607 100644 --- a/src/OpenFOAM/db/error/error.C +++ b/src/OpenFOAM/db/error/error.C @@ -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) diff --git a/src/OpenFOAM/db/error/error.H b/src/OpenFOAM/db/error/error.H index 79dd7f3166832bc5fe2dbdea60997be63aceda16..39d45d72e8b3f63835ce81444848652431241f91 100644 --- a/src/OpenFOAM/db/error/error.H +++ b/src/OpenFOAM/db/error/error.H @@ -114,6 +114,9 @@ public: string message() const; + //- Clear any messages + void clear() const; + inline const string& functionName() const { return functionName_; diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.C b/src/dynamicMesh/motionSmoother/motionSmoother.C index d30ee35732cfe4e8032eb6035948eea2b63f5a56..9c4a300e1644e89608fa2d875e9596a98d7e32d4 100644 --- a/src/dynamicMesh/motionSmoother/motionSmoother.C +++ b/src/dynamicMesh/motionSmoother/motionSmoother.C @@ -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 ) {} diff --git a/src/dynamicMesh/motionSmoother/motionSmoother.H b/src/dynamicMesh/motionSmoother/motionSmoother.H index 2b45646004961b28e9ecfe53bfaa588de0bc816f..31c29bb8b54e003a460fb8f38b99ff15b50bfd20 100644 --- a/src/dynamicMesh/motionSmoother/motionSmoother.H +++ b/src/dynamicMesh/motionSmoother/motionSmoother.H @@ -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 ); }; diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C index ee8a1137e4dd7eb8abc8fd45f2671ee80893ea90..35b5c1cdf7e9d0ec85973dd57901ba8389383b52 100644 --- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C +++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.C @@ -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) { diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H index 742962e5bb0e1b84d45ecbd9c2de89398f0a31ea..cabc041af2747a7ba18202c6c029f20815a2895b 100644 --- a/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H +++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgo.H @@ -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 +// ); }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C index 361fbe216e106ecb567c0907283538a8e44821cc..003067656f4808e52763ba4841609b3cc26a2465 100644 --- a/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C +++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C @@ -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; diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C index a228c61bb79663030fae0fc8baa5cf6bc8e91bb0..eea46f3ee0b002c55383972470c4ec039a4b319c 100644 --- a/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C +++ b/src/dynamicMesh/motionSmoother/motionSmootherAlgoTemplates.C @@ -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; +} + + // ************************************************************************* // diff --git a/src/mesh/snappyHexMesh/Make/files b/src/mesh/snappyHexMesh/Make/files index 49a981e4afdba010b23b2b2969ead66be61026c1..00c46b066eac113306e1aa34d353bea023460d9d 100644 --- a/src/mesh/snappyHexMesh/Make/files +++ b/src/mesh/snappyHexMesh/Make/files @@ -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 diff --git a/src/mesh/snappyHexMesh/Make/options b/src/mesh/snappyHexMesh/Make/options index bff1be5f77e38ed536d5711be5365565442813a7..2cf98acf6b29ca68bdb02f18d3c558b356bc7f99 100644 --- a/src/mesh/snappyHexMesh/Make/options +++ b/src/mesh/snappyHexMesh/Make/options @@ -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 diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C index 5fd4aa2af6f4d686d7ea90f28db697c60582765a..e4fa24a87483b54e3fc5f1fe2b5dacff33480c93 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.C @@ -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; diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H index b4df2ef5c4c4a2b52a744bf0221533e56a88a8bb..6462958f99588a8d35855f61bcafaa24adf48080 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/displacementMotionSolverMeshMover.H @@ -112,7 +112,8 @@ public: ( const dictionary& dict, const List<labelPair>& baffles, - pointVectorField& pointDisplacemen + pointVectorField& pointDisplacement, + const bool dryRun ); diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C index 43ee2771ad44aacbaf1ded79e9a233c8b3b5b67a..f7fa6213ee2ee313c1092204a3e0e2f00e59d617 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.C @@ -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) ); } diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H index f9d65ef82517ac2251a3b858b3824cdbe8d799b5..9656145f64039dae8ed32d95a08b90fa14f07f21 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/externalDisplacementMeshMover.H @@ -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 ); diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C index 523b2a8bc50799ecfd212a1603a18f0310d95e06..6b044a87001cc8dcf2c20f8c12067e92e9acc34f 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C @@ -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> diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H index 8f3ea990f4cc8d08febdf74099c5f6b52e8bdaa0..f47a99c19a343500ba16f9da8245a3b58e0c58f7 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H @@ -207,7 +207,8 @@ public: ( const dictionary& dict, const List<labelPair>& baffles, - pointVectorField& pointDisplacement + pointVectorField& pointDisplacement, + const bool dryRun ); diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C index 404042246b252a82594d77f692aa52f007c288eb..effc9cb93549a91adb994715934f542757eed157 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C @@ -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); + } +} + + // ************************************************************************* // diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H index f67d0722ecde6c0e9b2c5f0bae826859dc7db5c9..30dcbd8878301b824ad94e8214524094ca4bb873 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H @@ -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 + ); }; diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C index 01da7090da445d7596a720b1015f1d0b01f325bc..4982bba7f9a460f3df49375b823c8b4e73552811 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C @@ -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 ( diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C index a7cc9e377e14d350d418a6ddb9ff2f283b8bb271..14c293e2fbdc392a29e80edbad6abd03463c6902 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C @@ -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) //{ diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C index 2476b0922e7786934c4dd9a0ce2c7662df88af6b..b0eb7cc6eb6b6a88d600b39402dcd6471be8876c 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C @@ -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 diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C index 1315caa1c74e0260f017866597e89372ac27bf3c..5ed47ac1ef2b41a4b2b9a8fe355a5af1548e6958 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C @@ -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]; diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C index ff7c829fcf198c4c8ff5ab603302af0513faa527..2d7c6a3680235097b8f05c958bd01ee8eb010c8c 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C @@ -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; +} + + // ************************************************************************* // diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C index 86d322c192d62dae59c1cfa7cf9da71ff468a32e..d16b9eb0df3f3a9881a5b541dc71b9adc49fcecb 100644 --- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C +++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C @@ -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, diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H index 6147c81bcf311dbcac00a69ea6550a0f2279dfd8..b8b46ac8f3b4471dfa6a50a409e584b96ee5e96d 100644 --- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H +++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.H @@ -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; }; diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C index c300ceaa408f6f4882f2d1b06d0b4f883826affc..61d5c9607cd6808e9d6ccc5a0c26f5839b828707 100644 --- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C +++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C @@ -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) { diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H index d3939c763e46082a771399dd34db2fce44d9be42..c716f7ef5ad80248efe7fab786920f8c09f89091 100644 --- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H +++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H @@ -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 ); diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C index 5739217fe2b528c72d1dd6b4b3b8914cafe5665a..35c83377c512f81114fdc131c5be39cedb9cf58b 100644 --- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C +++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C @@ -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. diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H index 2390a98e95566c9d4e128812dc20c2d31dbc0b2b..f5bc6f7b636078082c21671d3cfb3c0205a42b82 100644 --- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H +++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H @@ -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 ); diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C index 4e5f08bfb5cf32e5bf9c3197d52e7e98699ed68d..00a76b59aa28557fb31b32ce7772b38be569730b 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.C @@ -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) { diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H index 34a9e813ec682dcdc5ea906c9c50c6b0ffa7ae9e..e5d586e6a6211d548488ab314b496f922f43964b 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/layerParameters/layerParameters.H @@ -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 diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C index c63345f9939a2dd0152e290cf15374e0c2a713dd..ab84e69831157ce90557e6004879c4726002317b 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.C @@ -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) { diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H index a6475c00d4b8642a67f687a46c2513007cc0c912..29112d02efddb1b53ae139e32c2c43bf6c8886be 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/refinementParameters/refinementParameters.H @@ -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 diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.C index 3f19065202193849ae9dae92b91eabf8f0248203..7c9693dc598e2f6a86ec30649b2b1cdf2502a83f 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.C @@ -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) {} diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.H index 18aca4d8c790756ed2b92b1fe3697af26195f2d8..1f72e3b8bd7be49f8730efb3f0cc194bfb086f7c 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.H +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snapParameters/snapParameters.H @@ -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 diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C index ebae51274f4f68697f9a4a011c82dcad239bc8b1..f794e54437142d5efaf25b2e4b858cd79ff8f529 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C @@ -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(), diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H index 9c663cf70b9f0cb7a2e2056ebddfd908737659c4..c0df3db8a2d0cfd3158879a8ddc4ab7960e2755b 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H @@ -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 ); diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C index 8bd94a084e36f2317f52d24e3c89db14266eab80..f05b6fe65978f4781acdffc17f9229db1354e83d 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C @@ -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 diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H index bab0d95ed7ce1346ad0946be6a9ac8607659f8e4..fdb9ca7dc0eb859c9a346fde900d25403dd342d7 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H @@ -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 ); diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C index b90d4b33272c28a38487f720ab09adc406d75d6f..c7a7f043b93261551e8e6f5b6bc67b98e72a4b80 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C @@ -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_ ) ); diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H index 7d7f9c4056599a5edce5f0a51cd17aa31f3b77b4..42783d354cdcdd579dbde93192cc72db3502be04 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H @@ -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 ); diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.C new file mode 100644 index 0000000000000000000000000000000000000000..126444be2b935840e0680e6d0b7a8479619c4a5d --- /dev/null +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.C @@ -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; + } +} + + +// ************************************************************************* // diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.H new file mode 100644 index 0000000000000000000000000000000000000000..e241bfcd85b2ad04833173104a7da9b9ddc8cb82 --- /dev/null +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyVoxelMeshDriver.H @@ -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 + +// ************************************************************************* // diff --git a/src/meshTools/searchableSurfaces/searchableSurfaces/searchableSurfaces.C b/src/meshTools/searchableSurfaces/searchableSurfaces/searchableSurfaces.C index 144373d15dea4b99e2717cc70a6e9bab6467535c..91ccdc8d3d037d449bbf3116d848da36988cddaa 100644 --- a/src/meshTools/searchableSurfaces/searchableSurfaces/searchableSurfaces.C +++ b/src/meshTools/searchableSurfaces/searchableSurfaces/searchableSurfaces.C @@ -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) diff --git a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C index d1e23e278e531b6b6070f09004cdff5c7a5db375..2b6e23baf83e3c14513d4102133a9db779a2ac54 100644 --- a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C +++ b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.C @@ -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); } diff --git a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H index d7ce1e543f87244d4532658d7a758d42fa607c24..691c3c73cc9b09226edeafa3e1d831216b9a5d67 100644 --- a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H +++ b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearch.H @@ -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,23 +168,27 @@ public: const bool clip ); - //- Voxel index to voxel centre - static point centre + //- Fill voxels indicated by bounding box + template<class Container, class Type> + static void fill ( + Container& elems, const boundBox& bb, const labelVector& nDivs, - const labelVector& voxel + const boundBox& subBb, + const Type val ); //- Fill voxels indicated by bounding box - template<class Container, class Type> + 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 Type val, + const CombineOp& cop = eqOp<Type>() ); //- Check if any voxel inside bounding box is set to val or diff --git a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearchTemplates.C b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearchTemplates.C index e748d23281d415c0a2cb14ccbed44dd2deb55b3f..adcac1dad3923b9f193537d863a5848b8f9f9a14 100644 --- a/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearchTemplates.C +++ b/src/overset/cellCellStencil/trackingInverseDistance/voxelMeshSearchTemplates.C @@ -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]; } }