diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C index afdba50767..5fc097dcc7 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C @@ -34,9 +34,6 @@ Description #include "Time.H" #include "fvMesh.H" #include "autoHexMeshDriver.H" -#include "pointMesh.H" -#include "motionSmoother.H" -#include "mapDistributePolyMesh.H" using namespace Foam; @@ -52,6 +49,18 @@ int main(int argc, char *argv[]) Info<< "Read mesh in = " << runTime.cpuTimeIncrement() << " s" << endl; + // Read decomposePar dictionary + IOdictionary decomposeDict + ( + IOobject + ( + "decomposeParDict", + runTime.system(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); // Read meshing dictionary IOdictionary meshDict @@ -66,24 +75,46 @@ int main(int argc, char *argv[]) ) ); - // Read decomposePar dictionary - IOdictionary decomposeDict + // refinement parameters + const dictionary& refineDict = meshDict.subDict("refineDict"); + + // snap-to-surface parameters + const dictionary& snapDict = meshDict.subDict("snapDict"); + + // mesh motion and mesh quality parameters + const dictionary& motionDict = meshDict.subDict("motionDict"); + + // layer addition parameters + const dictionary& layerDict = meshDict.subDict("layerDict"); + + + // Main meshing driver. Read surfaces. Determine initial intersections. + autoHexMeshDriver meshEngine ( - IOobject - ( - "decomposeParDict", - runTime.system(), - mesh, - IOobject::MUST_READ, - IOobject::NO_WRITE - ) + mesh, + meshDict, // global control parameters + refineDict, // refinement parameters + decomposeDict ); - // Main meshing driver. Read surfaces. Determine intersections. - autoHexMeshDriver meshEngine(mesh, meshDict, decomposeDict); + Switch wantRefine(meshDict.lookup("doRefine")); + Switch wantSnap(meshDict.lookup("doSnap")); + Switch wantLayers(meshDict.lookup("doLayers")); - // Do all: refine, snap, add layers - meshEngine.doMesh(); + if (wantRefine) + { + meshEngine.doRefine(refineDict, wantSnap); + } + + if (wantSnap) + { + meshEngine.doSnap(snapDict, motionDict); + } + + if (wantLayers) + { + meshEngine.doLayers(layerDict, motionDict); + } Info<< "Finished meshing in = " << runTime.elapsedCpuTime() << " s." << endl; diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict new file mode 100644 index 0000000000..e9aa8f189d --- /dev/null +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -0,0 +1,314 @@ +/*---------------------------------------------------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: 1.0 | +| \\ / A nd | Web: http://www.openfoam.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ + +FoamFile +{ + version 2.0; + format ascii; + + root "/home/penfold/mattijs/foam/mattijs2.1/run/icoFoam"; + case "cavity"; + instance "system"; + local ""; + + class dictionary; + object autoHexMeshDict; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// Which phases to run. +doRefine true; +doSnap true; +doLayers true; // includes autoMergeFaces + + +// Whether to dump intermediate meshes and print lots +// 1 : write mesh +// 2 : write volScalarField with cellLevel for postprocessing +// 4 : write current intersections as .obj files +debug 0; + +refineDict +{ + // Which part to keep. + // NOTE: This point should never be on a face, always inside a cell, even + // after refinement. + keepPoints ((3 0.28 0.43)); + + // Whether to remove/split cells likely to give problems when snapping + handleSnapProblems on; + + // Merge tolerance. Is fraction of overall bounding box of initial mesh + mergeTolerance 1E-6; + + // While refining maximum number of cells per processor. This is basically + // the number of cells that fit on a processor. If you choose this too small + // it will do just more refinement iterations to obtain a similar mesh. + procCellLimit 1000000; + + + // Overall cell limit (approximately). Refinement will stop immediately + // upon reaching this number so a refinement level might not complete. + // Note that this is the number of cells before removing the part which + // is not 'visible' from the keepPoint. The final number of cells might actually + // be a lot less. + cellLimit 2000000; + + + // The surface refinement loop might spend lots of iterations refining just a + // few cells. This setting will cause refinement to stop if <= minimumRefine + // are selected for refinement. Note: it will at least do one iteration + // (unless the number of cells to refine is 0) + minimumRefine 0; + + + + + // Number of buffer layers between different levels. + // 1 means normal 2:1 refinement restriction, larger means slower + // refinement. + nBufferLayers 1; + + + // Feature Edge Refinement + // ~~~~~~~~~~~~~~~~~~~~~~~ + + // External feature file. Read from constant/triSurface for now. + // Limitations: + // - either start or edge of any feature line has to be inside domain + // and be visible from keepPoint on starting mesh. + // - refining cells containing features is separate phase before refining + // based on surface. + // - no load balancing, no check for cellLimit is done while doing this. + features + ( + // { + // file "someLine.eMesh"; + // level 2; + // } + ); + + + // Internal Mesh Refinement + // ~~~~~~~~~~~~~~~~~~~~~~~~ + + // Specifies the areas where the refinement has to be a certain level. + // These surfaces have to be closed. Refinement is either inside + // (refineInside = true) or outside. (note:insideness or outsideness + // is with respect to far away) + // Note:that even using these the transition can never be faster than + // nBufferLayers! + // Note:the refinement level can never be higher than any of the surfaces + // they overlap with. See below for the surface refinement level specification. + refinementShells + ( + { + //type triSurfaceMesh; + //name cube1x1x1.stl; + type searchableBox; + name box1x1x1; + min (2.5 -0.5 -0.5); + max (3.5 0.5 0.5); + + level 4; + refineInside true; + } + ); + + + // Surface based refinement + // ~~~~~~~~~~~~~~~~~~~~~~~~ + + // Curvature. Cosine of angle between two neighbouring surface intersections. + // Only used if cell level > minLevel and < maxLevel. + curvature 0.5; + + + // Overall the refinement is according to the following rules: + // - if the cell-cell vector intersects a surface any cell that + // is less refined than the minRefinementLevel of the surface gets refined. + // - else if the refinement level of the cell is between the + // minRefinementLevel and maxRefinementLevel the cell gets refined if + // - the normal of neighbouring surface intersections differ by more + // than above curvature + // - or if neighbouring surface intersections are on different surfaces or + // different surface regions. + + + // surfaces + surfaces + ( + { + name sphere; + file "sphere.stl"; + + // Surface wide refinement level + minRefinementLevel 1; + maxRefinementLevel 1; + + // Layers + surfaceLayers 1; + + // Region specific refinement level + regions + ( + { + name firstSolid; + minRefinementLevel 3; + maxRefinementLevel 3; + surfaceLayers 2; + } + { + name secondSolid; + minRefinementLevel 1; + maxRefinementLevel 1; + surfaceLayers 1; + } + ); + } + ); +} + +// For snapping +snapDict +{ + //- Number of patch smoothing iterations before finding correspondence + // to surface + nSmoothPatch 3; + + //- Relative distance for points to be attracted by surface feature point + // or edge. True distance is this factor times local + // maximum edge length. + snapTol 4.0; + + //- Whether to move internal mesh as well as boundary + smoothMesh true; + + //- Number of mesh displacement smoothing iterations. + nSmoothDispl 30; + + //- Maximum number of snapping relaxation iterations. Should stop + // before upon reaching a correct mesh. + nSnap 5; +} + + +// For cell layers +layerDict +{ + //- When not to extrude surface. 0 is flat surface, 90 is when two faces + // make straight angle. + featureAngle 60; + + //- Maximum number of snapping relaxation iterations. Should stop + // before upon reaching a correct mesh. + nSnap 5; + + + //- Minimum thickness of cell layer. If for any reason layer cannot be + // above minThickness do not add layer if thickness below minThickNess. + // Relative to undistorted cell size + minThickness 0.25; + + //- If points get not extruded do nGrow layers of connected faces that are + // not grown. Is used to not do layers at all close to features. + nGrow 1; + + // Expansion factor for layer mesh + expansionRatio 1.3; + + // Ratio of cell size in final added cell layer to cell size + // outside layer + finalLayerRatio 0.3; + + // Number of smoothing iterations of surface normals + nSmoothSurfaceNormals 1; + + // Number of smoothing iterations of interior mesh movement direction + nSmoothNormals 3; + + // Smooth layer thickness over surface patches + nSmoothThickness 10; + + // Stop layer growth on highly warped cells + maxFaceThicknessRatio 0.5; + + // Reduce layer growth where ratio thickness to medial + // distance is large + maxThicknessToMedialRatio 0.3; + + // Angle used to pick up medial axis points + minMedianAxisAngle 130; + + // Create buffer region for new layer terminations + nBufferCellsNoExtrude 0; + + thickness 0.5; + nSmoothDispl 4; +} + + +// For mesh motion +motionDict +{ + // + // Mesh Quality Parameters. Decide when mesh is good enough to stop + // smoothing. + // + + //- Maximum non-orthogonality allowed. Set to 180 to disable. + maxNonOrtho 65; + + //- Max skewness allowed. Set to <0 to disable. + maxBoundarySkewness 20; + maxInternalSkewness 4; + + //- Max concaveness allowed. Is angle (in degrees) below which concavity + // is allowed. 0 is straight face, <0 would be convex face. + // Set to 180 to disable. + maxConcave 80; + + //- Minimum projected area v.s. actual area. Set to -1 to disable. + minFlatness 0.5; + + //- Minimum pyramid volume. Is absolute volume of cell pyramid. + // Set to very negative number (e.g. -1E30) to disable. + minVol 1e-13; + + //- Minimum face area. Set to <0 to disable. + minArea -1; + + //- Minimum face twist. Set to <-1 to disable. dot product of face normal + //- and face centre triangles normal + minTwist 0.05; + + //- minimum normalised cell determinant + //- 1 = hex, <= 0 = folded or flattened illegal cell + minDeterminant 0.001; + + //- minFaceWeight (0 -> 0.5) + minFaceWeight 0.05; + + //- minVolRatio (0 -> 1) + minVolRatio 0.01; + + + //must be >0 for Fluent compatibility + minTriangleTwist -1; + + // Advanced + + //- Number of error distribution iterations + nSmoothScale 4; + //- amount to scale back displacement at error points + errorReduction 0.75; +} + + +// ************************************************************************* // diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C index f01b31d70d..6435d491d9 100644 --- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C +++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoHexMeshDriver.C @@ -52,9 +52,9 @@ defineTypeNameAndDebug(autoHexMeshDriver, 0); // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // Check writing tolerance before doing any serious work -Foam::scalar Foam::autoHexMeshDriver::getMergeDistance() const +Foam::scalar Foam::autoHexMeshDriver::getMergeDistance(const scalar mergeTol) + const { - const scalar mergeTol = readScalar(dict_.lookup("mergeTolerance")); const boundBox& meshBb = mesh_.bounds(); scalar mergeDist = mergeTol*mag(meshBb.max() - meshBb.min()); scalar writeTol = std::pow @@ -70,7 +70,7 @@ Foam::scalar Foam::autoHexMeshDriver::getMergeDistance() const if (mesh_.time().writeFormat() == IOstream::ASCII && mergeTol < writeTol) { - FatalErrorIn("autoHexMeshDriver::getMergeDistance() const") + FatalErrorIn("autoHexMeshDriver::getMergeDistance(const scalar) const") << "Your current settings specify ASCII writing with " << IOstream::defaultPrecision() << " digits precision." << endl << "Your merging tolerance (" << mergeTol << ") is finer than this." @@ -306,9 +306,8 @@ Foam::autoHexMeshDriver::autoHexMeshDriver curvature_(readScalar(dict_.lookup("curvature"))), nBufferLayers_(readLabel(dict_.lookup("nBufferLayers"))), keepPoints_(dict_.lookup("keepPoints")), - mergeDist_(getMergeDistance()) + mergeDist_(getMergeDistance(readScalar(dict_.lookup("mergeTolerance")))) { - if (debug_ > 0) { meshRefinement::debug = debug_; @@ -471,17 +470,34 @@ Foam::autoHexMeshDriver::autoHexMeshDriver { if (nTrisPerRegion[i] > 0) { + label globalRegionI = surfaces().globalRegion(surfI, i); + + // Use optionally specified patch type and name + word patchType = surfaces().patchType()[globalRegionI]; + if (patchType == "") + { + patchType = wallPolyPatch::typeName; + } + + word patchName = surfaces().patchName()[globalRegionI]; + if (patchName == "") + { + patchName = + surfaces().names()[surfI] + + '_' + + regions[i].name(); + } + label patchI = meshRefinement::addPatch ( mesh, - //s.searchableSurface::name() + '_' + regions[i].name(), - surfaces().names()[surfI] + '_' + regions[i].name(), - wallPolyPatch::typeName + patchName, + patchType ); Info<< patchI << '\t' << regions[i].name() << nl; - globalToPatch_[surfaces().globalRegion(surfI, i)] = patchI; + globalToPatch_[globalRegionI] = patchI; } } @@ -593,15 +609,283 @@ Foam::autoHexMeshDriver::autoHexMeshDriver } +// Construct from separate dictionaries. +Foam::autoHexMeshDriver::autoHexMeshDriver +( + fvMesh& mesh, + const dictionary& controlDict, + const dictionary& refineDict, + const dictionary& decomposeDict +) +: + mesh_(mesh), + dict_(controlDict), + debug_(readLabel(controlDict.lookup("debug"))), + maxGlobalCells_(readLabel(refineDict.lookup("cellLimit"))), + maxLocalCells_(readLabel(refineDict.lookup("procCellLimit"))), + minRefineCells_(readLabel(refineDict.lookup("minimumRefine"))), + curvature_(readScalar(refineDict.lookup("curvature"))), + nBufferLayers_(readLabel(refineDict.lookup("nBufferLayers"))), + keepPoints_(refineDict.lookup("keepPoints")), + mergeDist_ + ( + getMergeDistance(readScalar(refineDict.lookup("mergeTolerance"))) + ) +{ + if (debug_ > 0) + { + meshRefinement::debug = debug_; + autoHexMeshDriver::debug = debug_; + } + + Info<< "Overall cell limit : " << maxGlobalCells_ + << endl; + Info<< "Per processor cell limit : " << maxLocalCells_ + << endl; + Info<< "Minimum number of cells to refine : " << minRefineCells_ + << endl; + Info<< "Curvature : " << curvature_ + << nl << endl; + Info<< "Layers between different refinement levels : " << nBufferLayers_ + << endl; + + // Check keepPoints are sensible + findCells(keepPoints_); + + + // Read refinement shells + // ~~~~~~~~~~~~~~~~~~~~~~ + + { + Info<< "Reading refinement shells." << endl; + + PtrList shellDicts(refineDict.lookup("refinementShells")); + + shells_.setSize(shellDicts.size()); + shellLevels_.setSize(shellDicts.size()); + shellRefineInside_.setSize(shellDicts.size()); + + forAll(shellDicts, i) + { + const dictionary& dict = shellDicts[i]; + + shells_.set + ( + i, + searchableSurface::New + ( + dict.lookup("type"), + dict.lookup("name"), + mesh_.time(), + dict + ) + ); + shellLevels_[i] = readLabel(dict.lookup("level")); + shellRefineInside_[i] = Switch(dict.lookup("refineInside")); + + if (shellRefineInside_[i]) + { + Info<< "Refinement level " << shellLevels_[i] + << " for all cells inside " << shells_[i].name() << endl; + } + else + { + Info<< "Refinement level " << shellLevels_[i] + << " for all cells outside " << shells_[i].name() << endl; + } + } + + Info<< "Read refinement shells in = " + << mesh_.time().cpuTimeIncrement() << " s" << endl; + + // Orient shell surfaces before any searching is done. + Info<< "Orienting triSurface shells so point far away is outside." + << endl; + orientOutside(shells_); + Info<< "Oriented shells in = " + << mesh_.time().cpuTimeIncrement() << " s" << endl; + } + + + // Read refinement surfaces + // ~~~~~~~~~~~~~~~~~~~~~~~~ + + { + Info<< "Reading surfaces and constructing search trees." << endl; + + surfacesPtr_.reset + ( + new refinementSurfaces + ( + IOobject + ( + "", // dummy name + mesh_.time().constant(), // directory + "triSurface", // instance + mesh_.time(), // registry + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + refineDict.lookup("surfaces") + ) + ); + Info<< "Read surfaces in = " + << mesh_.time().cpuTimeIncrement() << " s" << endl; + + // Orient surfaces (if they're closed) before any searching is done. + Info<< "Orienting (closed) surfaces so keepPoint is outside." << endl; + forAll(surfaces(), i) + { + if (refinementSurfaces::isSurfaceClosed(surfaces()[i])) + { + refinementSurfaces::orientSurface + ( + keepPoints_[0], + surfacesPtr_()[i] + ); + } + } + Info<< "Oriented closed surfaces in = " + << mesh_.time().cpuTimeIncrement() << " s" << endl; + + Info<< "Setting refinement level of surface to be consistent" + << " with shells." << endl; + surfacesPtr_().setMinLevelFields + ( + shells_, + shellLevels_, + shellRefineInside_ + ); + Info<< "Checked shell refinement in = " + << mesh_.time().cpuTimeIncrement() << " s" << endl; + } + + // Check faceZones are synchronised + checkCoupledFaceZones(); + + + // Add all the surface regions as patches + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + { + Info<< nl + << "Adding patches for surface regions" << nl + << "----------------------------------" << nl + << endl; + + // From global region number to mesh patch. + globalToPatch_.setSize(surfaces().nRegions(), -1); + + Info<< "Patch\tRegion" << nl + << "-----\t------" + << endl; + + forAll(surfaces(), surfI) + { + const triSurfaceMesh& s = surfaces()[surfI]; + + Info<< surfaces().names()[surfI] << ':' << nl << nl; + + const geometricSurfacePatchList& regions = s.patches(); + + labelList nTrisPerRegion(surfaces().countRegions(s)); + + forAll(regions, i) + { + if (nTrisPerRegion[i] > 0) + { + label patchI = meshRefinement::addPatch + ( + mesh, + //s.searchableSurface::name() + '_' + regions[i].name(), + surfaces().names()[surfI] + '_' + regions[i].name(), + wallPolyPatch::typeName + ); + + Info<< patchI << '\t' << regions[i].name() << nl; + + globalToPatch_[surfaces().globalRegion(surfI, i)] = patchI; + } + } + + Info<< nl; + } + Info<< "Added patches in = " + << mesh_.time().cpuTimeIncrement() << " s" << endl; + } + + // Parallel + // ~~~~~~~~ + + { + // Decomposition + decomposerPtr_ = decompositionMethod::New + ( + decomposeDict, + mesh_ + ); + decompositionMethod& decomposer = decomposerPtr_(); + + + if (Pstream::parRun() && !decomposer.parallelAware()) + { + FatalErrorIn("autoHexMeshDriver::autoHexMeshDriver(const IOobject&, fvMesh&)") + << "You have selected decomposition method " + << decomposer.typeName + << " which is not parallel aware." << endl + << "Please select one that is (parMetis, hierarchical)" + << exit(FatalError); + } + + // Mesh distribution engine (uses tolerance to reconstruct meshes) + distributorPtr_.reset(new fvMeshDistribute(mesh_, mergeDist_)); + } + + + // Refinement engine + // ~~~~~~~~~~~~~~~~~ + + { + Info<< nl + << "Determining initial surface intersections" << nl + << "-----------------------------------------" << nl + << endl; + + // Main refinement engine + meshRefinerPtr_.reset + ( + new meshRefinement + ( + mesh, + mergeDist_, // tolerance used in sorting coordinates + surfaces() + ) + ); + Info<< "Calculated surface intersections in = " + << mesh_.time().cpuTimeIncrement() << " s" << endl; + + // Some stats + meshRefinerPtr_().printMeshInfo(debug_, "Initial mesh"); + + meshRefinerPtr_().write + ( + debug_&meshRefinement::OBJINTERSECTIONS, + mesh_.time().path()/mesh_.time().timeName() + ); + } +} + + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // Read explicit feature edges -Foam::label Foam::autoHexMeshDriver::readFeatureEdges() +Foam::label Foam::autoHexMeshDriver::readFeatureEdges +( + const PtrList& featDicts +) { Info<< "Reading external feature lines." << endl; - PtrList featDicts(dict_.lookup("features")); - featureMeshes_.setSize(featDicts.size()); featureLevels_.setSize(featDicts.size()); @@ -647,13 +931,13 @@ Foam::label Foam::autoHexMeshDriver::readFeatureEdges() Foam::label Foam::autoHexMeshDriver::featureEdgeRefine ( + const PtrList& featDicts, const label maxIter, const label minRefine ) { // Read explicit feature edges - readFeatureEdges(); - + readFeatureEdges(featDicts); meshRefinement& meshRefiner = meshRefinerPtr_(); @@ -1359,122 +1643,232 @@ void Foam::autoHexMeshDriver::writeMesh(const string& msg) const } -void Foam::autoHexMeshDriver::doMesh() -{ - Switch doRefine(dict_.lookup("doRefine")); - Switch doSnap(dict_.lookup("doSnap")); - Switch doLayers(dict_.lookup("doLayers")); - Info<< "Do refinement : " << doRefine << nl - << "Do snapping : " << doSnap << nl - << "Do layers : " << doLayers << nl + +void Foam::autoHexMeshDriver::doRefine +( + const dictionary& refineDict, + const bool prepareForSnapping +) +{ + Info<< nl + << "Refinement phase" << nl + << "----------------" << nl << endl; - if (doRefine) + const_cast(mesh_.time())++; + + PtrList featDicts(refineDict.lookup("features")); + + // Refine around feature edges + featureEdgeRefine + ( + featDicts, + 100, // maxIter + 0 // min cells to refine + ); + + // Refine based on surface + surfaceOnlyRefine + ( + 100 // maxIter + ); + + // Remove cells (a certain distance) beyond surface intersections + removeInsideCells + ( + 1 // nBufferLayers + ); + + // Internal mesh refinement + shellRefine + ( + 100 // maxIter + ); + + // Introduce baffles at surface intersections + baffleAndSplitMesh(prepareForSnapping); + + // Mesh is at its finest. Do optional zoning. + zonify(); + + // Pull baffles apart + splitAndMergeBaffles(prepareForSnapping); + + // Do something about cells with refined faces on the boundary + if (prepareForSnapping) { - Info<< nl - << "Refinement phase" << nl - << "----------------" << nl - << endl; - - const_cast(mesh_.time())++; - - // Refine around feature edges - featureEdgeRefine - ( - 100, // maxIter - 0 // min cells to refine - ); - - // Refine based on surface - surfaceOnlyRefine - ( - 100 // maxIter - ); - - // Remove cells (a certain distance) beyond surface intersections - removeInsideCells - ( - 1 // nBufferLayers - ); - - // Internal mesh refinement - shellRefine - ( - 100 // maxIter - ); - - // Introduce baffles at surface intersections - baffleAndSplitMesh(doSnap); - - // Mesh is at its finest. Do optional zoning. - zonify(); - - // Pull baffles apart - splitAndMergeBaffles(doSnap); - - // Do something about cells with refined faces on the boundary - if (doSnap) - { - mergePatchFaces(); - } - - // Do final balancing. Keep zoned faces on one processor. - balance(true, false); - - // Write mesh - writeMesh("Refined mesh"); + mergePatchFaces(); } + // Do final balancing. Keep zoned faces on one processor. + balance(true, false); + + // Write mesh + writeMesh("Refined mesh"); +} + + +void Foam::autoHexMeshDriver::doSnap +( + const dictionary& snapDict, + const dictionary& motionDict +) +{ + Info<< nl + << "Morphing phase" << nl + << "--------------" << nl + << endl; + + const_cast(mesh_.time())++; + + // Get the labels of added patches. + labelList adaptPatchIDs(getSurfacePatches()); + + // Create baffles (pairs of faces that share the same points) + // Baffles stored as owner and neighbour face that have been created. + List baffles; + createZoneBaffles(baffles); - if (doSnap) { - Info<< nl - << "Morphing phase" << nl - << "--------------" << nl + autoPtr ppPtr + ( + meshRefinement::makePatch + ( + mesh_, + adaptPatchIDs + ) + ); + indirectPrimitivePatch& pp = ppPtr(); + + // Distance to attact to nearest feature on surface + const scalarField snapDist(calcSnapDistance(snapDict, pp)); + + + // Construct iterative mesh mover. + Info<< "Constructing mesh displacer ..." << endl; + Info<< "Using mesh parameters " << motionDict << nl << endl; + + pointMesh pMesh(mesh_); + + motionSmoother meshMover + ( + mesh_, + pp, + adaptPatchIDs, + meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs), + motionDict + ); + + + // Check initial mesh + Info<< "Checking initial mesh ..." << endl; + labelHashSet wrongFaces(mesh_.nFaces()/100); + motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces); + const label nInitErrors = returnReduce + ( + wrongFaces.size(), + sumOp