diff --git a/applications/utilities/mesh/generation/PDRblockMesh/PDRblockMesh.C b/applications/utilities/mesh/generation/PDRblockMesh/PDRblockMesh.C index 1b795379836547ed04d0784a008273c5775a81d9..cba7c7d968205ec5ac8e44641714a61cf570f903 100644 --- a/applications/utilities/mesh/generation/PDRblockMesh/PDRblockMesh.C +++ b/applications/utilities/mesh/generation/PDRblockMesh/PDRblockMesh.C @@ -57,7 +57,6 @@ Usage #include "Time.H" #include "IOdictionary.H" #include "OSspecific.H" -#include "OFstream.H" using namespace Foam; @@ -97,6 +96,22 @@ int main(int argc, char *argv[]) ); argList::addOptionCompat("no-clean", {"noClean", -2006}); + argList::addBoolOption + ( + "no-outer", + "Create without any other region" + ); + argList::addBoolOption + ( + "print-dict", + "Print blockMeshDict equivalent and exit" + ); + argList::addBoolOption + ( + "write-dict", + "Write system/blockMeshDict.PDRblockMesh and exit" + ); + argList::addOption("dict", "file", "Alternative PDRblockMeshDict"); argList::addOption ( @@ -111,6 +126,9 @@ int main(int argc, char *argv[]) // Remove old files, unless disabled const bool removeOldFiles = !args.found("no-clean"); + // Suppress creation of the outer region + const bool noOuterRegion = args.found("no-outer"); + const word regionName(polyMesh::defaultRegion); const word regionPath; @@ -149,8 +167,38 @@ int main(int argc, char *argv[]) Info<< "Creating PDRblockMesh from " << runTime.relativePath(dictIO.objectPath()) << endl; + // Always start from a PDRblock PDRblock blkMesh(meshDict, true); + if (args.found("print-dict")) + { + Info<< nl << "Equivalent blockMeshDict" << nl << nl; + + blkMesh.blockMeshDict(Info, true); + + Info<< "\nEnd\n" << endl; + return 0; + } + + if (args.found("write-dict")) + { + // Generate system/blockMeshDict and exit + blkMesh.writeBlockMeshDict + ( + IOobject + ( + "blockMeshDict.PDRblockMesh", + runTime.system(), // instance + runTime, // registry + IOobject::NO_READ, + IOobject::NO_WRITE, + false // Do not register + ) + ); + + Info<< "\nEnd\n" << endl; + return 0; + } // Instance for resulting mesh if (useTime) @@ -170,14 +218,19 @@ int main(int argc, char *argv[]) Info<< nl << "Creating polyMesh from PDRblockMesh" << endl; + if (noOuterRegion) + { + Info<< "Outer region disabled, using ijk generation" << nl; + } autoPtr<polyMesh> meshPtr = - blkMesh.mesh - ( - IOobject(regionName, meshInstance, runTime) - ); + ( + args.found("no-outer") + ? blkMesh.innerMesh(IOobject(regionName, meshInstance, runTime)) + : blkMesh.mesh(IOobject(regionName, meshInstance, runTime)) + ); - const polyMesh& mesh = *meshPtr; + polyMesh& mesh = *meshPtr; // Set the precision of the points data to 10 IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision())); diff --git a/applications/utilities/preProcessing/PDRsetFields/Allwclean b/applications/utilities/preProcessing/PDR/Allwclean similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/Allwclean rename to applications/utilities/preProcessing/PDR/Allwclean diff --git a/applications/utilities/preProcessing/PDRsetFields/Allwmake b/applications/utilities/preProcessing/PDR/Allwmake similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/Allwmake rename to applications/utilities/preProcessing/PDR/Allwmake diff --git a/applications/utilities/preProcessing/PDRsetFields/PDRsetFields/Make/files b/applications/utilities/preProcessing/PDR/PDRsetFields/Make/files similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/PDRsetFields/Make/files rename to applications/utilities/preProcessing/PDR/PDRsetFields/Make/files diff --git a/applications/utilities/preProcessing/PDRsetFields/PDRsetFields/Make/options b/applications/utilities/preProcessing/PDR/PDRsetFields/Make/options similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/PDRsetFields/Make/options rename to applications/utilities/preProcessing/PDR/PDRsetFields/Make/options diff --git a/applications/utilities/preProcessing/PDRsetFields/PDRsetFields/PDRsetFields.C b/applications/utilities/preProcessing/PDR/PDRsetFields/PDRsetFields.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/PDRsetFields/PDRsetFields.C rename to applications/utilities/preProcessing/PDR/PDRsetFields/PDRsetFields.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/Make/files b/applications/utilities/preProcessing/PDR/pdrFields/Make/files similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/Make/files rename to applications/utilities/preProcessing/PDR/pdrFields/Make/files diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/Make/options b/applications/utilities/preProcessing/PDR/pdrFields/Make/options similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/Make/options rename to applications/utilities/preProcessing/PDR/pdrFields/Make/options diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarrays.C b/applications/utilities/preProcessing/PDR/pdrFields/PDRarrays.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarrays.C rename to applications/utilities/preProcessing/PDR/pdrFields/PDRarrays.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarrays.H b/applications/utilities/preProcessing/PDR/pdrFields/PDRarrays.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarrays.H rename to applications/utilities/preProcessing/PDR/pdrFields/PDRarrays.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarraysAnalyse.C b/applications/utilities/preProcessing/PDR/pdrFields/PDRarraysAnalyse.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarraysAnalyse.C rename to applications/utilities/preProcessing/PDR/pdrFields/PDRarraysAnalyse.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarraysCalc.C b/applications/utilities/preProcessing/PDR/pdrFields/PDRarraysCalc.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRarraysCalc.C rename to applications/utilities/preProcessing/PDR/pdrFields/PDRarraysCalc.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRlegacy.H b/applications/utilities/preProcessing/PDR/pdrFields/PDRlegacy.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRlegacy.H rename to applications/utilities/preProcessing/PDR/pdrFields/PDRlegacy.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRlegacyMeshSpec.C b/applications/utilities/preProcessing/PDR/pdrFields/PDRlegacyMeshSpec.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRlegacyMeshSpec.C rename to applications/utilities/preProcessing/PDR/pdrFields/PDRlegacyMeshSpec.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.C b/applications/utilities/preProcessing/PDR/pdrFields/PDRmeshArrays.C similarity index 99% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.C rename to applications/utilities/preProcessing/PDR/pdrFields/PDRmeshArrays.C index 230638132291cbafc15377df04ccaf8f880b5f78..d09aa71f19e81c9e9bb3a952f41a0f92d66f558b 100644 --- a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.C +++ b/applications/utilities/preProcessing/PDR/pdrFields/PDRmeshArrays.C @@ -94,7 +94,7 @@ void Foam::PDRmeshArrays::classify << " nPoints:" << pdrBlock.nPoints() << " nCells:" << pdrBlock.nCells() << " nFaces:" << pdrBlock.nFaces() << nl - << " min-edge:" << pdrBlock.minEdgeLen() << nl; + << " min-edge:" << pdrBlock.edgeLimits().min() << nl; Info<< "Classifying ijk indexing... " << nl; diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.H b/applications/utilities/preProcessing/PDR/pdrFields/PDRmeshArrays.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRmeshArrays.H rename to applications/utilities/preProcessing/PDR/pdrFields/PDRmeshArrays.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRparams.C b/applications/utilities/preProcessing/PDR/pdrFields/PDRparams.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRparams.C rename to applications/utilities/preProcessing/PDR/pdrFields/PDRparams.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRparams.H b/applications/utilities/preProcessing/PDR/pdrFields/PDRparams.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRparams.H rename to applications/utilities/preProcessing/PDR/pdrFields/PDRparams.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRpatchDef.C b/applications/utilities/preProcessing/PDR/pdrFields/PDRpatchDef.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRpatchDef.C rename to applications/utilities/preProcessing/PDR/pdrFields/PDRpatchDef.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRpatchDef.H b/applications/utilities/preProcessing/PDR/pdrFields/PDRpatchDef.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRpatchDef.H rename to applications/utilities/preProcessing/PDR/pdrFields/PDRpatchDef.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRsetFields.H b/applications/utilities/preProcessing/PDR/pdrFields/PDRsetFields.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRsetFields.H rename to applications/utilities/preProcessing/PDR/pdrFields/PDRsetFields.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRutils.H b/applications/utilities/preProcessing/PDR/pdrFields/PDRutils.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRutils.H rename to applications/utilities/preProcessing/PDR/pdrFields/PDRutils.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRutilsInternal.H b/applications/utilities/preProcessing/PDR/pdrFields/PDRutilsInternal.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRutilsInternal.H rename to applications/utilities/preProcessing/PDR/pdrFields/PDRutilsInternal.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRutilsIntersect.C b/applications/utilities/preProcessing/PDR/pdrFields/PDRutilsIntersect.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRutilsIntersect.C rename to applications/utilities/preProcessing/PDR/pdrFields/PDRutilsIntersect.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRutilsOverlap.C b/applications/utilities/preProcessing/PDR/pdrFields/PDRutilsOverlap.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/PDRutilsOverlap.C rename to applications/utilities/preProcessing/PDR/pdrFields/PDRutilsOverlap.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/ObstaclesDict b/applications/utilities/preProcessing/PDR/pdrFields/obstacles/ObstaclesDict similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/ObstaclesDict rename to applications/utilities/preProcessing/PDR/pdrFields/obstacles/ObstaclesDict diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacle.C b/applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacle.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacle.C rename to applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacle.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacle.H b/applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacle.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacle.H rename to applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacle.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleI.H b/applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleI.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleI.H rename to applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleI.H diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleIO.C b/applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleIO.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleIO.C rename to applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleIO.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleLegacyIO.C b/applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleLegacyIO.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleLegacyIO.C rename to applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleLegacyIO.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleLegacyRead.C b/applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleLegacyRead.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleLegacyRead.C rename to applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleLegacyRead.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleTypes.C b/applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleTypes.C similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleTypes.C rename to applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleTypes.C diff --git a/applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleTypes.H b/applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleTypes.H similarity index 100% rename from applications/utilities/preProcessing/PDRsetFields/pdrFields/obstacles/PDRobstacleTypes.H rename to applications/utilities/preProcessing/PDR/pdrFields/obstacles/PDRobstacleTypes.H diff --git a/etc/caseDicts/annotated/PDRblockMeshDict b/etc/caseDicts/annotated/PDRblockMeshDict new file mode 100644 index 0000000000000000000000000000000000000000..c3eb54a3774de1c709d10d825228d90de7975fb2 --- /dev/null +++ b/etc/caseDicts/annotated/PDRblockMeshDict @@ -0,0 +1,97 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: v2006 | +| \\ / A nd | Website: www.openfoam.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object PDRblockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Point scaling (optional) +scale 1; + +x +{ + points ( -13.28 -0.10 6.0 19.19 ); + nCells ( 10 12 10 ); + ratios ( -5.16 1 5.16 ); +} + +y +{ + points ( -12.98 0 5.50 18.48 ); + nCells ( 10 11 10 ); + ratios ( -5.16 1 5.16 ); +} + +z +{ + points ( 0.00 4.80 17.26 ); + nCells ( 10 10 ); + ratios ( 1 5.16 ); +} + + +defaultPatch +{ + name walls; + type wall; +} + + +// Faces: 0=x-min, 1=x-max, 2=y-min, 3=y-max, 4=z-min, 5=z-max + +boundary +( + outer + { + type patch; + faces ( 0 1 2 3 5 ); + } + + mergingFaces + { + type wall; + faces (); + } + + blockedFaces + { + type wall; + faces (); + } + + ground + { + type wall; + faces ( 4 ); + } +); + + +// Treatment of the outer region + +outer +{ + type sphere; // (none | extend | box | sphere) [default: none] + onGround true; // Module on the ground? [default: false] + expansion relative; // (uniform | ratio | relative) [default: ratio] + + ratios 1.1; + + size 3; // Overall outer/inner size + nCells 10; // Number of cells for outer region + + // size (3 4); + // nCells (10 12); +} + + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files index 2f88d0e1e3aaa8324a2e29ae3dcc57d24a47ea1d..ab22fd286d51796f14bcefd33f9a2df77ebaea26 100644 --- a/src/mesh/blockMesh/Make/files +++ b/src/mesh/blockMesh/Make/files @@ -40,7 +40,9 @@ blockMesh/blockMeshMergeTopological.C blockMeshTools/blockMeshTools.C PDRblockMesh/PDRblock.C +PDRblockMesh/PDRblockBlockMesh.C PDRblockMesh/PDRblockCreate.C PDRblockMesh/PDRblockLocation.C +PDRblockMesh/PDRblockOuter.C LIB = $(FOAM_LIBBIN)/libblockMesh diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.C b/src/mesh/blockMesh/PDRblockMesh/PDRblock.C index aa6170b92467366b435de88c7e5598f0374fd434..7a97637c5743a536321ca0bd37871976a231a4da 100644 --- a/src/mesh/blockMesh/PDRblockMesh/PDRblock.C +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblock.C @@ -28,6 +28,7 @@ License #include "PDRblock.H" #include "ListOps.H" #include "emptyPolyPatch.H" +#include "gradingDescriptors.H" // Output when verbosity is enabled #undef Log @@ -52,21 +53,22 @@ Foam::PDRblock::expansionNames_ namespace Foam { - //- Calculate geometric expansion factor from expansion ratio - inline scalar calcGexp(const scalar expRatio, const label nDiv) - { - return nDiv > 1 ? pow(expRatio, 1.0/(nDiv - 1)) : 0.0; - } - //- Calculate geometric ratio from relative ratio - inline scalar relativeToGeometricRatio - ( - const scalar expRatio, - const label nDiv - ) - { - return nDiv > 1 ? pow(expRatio, (nDiv - 1)) : 1.0; - } +//- Calculate geometric expansion factor from expansion ratio +inline scalar calcGexp(const scalar expRatio, const label nDiv) +{ + return nDiv > 1 ? pow(expRatio, 1.0/(nDiv - 1)) : 0.0; +} + +//- Calculate geometric ratio from relative ratio +inline scalar relativeToGeometricRatio +( + const scalar expRatio, + const label nDiv +) +{ + return nDiv > 1 ? pow(expRatio, (nDiv - 1)) : 1.0; +} } // End namespace Foam @@ -129,32 +131,20 @@ void Foam::PDRblock::adjustSizes() grid_.z().clear(); bounds_ = boundBox::invertedBox; - minEdgeLen_ = Zero; + edgeLimits_.min() = 0; + edgeLimits_.max() = 0; return; } // Adjust boundBox - bounds_.min().x() = grid_.x().first(); - bounds_.min().y() = grid_.y().first(); - bounds_.min().z() = grid_.z().first(); + bounds_ = bounds(grid_.x(), grid_.y(), grid_.z()); - bounds_.max().x() = grid_.x().last(); - bounds_.max().y() = grid_.y().last(); - bounds_.max().z() = grid_.z().last(); + // Min/max edge lengths + edgeLimits_.clear(); - // Min edge length - minEdgeLen_ = GREAT; - - for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt) - { - const label nEdge = grid_[cmpt].nCells(); - - for (label edgei=0; edgei < nEdge; ++edgei) - { - const scalar len = grid_[cmpt].width(edgei); - minEdgeLen_ = min(minEdgeLen_, len); - } - } + edgeLimits_.add(grid_.x().edgeLimits()); + edgeLimits_.add(grid_.y().edgeLimits()); + edgeLimits_.add(grid_.z().edgeLimits()); } @@ -166,14 +156,18 @@ void Foam::PDRblock::readGridControl expansionType expandType ) { - //- The begin/end nodes for each segment - scalarList knots; + gridControl& ctrl = control_[cmpt]; + + // Begin/end nodes for each segment + scalarList& knots = static_cast<scalarList&>(ctrl); // The number of division per segment - labelList divisions; + labelList& divisions = ctrl.divisions_; // The expansion ratio per segment - scalarList expansion; + scalarList& expansion = ctrl.expansion_; + + expansion.clear(); // expansion is optional Log << "Reading grid control for " << vector::componentNames[cmpt] << " direction" << nl; @@ -514,6 +508,12 @@ void Foam::PDRblock::readBoundary(const dictionary& dict) // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +Foam::PDRblock::PDRblock() +: + PDRblock(dictionary::null, false) +{} + + Foam::PDRblock::PDRblock ( const UList<scalar>& xgrid, @@ -521,7 +521,7 @@ Foam::PDRblock::PDRblock const UList<scalar>& zgrid ) : - PDRblock() + PDRblock(dictionary::null, false) { // Default boundaries with patchi == shapeFacei patches_.resize(6); @@ -541,13 +541,21 @@ Foam::PDRblock::PDRblock } -Foam::PDRblock::PDRblock(const dictionary& dict, bool verbose) +Foam::PDRblock::PDRblock(const dictionary& dict, bool verboseOutput) : - PDRblock() + ijkMesh(), + meshDict_(dict), + grid_(), + outer_(), + bounds_(), + patches_(), + edgeLimits_(0,0), + verbose_(verboseOutput) { - verbose_ = verbose; - - read(dict); + if (&dict != &dictionary::null) + { + read(dict); + } } @@ -575,6 +583,16 @@ bool Foam::PDRblock::read(const dictionary& dict) readBoundary(dict); + // Outer treatment: (none | extend | box | sphere) + outer_.clear(); + + const dictionary* outerDictPtr = dict.findDict("outer"); + if (outerDictPtr) + { + outer_.read(*outerDictPtr); + } + outer_.report(Info); + return true; } @@ -644,7 +662,7 @@ bool Foam::PDRblock::gridIndex const scalar relTol ) const { - const scalar tol = relTol * minEdgeLen_; + const scalar tol = relTol * edgeLimits_.min(); for (direction cmpt=0; cmpt < labelVector::nComponents; ++cmpt) { @@ -688,4 +706,33 @@ Foam::labelVector Foam::PDRblock::gridIndex } +Foam::Vector<Foam::gradingDescriptors> Foam::PDRblock::grading() const +{ + return grading(control_); +} + + +Foam::gradingDescriptors Foam::PDRblock::grading(const direction cmpt) const +{ + switch (cmpt) + { + case vector::X : + case vector::Y : + case vector::Z : + { + return control_[cmpt].grading(); + break; + } + + default : + FatalErrorInFunction + << "Not gridControl for direction " << label(cmpt) << endl + << exit(FatalError); + break; + } + + return gradingDescriptors(); +} + + // ************************************************************************* // diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.H b/src/mesh/blockMesh/PDRblockMesh/PDRblock.H index 94a841e7c4484280e6730f9d8e1a377ad222ea56..0072cd01ea02a84330f2b7ce3697bd3c87f8e70a 100644 --- a/src/mesh/blockMesh/PDRblockMesh/PDRblock.H +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblock.H @@ -73,6 +73,8 @@ SourceFiles #include "pointField.H" #include "faceList.H" #include "Enum.H" +#include "vector2D.H" +#include "labelVector2D.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -81,7 +83,9 @@ namespace Foam // Forward Declarations class IOobject; +class blockMesh; class polyMesh; +class gradingDescriptors; /*---------------------------------------------------------------------------*\ Class PDRblock Declaration @@ -96,7 +100,7 @@ public: // Data Types //- The expansion type - enum expansionType + enum expansionType : uint8_t { EXPAND_UNIFORM, //!< Uniform expansion (ie, no expansion) EXPAND_RATIO, //!< End/start ratio @@ -115,7 +119,7 @@ public: { public: - //- The locations are valid if they contain 2 or more points + //- The location list is valid if it contains 2 or more points inline bool valid() const; //- The number of cells in this direction. @@ -136,6 +140,9 @@ public: //- Mid-point location, zero for an empty list. inline scalar centre() const; + //- The difference between min/max values, zero for an empty list. + inline scalar length() const; + //- Check that element index is within valid range. inline void checkIndex(const label i) const; @@ -146,6 +153,13 @@ public: // Treats -1 and nCells positions like a halo cell. inline scalar C(const label i) const; + //- Return min/max edge lengths + scalarMinMax edgeLimits() const; + + //- Return edge grading descriptors for the locations + // \see Foam::gradingDescriptor + gradingDescriptors grading() const; + //- Find the cell index enclosing this location // \return -1 for out-of-bounds label findCell(const scalar p) const; @@ -172,6 +186,33 @@ private: // Private Classes + //- The begin/end nodes for each segment, + //- with divisions and expansion for each segment + struct gridControl + : + public scalarList + { + //- The number of division per segment + labelList divisions_; + + //- The expansion ratio per segment + scalarList expansion_; + + //- Total number of cells in this direction + label nCells() const; + + //- Return edge grading descriptors for the locations + // \see Foam::gradingDescriptor + gradingDescriptors grading() const; + + //- Add point/divisions/expand to end of list (push_back) + void append(const scalar p, label nDiv, scalar expRatio=1); + + //- Add point/divisions/expand to front of list (push_front) + void prepend(const scalar p, label nDiv, scalar expRatio=1); + }; + + //- Extracted patch settings struct boundaryEntry { @@ -188,23 +229,100 @@ private: labelList faces_; }; + //- The begin/end nodes for each segment, + //- with divisions and expansion for each segment + struct outerControl + { + //- The control type + enum controlType : uint8_t + { + OUTER_NONE = 0, //!< No outer region + OUTER_EXTEND, //!< Extend inner region (orthogonal) + OUTER_BOX, //!< Cuboid + OUTER_SPHERE //!< Spherical + }; + + //- Named enumerations for the control type + const static Enum<controlType> controlNames_; + + //- The control type + controlType type_; + + //- The expansion type + expansionType expandType_; + + //- True if on the ground + bool onGround_; + + //- Relative size(s) for the outer region + vector2D relSize_; + + //- Number of cells in outer region + // Generally only single component is used + labelVector2D nCells_; + + //- Expansion ratio(s) for the outer region + vector2D expansion_; + + + // Constructors + + //- Default construct. NONE + outerControl(); + + + // Member Functions + + //- Reset to default (NONE) values + void clear(); + + //- Is enabled (not NONE) + bool active() const; + + //- Project on to sphere (is SPHERE) + bool isSphere() const; + + //- Is the outer region on the ground? + bool onGround() const; + + //- Define that the outer region is on the ground or not + // \return the old value + bool onGround(const bool on); + + //- Read content from dictionary + void read(const dictionary& dict); + + //- Report information about outer region + void report(Ostream& os) const; + }; + // Private Data - //- Verbosity - bool verbose_; + //- Reference to mesh dictionary + const dictionary& meshDict_; + + //- The grid controls in (i,j,k / x,y,z) directions. + Vector<gridControl> control_; - //- The grid points in all (i,j,k / x,y,z) directions. + //- The grid points in all (i,j,k / x,y,z) directions, + //- after applying the internal subdivisions. Vector<location> grid_; + //- Control for the outer-region (if any) + outerControl outer_; + //- The mesh bounding box boundBox bounds_; //- The boundary patch information PtrList<boundaryEntry> patches_; - //- The min edge length - scalar minEdgeLen_; + //- The min/max edge lengths + scalarMinMax edgeLimits_; + + //- Verbosity + bool verbose_; // Private Member Functions @@ -267,6 +385,35 @@ private: const scalar tol ) const; + //- The bounding box of the grid points + static boundBox bounds + ( + const scalarList& x, //!< X-points, monotonically increasing + const scalarList& y, //!< Y-points, monotonically increasing + const scalarList& z //!< T-points, monotonically increasing + ); + + //- Equivalent edge grading descriptors in (x,y,z) directions. + static Vector<gradingDescriptors> grading + ( + const Vector<gridControl>& ctrl + ); + + //- Mesh sizes based on the controls + static labelVector sizes + ( + const Vector<gridControl>& ctrl + ); + + + // Mesh Generation + + //- Create a blockMesh + autoPtr<blockMesh> createBlockMesh(const IOobject& io) const; + + //- Create polyMesh via blockMesh + autoPtr<polyMesh> meshBlockMesh(const IOobject& io) const; + public: @@ -278,8 +425,8 @@ public: // Constructors - //- Construct zero-size - inline PDRblock(); + //- Default construct, zero-size, inverted bounds etc + PDRblock(); //- Construct from components PDRblock @@ -290,7 +437,7 @@ public: ); //- Construct from dictionary - explicit PDRblock(const dictionary& dict, bool verbose=false); + explicit PDRblock(const dictionary& dict, bool verboseOutput=false); // Member Functions @@ -312,14 +459,23 @@ public: //- The grid point locations in the i,j,k (x,y,z) directions. inline const Vector<location>& grid() const; + //- Equivalent edge grading descriptors in (x,y,z) directions. + Vector<gradingDescriptors> grading() const; + + //- Equivalent edge grading descriptors in specified (x,y,z) direction. + gradingDescriptors grading(const direction cmpt) const; + + + // Mesh Information - // Mesh information + //- Mesh sizing as per ijkMesh + using ijkMesh::sizes; //- The mesh bounding box inline const boundBox& bounds() const; - //- The min edge length - inline const scalar& minEdgeLen() const; + //- The min/max edge length + inline const scalarMinMax& edgeLimits() const; //- Cell size in x-direction at i position. inline scalar dx(const label i) const; @@ -379,17 +535,30 @@ public: labelVector findCell(const point& pt) const; //- Obtain i,j,k grid index for point location within specified - // relative tolerance of the minEdgeLen. + // relative tolerance of the min edge length // The value (-1,-1,-1) is returned for out-of-bounds (not found). // and off-grid labelVector gridIndex(const point& pt, const scalar relTol=0.01) const; - // Mesh generation + // Mesh Generation + + //- Output content for an equivalent blockMeshDict + // Optionally generate header/footer content + Ostream& blockMeshDict(Ostream& os, const bool withHeader=false) const; + + //- Content for an equivalent blockMeshDict + dictionary blockMeshDict() const; + + //- Write an equivalent blockMeshDict + void writeBlockMeshDict(const IOobject& io) const; //- Create polyMesh for grid definition and patch information autoPtr<polyMesh> mesh(const IOobject& io) const; + //- Create polyMesh for inner-mesh only, + //- ignore any outer block definitions + autoPtr<polyMesh> innerMesh(const IOobject& io) const; }; diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockBlockMesh.C b/src/mesh/blockMesh/PDRblockMesh/PDRblockBlockMesh.C new file mode 100644 index 0000000000000000000000000000000000000000..98058904d38c5c55c3c0a96c50295a92bc407feb --- /dev/null +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockBlockMesh.C @@ -0,0 +1,920 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2020 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 "PDRblock.H" +#include "ListOps.H" +#include "cellModeller.H" +#include "gradingDescriptors.H" +#include "objectRegistry.H" +#include "Time.H" +#include "IOdictionary.H" +#include "Fstream.H" +#include "OTstream.H" +#include "edgeHashes.H" + +#include "cellModel.H" +#include "blockMesh.H" +#include "polyMesh.H" +#include "searchableSphere.H" + +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Name for the projection geometry +static const word projKeyword("project"); + +// Name for the projection geometry +static const word projGeomName("sphere"); + + +//- Calculate geometric ratio from relative ratio +inline scalar relativeToGeometricRatio +( + const scalar expRatio, + const label nDiv +) +{ + return nDiv > 1 ? pow(expRatio, (nDiv - 1)) : 1.0; +} + + +// Output space-separated flat list. No size prefix. +template<class T> +static Ostream& outputFlatList(Ostream& os, const UList<T>& list) +{ + os << token::BEGIN_LIST; + label i = 0; + for (const label val : list) + { + if (i++) os << token::SPACE; + os << val; + } + os << token::END_LIST; + + return os; +} + + +// Begin indent list +static inline Ostream& begIndentList(Ostream& os) +{ + os << indent << incrIndent << token::BEGIN_LIST << nl; + return os; +} + +// End indent list +static inline Ostream& endIndentList(Ostream& os) +{ + os << decrIndent << indent << token::END_LIST; + return os; +} + + +// Output list contents (newline separated) indented. +template<class T> +static Ostream& outputIndent(Ostream& os, const UList<T>& list) +{ + for (const T& val : list) + { + os << indent << val << nl; + } + return os; +} + + +// +static Ostream& serializeHex +( + Ostream& os, + const labelUList& hexVerts, + const labelVector& hexCount, + const Vector<gradingDescriptors> hexGrade, + const word& zoneName = word::null +) +{ + os << indent << cellModel::modelNames[cellModel::HEX] << token::SPACE; + outputFlatList(os, hexVerts); + + if (!zoneName.empty()) + { + os << token::SPACE << zoneName; + } + + os << token::SPACE << hexCount << nl + << indent << word("edgeGrading") << nl; + + begIndentList(os); + + // Grading (x/y/z) + for (const gradingDescriptors& gds : hexGrade) + { + begIndentList(os); + outputIndent(os, gds); + endIndentList(os) << nl; + } + + endIndentList(os) << nl; + return os; +} + + +// Generate list with entries: +// +// project (x y z) (geometry) +// +static Ostream& serializeProjectPoints +( + Ostream& os, + const UList<point>& list +) +{ + for (const point& p : list) + { + os << indent << projKeyword << token::SPACE + << p + << token::SPACE + << token::BEGIN_LIST << projGeomName << token::END_LIST << nl; + } + + return os; +} + + +// Generate entry: +// +// project (beg end) (geometry) +// +static Ostream& serializeProjectEdge +( + Ostream& os, + const edge& e +) +{ + os << indent << projKeyword << token::SPACE; + + if (e.sorted()) + { + os << e.first() << token::SPACE << e.second(); + } + else + { + os << e.second() << token::SPACE << e.first(); + } + + os << token::SPACE + << token::BEGIN_LIST << projGeomName << token::END_LIST << nl; + + return os; +} + + +// Generate entry: +// +// (0 1 2 ..) +// +static Ostream& serializeFace +( + Ostream& os, + const face& list +) +{ + os << indent; + outputFlatList(os, list); + os << nl; + return os; +} + + +// Generate entry: +// +// project (0 1 2 ..) geometry +// +static Ostream& serializeProjectFace +( + Ostream& os, + const face& list +) +{ + os << indent << projKeyword << token::SPACE; + outputFlatList(os, list); + os << token::SPACE << projGeomName << nl; + + return os; +} + +} // namespace Foam + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::Ostream& Foam::PDRblock::blockMeshDict +( + Ostream& os, + const bool withHeader +) const +{ + if (withHeader) + { + // Use dummy time for fake objectRegistry + autoPtr<Time> dummyTimePtr(Time::New()); + + IOdictionary iodict + ( + IOobject + ( + "blockMeshDict", + dummyTimePtr->system(), + *dummyTimePtr, + IOobject::NO_READ, + IOobject::NO_WRITE, + false // no register + ) + ); + + iodict.writeHeader(os); + } + + const cellModel& hex = cellModel::ref(cellModel::HEX); + + // The mesh topology will normally be an O-grid with a central (inner) + // block and 6 outer blocks. + // The inner block is described by (0 1 2 3 4 5 6 7). + // The additional points for the outer region: (8 9 19 11 12 13 14 15) + + // The outer blocks will be addressed according to their + // placement w.r.t. the inner block, and defined such that they retain + // the same global orientation as the inner block. + // For example, face 0 of all blocks will be on the logical x-min side + // of the mesh. + + // List of hex vertices, ordered with outer blocks first to allow + // direct addressing by their logical position. + const FixedList<labelList, 8> + hexVerts + ({ + {8, 0, 3, 11, 12, 4, 7, 15}, // x-min block (face 0) + {1, 9, 10, 2, 5, 13, 14, 6}, // x-max block (face 1) + {8, 9, 1, 0, 12, 13, 5, 4}, // y-min block (face 2) + {3, 2, 10, 11, 7, 6, 14, 15}, // y-max block (face 3) + {8, 9, 10, 11, 0, 1, 2, 3}, // z-min block (face 4) + {4, 5, 6, 7, 12, 13, 14, 15}, // z-max block (face 5) + {0, 1, 2, 3, 4, 5, 6, 7}, // Inner box description + {8, 9, 10, 11, 12, 13, 14, 15}, // Outer box description + }); + + // The face or the logical block index (for the O-grid) + enum faceIndex + { + X_Min = 0, + X_Max = 1, + Y_Min = 2, + Y_Max = 3, + Z_Min = 4, + Z_Max = 5, + Inner_Block = 6, // Inner block description + Outer_Block = 7, // Outer bounding box description + }; + + // Lists of block/face for outside and ground faces + DynamicList<labelPair> outerFaces(8); + DynamicList<labelPair> groundFaces(8); + + // We handle a few fixed topology configurations + + enum blockTopologyType + { + INNER_ONLY = 0, // No outer region + EXTENDED = 1, // Outer created by extending inner region + CLIP_BOTTOM = 5, // Outer O-grid on 5 sides + FULL_OUTER = 6 // Outer O-grid on 6 sides + }; + + + + // Expansion ratios need conversion from relative to geometric + const bool useRelToGeom = + (expansionType::EXPAND_RATIO == outer_.expandType_); + + + // Physical dimensions + + Vector<gridControl> ctrl(control_); + boundBox innerCorners(bounds(ctrl.x(), ctrl.y(), ctrl.z())); + + boundBox outerCorners; + + + point radialCentre(innerCorners.centre()); + vector radialSizes(0.5*innerCorners.span()); + + blockTopologyType outerTopology = INNER_ONLY; + + + if (outer_.active()) + { + outerTopology = FULL_OUTER; + + // Convert from relative size + radialSizes.x() *= outer_.relSize_.x(); + radialSizes.y() *= outer_.relSize_.y(); + radialSizes.z() *= min(outer_.relSize_.x(), outer_.relSize_.y()); + + if (outer_.onGround()) + { + outerTopology = CLIP_BOTTOM; + radialCentre.z() = innerCorners.min().z(); + radialSizes.z() *= 2; + } + + // Corners for a box + outerCorners.min() = radialCentre - radialSizes; + outerCorners.max() = radialCentre + radialSizes; + + if (outer_.onGround()) + { + outerCorners.min().z() = innerCorners.min().z(); + } + + + if (outer_.isSphere()) + { + // For spheroid projection, don't trust that blockMesh does it + // properly. Give some reasonable estimates of the corners. + + // Use dummy Time for objectRegistry + autoPtr<Time> dummyTimePtr(Time::New()); + + const IOobject io + ( + "sphere", + *dummyTimePtr, + IOobject::NO_READ, + IOobject::NO_WRITE, + false // do not register + ); + + searchableSphere sphere(io, radialCentre, radialSizes); + + pointField queries(2); + queries[0] = outerCorners.min(); + queries[1] = outerCorners.max(); + + List<pointIndexHit> hits; + sphere.findNearest + ( + queries, + scalarField(2, GREAT), + hits + ); + + outerCorners.min() = hits[0].hitPoint(); + outerCorners.max() = hits[1].hitPoint(); + } + else if (outerControl::OUTER_EXTEND == outer_.type_) + { + outerTopology = EXTENDED; + + // Extend the inner block + label outerCount; + scalar expRatio; + + outerCount = outer_.nCells_.x(); + expRatio = outer_.expansion_.x(); + if (useRelToGeom) + { + expRatio = relativeToGeometricRatio(expRatio, outerCount); + } + + ctrl.x().prepend(outerCorners.min().x(), outerCount, -expRatio); + ctrl.x().append(outerCorners.max().x(), outerCount, expRatio); + + + outerCount = outer_.nCells_.y(); + expRatio = outer_.expansion_.y(); + if (useRelToGeom) + { + expRatio = relativeToGeometricRatio(expRatio, outerCount); + } + + ctrl.y().prepend(outerCorners.min().y(), outerCount, -expRatio); + ctrl.y().append(outerCorners.max().y(), outerCount, expRatio); + + outerCount = max(outer_.nCells_.x(), outer_.nCells_.y()); + expRatio = min(outer_.expansion_.x(), outer_.expansion_.y()); + if (useRelToGeom) + { + expRatio = relativeToGeometricRatio(expRatio, outerCount); + } + + if (!outer_.onGround()) + { + ctrl.z().prepend(outerCorners.min().z(), outerCount, -expRatio); + } + ctrl.z().append(outerCorners.max().z(), outerCount, expRatio); + + // Update corners + innerCorners = bounds(ctrl.x(), ctrl.y(), ctrl.z()); + outerCorners = innerCorners; + } + } + + + const Vector<gradingDescriptors> innerGrading(grading(ctrl)); + const labelVector innerCount(sizes(ctrl)); + + labelVector hexCount; + Vector<gradingDescriptors> hexGrade; + + + const label radialCount = outer_.nCells_.x(); + scalar expRatio = outer_.expansion_.x(); + + if (useRelToGeom) + { + expRatio = relativeToGeometricRatio(expRatio, radialCount); + } + + const gradingDescriptors radialInward + ( + gradingDescriptor{-expRatio} + ); + + const gradingDescriptors radialOutward + ( + gradingDescriptor{expRatio} + ); + + + if (EXTENDED == outerTopology) + { + // The inner block is extended to become the outer faces + outerFaces.append + ({ + labelPair(Inner_Block, X_Min), + labelPair(Inner_Block, X_Max), + labelPair(Inner_Block, Y_Min), + labelPair(Inner_Block, Y_Max), + labelPair(Inner_Block, Z_Max) + }); + + // The ground faces vs outside faces + if (outer_.onGround()) + { + groundFaces.append + ( + labelPair(Inner_Block, Z_Min) + ); + } + else + { + outerFaces.append + ( + labelPair(Inner_Block, Z_Min) + ); + } + } + else if (CLIP_BOTTOM == outerTopology || FULL_OUTER == outerTopology) + { + // The outside faces + outerFaces.append + ({ + labelPair(X_Min, X_Min), + labelPair(X_Max, X_Max), + labelPair(Y_Min, Y_Min), + labelPair(Y_Max, Y_Max), + labelPair(Z_Max, Z_Max) + }); + + // The ground faces + if (CLIP_BOTTOM == outerTopology) + { + groundFaces.append + ({ + labelPair(X_Min, Z_Min), + labelPair(X_Max, Z_Min), + labelPair(Y_Min, Z_Min), + labelPair(Y_Max, Z_Min), + // Note: {Z_Min, Z_Min} will not exist + labelPair(Inner_Block, Z_Min) + }); + } + else + { + outerFaces.append + ( + labelPair(Z_Min, Z_Min) + ); + } + } + + + if (outer_.isSphere()) + { + os.beginBlock("geometry"); + { + os.beginBlock(projGeomName); + { + os.writeEntry("type", "sphere"); + os.writeEntry("origin", radialCentre); + os.writeEntry("radius", radialSizes); + } + os.endBlock(); + } + os.endBlock(); + } + + // vertices + { + os << nl << word("vertices") << nl; + begIndentList(os); + + pointField corners(innerCorners.points()); + + // inner + outputIndent(os, corners); + + // outer + if (CLIP_BOTTOM == outerTopology || FULL_OUTER == outerTopology) + { + corners = outerCorners.points(); + + if (outer_.isSphere()) + { + serializeProjectPoints(os, corners); + } + else + { + outputIndent(os, corners); + } + } + + endIndentList(os) << token::END_STATEMENT << nl; + } + + + // blocks + { + word innerZoneName = "inner"; + if (INNER_ONLY == outerTopology || EXTENDED == outerTopology) + { + innerZoneName.clear(); + } + + os << nl << word("blocks") << nl; + begIndentList(os); + + // Inner block + hexCount = innerCount; + + serializeHex + ( + os, + hexVerts[Inner_Block], + hexCount, + innerGrading, + innerZoneName + ); + + // outer + if (CLIP_BOTTOM == outerTopology || FULL_OUTER == outerTopology) + { + // Radial direction = X + hexCount = innerCount; + hexGrade = innerGrading; + + hexCount.x() = radialCount; + + // Face 0: x-min + { + hexGrade.x() = radialInward; + serializeHex(os, hexVerts[X_Min], hexCount, hexGrade); + } + // Face 1: x-max + { + hexGrade.x() = radialOutward; + serializeHex(os, hexVerts[X_Max], hexCount, hexGrade); + } + + + // Radial direction = Y + hexCount = innerCount; + hexGrade = innerGrading; + + hexCount.y() = radialCount; + + // Face 2: y-min + { + hexGrade.y() = radialInward; + serializeHex(os, hexVerts[Y_Min], hexCount, hexGrade); + } + // Face 3: y-max + { + hexGrade.y() = radialOutward; + serializeHex(os, hexVerts[Y_Max], hexCount, hexGrade); + } + + + // Radial direction = Z + hexCount = innerCount; + hexGrade = innerGrading; + + hexCount.z() = radialCount; + + // Face 4: z-min + if (!outer_.onGround()) + { + hexGrade.z() = radialInward; + serializeHex(os, hexVerts[Z_Min], hexCount, hexGrade); + } + // Face 5: z-max + { + hexGrade.z() = radialOutward; + serializeHex(os, hexVerts[Z_Max], hexCount, hexGrade); + } + } + + endIndentList(os) << token::END_STATEMENT << nl; + } + + + // edges + { + os << nl << word("edges") << nl; + begIndentList(os); + + if (outer_.isSphere() && outerFaces.size()) + { + // Edges for the outer face of the block + edgeHashSet projEdges(32); + + for (const labelPair& pr : outerFaces) + { + projEdges.insert + ( + hex.face(pr.second(), hexVerts[pr.first()]).edges() + ); + } + + for (const edge& e : projEdges.sortedToc()) + { + serializeProjectEdge(os, e); + } + } + + endIndentList(os) << token::END_STATEMENT << nl; + } + + + // faces + { + os << nl << word("faces") << nl; + begIndentList(os); + + if (outer_.isSphere() && outerFaces.size()) + { + for (const labelPair& pr : outerFaces) + { + serializeProjectFace + ( + os, + hex.face(pr.second(), hexVerts[pr.first()]) + ); + } + } + + endIndentList(os) << token::END_STATEMENT << nl; + } + + + // boundary + { + os << nl << word("boundary") << nl; + begIndentList(os); + + // outer + { + os.beginBlock("outer"); + os.writeEntry("type", word("patch")); + + os << indent << word("faces") << nl; + begIndentList(os); + + for (const labelPair& pr : outerFaces) + { + serializeFace + ( + os, + hex.face(pr.second(), hexVerts[pr.first()]) + ); + } + + endIndentList(os) << token::END_STATEMENT << nl; + + os.endBlock(); + } + + if (outer_.onGround()) + { + os.beginBlock("ground"); + os.writeEntry("type", word("wall")); + + os << indent << word("faces") << nl; + begIndentList(os); + + for (const labelPair& pr : groundFaces) + { + serializeFace + ( + os, + hex.face(pr.second(), hexVerts[pr.first()]) + ); + } + + endIndentList(os) << token::END_STATEMENT << nl; + os.endBlock(); + } + + endIndentList(os) << token::END_STATEMENT << nl; + } + + + if (withHeader) + { + IOobject::writeEndDivider(os); + } + + return os; +} + + +Foam::dictionary Foam::PDRblock::blockMeshDict() const +{ + OTstream os; + blockMeshDict(os); + + ITstream is("blockMeshDict", tokenList()); + is.transfer(os.tokens()); + + return dictionary(is); +} + + +void Foam::PDRblock::writeBlockMeshDict(const IOobject& io) const +{ + IOdictionary iodict + ( + IOobject + ( + io.name(), + io.db().time().system(), + io.local(), + io.db(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false // no register + ) + ); + + OFstream os(iodict.objectPath()); + + Info<< nl + << "Generate blockMeshDict: " + << iodict.db().time().relativePath(os.name()) << endl; + + // Set precision for points to 10 + os.precision(max(10u, IOstream::defaultPrecision())); + + iodict.writeHeader(os); + + // Just like writeData, but without copying beforehand + this->blockMeshDict(os); + + iodict.writeEndDivider(os); +} + + +Foam::autoPtr<Foam::blockMesh> +Foam::PDRblock::createBlockMesh(const IOobject& io) const +{ + IOdictionary iodict + ( + IOobject + ( + "blockMeshDict.PDRblockMesh", + io.db().time().system(), + io.local(), + io.db(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false // no register + ), + blockMeshDict() + ); + + return autoPtr<blockMesh>::New(iodict); +} + + +Foam::autoPtr<Foam::polyMesh> +Foam::PDRblock::meshBlockMesh(const IOobject& io) const +{ + const bool oldVerbose = blockMesh::verboseOutput; + blockMesh::verboseOutput = false; + + autoPtr<polyMesh> meshPtr(createBlockMesh(io)->mesh(io)); + + blockMesh::verboseOutput = oldVerbose; + + // This is a bit ugly. + // For extend, we still wish to have an 'inner' cellZone, + // but we meshed the entirety. + + if + ( + outerControl::OUTER_EXTEND == outer_.type_ + && meshPtr->cellZones().empty() + ) + { + const boundBox innerBox + ( + bounds(control_.x(), control_.y(), control_.z()) + ); + + const label nZoneCellsMax = + ( + control_.x().nCells() + * control_.y().nCells() + * control_.z().nCells() + ); + + + polyMesh& pmesh = *meshPtr; + + List<cellZone*> cz(1); + cz[0] = new cellZone + ( + "inner", + labelList(nZoneCellsMax), + 0, // zonei + pmesh.cellZones() + ); + + cellZone& innerZone = *(cz[0]); + + const vectorField& cc = pmesh.cellCentres(); + + label nZoneCells = 0; + + for + ( + label celli = 0; + celli < cc.size() && nZoneCells < nZoneCellsMax; + ++celli + ) + { + if (innerBox.contains(cc[celli])) + { + innerZone[nZoneCells] = celli; + ++nZoneCells; + } + } + + innerZone.resize(nZoneCells); + + pmesh.pointZones().clear(); + pmesh.faceZones().clear(); + pmesh.cellZones().clear(); + pmesh.addZones(List<pointZone*>(), List<faceZone*>(), cz); + } + + return meshPtr; +} + + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C b/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C index 61f07e7582d8ee56d3a8455e65090b66efd6def7..df1e3c9b5efa831ae7923f539ad7d4bfc730cbbd 100644 --- a/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockCreate.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -300,7 +300,8 @@ Foam::label Foam::PDRblock::addBoundaryFaces } -Foam::autoPtr<Foam::polyMesh> Foam::PDRblock::mesh(const IOobject& io) const +Foam::autoPtr<Foam::polyMesh> +Foam::PDRblock::innerMesh(const IOobject& io) const { pointField pts(nPoints()); @@ -376,4 +377,20 @@ Foam::autoPtr<Foam::polyMesh> Foam::PDRblock::mesh(const IOobject& io) const } +Foam::autoPtr<Foam::polyMesh> +Foam::PDRblock::mesh(const IOobject& io) const +{ + if (outer_.active()) + { + Info<< "Outer region is active, using blockMesh generation" << nl; + return meshBlockMesh(io); + } + else + { + Info<< "Outer region is inactive, using ijk generation" << nl; + return innerMesh(io); + } +} + + // ************************************************************************* // diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H b/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H index 42ada92b255fdf3267501d88d52dfb6e6873a322..d05c279dd6c51f1d0b60bf394be1ce4e9ab02ad0 100644 --- a/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -25,20 +25,8 @@ License \*---------------------------------------------------------------------------*/ -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -inline Foam::PDRblock::PDRblock() -: - ijkMesh(), - verbose_(false), - grid_(), - bounds_(), - patches_(), - minEdgeLen_(Zero) -{} - - // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +// Location inline bool Foam::PDRblock::location::valid() const { @@ -82,6 +70,12 @@ inline Foam::scalar Foam::PDRblock::location::centre() const } +inline Foam::scalar Foam::PDRblock::location::length() const +{ + return scalarList::empty() ? 0 : mag(last() - first()); +} + + inline void Foam::PDRblock::location::checkIndex(const label i) const { if (i < 0 || i >= nCells()) @@ -148,6 +142,8 @@ Foam::PDRblock::location::clip(const scalar& val) const } +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + inline const Foam::Vector<Foam::PDRblock::location>& Foam::PDRblock::grid() const { @@ -155,9 +151,9 @@ Foam::PDRblock::grid() const } -inline const Foam::scalar& Foam::PDRblock::minEdgeLen() const +inline const Foam::scalarMinMax& Foam::PDRblock::edgeLimits() const { - return minEdgeLen_; + return edgeLimits_; } diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C b/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C index 1f8311ac1acba42f0cd8311296e7aa74455017fe..6003c08a798a8034dc0d411afa5ac79f27546305 100644 --- a/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2019 OpenCFD Ltd. + Copyright (C) 2019-2020 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -26,13 +26,211 @@ License \*---------------------------------------------------------------------------*/ #include "PDRblock.H" -#include "ListOps.H" +#include "gradingDescriptors.H" + +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Prepend a value by shifting contents +template<class T> +static void prependList(List<T>& list, const T& val) +{ + const label oldLen = list.size(); + list.resize(oldLen + 1); + + for (label i = oldLen; i > 0; --i) + { + list[i] = std::move(list[i-1]); + } + + list[0] = val; +} + +} // End namespace Foam + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + +Foam::boundBox Foam::PDRblock::bounds +( + const scalarList& x, + const scalarList& y, + const scalarList& z +) +{ + return boundBox + ( + point(x.first(), y.first(), z.first()), + point(x.last(), y.last(), z.last()) + ); +} + + +Foam::Vector<Foam::gradingDescriptors> +Foam::PDRblock::grading(const Vector<gridControl>& ctrl) +{ + return Vector<gradingDescriptors> + ( + ctrl.x().grading(), + ctrl.y().grading(), + ctrl.z().grading() + ); +} + +Foam::labelVector +Foam::PDRblock::sizes(const Vector<gridControl>& ctrl) +{ + return labelVector + ( + ctrl.x().nCells(), + ctrl.y().nCells(), + ctrl.z().nCells() + ); +} + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +Foam::label Foam::PDRblock::gridControl::nCells() const +{ + label nTotal = 0; + for (const label nDiv : divisions_) + { + nTotal += nDiv; + } + + return nTotal; +} + + +Foam::gradingDescriptors Foam::PDRblock::gridControl::grading() const +{ + // Begin/end nodes for each segment + const scalarList& knots = *this; + + gradingDescriptors gds(divisions_.size()); + + forAll(gds, i) + { + //- Construct from components + gds[i] = gradingDescriptor + ( + knots[i+1] - knots[i], // blockFraction from delta + divisions_[i], // nDivFraction from nDivs + expansion_[i] + ); + } + + gds.normalise(); + + return gds; +} + + +void Foam::PDRblock::gridControl::append +( + const scalar p, + const label nDiv, + scalar expRatio +) +{ + // Begin/end nodes for each segment + scalarList& knots = *this; + + // Is monotonic? + if (knots.size() && (p <= knots.last())) + { + WarningInFunction + << "Cannot append point " << p + << " which is <= last value " << knots.last() << endl; + return; + } + + if (nDiv < 1) + { + WarningInFunction + << "Negative or zero divisions " << nDiv << endl; + return; + } + + // Correct expansion ratios - negative is the same as inverse. + if (expRatio < 0) + { + expRatio = 1.0/(-expRatio); + } + else if (equal(expRatio, 0)) + { + expRatio = 1; + } + + // Now append (push_back) + knots.append(p); + divisions_.append(nDiv); + expansion_.append(expRatio); +} + + +void Foam::PDRblock::gridControl::prepend +( + const scalar p, + const label nDiv, + scalar expRatio +) +{ + // Begin/end nodes for each segment + scalarList& knots = static_cast<scalarList&>(*this); + + // Is monotonic? + if (knots.size() && (p >= knots.first())) + { + WarningInFunction + << "Cannot prepend point " << p + << " which is >= first value " << knots.first() << endl; + return; + } + + if (nDiv < 1) + { + WarningInFunction + << "Negative or zero divisions " << nDiv << endl; + return; + } + + // Correct expansion ratios - negative is the same as inverse. + if (expRatio < 0) + { + expRatio = 1.0/(-expRatio); + } + else if (equal(expRatio, 0)) + { + expRatio = 1; + } + + // Now prepend (push_front) + prependList(knots, p); + prependList(divisions_, nDiv); + prependList(expansion_, expRatio); +} + + +Foam::scalarMinMax Foam::PDRblock::location::edgeLimits() const +{ + scalarMinMax limits; + + for (label edgei = 0; edgei < this->nCells(); ++edgei) + { + limits.add(width(edgei)); + } + + return limits; +} + + Foam::label Foam::PDRblock::location::findCell(const scalar p) const { - if (scalarList::empty()) + if (scalarList::empty() || p < first() || p > last()) { return -1; } diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockOuter.C b/src/mesh/blockMesh/PDRblockMesh/PDRblockOuter.C new file mode 100644 index 0000000000000000000000000000000000000000..46692a28173324e9480199c8d8c8aa26e500bd57 --- /dev/null +++ b/src/mesh/blockMesh/PDRblockMesh/PDRblockOuter.C @@ -0,0 +1,245 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2020 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 "PDRblock.H" +#include "dictionary.H" +#include "Switch.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +const Foam::Enum +< + Foam::PDRblock::outerControl::controlType +> +Foam::PDRblock::outerControl::controlNames_ +({ + { controlType::OUTER_NONE, "none" }, + { controlType::OUTER_EXTEND, "extend" }, + { controlType::OUTER_BOX, "box" }, + { controlType::OUTER_SPHERE, "sphere" }, +}); + + +// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Get a single or a pair of values +template<class T> +static Vector2D<T> getLazyPair(const word& name, const dictionary& dict) +{ + if (token(dict.lookup(name)).isNumber()) + { + return Vector2D<T>::uniform(dict.get<T>(name)); + } + + return dict.get<Vector2D<T>>(name); +} + +} // End namespace Foam + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::PDRblock::outerControl::outerControl() +: + type_(controlType::OUTER_NONE), + expandType_(expansionType::EXPAND_RATIO), + onGround_(false), + relSize_(0,0), + nCells_(0,0), + expansion_(1,1) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::PDRblock::outerControl::clear() +{ + type_ = controlType::OUTER_NONE; + expandType_ = expansionType::EXPAND_RATIO; + onGround_ = false; + relSize_ = Zero; + nCells_ = Zero; + expansion_ = vector2D::uniform(1); +} + + +void Foam::PDRblock::outerControl::report(Ostream& os) const +{ + if (active()) + { + os << "Has outer region: " << controlNames_[type_] << nl + << " onGround : " << Switch::name(onGround_) << nl + << " sizes : " << relSize_ << nl + << " nCells : " << nCells_ << nl; + } + else + { + os << "No outer region" << nl; + } +} + + +bool Foam::PDRblock::outerControl::active() const +{ + return (controlType::OUTER_NONE != type_); +} + + +bool Foam::PDRblock::outerControl::isSphere() const +{ + return (controlType::OUTER_SPHERE == type_); +} + + +bool Foam::PDRblock::outerControl::onGround() const +{ + return onGround_; +} + + +bool Foam::PDRblock::outerControl::onGround(const bool on) +{ + bool old(onGround_); + onGround_ = on; + return old; +} + + +void Foam::PDRblock::outerControl::read(const dictionary& dict) +{ + clear(); + + type_ = controlNames_.getOrDefault("type", dict, controlType::OUTER_NONE); + onGround_ = dict.getOrDefault("onGround", false); + + if (controlType::OUTER_NONE == type_) + { + return; + } + + // Everything else + + nCells_ = getLazyPair<label>("nCells", dict); + relSize_ = getLazyPair<scalar>("size", dict); + + expandType_ = + expansionNames_.getOrDefault + ( + "expansion", + dict, + expansionType::EXPAND_RATIO + ); + + + if (dict.found("ratios")) + { + expansion_ = getLazyPair<scalar>("ratios", dict); + } + else + { + if (expandType_ != expansionType::EXPAND_UNIFORM) + { + expandType_ = expansionType::EXPAND_UNIFORM; + // Info << "Warning: no 'ratios', use uniform spacing" << nl; + } + } + + if (expandType_ == expansionType::EXPAND_UNIFORM) + { + expansion_ = vector2D::uniform(1); + } + + + // Errors + int die = 0; + + if (nCells_.x() <= 1 || nCells_.y() <= 1) + { + if (!die++) + { + FatalIOErrorInFunction(dict); + } + FatalIOError + << "Too few outer cells: " << nCells_ << nl; + } + + if (relSize_.x() <= 1 || relSize_.y() <= 1) + { + if (!die++) + { + FatalIOErrorInFunction(dict); + } + FatalIOError + << "Outer dimensions must be > 1. Had " << relSize_ << nl; + } + + if (die) + { + FatalIOError << nl << exit(FatalIOError); + } + + + // Warnings + + if + ( + controlType::OUTER_BOX == type_ + || controlType::OUTER_SPHERE == type_ + ) + { + if (relSize_.x() < 2 || relSize_.y() < 2) + { + WarningInFunction + << "Outer dimensions " + << relSize_ << " too small for " + << controlNames_[type_] << " - switching to " + << controlNames_[controlType::OUTER_EXTEND] << nl; + + type_ = controlType::OUTER_EXTEND; + } + } + + if (controlType::OUTER_SPHERE == type_) + { + if (relSize_.x() < 3 || relSize_.y() < 3) + { + WarningInFunction + << "Outer dimensions " + << relSize_ << " too small for " + << controlNames_[type_] << " - switching to " + << controlNames_[controlType::OUTER_BOX] << nl; + + type_ = controlType::OUTER_EXTEND; + } + } +} + + +// ************************************************************************* // diff --git a/tutorials/mesh/PDRblockMesh/box0/system/PDRblockMeshDict b/tutorials/mesh/PDRblockMesh/box0/system/PDRblockMeshDict index eb8065dfc94df3d2e15ae0190d222f7ab73efb5c..c3eb54a3774de1c709d10d825228d90de7975fb2 100644 --- a/tutorials/mesh/PDRblockMesh/box0/system/PDRblockMeshDict +++ b/tutorials/mesh/PDRblockMesh/box0/system/PDRblockMeshDict @@ -75,4 +75,23 @@ boundary } ); + +// Treatment of the outer region + +outer +{ + type sphere; // (none | extend | box | sphere) [default: none] + onGround true; // Module on the ground? [default: false] + expansion relative; // (uniform | ratio | relative) [default: ratio] + + ratios 1.1; + + size 3; // Overall outer/inner size + nCells 10; // Number of cells for outer region + + // size (3 4); + // nCells (10 12); +} + + // ************************************************************************* // diff --git a/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRblockMeshDict b/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRblockMeshDict index 7232300ed7af4f7bb15c8f9cac0e89b729df373d..8fd2622338be981c2ec738a061bf06024ced9369 100644 --- a/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRblockMeshDict +++ b/tutorials/preProcessing/PDRsetFields/simplePipeCage/system/PDRblockMeshDict @@ -74,4 +74,23 @@ boundary } ); + +// Treatment of the outer region + +outer +{ + type sphere; // (none | extend | box | sphere) [default: none] + onGround true; // Module on the ground? [default: false] + expansion relative; // (uniform | ratio | relative) [default: ratio] + + ratios 1.1; + + size 3; // Overall outer/inner size + nCells 10; // Number of cells for outer region + + // size (3 4); + // nCells (10 12); +} + + //***************************************************************************//