diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C index f1ce1aa996..b34172130d 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C @@ -65,7 +65,8 @@ autoPtr createRefinementSurfaces ( const searchableSurfaces& allGeometry, const dictionary& surfacesDict, - const dictionary& shapeControlDict + const dictionary& shapeControlDict, + const label gapLevelIncrement ) { autoPtr surfacePtr; @@ -145,7 +146,7 @@ autoPtr createRefinementSurfaces globalLevelIncr[surfI] = shapeDict.lookupOrDefault ( "gapLevelIncrement", - 0 + gapLevelIncrement ); } else @@ -356,7 +357,7 @@ autoPtr createRefinementSurfaces label levelIncr = regionDict.lookupOrDefault ( "gapLevelIncrement", - 0 + gapLevelIncrement ); regionLevelIncr[surfI].insert(regionI, levelIncr); @@ -1031,7 +1032,8 @@ int main(int argc, char *argv[]) ( allGeometry, conformationDict, - shapeControlDict + shapeControlDict, + refineDict.lookupOrDefault("gapLevelIncrement", 0) ); } else @@ -1041,7 +1043,8 @@ int main(int argc, char *argv[]) new refinementSurfaces ( allGeometry, - refineDict.subDict("refinementSurfaces") + refineDict.subDict("refinementSurfaces"), + refineDict.lookupOrDefault("gapLevelIncrement", 0) ) ); @@ -1052,7 +1055,6 @@ int main(int argc, char *argv[]) refinementSurfaces& surfaces = surfacesPtr(); - // Checking only? if (checkGeometry) @@ -1351,6 +1353,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; @@ -1364,15 +1373,20 @@ int main(int argc, char *argv[]) globalToSlavePatch ); - // Refinement parameters - refinementParameters refineParams(refineDict); if (!overwrite && !debug) { const_cast(mesh.time())++; } - refineDriver.doRefine(refineDict, refineParams, wantSnap, motionDict); + refineDriver.doRefine + ( + refineDict, + refineParams, + snapParams, + wantSnap, + motionDict + ); writeMesh ( @@ -1396,20 +1410,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(mesh.time())++; } + // Use the resolveFeatureAngle from the refinement parameters + scalar curvature = refineParams.curvature(); + snapDriver.doSnap(snapDict, motionDict, curvature, snapParams); writeMesh @@ -1437,17 +1445,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() - ); - } + // Use the maxLocalCells from the refinement parameters + bool preBalance = returnReduce + ( + (mesh.nCells() >= refineParams.maxLocalCells()), + orOp() + ); if (!overwrite && !debug) diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index ca4a9dc6b5..ba16b1c4aa 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -152,6 +152,10 @@ castellatedMeshControls inGroups (meshedPatches); } + + //- Optional increment (on top of max level) in small gaps + //gapLevelIncrement 2; + //- Optional angle to detect small-large cell situation // perpendicular to the surface. Is the angle of face w.r.t. // the local surface normal. Use on flat(ish) surfaces only. @@ -180,11 +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 neighbouring surface normals - // are roughly the same so e.g. free-standing baffles can be merged + // - used to determine if surface normals + // 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 30; + planarAngle 30; // Region-wise refinement @@ -229,10 +239,8 @@ castellatedMeshControls // free-standing zone faces. Not used if there are no faceZones. allowFreeStandingZoneFaces true; - - // Optional: whether all baffles get eroded away. WIP. Used for - // surface simplification. - //allowFreeStandingBaffles false; + // + //useTopologicalSnapDetection false; } // Settings for the snapping. @@ -462,20 +470,21 @@ meshQualityControls // 1 = hex, <= 0 = folded or flattened illegal cell minDeterminant 0.001; - // minFaceWeight (0 -> 0.5) + // Relative position of face in relation to cell centres (0.5 for orthogonal + // mesh) (0 -> 0.5) minFaceWeight 0.05; - // minVolRatio (0 -> 1) + // Volume ratio of neighbouring cells (0 -> 1) minVolRatio 0.01; // must be >0 for Fluent compatibility minTriangleTwist -1; - //- if >0 : preserve single cells with all points on the surface if the + //- if >0 : preserve cells with all points on the surface if the // resulting volume after snapping (by approximation) is larger than // minVolCollapseRatio times old volume (i.e. not collapsed to flat cell). // If <0 : delete always. - //minVolCollapseRatio 0.5; + //minVolCollapseRatio 0.1; // Advanced diff --git a/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C b/applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C index 74c5fe0bf6..c7ad881527 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 75af8eb20a..dd51eb645b 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/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C index 3af85508be..4945eb0b09 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C @@ -97,7 +97,6 @@ Foam::tmp Foam::polyMeshTools::faceSkewness { const labelList& own = mesh.faceOwner(); const labelList& nei = mesh.faceNeighbour(); - const faceList& fcs = mesh.faces(); const polyBoundaryMesh& pbm = mesh.boundaryMesh(); tmp tskew(new scalarField(mesh.nFaces())); @@ -154,31 +153,16 @@ Foam::tmp Foam::polyMeshTools::faceSkewness { label faceI = pp.start() + i; - vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]]; + skew[faceI] = primitiveMeshTools::boundaryFaceSkewness + ( + mesh, + p, + fCtrs, + fAreas, - vector normal = fAreas[faceI]; - normal /= mag(normal) + VSMALL; - vector d = normal*(normal & Cpf); - - - // Skewness vector - vector sv = - Cpf - - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d; - vector svHat = sv/(mag(sv) + VSMALL); - - // Normalisation distance calculated as the approximate distance - // from the face centre to the edge of the face in the direction - // of the skewness - scalar fd = 0.4*mag(d) + VSMALL; - const face& f = fcs[faceI]; - forAll(f, pi) - { - fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI]))); - } - - // Normalised skewness - skew[faceI] = mag(sv)/fd; + faceI, + cellCtrs[own[faceI]] + ); } } } diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C index e94c8edbf7..241f4eef8e 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -64,6 +64,44 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness return mag(sv)/fd; } +Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness +( + const primitiveMesh& mesh, + const pointField& p, + const vectorField& fCtrs, + const vectorField& fAreas, + + const label faceI, + const point& ownCc +) +{ + vector Cpf = fCtrs[faceI] - ownCc; + + vector normal = fAreas[faceI]; + normal /= mag(normal) + VSMALL; + vector d = normal*(normal & Cpf); + + + // Skewness vector + vector sv = + Cpf + - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d; + vector svHat = sv/(mag(sv) + VSMALL); + + // Normalisation distance calculated as the approximate distance + // from the face centre to the edge of the face in the direction + // of the skewness + scalar fd = 0.4*mag(d) + VSMALL; + const face& f = mesh.faces()[faceI]; + forAll(f, pi) + { + fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI]))); + } + + // Normalised skewness + return mag(sv)/fd; +} + Foam::scalar Foam::primitiveMeshTools::faceOrthogonality ( @@ -119,7 +157,6 @@ Foam::tmp Foam::primitiveMeshTools::faceSkewness { const labelList& own = mesh.faceOwner(); const labelList& nei = mesh.faceNeighbour(); - const faceList& fcs = mesh.faces(); tmp tskew(new scalarField(mesh.nFaces())); scalarField& skew = tskew(); @@ -145,31 +182,15 @@ Foam::tmp Foam::primitiveMeshTools::faceSkewness for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++) { - vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]]; - - vector normal = fAreas[faceI]; - normal /= mag(normal) + VSMALL; - vector d = normal*(normal & Cpf); - - - // Skewness vector - vector sv = - Cpf - - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d; - vector svHat = sv/(mag(sv) + VSMALL); - - // Normalisation distance calculated as the approximate distance - // from the face centre to the edge of the face in the direction - // of the skewness - scalar fd = 0.4*mag(d) + VSMALL; - const face& f = fcs[faceI]; - forAll(f, pi) - { - fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI]))); - } - - // Normalised skewness - skew[faceI] = mag(sv)/fd; + skew[faceI] = boundaryFaceSkewness + ( + mesh, + p, + fCtrs, + fAreas, + faceI, + cellCtrs[own[faceI]] + ); } return tskew; diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H index 5e3b2c07fd..170d089e11 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -48,32 +48,6 @@ namespace Foam class primitiveMeshTools { - -protected: - - // Protected Member Functions - - //- Skewness of single face - static scalar faceSkewness - ( - const primitiveMesh& mesh, - const pointField& p, - const vectorField& fCtrs, - const vectorField& fAreas, - - const label faceI, - const point& ownCc, - const point& neiCc - ); - - //- Orthogonality of single face - static scalar faceOrthogonality - ( - const point& ownCc, - const point& neiCc, - const vector& s - ); - public: //- Generate non-orthogonality field (internal faces only) @@ -154,6 +128,41 @@ public: const PackedBoolList& internalOrCoupledFace ); + + // Helpers: single face check + + //- Skewness of single face + static scalar faceSkewness + ( + const primitiveMesh& mesh, + const pointField& p, + const vectorField& fCtrs, + const vectorField& fAreas, + + const label faceI, + const point& ownCc, + const point& neiCc + ); + + //- Skewness of single boundary face + static scalar boundaryFaceSkewness + ( + const primitiveMesh& mesh, + const pointField& p, + const vectorField& fCtrs, + const vectorField& fAreas, + + const label faceI, + const point& ownCc + ); + + //- Orthogonality of single face + static scalar faceOrthogonality + ( + const point& ownCc, + const point& neiCc, + const vector& s + ); }; diff --git a/src/dynamicMesh/motionSmoother/motionSmootherCheck.C b/src/dynamicMesh/motionSmoother/motionSmootherCheck.C index 77ce89ca71..3efb62a995 100644 --- a/src/dynamicMesh/motionSmoother/motionSmootherCheck.C +++ b/src/dynamicMesh/motionSmoother/motionSmootherCheck.C @@ -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 @@ -239,6 +239,7 @@ bool Foam::motionSmoother::checkMesh maxIntSkew, maxBounSkew, mesh, + mesh.points(), mesh.cellCentres(), mesh.faceCentres(), mesh.faceAreas(), @@ -481,14 +482,14 @@ bool Foam::motionSmoother::checkMesh ( readScalar(dict.lookup("minArea", true)) ); - const scalar maxIntSkew - ( - readScalar(dict.lookup("maxInternalSkewness", true)) - ); - const scalar maxBounSkew - ( - readScalar(dict.lookup("maxBoundarySkewness", true)) - ); + //const scalar maxIntSkew + //( + // readScalar(dict.lookup("maxInternalSkewness", true)) + //); + //const scalar maxBounSkew + //( + // readScalar(dict.lookup("maxBoundarySkewness", true)) + //); const scalar minWeight ( readScalar(dict.lookup("minFaceWeight", true)) @@ -615,27 +616,30 @@ bool Foam::motionSmoother::checkMesh nWrongFaces = nNewWrongFaces; } - if (maxIntSkew > 0 || maxBounSkew > 0) - { - meshGeom.checkFaceSkewness - ( - report, - maxIntSkew, - maxBounSkew, - checkFaces, - baffles, - &wrongFaces - ); - label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp