diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C index 14bf0ee9eb9d5eec0442060b3ca8a308f2bd1e74..f177e7de8187f4d9deda4e71de7e77f4246a2b0c 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C @@ -247,7 +247,8 @@ int main(int argc, char *argv[]) refinementSurfaces surfaces ( allGeometry, - refineDict.subDict("refinementSurfaces") + refineDict.subDict("refinementSurfaces"), + refineDict.lookupOrDefault("gapLevelIncrement", 0) ); Info<< "Read refinement surfaces in = " << mesh.time().cpuTimeIncrement() << " s" << nl << endl; @@ -551,6 +552,13 @@ int main(int argc, char *argv[]) const Switch wantSnap(meshDict.lookup("snap")); const Switch wantLayers(meshDict.lookup("addLayers")); + // Refinement parameters + const refinementParameters refineParams(refineDict); + + // Snap parameters + const snapParameters snapParams(snapDict); + + if (wantRefine) { cpuTime timer; @@ -564,15 +572,20 @@ int main(int argc, char *argv[]) globalToSlavePatch ); - // Refinement parameters - refinementParameters refineParams(refineDict); if (!overwrite && !debug) { const_cast<Time&>(mesh.time())++; } - refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict); + refineDriver.doRefine + ( + refineDict, + refineParams, + snapParams, + wantSnap, + motionDict + ); writeMesh ( @@ -596,20 +609,14 @@ int main(int argc, char *argv[]) globalToSlavePatch ); - // Snap parameters - snapParameters snapParams(snapDict); - // Temporary hack to get access to resolveFeatureAngle - scalar curvature; - { - refinementParameters refineParams(refineDict); - curvature = refineParams.curvature(); - } - if (!overwrite && !debug) { const_cast<Time&>(mesh.time())++; } + // Use the resolveFeatureAngle from the refinement parameters + scalar curvature = refineParams.curvature(); + snapDriver.doSnap(snapDict, motionDict, curvature, snapParams); writeMesh @@ -637,17 +644,12 @@ int main(int argc, char *argv[]) // Layer addition parameters layerParameters layerParams(layerDict, mesh.boundaryMesh()); - //!!! Temporary hack to get access to maxLocalCells - bool preBalance; - { - refinementParameters refineParams(refineDict); - - preBalance = returnReduce - ( - (mesh.nCells() >= refineParams.maxLocalCells()), - orOp<bool>() - ); - } + // Use the maxLocalCells from the refinement parameters + bool preBalance = returnReduce + ( + (mesh.nCells() >= refineParams.maxLocalCells()), + orOp<bool>() + ); if (!overwrite && !debug) diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index 7d1ad76db3007982d4977face4a0b095f9697d57..ba16b1c4aadcbbcc05b83075b4cc7b447aca471f 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -153,7 +153,7 @@ castellatedMeshControls } - //- Optional increment (on top of max level) in small gaps. + //- Optional increment (on top of max level) in small gaps //gapLevelIncrement 2; //- Optional angle to detect small-large cell situation @@ -184,14 +184,17 @@ castellatedMeshControls // - used if feature snapping (see snapControls below) is used resolveFeatureAngle 30; + //- Optional increment (on top of max level) in small gaps + //gapLevelIncrement 2; + // Planar angle: // - used to determine if surface normals - // are roughly the same or opposite. Used e.g. in gap refinement + // are roughly the same or opposite. Used e.g. in proximity refinement // and to decide when to merge free-standing baffles // // If not specified same as resolveFeatureAngle - planarAngle 15; + planarAngle 30; // Region-wise refinement @@ -235,6 +238,9 @@ castellatedMeshControls // are only on the boundary of corresponding cellZones or also allow // free-standing zone faces. Not used if there are no faceZones. allowFreeStandingZoneFaces true; + + // + //useTopologicalSnapDetection false; } // Settings for the snapping. diff --git a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C index 74c5fe0bf635e34531220feb05165217df381a9a..c7ad881527ae51d82f42eedbe579412fc80abd7e 100644 --- a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C +++ b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C @@ -168,6 +168,12 @@ int main(int argc, char *argv[]) "add exposed internal faces to specified patch instead of to " "'oldInternalFaces'" ); + argList::addOption + ( + "resultTime", + "time", + "specify a time for the resulting mesh" + ); #include "setRootCase.H" #include "createTime.H" runTime.functionObjects().off(); @@ -178,10 +184,12 @@ int main(int argc, char *argv[]) #include "createNamedMesh.H" - const word oldInstance = mesh.pointsInstance(); - const word setName = args[1]; - const bool overwrite = args.optionFound("overwrite"); + + word resultInstance = mesh.pointsInstance(); + const bool specifiedInstance = + args.optionFound("overwrite") + || args.optionReadIfPresent("resultTime", resultInstance); Info<< "Reading cell set from " << setName << endl << endl; @@ -347,13 +355,14 @@ int main(int argc, char *argv[]) // Write mesh and fields to new time // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - if (!overwrite) + if (!specifiedInstance) { runTime++; } else { - subsetter.subMesh().setInstance(oldInstance); + runTime.setTime(instant(resultInstance), 0); + subsetter.subMesh().setInstance(runTime.timeName()); } Info<< "Writing subsetted mesh and fields to time " << runTime.timeName() diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C index 75af8eb20a2744a24952b9e44fcfc91a570745c1..dd51eb645b8daf05e192f4d88eb393dfff5d5968 100644 --- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C +++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGAgglomerations/GAMGAgglomeration/GAMGAgglomeration.C @@ -165,7 +165,7 @@ void Foam::GAMGAgglomeration::compactLevels(const label nCreatedLevels) << setw(8) << setprecision(4) << totFaceCellRatio/totNprocs << setw(8) << setprecision(4) << maxFaceCellRatio << " " - << setw(8) << totNInt/totNprocs + << setw(8) << scalar(totNInt)/totNprocs << setw(8) << maxNInt << " " << setw(8) << setprecision(4) << totRatio/totNprocs diff --git a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrads.C b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrads.C index 137850e65fe856fa4839fcedaaf2596e22fffa1d..4006fd5e706204c2495defa260f066525148caf6 100644 --- a/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrads.C +++ b/src/finiteVolume/finiteVolume/gradSchemes/limitedGradSchemes/faceMDLimitedGrad/faceMDLimitedGrads.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -235,7 +235,7 @@ Foam::fv::faceMDLimitedGrad<Foam::vector>::calcGrad g[own], maxFace - vvfOwn, minFace - vvfOwn, - Cf[facei] - C[nei] + Cf[facei] - C[own] ); diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.H index 453603d287a78950353c577b3151d6974d71b44e..1ef0c000f7840e31c7b8cf9be542afcf8d11f83e 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwind.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -111,16 +111,7 @@ public: mesh.gradScheme(gradSchemeName_) ) ) - { - if (!schemeData.eof()) - { - IOWarningIn("linearUpwind(const fvMesh&, Istream&)", schemeData) - << "unexpected additional entries in stream." << nl - << " Only the name of the gradient scheme in the" - " 'gradSchemes' dictionary should be specified." - << endl; - } - } + {} //- Construct from faceFlux and Istream linearUpwind @@ -140,21 +131,7 @@ public: mesh.gradScheme(gradSchemeName_) ) ) - { - - if (!schemeData.eof()) - { - IOWarningIn - ( - "linearUpwind(const fvMesh&, " - "const surfaceScalarField& faceFlux, Istream&)", - schemeData - ) << "unexpected additional entries in stream." << nl - << " Only the name of the gradient scheme in the" - " 'gradSchemes' dictionary should be specified." - << endl; - } - } + {} // Member Functions diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H index 862e79b812bcf8c091af5add10419d0071e7d7f9..88034a7de9cea44b04d3136ec251f780321dad5d 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/linearUpwind/linearUpwindV.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -110,19 +110,7 @@ public: mesh.gradScheme(gradSchemeName_) ) ) - { - if (!schemeData.eof()) - { - IOWarningIn - ( - "linearUpwindV(const fvMesh&, Istream&)", - schemeData - ) << "unexpected additional entries in stream." << nl - << " Only the name of the gradient scheme in the" - " 'gradSchemes' dictionary should be specified." - << endl; - } - } + {} //- Construct from faceFlux and Istream linearUpwindV @@ -142,20 +130,7 @@ public: mesh.gradScheme(gradSchemeName_) ) ) - { - if (!schemeData.eof()) - { - IOWarningIn - ( - "linearUpwindV(const fvMesh&, " - "const surfaceScalarField& faceFlux, Istream&)", - schemeData - ) << "unexpected additional entries in stream." << nl - << " Only the name of the gradient scheme in the" - " 'gradSchemes' dictionary should be specified." - << endl; - } - } + {} // Member Functions diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C index d78866bcd78b781a5814552394fbdae5548dbe4b..82bfa13522c573e57f274d02e34174f8b0179251 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C @@ -35,6 +35,7 @@ License #include "shellSurfaces.H" #include "mapDistributePolyMesh.H" #include "unitConversion.H" +#include "snapParameters.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -610,9 +611,16 @@ Foam::label Foam::autoRefineDriver::danglingCellRefine << " cells (out of " << mesh.globalData().nTotalCells() << ')' << endl; - // Stop when no cells to refine. No checking of minRefineCells since - // too few cells - if (nCellsToRefine == 0) + // Stop when no cells to refine. After a few iterations check if too + // few cells + if + ( + nCellsToRefine == 0 + || ( + iter >= 1 + && nCellsToRefine <= refineParams.minRefineCells() + ) + ) { Info<< "Stopping refining since too few cells selected." << nl << endl; @@ -870,6 +878,7 @@ Foam::label Foam::autoRefineDriver::shellRefine void Foam::autoRefineDriver::baffleAndSplitMesh ( const refinementParameters& refineParams, + const snapParameters& snapParams, const bool handleSnapProblems, const dictionary& motionDict ) @@ -887,10 +896,17 @@ void Foam::autoRefineDriver::baffleAndSplitMesh meshRefiner_.baffleAndSplitMesh ( handleSnapProblems, // detect&remove potential snap problem + + // Snap problem cell detection + snapParams, + refineParams.useTopologicalSnapDetection(), false, // perpendicular edge connected cells scalarField(0), // per region perpendicular angle + + // Free standing baffles !handleSnapProblems, // merge free standing baffles? refineParams.planarAngle(), + motionDict, const_cast<Time&>(mesh.time()), globalToMasterPatch_, @@ -951,6 +967,7 @@ void Foam::autoRefineDriver::zonify void Foam::autoRefineDriver::splitAndMergeBaffles ( const refinementParameters& refineParams, + const snapParameters& snapParams, const bool handleSnapProblems, const dictionary& motionDict ) @@ -973,10 +990,17 @@ void Foam::autoRefineDriver::splitAndMergeBaffles meshRefiner_.baffleAndSplitMesh ( handleSnapProblems, + + // Snap problem cell detection + snapParams, + refineParams.useTopologicalSnapDetection(), handleSnapProblems, // remove perp edge connected cells perpAngle, // perp angle + + // Free standing baffles true, // merge free standing baffles? refineParams.planarAngle(), // planar angle + motionDict, const_cast<Time&>(mesh.time()), globalToMasterPatch_, @@ -1081,6 +1105,7 @@ void Foam::autoRefineDriver::doRefine ( const dictionary& refineDict, const refinementParameters& refineParams, + const snapParameters& snapParams, const bool prepareForSnapping, const dictionary& motionDict ) @@ -1146,13 +1171,25 @@ void Foam::autoRefineDriver::doRefine // Introduce baffles at surface intersections. Remove sections unreachable // from keepPoint. - baffleAndSplitMesh(refineParams, prepareForSnapping, motionDict); + baffleAndSplitMesh + ( + refineParams, + snapParams, + prepareForSnapping, + motionDict + ); // Mesh is at its finest. Do optional zoning. zonify(refineParams); // Pull baffles apart - splitAndMergeBaffles(refineParams, prepareForSnapping, motionDict); + splitAndMergeBaffles + ( + refineParams, + snapParams, + prepareForSnapping, + motionDict + ); // Do something about cells with refined faces on the boundary if (prepareForSnapping) diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H index c44a14a79457e33cdb82643d24c0701781fea566..a72f3c0796d291f8d16409cefdfd09bec479f9c0 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.H @@ -43,6 +43,8 @@ namespace Foam // Forward declaration of classes class refinementParameters; +class snapParameters; + class meshRefinement; class decompositionMethod; class fvMeshDistribute; @@ -121,6 +123,7 @@ class autoRefineDriver void baffleAndSplitMesh ( const refinementParameters& refineParams, + const snapParameters& snapParams, const bool handleSnapProblems, const dictionary& motionDict ); @@ -131,6 +134,7 @@ class autoRefineDriver void splitAndMergeBaffles ( const refinementParameters& refineParams, + const snapParameters& snapParams, const bool handleSnapProblems, const dictionary& motionDict ); @@ -176,6 +180,7 @@ public: ( const dictionary& refineDict, const refinementParameters& refineParams, + const snapParameters& snapParams, const bool prepareForSnapping, const dictionary& motionDict ); diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C index 4c262bd132049c99532332526fc94d9e4c059e80..39273051772b76f7d4fccca6a001a8b88c4a80a8 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C @@ -121,7 +121,7 @@ Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement ( const motionSmoother& meshMover, const List<labelPair>& baffles -) const +) { const indirectPrimitivePatch& pp = meshMover.patch(); @@ -615,14 +615,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles Foam::scalarField Foam::autoSnapDriver::calcSnapDistance ( + const fvMesh& mesh, const snapParameters& snapParams, const indirectPrimitivePatch& pp -) const +) { const edgeList& edges = pp.edges(); const labelListList& pointEdges = pp.pointEdges(); const pointField& localPoints = pp.localPoints(); - const fvMesh& mesh = meshRefiner_.mesh(); scalarField maxEdgeLen(localPoints.size(), -GREAT); @@ -655,13 +655,14 @@ Foam::scalarField Foam::autoSnapDriver::calcSnapDistance void Foam::autoSnapDriver::preSmoothPatch ( + const meshRefinement& meshRefiner, const snapParameters& snapParams, const label nInitErrors, const List<labelPair>& baffles, motionSmoother& meshMover -) const +) { - const fvMesh& mesh = meshRefiner_.mesh(); + const fvMesh& mesh = meshRefiner.mesh(); labelList checkFaces; @@ -724,11 +725,11 @@ void Foam::autoSnapDriver::preSmoothPatch { const_cast<Time&>(mesh.time())++; Info<< "Writing patch smoothed mesh to time " - << meshRefiner_.timeName() << '.' << endl; - meshRefiner_.write + << meshRefiner.timeName() << '.' << endl; + meshRefiner.write ( debug, - mesh.time().path()/meshRefiner_.timeName() + mesh.time().path()/meshRefiner.timeName() ); Info<< "Dumped mesh in = " << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl; @@ -742,12 +743,11 @@ void Foam::autoSnapDriver::preSmoothPatch // Get (pp-local) indices of points that are both on zone and on patched surface Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints ( + const fvMesh& mesh, const indirectPrimitivePatch& pp, const word& zoneName -) const +) { - const fvMesh& mesh = meshRefiner_.mesh(); - label zoneI = mesh.faceZones().findZoneID(zoneName); if (zoneI == -1) @@ -755,7 +755,7 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints FatalErrorIn ( "autoSnapDriver::getZoneSurfacePoints" - "(const indirectPrimitivePatch&, const word&)" + "(const fvMesh&, const indirectPrimitivePatch&, const word&)" ) << "Cannot find zone " << zoneName << exit(FatalError); } @@ -793,17 +793,17 @@ Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints Foam::vectorField Foam::autoSnapDriver::calcNearestSurface ( + const meshRefinement& meshRefiner, const scalarField& snapDist, - motionSmoother& meshMover -) const + const indirectPrimitivePatch& pp +) { Info<< "Calculating patchDisplacement as distance to nearest surface" << " point ..." << endl; - const indirectPrimitivePatch& pp = meshMover.patch(); const pointField& localPoints = pp.localPoints(); - const refinementSurfaces& surfaces = meshRefiner_.surfaces(); - const fvMesh& mesh = meshRefiner_.mesh(); + const refinementSurfaces& surfaces = meshRefiner.surfaces(); + const fvMesh& mesh = meshRefiner.mesh(); // Displacement per patch point vectorField patchDisp(localPoints.size(), vector::zero); @@ -815,9 +815,9 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface // Divide surfaces into zoned and unzoned labelList zonedSurfaces = - meshRefiner_.surfaces().getNamedSurfaces(); + meshRefiner.surfaces().getNamedSurfaces(); labelList unzonedSurfaces = - meshRefiner_.surfaces().getUnnamedSurfaces(); + meshRefiner.surfaces().getUnnamedSurfaces(); // 1. All points to non-interface surfaces @@ -870,6 +870,7 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface ( getZoneSurfacePoints ( + mesh, pp, faceZoneNames[zoneSurfI] ) @@ -966,6 +967,254 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurface } +////XXXXXXXXX +//// Get (pp-local) indices of points that are on both patches +//Foam::labelList Foam::autoSnapDriver::getPatchSurfacePoints +//( +// const fvMesh& mesh, +// const indirectPrimitivePatch& allPp, +// const polyPatch& pp +//) +//{ +// // Could use PrimitivePatch & localFaces to extract points but might just +// // as well do it ourselves. +// +// boolList pointOnZone(allPp.nPoints(), false); +// +// forAll(pp, i) +// { +// const face& f = pp[i]; +// +// forAll(f, fp) +// { +// label meshPointI = f[fp]; +// +// Map<label>::const_iterator iter = +// allPp.meshPointMap().find(meshPointI); +// +// if (iter != allPp.meshPointMap().end()) +// { +// label pointI = iter(); +// pointOnZone[pointI] = true; +// } +// } +// } +// +// return findIndices(pointOnZone, true); +//} +//Foam::vectorField Foam::autoSnapDriver::calcNearestLocalSurface +//( +// const meshRefinement& meshRefiner, +// const scalarField& snapDist, +// const indirectPrimitivePatch& pp +//) +//{ +// Info<< "Calculating patchDisplacement as distance to nearest" +// << " local surface point ..." << endl; +// +// const pointField& localPoints = pp.localPoints(); +// const refinementSurfaces& surfaces = meshRefiner.surfaces(); +// const fvMesh& mesh = meshRefiner.mesh(); +// const polyBoundaryMesh& pbm = mesh.boundaryMesh(); +// +// +//// // Assume that all patch-internal points get attracted to their surface +//// // only. So we want to know if point is on multiple regions +//// +//// labelList minPatch(mesh.nPoints(), labelMax); +//// labelList maxPatch(mesh.nPoints(), labelMin); +//// +//// forAll(meshMover.adaptPatchIDs(), i) +//// { +//// label patchI = meshMover.adaptPatchIDs()[i]; +//// const labelList& meshPoints = pbm[patchI].meshPoints(); +//// +//// forAll(meshPoints, meshPointI) +//// { +//// label meshPointI = meshPoints[meshPointI]; +//// minPatch[meshPointI] = min(minPatch[meshPointI], patchI); +//// maxPatch[meshPointI] = max(maxPatch[meshPointI], patchI); +//// } +//// } +//// +//// syncTools::syncPointList +//// ( +//// mesh, +//// minPatch, +//// minEqOp<label>(), // combine op +//// labelMax // null value +//// ); +//// syncTools::syncPointList +//// ( +//// mesh, +//// maxPatch, +//// maxEqOp<label>(), // combine op +//// labelMin // null value +//// ); +// +// // Now all points with minPatch != maxPatch will be on the outside of +// // the patch. +// +// // Displacement per patch point +// vectorField patchDisp(localPoints.size(), vector::zero); +// // Current best snap distance +// scalarField minSnapDist(snapDist); +// // Current surface snapped to +// labelList snapSurf(localPoints.size(), -1); +// +// const labelList& surfaceGeometry = surfaces.surfaces(); +// forAll(surfaceGeometry, surfI) +// { +// label geomI = surfaceGeometry[surfI]; +// const wordList& regNames = allGeometry.regionNames()[geomI]; +// forAll(regNames, regionI) +// { +// label globalRegionI = surfaces.globalRegion(surfI, regionI); +// // Collect master patch points +// label masterPatchI = globalToMasterPatch_[globalRegionI]; +// label slavePatchI = globalToSlavePatch_[globalRegionI]; +// +// labelList patchPointIndices +// ( +// getPatchSurfacePoints +// ( +// mesh, +// pp, +// pbm[masterPatchI] +// ) +// ); +// +// // Find nearest for points both on faceZone and pp. +// List<pointIndexHit> hitInfo; +// surfaces.findNearest +// ( +// surfI, +// regionI, +// pointField(localPoints, patchPointIndices), +// sqr(scalarField(minSnapDist, patchPointIndices)), +// hitInfo +// ); +// +// forAll(hitInfo, i) +// { +// label pointI = patchPointIndices[i]; +// +// if (hitInfo[i].hit()) +// { +// const point& pt = hitInfo[i].hitPoint(); +// patchDisp[pointI] = pt-localPoints[pointI]; +// minSnapDist[pointI] = min +// ( +// minSnapDist[pointI], +// mag(patchDisp[pointI]) +// ); +// snapSurf[pointI] = surfI; +// } +// } +// +// // Slave patch +// if (slavePatchI != masterPatchI) +// { +// labelList patchPointIndices +// ( +// getPatchSurfacePoints +// ( +// mesh, +// pp, +// pbm[slavePatchI] +// ) +// ); +// +// // Find nearest for points both on faceZone and pp. +// List<pointIndexHit> hitInfo; +// surfaces.findNearest +// ( +// surfI, +// regionI, +// pointField(localPoints, patchPointIndices), +// sqr(scalarField(minSnapDist, patchPointIndices)), +// hitInfo +// ); +// +// forAll(hitInfo, i) +// { +// label pointI = patchPointIndices[i]; +// +// if (hitInfo[i].hit()) +// { +// const point& pt = hitInfo[i].hitPoint(); +// patchDisp[pointI] = pt-localPoints[pointI]; +// minSnapDist[pointI] = min +// ( +// minSnapDist[pointI], +// mag(patchDisp[pointI]) +// ); +// snapSurf[pointI] = surfI; +// } +// } +// } +// } +// } +// +// +// // Check if all points are being snapped +// forAll(snapSurf, pointI) +// { +// if (snapSurf[pointI] == -1) +// { +// WarningIn("autoSnapDriver::calcNearestLocalSurface(..)") +// << "For point:" << pointI +// << " coordinate:" << localPoints[pointI] +// << " did not find any surface within:" +// << minSnapDist[pointI] +// << " metre." << endl; +// } +// } +// +// { +// scalarField magDisp(mag(patchDisp)); +// +// Info<< "Wanted displacement : average:" +// << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>()) +// << " min:" << gMin(magDisp) +// << " max:" << gMax(magDisp) << endl; +// } +// +// +// Info<< "Calculated surface displacement in = " +// << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl; +// +// +// // Limit amount of movement. +// forAll(patchDisp, patchPointI) +// { +// scalar magDisp = mag(patchDisp[patchPointI]); +// +// if (magDisp > snapDist[patchPointI]) +// { +// patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp; +// +// Pout<< "Limiting displacement for " << patchPointI +// << " from " << magDisp << " to " << snapDist[patchPointI] +// << endl; +// } +// } +// +// // Points on zones in one domain but only present as point on other +// // will not do condition 2 on all. Sync explicitly. +// syncTools::syncPointList +// ( +// mesh, +// pp.meshPoints(), +// patchDisp, +// minMagSqrEqOp<point>(), // combine op +// vector(GREAT, GREAT, GREAT)// null value (note: cannot use VGREAT) +// ); +// +// return patchDisp; +//} +////XXXXXXXXX + void Foam::autoSnapDriver::smoothDisplacement ( const snapParameters& snapParams, @@ -1153,7 +1402,15 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::repatchToSurface scalarField faceSnapDist(pp.size(), -GREAT); { // Distance to attract to nearest feature on surface - const scalarField snapDist(calcSnapDistance(snapParams, pp)); + const scalarField snapDist + ( + calcSnapDistance + ( + mesh, + snapParams, + pp + ) + ); const faceList& localFaces = pp.localFaces(); @@ -1429,7 +1686,7 @@ void Foam::autoSnapDriver::doSnap indirectPrimitivePatch& pp = ppPtr(); // Distance to attract to nearest feature on surface - const scalarField snapDist(calcSnapDistance(snapParams, pp)); + const scalarField snapDist(calcSnapDistance(mesh, snapParams, pp)); // Construct iterative mesh mover. @@ -1467,7 +1724,14 @@ void Foam::autoSnapDriver::doSnap << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl; // Pre-smooth patch vertices (so before determining nearest) - preSmoothPatch(snapParams, nInitErrors, baffles, meshMover); + preSmoothPatch + ( + meshRefiner_, + snapParams, + nInitErrors, + baffles, + meshMover + ); for (label iter = 0; iter < nFeatIter; iter++) @@ -1478,7 +1742,12 @@ void Foam::autoSnapDriver::doSnap // Calculate displacement at every patch point. Insert into // meshMover. - vectorField disp = calcNearestSurface(snapDist, meshMover); + vectorField disp = calcNearestSurface + ( + meshRefiner_, + snapDist, + meshMover.patch() + ); // Override displacement with feature edge attempt if (doFeatures) diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H index e6763bc084850a9150d726e7a880544093755f59..f93ef62e081067520a252fab12810736219b86d4 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -81,11 +81,11 @@ class autoSnapDriver //- Calculate displacement per patch point to smooth out patch. // Quite complicated in determining which points to move where. - pointField smoothPatchDisplacement + static pointField smoothPatchDisplacement ( const motionSmoother&, const List<labelPair>& - ) const; + ); //- Check that face zones are synced void checkCoupledFaceZones() const; @@ -406,37 +406,51 @@ public: autoPtr<mapPolyMesh> mergeZoneBaffles(const List<labelPair>&); //- Calculate edge length per patch point. - scalarField calcSnapDistance + static scalarField calcSnapDistance ( + const fvMesh& mesh, const snapParameters& snapParams, const indirectPrimitivePatch& - ) const; + ); //- Smooth the mesh (patch and internal) to increase visibility // of surface points (on castellated mesh) w.r.t. surface. - void preSmoothPatch + static void preSmoothPatch ( + const meshRefinement& meshRefiner, const snapParameters& snapParams, const label nInitErrors, const List<labelPair>& baffles, motionSmoother& - ) const; + ); //- Get points both on patch and facezone. - labelList getZoneSurfacePoints + static labelList getZoneSurfacePoints ( + const fvMesh& mesh, const indirectPrimitivePatch&, const word& zoneName - ) const; + ); //- Per patch point calculate point on nearest surface. Set as // boundary conditions of motionSmoother displacement field. Return // displacement of patch points. - vectorField calcNearestSurface + static vectorField calcNearestSurface ( + const meshRefinement& meshRefiner, const scalarField& snapDist, - motionSmoother& meshMover - ) const; + const indirectPrimitivePatch& + ); + + ////- Per patch point calculate point on nearest surface. Set as + //// boundary conditions of motionSmoother displacement field. + //// Return displacement of patch points. + //static vectorField calcNearestLocalSurface + //( + // const meshRefinement& meshRefiner, + // const scalarField& snapDist, + // const indirectPrimitivePatch& + //); //- Smooth the displacement field to the internal. void smoothDisplacement diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C index 7bf8d8cf6c40f65f1bd2a6eea082d551e3cecb38..0e2fcd07ce7b088dba1a648903f38474869d293c 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C @@ -45,7 +45,11 @@ Foam::refinementParameters::refinementParameters nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))), keepPoints_(dict.lookup("keepPoints")), allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")), - maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0)) + useTopologicalSnapDetection_ + ( + dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true) + ), + maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0)) {} @@ -65,7 +69,11 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict) nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))), keepPoints_(pointField(1, dict.lookup("locationInMesh"))), allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")), - maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance",0)) + useTopologicalSnapDetection_ + ( + dict.lookupOrDefault<bool>("useTopologicalSnapDetection", true) + ), + maxLoadUnbalance_(dict.lookupOrDefault<scalar>("maxLoadUnbalance", 0)) { scalar featAngle(readScalar(dict.lookup("resolveFeatureAngle"))); @@ -83,7 +91,7 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // Foam::labelList Foam::refinementParameters::findCells(const polyMesh& mesh) - const +const { // Force calculation of tet-diag decomposition (for use in findCell) (void)mesh.tetBasePtIs(); diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H index 4fb05644e80b1748c7a39c59dfd9f76112943ba0..427fbd64b1ae6ffde6abe74caa91a7684c87a98e 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H @@ -80,6 +80,10 @@ class refinementParameters // cellZone? Switch allowFreeStandingZoneFaces_; + //- Use old topology based problem-cell removal (cells with 8 points + // on surface) + Switch useTopologicalSnapDetection_; + //- Allowed load unbalance scalar maxLoadUnbalance_; @@ -157,6 +161,13 @@ public: return allowFreeStandingZoneFaces_; } + //- Use old topology based problem-cell removal + // (cells with 8 points on surface) + bool useTopologicalSnapDetection() const + { + return useTopologicalSnapDetection_; + } + //- Allowed load unbalance scalar maxLoadUnbalance() const { diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H index 7487e8a2bf12b00072e8913133d01af81b5a0e51..09708a3d1a981672e7357f308c28d13a1d6c5cbf 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H @@ -72,6 +72,8 @@ class globalIndex; class removePoints; class localPointRegion; +class snapParameters; + /*---------------------------------------------------------------------------*\ Class meshRefinement Declaration \*---------------------------------------------------------------------------*/ @@ -470,6 +472,14 @@ private: const labelList& globalToMasterPatch ) const; + //- Returns list with for every internal face -1 or the patch + // they should be baffled into. + labelList markFacesOnProblemCellsGeometric + ( + const snapParameters& snapParams, + const dictionary& motionDict + ) const; + // Baffle merging @@ -536,6 +546,8 @@ private: //- Remove any loose standing cells void handleSnapProblems ( + const snapParameters& snapParams, + const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField& perpendicularAngle, const dictionary& motionDict, @@ -760,10 +772,17 @@ public: void baffleAndSplitMesh ( const bool handleSnapProblems, + + // How to remove problem snaps + const snapParameters& snapParams, + const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField& perpendicularAngle, + + // How to handle free-standing baffles const bool mergeFreeStanding, const scalar freeStandingAngle, + const dictionary& motionDict, Time& runTime, const labelList& globalToMasterPatch, diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C index 628dd7898ac9fb7327eb45f8df981152601b2794..6b9d1074c55438384b3e5af61b06033fe576df6c 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C @@ -1872,6 +1872,8 @@ void Foam::meshRefinement::makeConsistentFaceIndex void Foam::meshRefinement::handleSnapProblems ( + const snapParameters& snapParams, + const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField& perpendicularAngle, const dictionary& motionDict, @@ -1885,44 +1887,39 @@ void Foam::meshRefinement::handleSnapProblems << "----------------------------------------------" << nl << endl; - labelList facePatch - ( - markFacesOnProblemCells + labelList facePatch; + if (useTopologicalSnapDetection) + { + facePatch = markFacesOnProblemCells ( motionDict, removeEdgeConnectedCells, perpendicularAngle, globalToMasterPatch - ) - ); + ); + } + else + { + facePatch = markFacesOnProblemCellsGeometric(snapParams, motionDict); + } Info<< "Analyzed problem cells in = " << runTime.cpuTimeIncrement() << " s\n" << nl << endl; if (debug&meshRefinement::MESH) { - faceSet problemTopo(mesh_, "problemFacesTopo", 100); + faceSet problemFaces(mesh_, "problemFaces", 100); - const labelList facePatchTopo - ( - markFacesOnProblemCells - ( - motionDict, - removeEdgeConnectedCells, - perpendicularAngle, - globalToMasterPatch - ) - ); - forAll(facePatchTopo, faceI) + forAll(facePatch, faceI) { - if (facePatchTopo[faceI] != -1) + if (facePatch[faceI] != -1) { - problemTopo.insert(faceI); + problemFaces.insert(faceI); } } - problemTopo.instance() = timeName(); - Pout<< "Dumping " << problemTopo.size() - << " problem faces to " << problemTopo.objectPath() << endl; - problemTopo.write(); + problemFaces.instance() = timeName(); + Pout<< "Dumping " << problemFaces.size() + << " problem faces to " << problemFaces.objectPath() << endl; + problemFaces.write(); } Info<< "Introducing baffles to delete problem cells." << nl << endl; @@ -1962,6 +1959,8 @@ void Foam::meshRefinement::handleSnapProblems void Foam::meshRefinement::baffleAndSplitMesh ( const bool doHandleSnapProblems, + const snapParameters& snapParams, + const bool useTopologicalSnapDetection, const bool removeEdgeConnectedCells, const scalarField& perpendicularAngle, const bool mergeFreeStanding, @@ -2029,83 +2028,10 @@ void Foam::meshRefinement::baffleAndSplitMesh if (doHandleSnapProblems) { - //Info<< nl - // << "Introducing baffles to block off problem cells" << nl - // << "----------------------------------------------" << nl - // << endl; - // - //labelList facePatch - //( - // markFacesOnProblemCells - // ( - // motionDict, - // removeEdgeConnectedCells, - // perpendicularAngle, - // globalToMasterPatch - // ) - // //markFacesOnProblemCellsGeometric(motionDict) - //); - //Info<< "Analyzed problem cells in = " - // << runTime.cpuTimeIncrement() << " s\n" << nl << endl; - // - //if (debug&meshRefinement::MESH) - //{ - // faceSet problemTopo(mesh_, "problemFacesTopo", 100); - // - // const labelList facePatchTopo - // ( - // markFacesOnProblemCells - // ( - // motionDict, - // removeEdgeConnectedCells, - // perpendicularAngle, - // globalToMasterPatch - // ) - // ); - // forAll(facePatchTopo, faceI) - // { - // if (facePatchTopo[faceI] != -1) - // { - // problemTopo.insert(faceI); - // } - // } - // problemTopo.instance() = timeName(); - // Pout<< "Dumping " << problemTopo.size() - // << " problem faces to " << problemTopo.objectPath() << endl; - // problemTopo.write(); - //} - // - //Info<< "Introducing baffles to delete problem cells." << nl << endl; - // - //if (debug) - //{ - // runTime++; - //} - // - //// Create baffles with same owner and neighbour for now. - //createBaffles(facePatch, facePatch); - // - //if (debug) - //{ - // // Debug:test all is still synced across proc patches - // checkData(); - //} - //Info<< "Created baffles in = " - // << runTime.cpuTimeIncrement() << " s\n" << nl << endl; - // - //printMeshInfo(debug, "After introducing baffles"); - // - //if (debug&meshRefinement::MESH) - //{ - // Pout<< "Writing extra baffled mesh to time " - // << timeName() << endl; - // write(debug, runTime.path()/"extraBaffles"); - // Pout<< "Dumped debug data in = " - // << runTime.cpuTimeIncrement() << " s\n" << nl << endl; - //} - handleSnapProblems ( + snapParams, + useTopologicalSnapDetection, removeEdgeConnectedCells, perpendicularAngle, motionDict, @@ -2192,6 +2118,8 @@ void Foam::meshRefinement::baffleAndSplitMesh // and delete them handleSnapProblems ( + snapParams, + useTopologicalSnapDetection, removeEdgeConnectedCells, perpendicularAngle, motionDict, diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C index 8676847fd4124c4254747b35cf57d85bdd9e2f19..dcc268c238cbef33baac0706b021a64ea75374a8 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C @@ -36,6 +36,10 @@ License #include "polyMeshGeometry.H" #include "IOmanip.H" #include "unitConversion.H" +#include "autoSnapDriver.H" + +#include "snapParameters.H" +#include "motionSmoother.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -968,164 +972,241 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells } -//// Mark faces to be baffled to prevent snapping problems. Does -//// test to find nearest surface and checks which faces would get squashed. -//Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric -//( -// const dictionary& motionDict -//) const -//{ -// // Construct addressing engine. -// autoPtr<indirectPrimitivePatch> ppPtr -// ( -// meshRefinement::makePatch -// ( -// mesh_, -// meshedPatches() -// ) -// ); -// const indirectPrimitivePatch& pp = ppPtr(); -// const pointField& localPoints = pp.localPoints(); -// const labelList& meshPoints = pp.meshPoints(); -// -// // Find nearest (non-baffle) surface -// pointField newPoints(mesh_.points()); -// { -// List<pointIndexHit> hitInfo; -// labelList hitSurface; -// surfaces_.findNearest -// ( -// surfaces_.getUnnamedSurfaces(), -// localPoints, -// scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction -// hitSurface, -// hitInfo -// ); -// -// forAll(hitInfo, i) -// { -// if (hitInfo[i].hit()) -// { -// //label pointI = meshPoints[i]; -// //Pout<< " " << pointI << " moved from " -// // << mesh_.points()[pointI] << " by " -// // << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI]) -// // << endl; -// newPoints[meshPoints[i]] = hitInfo[i].hitPoint(); -// } -// } -// } -// -// // Per face (internal or coupled!) the patch that the -// // baffle should get (or -1). -// labelList facePatch(mesh_.nFaces(), -1); -// // Count of baffled faces -// label nBaffleFaces = 0; -// -// { -// pointField oldPoints(mesh_.points()); -// mesh_.movePoints(newPoints); -// faceSet wrongFaces(mesh_, "wrongFaces", 100); -// { -// //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces); -// -// // Just check the errors from squashing -// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -// -// const labelList allFaces(identity(mesh_.nFaces())); -// label nWrongFaces = 0; -// -// scalar minArea(readScalar(motionDict.lookup("minArea"))); -// if (minArea > -SMALL) -// { -// polyMeshGeometry::checkFaceArea -// ( -// false, -// minArea, -// mesh_, -// mesh_.faceAreas(), -// allFaces, -// &wrongFaces -// ); -// -// label nNewWrongFaces = returnReduce -// ( -// wrongFaces.size(), -// sumOp<label>() -// ); -// -// Info<< " faces with area < " -// << setw(5) << minArea -// << " m^2 : " -// << nNewWrongFaces-nWrongFaces << endl; -// -// nWrongFaces = nNewWrongFaces; -// } -// -//// scalar minDet(readScalar(motionDict.lookup("minDeterminant"))); -// scalar minDet = 0.01; -// if (minDet > -1) -// { -// polyMeshGeometry::checkCellDeterminant -// ( -// false, -// minDet, -// mesh_, -// mesh_.faceAreas(), -// allFaces, -// polyMeshGeometry::affectedCells(mesh_, allFaces), -// &wrongFaces -// ); -// -// label nNewWrongFaces = returnReduce -// ( -// wrongFaces.size(), -// sumOp<label>() -// ); -// -// Info<< " faces on cells with determinant < " -// << setw(5) << minDet << " : " -// << nNewWrongFaces-nWrongFaces << endl; -// -// nWrongFaces = nNewWrongFaces; -// } -// } -// -// -// forAllConstIter(faceSet, wrongFaces, iter) -// { -// label patchI = mesh_.boundaryMesh().whichPatch(iter.key()); -// -// if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled()) -// { -// facePatch[iter.key()] = getBafflePatch(facePatch, iter.key()); -// nBaffleFaces++; -// -// //Pout<< " " << iter.key() -// // //<< " on patch " << mesh_.boundaryMesh()[patchI].name() -// // << " is destined for patch " << facePatch[iter.key()] -// // << endl; -// } -// } -// // Restore points. -// mesh_.movePoints(oldPoints); -// } -// -// -// Info<< "markFacesOnProblemCellsGeometric : marked " -// << returnReduce(nBaffleFaces, sumOp<label>()) -// << " additional internal and coupled faces" -// << " to be converted into baffles." << endl; -// -// syncTools::syncFaceList -// ( -// mesh_, -// facePatch, -// maxEqOp<label>() -// ); -// -// return facePatch; -//} +// Mark faces to be baffled to prevent snapping problems. Does +// test to find nearest surface and checks which faces would get squashed. +Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric +( + const snapParameters& snapParams, + const dictionary& motionDict +) const +{ + pointField oldPoints(mesh_.points()); + + // Repeat (most of) autoSnapDriver::doSnap + { + labelList adaptPatchIDs(meshedPatches()); + + // Construct addressing engine. + autoPtr<indirectPrimitivePatch> ppPtr + ( + meshRefinement::makePatch + ( + mesh_, + adaptPatchIDs + ) + ); + indirectPrimitivePatch& pp = ppPtr(); + + // Distance to attract to nearest feature on surface + const scalarField snapDist + ( + autoSnapDriver::calcSnapDistance(mesh_, snapParams, pp) + ); + + + // Construct iterative mesh mover. + Info<< "Constructing mesh displacer ..." << endl; + Info<< "Using mesh parameters " << motionDict << nl << endl; + + const pointMesh& pMesh = pointMesh::New(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<label>() + ); + + Info<< "Detected " << nInitErrors << " illegal faces" + << " (concave, zero area or negative cell pyramid volume)" + << endl; + + + Info<< "Checked initial mesh in = " + << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl; + + // Pre-smooth patch vertices (so before determining nearest) + autoSnapDriver::preSmoothPatch + ( + *this, + snapParams, + nInitErrors, + List<labelPair>(0), //baffles + meshMover + ); + + const vectorField disp + ( + autoSnapDriver::calcNearestSurface + ( + *this, + snapDist, // attraction + pp + ) + ); + + const labelList& meshPoints = pp.meshPoints(); + + pointField newPoints(mesh_.points()); + forAll(meshPoints, i) + { + newPoints[meshPoints[i]] += disp[i]; + } + mesh_.movePoints(newPoints); + } + + + // Per face (internal or coupled!) the patch that the + // baffle should get (or -1). + labelList facePatch(mesh_.nFaces(), -1); + // Count of baffled faces + label nBaffleFaces = 0; + + { + faceSet wrongFaces(mesh_, "wrongFaces", 100); + { + //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces); + + // Just check the errors from squashing + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + const labelList allFaces(identity(mesh_.nFaces())); + label nWrongFaces = 0; + + //const scalar minV(readScalar(motionDict.lookup("minVol", true))); + //if (minV > -GREAT) + //{ + // polyMeshGeometry::checkFacePyramids + // ( + // false, + // minV, + // mesh_, + // mesh_.cellCentres(), + // mesh_.points(), + // allFaces, + // List<labelPair>(0), + // &wrongFaces + // ); + // + // label nNewWrongFaces = returnReduce + // ( + // wrongFaces.size(), + // sumOp<label>() + // ); + // + // Info<< " faces with pyramid volume < " + // << setw(5) << minV + // << " m^3 : " + // << nNewWrongFaces-nWrongFaces << endl; + // + // nWrongFaces = nNewWrongFaces; + //} + + scalar minArea(readScalar(motionDict.lookup("minArea"))); + if (minArea > -SMALL) + { + polyMeshGeometry::checkFaceArea + ( + false, + minArea, + mesh_, + mesh_.faceAreas(), + allFaces, + &wrongFaces + ); + + label nNewWrongFaces = returnReduce + ( + wrongFaces.size(), + sumOp<label>() + ); + + Info<< " faces with area < " + << setw(5) << minArea + << " m^2 : " + << nNewWrongFaces-nWrongFaces << endl; + + nWrongFaces = nNewWrongFaces; + } + + scalar minDet(readScalar(motionDict.lookup("minDeterminant"))); + if (minDet > -1) + { + polyMeshGeometry::checkCellDeterminant + ( + false, + minDet, + mesh_, + mesh_.faceAreas(), + allFaces, + polyMeshGeometry::affectedCells(mesh_, allFaces), + &wrongFaces + ); + + label nNewWrongFaces = returnReduce + ( + wrongFaces.size(), + sumOp<label>() + ); + + Info<< " faces on cells with determinant < " + << setw(5) << minDet << " : " + << nNewWrongFaces-nWrongFaces << endl; + + nWrongFaces = nNewWrongFaces; + } + } + + + forAllConstIter(faceSet, wrongFaces, iter) + { + label patchI = mesh_.boundaryMesh().whichPatch(iter.key()); + + if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled()) + { + facePatch[iter.key()] = getBafflePatch(facePatch, iter.key()); + nBaffleFaces++; + + //Pout<< " " << iter.key() + // //<< " on patch " << mesh_.boundaryMesh()[patchI].name() + // << " is destined for patch " << facePatch[iter.key()] + // << endl; + } + } + } + + + // Restore points. + mesh_.movePoints(oldPoints); + + + Info<< "markFacesOnProblemCellsGeometric : marked " + << returnReduce(nBaffleFaces, sumOp<label>()) + << " additional internal and coupled faces" + << " to be converted into baffles." << endl; + + syncTools::syncFaceList + ( + mesh_, + facePatch, + maxEqOp<label>() + ); + + return facePatch; +} // ************************************************************************* // diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C index f1917ee6f94a14a5f8e9462319a0fbfdd12b4425..93c4600f9cb0b04fa3df4d7ea80d90932e28124f 100644 --- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C +++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C @@ -77,7 +77,8 @@ const Foam::NamedEnum<Foam::refinementSurfaces::faceZoneType, 3> Foam::refinementSurfaces::refinementSurfaces ( const searchableSurfaces& allGeometry, - const dictionary& surfacesDict + const dictionary& surfacesDict, + const label gapLevelIncrement ) : allGeometry_(allGeometry), @@ -143,7 +144,7 @@ Foam::refinementSurfaces::refinementSurfaces globalLevelIncr[surfI] = dict.lookupOrDefault ( "gapLevelIncrement", - 0 + gapLevelIncrement ); if @@ -274,7 +275,7 @@ Foam::refinementSurfaces::refinementSurfaces label levelIncr = regionDict.lookupOrDefault ( "gapLevelIncrement", - 0 + gapLevelIncrement ); regionLevelIncr[surfI].insert(regionI, levelIncr); diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H index c9773398e5f9734c31caf8644b32434f605ff91d..05ecdc8323ac0e10a985d4129ad4a29db900e1e5 100644 --- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H +++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H @@ -149,7 +149,8 @@ public: refinementSurfaces ( const searchableSurfaces& allGeometry, - const dictionary& + const dictionary&, + const label gapLevelIncrement ); diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/0.org/U b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/0.org/U index b82b7ce18aeae391934094cf3a54196467dce9dc..7ea3a0c32328d93fd01ba3c580c7d73b54473d22 100644 --- a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/0.org/U +++ b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phase/0.org/U @@ -38,7 +38,7 @@ boundaryField } atmosphere { - type fluxCorrectedVelocity; + type pressureInletOutletVelocity; value uniform (0 0 0); } defaultFaces diff --git a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/0.org/U b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/0.org/U index b82b7ce18aeae391934094cf3a54196467dce9dc..7ea3a0c32328d93fd01ba3c580c7d73b54473d22 100644 --- a/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/0.org/U +++ b/tutorials/multiphase/multiphaseInterFoam/laminar/damBreak4phaseFine/0.org/U @@ -38,7 +38,7 @@ boundaryField } atmosphere { - type fluxCorrectedVelocity; + type pressureInletOutletVelocity; value uniform (0 0 0); } defaultFaces