diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index 3a65e3491a..ebd3c3f849 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -299,6 +299,12 @@ castellatedMeshControls // // with e.g. dynamic refinement/unrefinement. // // - cellLevel will include any directional refinement // // (i.e. it will be the maximum of all three directions) + // + // // Optional directional expansion-ratio smoothing (after all + // // refinement). This will try to smooth out edge/cell size jumps + // // Specify smoothing direction and number of iterations + // nSmoothExpansion 100; + // smoothDirection (1 0 0); //} } diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C index d721b34170..81e52c58f6 100644 --- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C +++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C @@ -585,6 +585,8 @@ Foam::shellSurfaces::shellSurfaces distances_.setSize(shellI); levels_.setSize(shellI); dirLevels_.setSize(shellI); + smoothDirection_.setSize(shellI); + nSmoothExpansion_.setSize(shellI); extendedGapLevel_.setSize(shellI); extendedGapMode_.setSize(shellI); @@ -661,6 +663,13 @@ Foam::shellSurfaces::shellSurfaces } } + // Directional smoothing + // ~~~~~~~~~~~~~~~~~~~~~ + + smoothDirection_[shellI] = + dict.lookupOrDefault("smoothDirection", vector::zero); + nSmoothExpansion_[shellI] = + dict.lookupOrDefault("nSmoothExpansion", 0); // Gap specification @@ -803,6 +812,18 @@ Foam::labelPairList Foam::shellSurfaces::directionalSelectLevel() const } +const Foam::labelList& Foam::shellSurfaces::nSmoothExpansion() const +{ + return nSmoothExpansion_; +} + + +const Foam::vectorField& Foam::shellSurfaces::smoothDirection() const +{ + return smoothDirection_; +} + + void Foam::shellSurfaces::findHigherLevel ( const pointField& pt, diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H index de3edb7c02..8b4f8a8ef9 100644 --- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H +++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H @@ -86,8 +86,17 @@ private: //- Per shell per distance the refinement level labelListList levels_; - //- Per shell any additional directional refinement - List> dirLevels_; + + // Directional + + //- Per shell any additional directional refinement + List> dirLevels_; + + //- Per shell the smoothing direction + vectorField smoothDirection_; + + //- Per shell the directional smoothing iterations + labelList nSmoothExpansion_; // Gap level refinement @@ -231,6 +240,12 @@ public: const direction dir, labelList& shell ) const; + + //- Per shell the smoothing direction + const vectorField& smoothDirection() const; + + //- Per shell the directional smoothing iterations + const labelList& nSmoothExpansion() const; }; diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/PointIntegrateData/PointIntegrateData.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/PointIntegrateData/PointIntegrateData.H new file mode 100644 index 0000000000..a2d32ba9ed --- /dev/null +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/PointIntegrateData/PointIntegrateData.H @@ -0,0 +1,258 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::PointIntegrateData + +Description + Integrate along selected edges using PointEdgeWave. + +SourceFiles + PointIntegrateDataI.H + +\*---------------------------------------------------------------------------*/ + +#ifndef PointIntegrateData_H +#define PointIntegrateData_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class Istream; +class Ostream; +template +class PointIntegrateData; + +// Forward declaration of friend functions and operators +template +Ostream& operator<<(Ostream&, const PointIntegrateData&); +template +Istream& operator>>(Istream&, PointIntegrateData&); + +/*---------------------------------------------------------------------------*\ + Class PointIntegrateData declaration +\*---------------------------------------------------------------------------*/ + +template +class PointIntegrateData +{ +private: + + // Private data + + //- Valid flag + bool valid_; + + //- Integrated data + DataType data_; + + +public: + + //- Class used to pass extra data + class trackingData + { + public: + + UList& edgeData_; + + trackingData(UList& edgeData) + : + edgeData_(edgeData) + {} + }; + + + + // Constructors + + //- Construct null + inline PointIntegrateData(); + + //- Construct from data + inline PointIntegrateData(const DataType& data); + + + // Member Functions + + // Access + + //- Const access the data + inline const DataType& data() const + { + return data_; + }; + + + // Needed by PointEdgeWave + + //- Check whether original (invalid) value. + template + inline bool valid(TrackingData& td) const; + + //- Check for identical geometrical data. Used for cyclics checking. + template + inline bool sameGeometry + ( + const PointIntegrateData&, + const scalar tol, + TrackingData& td + ) const; + + //- Convert origin to relative vector to leaving point + // (= point coordinate) + template + inline void leaveDomain + ( + const polyPatch& patch, + const label patchPointi, + const point& pos, + TrackingData& td + ); + + //- Convert relative origin to absolute by adding entering point + template + inline void enterDomain + ( + const polyPatch& patch, + const label patchPointi, + const point& pos, + TrackingData& td + ); + + //- Apply rotation matrix to the data + template + inline void transform + ( + const tensor& rotTensor, + TrackingData& td + ); + + //- Influence of edge on point + template + inline bool updatePoint + ( + const polyMesh& mesh, + const label pointI, + const label edgeI, + const PointIntegrateData& edgeInfo, + const scalar tol, + TrackingData& td + ); + + //- Influence of different value on same point. + // Merge new and old info. + template + inline bool updatePoint + ( + const polyMesh& mesh, + const label pointI, + const PointIntegrateData& newPointInfo, + const scalar tol, + TrackingData& td + ); + + //- Influence of different value on same point. + // No information about current position whatsoever. + template + inline bool updatePoint + ( + const PointIntegrateData& newPointInfo, + const scalar tol, + TrackingData& td + ); + + //- Influence of point on edge. + template + inline bool updateEdge + ( + const polyMesh& mesh, + const label edgeI, + const label pointI, + const PointIntegrateData& pointInfo, + const scalar tol, + TrackingData& td + ); + + //- Same (like operator==) + template + inline bool equal + ( + const PointIntegrateData&, + TrackingData& td + ) const; + + + // Member Operators + + // Needed for List IO + inline bool operator==(const PointIntegrateData&) const; + inline bool operator!=(const PointIntegrateData&) const; + + + // IOstream Operators + + friend Ostream& operator<< + ( + Ostream&, + const PointIntegrateData& + ); + friend Istream& operator>> + ( + Istream&, + PointIntegrateData& + ); +}; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Data associated with PointIntegrateData types is contiguous + +template<> +inline bool contiguous>() +{ + return true; +} + +template<> +inline bool contiguous>() +{ + return true; +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "PointIntegrateDataI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/PointIntegrateData/PointIntegrateDataI.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/PointIntegrateData/PointIntegrateDataI.H new file mode 100644 index 0000000000..c22cdac898 --- /dev/null +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/PointIntegrateData/PointIntegrateDataI.H @@ -0,0 +1,304 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "polyMesh.H" +#include "transform.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +inline Foam::PointIntegrateData::PointIntegrateData() +: + valid_(false) +{} + + +template +inline Foam::PointIntegrateData::PointIntegrateData +( + const DataType& data +) +: + valid_(true), + data_(data) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +template +inline bool Foam::PointIntegrateData::valid(TrackingData& td) const +{ + return valid_; +} + + +template +template +inline bool Foam::PointIntegrateData::sameGeometry +( + const PointIntegrateData&, + const scalar tol, + TrackingData& td +) const +{ + return true; +} + + +template +template +inline void Foam::PointIntegrateData::leaveDomain +( + const polyPatch& patch, + const label patchPointi, + const point& pos, + TrackingData& td +) +{} + + +template +template +inline void Foam::PointIntegrateData::enterDomain +( + const polyPatch& patch, + const label patchPointi, + const point& pos, + TrackingData& td +) +{} + + +template +template +inline void Foam::PointIntegrateData::transform +( + const tensor& rotTensor, + TrackingData& td +) +{ + this->data_ = Foam::transform(rotTensor, this->data_); +} + + +template +template +inline bool Foam::PointIntegrateData::updatePoint +( + const polyMesh& mesh, + const label pointI, + const label edgeI, + const PointIntegrateData& edgeInfo, + const scalar tol, + TrackingData& td +) +{ + // Update point from an edge + if (!valid_) + { + if (!edgeInfo.valid_) + { + FatalErrorInFunction<< "edgeInfo:" << edgeInfo << exit(FatalError); + } + this->operator=(edgeInfo); + return true; + } + else + { + return false; + } +} + + +template +template +inline bool Foam::PointIntegrateData::updatePoint +( + const polyMesh& mesh, + const label pointI, + const PointIntegrateData& newPointInfo, + const scalar tol, + TrackingData& td +) +{ + // Update point from coupled point + if (!valid_) + { + if (!newPointInfo.valid_) + { + FatalErrorInFunction<< "newPointInfo:" << newPointInfo + << exit(FatalError); + } + this->operator=(newPointInfo); + return true; + } + else + { + return false; + } +} + + +template +template +inline bool Foam::PointIntegrateData::updatePoint +( + const PointIntegrateData& newPointInfo, + const scalar tol, + TrackingData& td +) +{ + if (!valid_) + { + if (!newPointInfo.valid_) + { + FatalErrorInFunction<< "newPointInfo:" << newPointInfo + << exit(FatalError); + } + this->operator=(newPointInfo); + return true; + } + else + { + return false; + } +} + + +template +template +inline bool Foam::PointIntegrateData::updateEdge +( + const polyMesh& mesh, + const label edgeI, + const label pointI, + const PointIntegrateData& pointInfo, + const scalar tol, + TrackingData& td +) +{ + if (!valid_) + { + if (!pointInfo.valid(td)) + { + FatalErrorInFunction<< "problem: invalid point:" + << mesh.points()[pointI] + << " data:" << pointInfo + << exit(FatalError); + } + this->data_ = pointInfo.data_ + td.edgeData_[edgeI]; + this->valid_ = true; + return true; + } + else + { + return false; + } +} + +//- Same (like operator==) +template +template +inline bool Foam::PointIntegrateData::equal +( + const PointIntegrateData& pi, + TrackingData& td +) const +{ + if (!valid_) + { + return false; + } + else if (!pi.valid_) + { + FatalErrorInFunction << "pi:" << pi + << exit(FatalError); + return false; + } + else + { + return this->data_ == pi.data_; + } +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +template +inline bool Foam::PointIntegrateData::operator== +( + const Foam::PointIntegrateData& rhs +) +const +{ + return this->data_ == rhs.data_; +} + + +template +inline bool Foam::PointIntegrateData::operator!= +( + const Foam::PointIntegrateData& rhs +) +const +{ + return !(*this == rhs); +} + + +// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // + +template +inline Foam::Ostream& Foam::operator<< +( + Ostream& os, + const PointIntegrateData& pd +) +{ + if (os.format() == IOstream::ASCII) + { + return os << pd.valid_ << token::SPACE << pd.data(); + } + else + { + return os << pd.valid_ << pd.data(); + } +} + + +template +inline Foam::Istream& Foam::operator>> +( + Istream& is, + PointIntegrateData& pd +) +{ + return is >> pd.valid_ >> pd.data_; +} + + +// ************************************************************************* // diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C index fd33ec7b6e..7b4db14fec 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C @@ -668,7 +668,7 @@ void Foam::snappyLayerDriver::handleFeatureAngle ( const indirectPrimitivePatch& pp, const labelList& meshEdges, - const scalar minCos, + const scalar minAngle, pointField& patchDisp, labelList& patchNLayers, List& extrudeStatus @@ -676,7 +676,10 @@ void Foam::snappyLayerDriver::handleFeatureAngle { const fvMesh& mesh = meshRefiner_.mesh(); - Info<< nl << "Handling feature edges ..." << endl; + const scalar minCos = Foam::cos(degToRad(minAngle)); + + Info<< nl << "Handling feature edges (angle < " << minAngle + << ") ..." << endl; if (minCos < 1-SMALL) { @@ -3647,7 +3650,7 @@ void Foam::snappyLayerDriver::addLayers ( pp, meshEdges, - degToRad(layerParams.featureAngle()), + layerParams.featureAngle(), patchDisp, patchNLayers, diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H index f015f01e0f..e7c95c953d 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H @@ -218,7 +218,7 @@ private: ( const indirectPrimitivePatch& pp, const labelList& meshEdges, - const scalar minCos, + const scalar minAngle, pointField& patchDisp, labelList& patchNLayers, List& extrudeStatus diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C index 0bc07420d2..21ca736aa3 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2015-2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2015-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -41,6 +41,9 @@ License #include "labelVector.H" #include "profiling.H" #include "searchableSurfaces.H" +#include "PointIntegrateData.H" +#include "PointEdgeWave.H" +#include "OBJstream.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -1627,7 +1630,6 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine for (direction dir = 0; dir < vector::nComponents; dir++) { // Select the cells that need to be refined in certain direction: - // // - cell inside/outside shell // - original cellLevel (using mapping) mentioned in levelIncrement // - dirCellLevel not yet up to cellLevel+levelIncrement @@ -1786,6 +1788,576 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine } +void Foam::snappyRefineDriver::getVirtualEdgeLength +( + const vector& userDirection, + const labelList& pointLabels, + const PackedBoolList& isXFace, + scalarField& maxUserSize +) const +{ + const fvMesh& mesh = meshRefiner_.mesh(); + const labelList& own = mesh.faceOwner(); + const labelList& nei = mesh.faceNeighbour(); + const vectorField& fA = mesh.faceAreas(); + const pointField& faceCentres = mesh.faceCentres(); + const pointField& cellCentres = mesh.cellCentres(); + + scalarField sumWeights(mesh.nCells(), 0.0); + scalarField sumArea(mesh.nCells(), 0.0); + + forAll(own, facei) + { + const scalar fc = (faceCentres[facei]&userDirection); + const scalar cc = (cellCentres[own[facei]]&userDirection); + // Calculate 'pyramid volume' + scalar pyr3Vol = (fA[facei]&userDirection)*(fc-cc); + sumWeights[own[facei]] += pyr3Vol; + sumArea[own[facei]] += mag(fA[facei]&userDirection); + } + + forAll(nei, facei) + { + const scalar fc = (faceCentres[facei]&userDirection); + const scalar cc = (cellCentres[nei[facei]]&userDirection); + // Calculate 'pyramid volume' + scalar pyr3Vol = (fA[facei]&userDirection)*(cc-fc); + sumWeights[nei[facei]] += pyr3Vol; + sumArea[nei[facei]] += mag(fA[facei]&userDirection); + } + + const scalarField cellSize(2*sumWeights/sumArea); + + maxUserSize.setSize(mesh.nPoints()); + maxUserSize = -GREAT; + + forAll(pointLabels, i) + { + label pointi = pointLabels[i]; + const labelList& pFaces = mesh.pointFaces()[pointi]; + scalar& maxSz = maxUserSize[pointi]; + forAll(pFaces, pFacei) + { + label facei = pFaces[pFacei]; + if (isXFace[facei]) + { + const point& fc = faceCentres[facei]; + if (((cellCentres[own[facei]]-fc)&userDirection) > 0.0) + { + maxSz = max(maxSz, cellSize[own[facei]]); + } + if (mesh.isInternalFace(facei)) + { + if (((cellCentres[nei[facei]]-fc)&userDirection) > 0.0) + { + maxSz = max(maxSz, cellSize[nei[facei]]); + } + } + } + } + } + syncTools::syncPointList + ( + mesh, + maxUserSize, + maxEqOp(), + -GREAT + ); +} + + +Foam::label Foam::snappyRefineDriver::directionalSmooth +( + const refinementParameters& refineParams +) +{ + addProfiling(split, "snappyHexMesh::refine::smooth"); + Info<< nl + << "Directional expansion ratio smoothing" << nl + << "-------------------------------------" << nl + << endl; + + fvMesh& mesh = meshRefiner_.mesh(); + const pointField& points = mesh.points(); + const searchableSurfaces& geometry = meshRefiner_.surfaces().geometry(); + const shellSurfaces& shells = meshRefiner_.shells(); + + label iter = 0; + + forAll(shells.nSmoothExpansion(), shellI) + { + if (shells.nSmoothExpansion()[shellI] > 0) + { + label surfi = shells.shells()[shellI]; + const vector& userDirection = shells.smoothDirection()[shellI]; + + + // Extract inside points + labelList pointLabels; + PackedBoolList isInsidePoint(mesh.nPoints()); + { + // Get inside points + List volType; + geometry[surfi].getVolumeType(points, volType); + + forAll(volType, pointi) + { + if (volType[pointi] == volumeType::INSIDE) + { + isInsidePoint.set(pointi); + } + } + pointLabels = isInsidePoint.used(); + } + + + const edgeList& edges = mesh.edges(); + const labelListList& pointEdges = mesh.pointEdges(); + + // Mark all directed edges + PackedBoolList isXEdge(edges.size()); + { + forAll(edges, edgei) + { + const edge& e = edges[edgei]; + vector eVec(e.vec(points)); + eVec /= mag(eVec); + if (mag(eVec&userDirection) > 0.9) + { + isXEdge.set(edgei); + } + } + } + + + // Mark all directed faces (with all points inside the region) + PackedBoolList isXFace(mesh.nFaces()); + { + const vectorField& faceAreas = mesh.faceAreas(); + const faceList& faces = mesh.faces(); + + forAll(faces, facei) + { + const face& f = faces[facei]; + + bool allInRegion = true; + for (label pointi : f) + { + if (!isInsidePoint[pointi]) + { + allInRegion = false; + break; + } + } + + if (allInRegion) + { + vector n(faceAreas[facei]); + n /= mag(n); + if (mag(n&userDirection) > 0.9) + { + isXFace[facei] = true; + } + } + } + } + + + // Get left and right local edge + labelList leftEdge(points.size(), -1); + labelList rightEdge(points.size(), -1); + + forAll(pointLabels, i) + { + label pointi = pointLabels[i]; + const labelList& pEdges = pointEdges[pointi]; + forAll(pEdges, pEdgei) + { + label edgei = pEdges[pEdgei]; + const edge& e = edges[edgei]; + label otherPointi = e.otherVertex(pointi); + vector eVec(points[otherPointi]-points[pointi]); + eVec /= mag(eVec); + scalar cosAngle = (eVec&userDirection); + + if (cosAngle < -0.9) + { + leftEdge[pointi] = edgei; + } + else if (cosAngle > 0.9) + { + rightEdge[pointi] = edgei; + } + } + } + + + pointField newPoints(points); + + Info<< "Smoothing " << geometry[surfi].name() << ':' << nl + << " Direction : " << userDirection << nl + << " Number of points : " + << returnReduce(pointLabels.size(), sumOp