BUG: autoLayerDriver: update local data

This commit is contained in:
mattijs
2013-12-02 11:08:08 +00:00
parent 1f35db6234
commit a4d50f98ec
10 changed files with 388 additions and 437 deletions

View File

@ -55,6 +55,10 @@ Description
#include "fixedValueFvPatchFields.H" #include "fixedValueFvPatchFields.H"
#include "localPointRegion.H" #include "localPointRegion.H"
#include "externalDisplacementMeshMover.H"
#include "medialAxisMeshMover.H"
#include "scalarIOField.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam namespace Foam
@ -1378,8 +1382,8 @@ void Foam::autoLayerDriver::syncPatchDisplacement
mesh, mesh,
meshPoints, meshPoints,
patchDisp, patchDisp,
minEqOp<vector>(), minMagSqrEqOp<vector>(),
point::max // null value point::rootMax // null value
); );
// Unmark if displacement too small // Unmark if displacement too small
@ -1487,7 +1491,7 @@ void Foam::autoLayerDriver::syncPatchDisplacement
// patch point. // patch point.
void Foam::autoLayerDriver::getPatchDisplacement void Foam::autoLayerDriver::getPatchDisplacement
( (
const motionSmoother& meshMover, const indirectPrimitivePatch& pp,
const scalarField& thickness, const scalarField& thickness,
const scalarField& minThickness, const scalarField& minThickness,
pointField& patchDisp, pointField& patchDisp,
@ -1499,7 +1503,6 @@ void Foam::autoLayerDriver::getPatchDisplacement
<< " according to pointNormal ..." << endl; << " according to pointNormal ..." << endl;
const fvMesh& mesh = meshRefiner_.mesh(); const fvMesh& mesh = meshRefiner_.mesh();
const indirectPrimitivePatch& pp = meshMover.patch();
const vectorField& faceNormals = pp.faceNormals(); const vectorField& faceNormals = pp.faceNormals();
const labelListList& pointFaces = pp.pointFaces(); const labelListList& pointFaces = pp.pointFaces();
const pointField& localPoints = pp.localPoints(); const pointField& localPoints = pp.localPoints();
@ -1594,6 +1597,11 @@ void Foam::autoLayerDriver::getPatchDisplacement
nExtrudeRemove++; nExtrudeRemove++;
} }
else
{
// All surrounding points are not extruded. Leave patchDisp
// intact.
}
} }
} }
@ -1720,7 +1728,7 @@ Foam::label Foam::autoLayerDriver::truncateDisplacement
( (
const globalIndex& globalFaces, const globalIndex& globalFaces,
const labelListList& edgeGlobalFaces, const labelListList& edgeGlobalFaces,
const motionSmoother& meshMover, const indirectPrimitivePatch& pp,
const scalarField& minThickness, const scalarField& minThickness,
const faceSet& illegalPatchFaces, const faceSet& illegalPatchFaces,
pointField& patchDisp, pointField& patchDisp,
@ -1728,8 +1736,7 @@ Foam::label Foam::autoLayerDriver::truncateDisplacement
List<extrudeMode>& extrudeStatus List<extrudeMode>& extrudeStatus
) const ) const
{ {
const polyMesh& mesh = meshMover.mesh(); const fvMesh& mesh = meshRefiner_.mesh();
const indirectPrimitivePatch& pp = meshMover.patch();
label nChanged = 0; label nChanged = 0;
@ -2034,7 +2041,7 @@ Foam::label Foam::autoLayerDriver::truncateDisplacement
// regions where layer mesh terminates. // regions where layer mesh terminates.
void Foam::autoLayerDriver::setupLayerInfoTruncation void Foam::autoLayerDriver::setupLayerInfoTruncation
( (
const motionSmoother& meshMover, const indirectPrimitivePatch& pp,
const labelList& patchNLayers, const labelList& patchNLayers,
const List<extrudeMode>& extrudeStatus, const List<extrudeMode>& extrudeStatus,
const label nBufferCellsNoExtrude, const label nBufferCellsNoExtrude,
@ -2044,8 +2051,7 @@ void Foam::autoLayerDriver::setupLayerInfoTruncation
{ {
Info<< nl << "Setting up information for layer truncation ..." << endl; Info<< nl << "Setting up information for layer truncation ..." << endl;
const indirectPrimitivePatch& pp = meshMover.patch(); const fvMesh& mesh = meshRefiner_.mesh();
const polyMesh& mesh = meshMover.mesh();
if (nBufferCellsNoExtrude < 0) if (nBufferCellsNoExtrude < 0)
{ {
@ -2917,26 +2923,6 @@ void Foam::autoLayerDriver::addLayers
); );
// Construct iterative mesh mover.
Info<< "Constructing mesh displacer ..." << endl;
autoPtr<motionSmoother> meshMover
(
new motionSmoother
(
mesh,
pp(),
patchIDs,
makeLayerDisplacementField
(
pointMesh::New(mesh),
layerParams.numLayers()
),
motionDict
)
);
// Point-wise extrusion data // Point-wise extrusion data
// ~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~
@ -2961,7 +2947,7 @@ void Foam::autoLayerDriver::addLayers
setNumLayers setNumLayers
( (
layerParams.numLayers(), // per patch the num layers layerParams.numLayers(), // per patch the num layers
meshMover().adaptPatchIDs(),// patches that are being moved patchIDs, // patches that are being moved
pp, // indirectpatch for all faces moving pp, // indirectpatch for all faces moving
patchDisp, patchDisp,
@ -3059,12 +3045,22 @@ void Foam::autoLayerDriver::addLayers
// Determine (wanted) point-wise overall layer thickness and expansion // Determine (wanted) point-wise overall layer thickness and expansion
// ratio // ratio
scalarField thickness(pp().nPoints()); scalarField thickness(pp().nPoints());
scalarField minThickness(pp().nPoints()); scalarIOField minThickness
(
IOobject
(
"minThickness",
meshRefiner_.timeName(),
mesh,
IOobject::NO_READ
),
pp().nPoints()
);
scalarField expansionRatio(pp().nPoints()); scalarField expansionRatio(pp().nPoints());
calculateLayerThickness calculateLayerThickness
( (
pp, pp,
meshMover().adaptPatchIDs(), patchIDs,
layerParams, layerParams,
cellLevel, cellLevel,
patchNLayers, patchNLayers,
@ -3077,89 +3073,36 @@ void Foam::autoLayerDriver::addLayers
// Calculate wall to medial axis distance for smoothing displacement // Overall displacement field
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ pointVectorField displacement
pointScalarField pointMedialDist
( (
IOobject makeLayerDisplacementField
( (
"pointMedialDist", pointMesh::New(mesh),
meshRefiner_.timeName(), layerParams.numLayers()
mesh, )
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
meshMover().pMesh(),
dimensionedScalar("pointMedialDist", dimLength, 0.0)
); );
pointVectorField dispVec // Allocate run-time selectable mesh mover
( autoPtr<externalDisplacementMeshMover> medialAxisMoverPtr;
IOobject {
// Set up controls for meshMover
dictionary combinedDict(layerParams.dict());
// Add mesh quality constraints
combinedDict.merge(motionDict);
// Where to get minThickness from
combinedDict.add("minThicknessName", minThickness.name());
// Take over patchDisp as boundary conditions on displacement
// pointVectorField
medialAxisMoverPtr = externalDisplacementMeshMover::New
( (
"dispVec", layerParams.meshShrinker(),
meshRefiner_.timeName(), combinedDict,
mesh, baffles,
IOobject::NO_READ, displacement
IOobject::NO_WRITE, );
false }
),
meshMover().pMesh(),
dimensionedVector("dispVec", dimLength, vector::zero)
);
pointScalarField medialRatio
(
IOobject
(
"medialRatio",
meshRefiner_.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
meshMover().pMesh(),
dimensionedScalar("medialRatio", dimless, 0.0)
);
pointVectorField medialVec
(
IOobject
(
"medialVec",
meshRefiner_.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
meshMover().pMesh(),
dimensionedVector("medialVec", dimLength, vector::zero)
);
// Setup information for medial axis smoothing. Calculates medial axis
// and a smoothed displacement direction.
// - pointMedialDist : distance to medial axis
// - dispVec : normalised direction of nearest displacement
// - medialRatio : ratio of medial distance to wall distance.
// (1 at wall, 0 at medial axis)
medialAxisSmoothingInfo
(
meshMover,
layerParams.nSmoothNormals(),
layerParams.nSmoothSurfaceNormals(),
layerParams.minMedianAxisAngleCos(),
layerParams.slipFeatureAngle(),
dispVec,
medialRatio,
pointMedialDist,
medialVec
);
// Saved old points // Saved old points
@ -3215,7 +3158,7 @@ void Foam::autoLayerDriver::addLayers
// Displacement acc. to pointnormals // Displacement acc. to pointnormals
getPatchDisplacement getPatchDisplacement
( (
meshMover, pp,
thickness, thickness,
minThickness, minThickness,
patchDisp, patchDisp,
@ -3229,43 +3172,39 @@ void Foam::autoLayerDriver::addLayers
{ {
pointField oldPatchPos(pp().localPoints()); pointField oldPatchPos(pp().localPoints());
//// Laplacian displacement shrinking. // Take over patchDisp into pointDisplacement field and
//shrinkMeshDistance // adjust both for multi-patch constraints
//( motionSmootherAlgo::setDisplacement
// debug,
// meshMover,
// -patchDisp, // Shrink in opposite direction of addedPoints
// layerParams.nSmoothDisp(),
// layerParams.nSnap()
//);
// Medial axis based shrinking
shrinkMeshMedialDistance
( (
meshMover(), patchIDs,
meshQualityDict, pp,
baffles,
layerParams.nSmoothThickness(),
layerParams.nSmoothDisplacement(),
layerParams.maxThicknessToMedialRatio(),
nAllowableErrors,
layerParams.nSnap(),
layerParams.layerTerminationCos(),
thickness,
minThickness,
dispVec,
medialRatio,
pointMedialDist,
medialVec,
extrudeStatus,
patchDisp, patchDisp,
patchNLayers displacement
); );
// Move mesh
// ~~~~~~~~~
// Set up controls for meshMover
dictionary combinedDict(layerParams.dict());
// Add standard quality constraints
combinedDict.merge(motionDict);
// Add relaxed constraints (overrides standard ones)
combinedDict.merge(meshQualityDict);
// Where to get minThickness from
combinedDict.add("minThicknessName", minThickness.name());
labelList checkFaces(identity(mesh.nFaces()));
medialAxisMoverPtr().move
(
combinedDict,
nAllowableErrors,
checkFaces
);
pp().movePoints(mesh.points());
// Update patchDisp (since not all might have been honoured) // Update patchDisp (since not all might have been honoured)
patchDisp = oldPatchPos - pp().localPoints(); patchDisp = oldPatchPos - pp().localPoints();
} }
@ -3278,7 +3217,7 @@ void Foam::autoLayerDriver::addLayers
( (
globalFaces, globalFaces,
edgeGlobalFaces, edgeGlobalFaces,
meshMover(), pp,
minThickness, minThickness,
dummySet, dummySet,
patchDisp, patchDisp,
@ -3331,7 +3270,7 @@ void Foam::autoLayerDriver::addLayers
labelList nPatchFaceLayers(pp().size(), -1); labelList nPatchFaceLayers(pp().size(), -1);
setupLayerInfoTruncation setupLayerInfoTruncation
( (
meshMover(), pp,
patchNLayers, patchNLayers,
extrudeStatus, extrudeStatus,
layerParams.nBufferCellsNoExtrude(), layerParams.nBufferCellsNoExtrude(),
@ -3456,7 +3395,7 @@ void Foam::autoLayerDriver::addLayers
<< endl; << endl;
newMesh.write(); newMesh.write();
cellSet addedCellSet(newMesh," addedCells", nAddedCells); cellSet addedCellSet(newMesh, "addedCells", nAddedCells);
forAll(cellNLayers, cellI) forAll(cellNLayers, cellI)
{ {
if (cellNLayers[cellI] > 0) if (cellNLayers[cellI] > 0)
@ -3521,10 +3460,7 @@ void Foam::autoLayerDriver::addLayers
// Reset mesh points and start again // Reset mesh points and start again
mesh.movePoints(oldPoints); mesh.movePoints(oldPoints);
// Update meshmover for change of mesh geometry pp().movePoints(mesh.points());
meshMover().movePoints();
meshMover().correct();
// Grow out region of non-extrusion // Grow out region of non-extrusion
for (label i = 0; i < layerParams.nGrow(); i++) for (label i = 0; i < layerParams.nGrow(); i++)

View File

@ -277,7 +277,7 @@ private:
//- Get nearest point on surface to snap to //- Get nearest point on surface to snap to
void getPatchDisplacement void getPatchDisplacement
( (
const motionSmoother& meshMover, const indirectPrimitivePatch& pp,
const scalarField& thickness, const scalarField& thickness,
const scalarField& minThickness, const scalarField& minThickness,
pointField& patchDisp, pointField& patchDisp,
@ -315,7 +315,7 @@ private:
( (
const globalIndex& globalFaces, const globalIndex& globalFaces,
const labelListList& edgeGlobalFaces, const labelListList& edgeGlobalFaces,
const motionSmoother& meshMover, const indirectPrimitivePatch& pp,
const scalarField& minThickness, const scalarField& minThickness,
const faceSet& illegalPatchFaces, const faceSet& illegalPatchFaces,
pointField& patchDisp, pointField& patchDisp,
@ -331,7 +331,7 @@ private:
// to go into the actual layer addition engine. // to go into the actual layer addition engine.
void setupLayerInfoTruncation void setupLayerInfoTruncation
( (
const motionSmoother& meshMover, const indirectPrimitivePatch& pp,
const labelList& patchNLayers, const labelList& patchNLayers,
const List<extrudeMode>& extrudeStatus, const List<extrudeMode>& extrudeStatus,
const label nBufferCellsNoExtrude, const label nBufferCellsNoExtrude,

View File

@ -752,15 +752,15 @@ void Foam::autoLayerDriver::findIsolatedRegions
nPointCounter++; nPointCounter++;
nChanged++; nChanged++;
} }
} }
} }
if (returnReduce(nChanged, sumOp<label>()) == 0) if (returnReduce(nChanged, sumOp<label>()) == 0)
{ {
break; break;
} }
} }
const edgeList& edges = pp.edges(); const edgeList& edges = pp.edges();

View File

@ -28,6 +28,7 @@ License
#include "unitConversion.H" #include "unitConversion.H"
#include "refinementSurfaces.H" #include "refinementSurfaces.H"
#include "searchableSurfaces.H" #include "searchableSurfaces.H"
#include "medialAxisMeshMover.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -109,6 +110,7 @@ Foam::layerParameters::layerParameters
const polyBoundaryMesh& boundaryMesh const polyBoundaryMesh& boundaryMesh
) )
: :
dict_(dict),
numLayers_(boundaryMesh.size(), -1), numLayers_(boundaryMesh.size(), -1),
relativeSizes_(dict.lookup("relativeSizes")), relativeSizes_(dict.lookup("relativeSizes")),
layerSpec_(ILLEGAL), layerSpec_(ILLEGAL),
@ -163,7 +165,15 @@ Foam::layerParameters::layerParameters
nSnap_(readLabel(dict.lookup("nRelaxIter"))), nSnap_(readLabel(dict.lookup("nRelaxIter"))),
nLayerIter_(readLabel(dict.lookup("nLayerIter"))), nLayerIter_(readLabel(dict.lookup("nLayerIter"))),
nRelaxedIter_(labelMax), nRelaxedIter_(labelMax),
additionalReporting_(dict.lookupOrDefault("additionalReporting", false)) additionalReporting_(dict.lookupOrDefault("additionalReporting", false)),
meshShrinker_
(
dict.lookupOrDefault
(
"meshShrinker",
medialAxisMeshMover::typeName
)
)
{ {
// Detect layer specification mode // Detect layer specification mode

View File

@ -86,6 +86,9 @@ private:
// Private data // Private data
const dictionary dict_;
// Per patch (not region!) information // Per patch (not region!) information
//- How many layers to add. //- How many layers to add.
@ -143,6 +146,8 @@ private:
Switch additionalReporting_; Switch additionalReporting_;
word meshShrinker_;
// Private Member Functions // Private Member Functions
//- Calculate expansion ratio from overall size v.s. thickness of //- Calculate expansion ratio from overall size v.s. thickness of
@ -170,6 +175,12 @@ public:
// Member Functions // Member Functions
const dictionary& dict() const
{
return dict_;
}
// Per patch information // Per patch information
//- How many layers to add. //- How many layers to add.
@ -230,98 +241,7 @@ public:
} }
scalar featureAngle() const // Control
{
return featureAngle_;
}
//- At non-patched sides allow mesh to slip if extrusion
// direction makes angle larger than slipFeatureAngle.
scalar slipFeatureAngle() const
{
return slipFeatureAngle_;
}
scalar concaveAngle() const
{
return concaveAngle_;
}
//- 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.
label nGrow() const
{
return nGrow_;
}
//- Number of smoothing iterations of surface normals
label nSmoothSurfaceNormals() const
{
return nSmoothSurfaceNormals_;
}
//- Number of smoothing iterations of interior mesh movement
// direction
label nSmoothNormals() const
{
return nSmoothNormals_;
}
//- Stop layer growth on highly warped cells
scalar maxFaceThicknessRatio() const
{
return maxFaceThicknessRatio_;
}
scalar layerTerminationCos() const
{
return layerTerminationCos_;
}
//- Smooth internal displacement
label nSmoothDisplacement() const
{
return nSmoothDisplacement_;
}
//- Smooth layer thickness over surface patches
label nSmoothThickness() const
{
return nSmoothThickness_;
}
//- Reduce layer growth where ratio thickness to medial
// distance is large
scalar maxThicknessToMedialRatio() const
{
return maxThicknessToMedialRatio_;
}
//- Angle used to pick up medial axis points
scalar minMedianAxisAngleCos() const
{
return minMedianAxisAngleCos_;
}
//- Create buffer region for new layer terminations
label nBufferCellsNoExtrude() const
{
return nBufferCellsNoExtrude_;
}
label nSnap() const
{
return nSnap_;
}
const Switch& additionalReporting() const
{
return additionalReporting_;
}
// Overall
//- Number of overall layer addition iterations //- Number of overall layer addition iterations
label nLayerIter() const label nLayerIter() const
@ -337,6 +257,109 @@ public:
} }
// Control
scalar featureAngle() const
{
return featureAngle_;
}
scalar concaveAngle() const
{
return concaveAngle_;
}
//- 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.
label nGrow() const
{
return nGrow_;
}
//- Stop layer growth on highly warped cells
scalar maxFaceThicknessRatio() const
{
return maxFaceThicknessRatio_;
}
//- Create buffer region for new layer terminations
label nBufferCellsNoExtrude() const
{
return nBufferCellsNoExtrude_;
}
const Switch& additionalReporting() const
{
return additionalReporting_;
}
//- Type of mesh shrinker
const word& meshShrinker() const
{
return meshShrinker_;
}
// Specific to medial axis mesh shrinking
//- At non-patched sides allow mesh to slip if extrusion
// direction makes angle larger than slipFeatureAngle.
scalar slipFeatureAngle() const
{
return slipFeatureAngle_;
}
//- Number of smoothing iterations of surface normals
label nSmoothSurfaceNormals() const
{
return nSmoothSurfaceNormals_;
}
//- Number of smoothing iterations of interior mesh movement
// direction
label nSmoothNormals() const
{
return nSmoothNormals_;
}
scalar layerTerminationCos() const
{
return layerTerminationCos_;
}
//- Smooth internal displacement
label nSmoothDisplacement() const
{
return nSmoothDisplacement_;
}
//- Smooth layer thickness over surface patches
label nSmoothThickness() const
{
return nSmoothThickness_;
}
//- Reduce layer growth where ratio thickness to medial
// distance is large
scalar maxThicknessToMedialRatio() const
{
return maxThicknessToMedialRatio_;
}
//- Angle used to pick up medial axis points
scalar minMedianAxisAngleCos() const
{
return minMedianAxisAngleCos_;
}
label nSnap() const
{
return nSnap_;
}
// Helper // Helper
//- Determine overall thickness. Uses two of the four parameters //- Determine overall thickness. Uses two of the four parameters

View File

@ -105,7 +105,8 @@ void Foam::displacementMeshMoverMotionSolver::solve()
pointDisplacement().boundaryField().updateCoeffs(); pointDisplacement().boundaryField().updateCoeffs();
label nAllowableErrors = 0; label nAllowableErrors = 0;
meshMover().move(nAllowableErrors); labelList checkFaces(identity(mesh().nFaces()));
meshMover().move(coeffDict(), nAllowableErrors, checkFaces);
} }

View File

@ -44,7 +44,6 @@ Foam::externalDisplacementMeshMover::externalDisplacementMeshMover
pointVectorField& pointDisplacement pointVectorField& pointDisplacement
) )
: :
dict_(dict),
baffles_(baffles), baffles_(baffles),
pointDisplacement_(pointDisplacement) pointDisplacement_(pointDisplacement)
{} {}

View File

@ -28,7 +28,10 @@ Description
Virtual base class for mesh movers with externally provided displacement Virtual base class for mesh movers with externally provided displacement
field giving the boundary conditions. Move the mesh from the current field giving the boundary conditions. Move the mesh from the current
location to a new location (so modify the mesh; v.s. motionSolver that location to a new location (so modify the mesh; v.s. motionSolver that
only returns the new location) only returns the new location).
All mesh movers are expected to read the dictionary settings at invocation
of move(), i.e. not cache any settings.
SourceFiles SourceFiles
externalDisplacementMeshMover.C externalDisplacementMeshMover.C
@ -57,9 +60,6 @@ protected:
// Protected data // Protected data
//- Settings
dictionary dict_;
//- Baffles in the mesh //- Baffles in the mesh
List<labelPair> baffles_; List<labelPair> baffles_;
@ -104,7 +104,8 @@ public:
// Constructors // Constructors
//- Construct from dictionary and displacement field //- Construct from dictionary and displacement field. Dictionary is
// allowed to go out of scope!
externalDisplacementMeshMover externalDisplacementMeshMover
( (
const dictionary& dict, const dictionary& dict,
@ -158,8 +159,16 @@ public:
// Mesh mover // Mesh mover
//- Move mesh. Return true if succesful //- Move mesh using current pointDisplacement boundary values
virtual bool move(const label nAllowableErrors) = 0; // and current dictionary settings. Return true if succesful
// (errors on checkFaces less than allowable). Updates
// pointDisplacement.
virtual bool move
(
const dictionary&,
const label nAllowableErrors,
labelList& checkFaces
) = 0;
//- Update local data for geometry changes //- Update local data for geometry changes
virtual void movePoints(const pointField&); virtual void movePoints(const pointField&);

View File

@ -52,51 +52,6 @@ namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
const Foam::dictionary& Foam::medialAxisMeshMover::coeffDict() const
{
return dict_;
}
const Foam::dictionary& Foam::medialAxisMeshMover::meshQualityDict() const
{
return dict_;
}
void Foam::medialAxisMeshMover::read()
{
nSmoothPatchThickness_ = readLabel(coeffDict().lookup("nSmoothThickness"));
nSmoothDisplacement_ = coeffDict().lookupOrDefault
(
"nSmoothDisplacement",
0
);
maxThicknessToMedialRatio_ = readScalar
(
coeffDict().lookup("maxThicknessToMedialRatio")
);
featureAngle_ = readScalar(coeffDict().lookup("featureAngle"));
minCosLayerTermination_ = Foam::cos(degToRad(0.5*featureAngle_));
nSmoothNormals_ = readLabel(coeffDict().lookup("nSmoothNormals"));
nSmoothSurfaceNormals_ = readLabel
(
coeffDict().lookup("nSmoothSurfaceNormals")
);
minMedianAxisAngleCos_ = Foam::cos
(
degToRad(readScalar(coeffDict().lookup("minMedianAxisAngle")))
);
slipFeatureAngle_ =
(
coeffDict().found("slipFeatureAngle")
? readScalar(coeffDict().lookup("slipFeatureAngle"))
: 0.5*featureAngle_
);
nSnap_ = readLabel(coeffDict().lookup("nRelaxIter"));
}
Foam::labelList Foam::medialAxisMeshMover::getFixedValueBCs Foam::labelList Foam::medialAxisMeshMover::getFixedValueBCs
( (
const pointVectorField& fld const pointVectorField& fld
@ -444,7 +399,7 @@ bool Foam::medialAxisMeshMover::isMaxEdge
} }
void Foam::medialAxisMeshMover::update() void Foam::medialAxisMeshMover::update(const dictionary& coeffDict)
{ {
Info<< typeName Info<< typeName
<< " : Calculate distance to Medial Axis ..." << endl; << " : Calculate distance to Medial Axis ..." << endl;
@ -455,29 +410,37 @@ void Foam::medialAxisMeshMover::update()
const labelList& meshPoints = pp.meshPoints(); const labelList& meshPoints = pp.meshPoints();
// Look up a few controls // Read a few parameters
// ~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~
// const label nSmoothNormals
// ( //- Smooth surface normals
// readLabel(coeffDict().lookup("nSmoothNormals")) const label nSmoothSurfaceNormals = readLabel
// ); (
// const label nSmoothSurfaceNormals coeffDict.lookup("nSmoothSurfaceNormals")
// ( );
// readLabel(coeffDict().lookup("nSmoothSurfaceNormals"))
// ); //- When is medial axis
// const scalar minMedianAxisAngleCos const scalar minMedianAxisAngleCos = Foam::cos
// ( (
// Foam::cos degToRad(readScalar(coeffDict.lookup("minMedianAxisAngle")))
// ( );
// degToRad(readScalar(coeffDict().lookup("minMedianAxisAngle")))
// ) //- Feature angle when to stop adding layers
// ); const scalar featureAngle = readScalar(coeffDict.lookup("featureAngle"));
// const scalar slipFeatureAngle
// ( //- When to slip along wall
// coeffDict().found("slipFeatureAngle") const scalar slipFeatureAngle =
// ? readScalar(coeffDict().lookup("slipFeatureAngle")) (
// : 0.5*readScalar(coeffDict().lookup("featureAngle")) coeffDict.found("slipFeatureAngle")
// ); ? readScalar(coeffDict.lookup("slipFeatureAngle"))
: 0.5*featureAngle
);
//- Smooth internal normals
const label nSmoothNormals = readLabel
(
coeffDict.lookup("nSmoothNormals")
);
// Predetermine mesh edges // Predetermine mesh edges
@ -505,7 +468,7 @@ void Foam::medialAxisMeshMover::update()
// Smooth patch normal vectors // Smooth patch normal vectors
smoothPatchNormals smoothPatchNormals
( (
nSmoothSurfaceNormals_, nSmoothSurfaceNormals,
isMasterPoint, isMasterPoint,
isMasterEdge, isMasterEdge,
meshEdges, meshEdges,
@ -602,7 +565,7 @@ void Foam::medialAxisMeshMover::update()
{ {
// Unvisited point. See above about nUnvisit warning // Unvisited point. See above about nUnvisit warning
} }
else if (isMaxEdge(pointWallDist, edgeI, minMedianAxisAngleCos_)) else if (isMaxEdge(pointWallDist, edgeI, minMedianAxisAngleCos))
{ {
// Both end points of edge have very different nearest wall // Both end points of edge have very different nearest wall
// point. Mark both points as medial axis points. // point. Mark both points as medial axis points.
@ -722,11 +685,11 @@ void Foam::medialAxisMeshMover::update()
// on this patch. // on this patch.
Info<< "Inserting points on patch " << pp.name() Info<< "Inserting points on patch " << pp.name()
<< " if angle to nearest layer patch > " << " if angle to nearest layer patch > "
<< slipFeatureAngle_ << " degrees." << endl; << slipFeatureAngle << " degrees." << endl;
scalar slipFeatureAngleCos = Foam::cos scalar slipFeatureAngleCos = Foam::cos
( (
degToRad(slipFeatureAngle_) degToRad(slipFeatureAngle)
); );
pointField pointNormals pointField pointNormals
( (
@ -818,7 +781,7 @@ void Foam::medialAxisMeshMover::update()
// Smooth normal vectors. Do not change normals on pp.meshPoints // Smooth normal vectors. Do not change normals on pp.meshPoints
smoothNormals smoothNormals
( (
nSmoothNormals_, nSmoothNormals,
isMasterPoint, isMasterPoint,
isMasterEdge, isMasterEdge,
meshPoints, meshPoints,
@ -854,7 +817,7 @@ void Foam::medialAxisMeshMover::update()
} }
//if (debug) if (debug)
{ {
Info<< typeName Info<< typeName
<< " : Writing medial axis fields:" << nl << " : Writing medial axis fields:" << nl
@ -926,8 +889,8 @@ void Foam::medialAxisMeshMover::syncPatchDisplacement
mesh(), mesh(),
meshPoints, meshPoints,
patchDisp, patchDisp,
minEqOp<vector>(), minMagSqrEqOp<vector>(),
point::max // null value point::rootMax // null value
); );
// Unmark if displacement too small // Unmark if displacement too small
@ -1323,15 +1286,14 @@ void Foam::medialAxisMeshMover::findIsolatedRegions
nPointCounter++; nPointCounter++;
nChanged++; nChanged++;
} }
} }
} }
if (returnReduce(nChanged, sumOp<label>()) == 0)
if (returnReduce(nChanged, sumOp<label>()) == 0) {
{ break;
break; }
} }
}
const edgeList& edges = pp.edges(); const edgeList& edges = pp.edges();
@ -1554,7 +1516,7 @@ Foam::medialAxisMeshMover::medialAxisMeshMover
scale_, scale_,
oldPoints_, oldPoints_,
adaptPatchIDs_, adaptPatchIDs_,
coeffDict() dict
), ),
dispVec_ dispVec_
( (
@ -1613,8 +1575,7 @@ Foam::medialAxisMeshMover::medialAxisMeshMover
dimensionedVector("medialVec", dimLength, vector::zero) dimensionedVector("medialVec", dimLength, vector::zero)
) )
{ {
read(); update(dict);
update();
} }
@ -1628,6 +1589,7 @@ Foam::medialAxisMeshMover::~medialAxisMeshMover()
void Foam::medialAxisMeshMover::calculateDisplacement void Foam::medialAxisMeshMover::calculateDisplacement
( (
const dictionary& coeffDict,
const scalarField& minThickness, const scalarField& minThickness,
List<autoLayerDriver::extrudeMode>& extrudeStatus, List<autoLayerDriver::extrudeMode>& extrudeStatus,
pointField& patchDisp pointField& patchDisp
@ -1639,29 +1601,36 @@ void Foam::medialAxisMeshMover::calculateDisplacement
const labelList& meshPoints = pp.meshPoints(); const labelList& meshPoints = pp.meshPoints();
// Look up a few controls // Read settings
// ~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~
// const label nSmoothPatchThickness //- (lambda-mu) smoothing of internal displacement
// ( const label nSmoothDisplacement = coeffDict.lookupOrDefault
// readLabel(coeffDict().lookup("nSmoothThickness")) (
// ); "nSmoothDisplacement",
// const label nSmoothDisplacement 0
// ( );
// coeffDict().lookupOrDefault("nSmoothDisplacement", 0)
// ); //- Layer thickness too big
// const scalar maxThicknessToMedialRatio const scalar maxThicknessToMedialRatio = readScalar
// ( (
// readScalar(coeffDict().lookup("maxThicknessToMedialRatio")) coeffDict.lookup("maxThicknessToMedialRatio")
// ); );
// const scalar featureAngle
// ( //- Feature angle when to stop adding layers
// readScalar(coeffDict().lookup("featureAngle")) const scalar featureAngle = readScalar(coeffDict.lookup("featureAngle"));
// );
// const scalar minCosLayerTermination //- Stop layer growth where mesh wraps around sharp edge
// ( const scalar minCosLayerTermination = Foam::cos
// Foam::cos(degToRad(0.5*featureAngle)) (
// ); degToRad(0.5*featureAngle)
);
//- Smoothing wanted patch thickness
const label nSmoothPatchThickness = readLabel
(
coeffDict.lookup("nSmoothThickness")
);
// Precalulate master points/edge (only relevant for shared points/edges) // Precalulate master points/edge (only relevant for shared points/edges)
@ -1749,7 +1718,7 @@ void Foam::medialAxisMeshMover::calculateDisplacement
mVec /= mag(mVec)+VSMALL; mVec /= mag(mVec)+VSMALL;
thicknessRatio *= (n&mVec); thicknessRatio *= (n&mVec);
if (thicknessRatio > maxThicknessToMedialRatio_) if (thicknessRatio > maxThicknessToMedialRatio)
{ {
// Truncate thickness. // Truncate thickness.
if (debug) if (debug)
@ -1810,7 +1779,7 @@ void Foam::medialAxisMeshMover::calculateDisplacement
findIsolatedRegions findIsolatedRegions
( (
minCosLayerTermination_, minCosLayerTermination,
isMasterPoint, isMasterPoint,
isMasterEdge, isMasterEdge,
meshEdges, meshEdges,
@ -1833,7 +1802,7 @@ void Foam::medialAxisMeshMover::calculateDisplacement
// smooth layer thickness on moving patch // smooth layer thickness on moving patch
minSmoothField minSmoothField
( (
nSmoothPatchThickness_, nSmoothPatchThickness,
isMasterPoint, isMasterPoint,
isMasterEdge, isMasterEdge,
meshEdges, meshEdges,
@ -1911,11 +1880,11 @@ void Foam::medialAxisMeshMover::calculateDisplacement
// Smear displacement away from fixed values (medialRatio=0 or 1) // Smear displacement away from fixed values (medialRatio=0 or 1)
if (nSmoothDisplacement_ > 0) if (nSmoothDisplacement > 0)
{ {
smoothLambdaMuDisplacement smoothLambdaMuDisplacement
( (
nSmoothDisplacement_, nSmoothDisplacement,
isMasterPoint, isMasterPoint,
isMasterEdge, isMasterEdge,
displacement displacement
@ -1927,25 +1896,29 @@ void Foam::medialAxisMeshMover::calculateDisplacement
bool Foam::medialAxisMeshMover::shrinkMesh bool Foam::medialAxisMeshMover::shrinkMesh
( (
const dictionary& meshQualityDict, const dictionary& meshQualityDict,
const label nAllowableErrors const label nAllowableErrors,
labelList& checkFaces
) )
{ {
//- Number of attempts shrinking the mesh
const label nSnap = readLabel(meshQualityDict.lookup("nRelaxIter"));
// Make sure displacement boundary conditions is uptodate with // Make sure displacement boundary conditions is uptodate with
// internal field // internal field
meshMover_.setDisplacementPatchFields(); meshMover_.setDisplacementPatchFields();
// Current faces to check. Gets modified in meshMover.scaleMesh
labelList checkFaces(identity(mesh().nFaces()));
Info<< typeName << " : Moving mesh ..." << endl; Info<< typeName << " : Moving mesh ..." << endl;
scalar oldErrorReduction = -1; scalar oldErrorReduction = -1;
bool meshOk = false; bool meshOk = false;
for (label iter = 0; iter < 2*nSnap_ ; iter++) for (label iter = 0; iter < 2*nSnap ; iter++)
{ {
Info<< "Iteration " << iter << endl; Info<< "Iteration " << iter << endl;
if (iter == nSnap_) if (iter == nSnap)
{ {
Info<< "Displacement scaling for error reduction set to 0." << endl; Info<< "Displacement scaling for error reduction set to 0." << endl;
oldErrorReduction = meshMover_.setErrorReduction(0.0); oldErrorReduction = meshMover_.setErrorReduction(0.0);
@ -1981,8 +1954,20 @@ bool Foam::medialAxisMeshMover::shrinkMesh
} }
bool Foam::medialAxisMeshMover::move(const label nAllowableErrors) bool Foam::medialAxisMeshMover::move
(
const dictionary& moveDict,
const label nAllowableErrors,
labelList& checkFaces
)
{ {
// Read a few settings
// ~~~~~~~~~~~~~~~~~~~
//- Name of field specifying min thickness
const word minThicknessName = word(moveDict.lookup("minThicknessName"));
// The points have moved so before calculation update // The points have moved so before calculation update
// the mesh and motionSolver accordingly // the mesh and motionSolver accordingly
movePoints(mesh().points()); movePoints(mesh().points());
@ -1994,6 +1979,19 @@ bool Foam::medialAxisMeshMover::move(const label nAllowableErrors)
// Extract out patch-wise displacement // Extract out patch-wise displacement
const indirectPrimitivePatch& pp = adaptPatchPtr_(); const indirectPrimitivePatch& pp = adaptPatchPtr_();
scalarField zeroMinThickness;
if (minThicknessName == "none")
{
zeroMinThickness = scalarField(pp.nPoints(), 0.0);
}
const scalarField& minThickness =
(
(minThicknessName == "none")
? zeroMinThickness
: mesh().lookupObject<scalarField>(minThicknessName)
);
pointField patchDisp pointField patchDisp
( (
pointDisplacement_.internalField(), pointDisplacement_.internalField(),
@ -2005,18 +2003,24 @@ bool Foam::medialAxisMeshMover::move(const label nAllowableErrors)
pp.nPoints(), pp.nPoints(),
autoLayerDriver::EXTRUDE autoLayerDriver::EXTRUDE
); );
forAll(extrudeStatus, pointI)
{
if (mag(patchDisp[pointI]) <= minThickness[pointI]+SMALL)
{
extrudeStatus[pointI] = autoLayerDriver::NOEXTRUDE;
}
}
const scalarField minThickness(patchDisp.size(), 0.0);
// Solve displacement // Solve displacement
calculateDisplacement(minThickness, extrudeStatus, patchDisp); calculateDisplacement(moveDict, minThickness, extrudeStatus, patchDisp);
//- Move mesh according to calculated displacement //- Move mesh according to calculated displacement
return shrinkMesh return shrinkMesh
( (
meshQualityDict(), // meshQualityDict, moveDict, // meshQualityDict,
nAllowableErrors // nAllowableErrors nAllowableErrors, // nAllowableErrors
checkFaces
); );
} }

View File

@ -73,39 +73,6 @@ class medialAxisMeshMover
motionSmootherAlgo meshMover_; motionSmootherAlgo meshMover_;
// Settings
//- Smoothing wanted patch thickness
label nSmoothPatchThickness_;
//- (lambda-mu) smoothing of internal displacement
label nSmoothDisplacement_;
//- Layer thickness too big
scalar maxThicknessToMedialRatio_;
//- Feature angle when to stop adding layers
scalar featureAngle_;
//- Stop layer growth where mesh wraps around sharp edge
scalar minCosLayerTermination_;
//- Smooth internal normals
label nSmoothNormals_;
//- Smooth surface normals
label nSmoothSurfaceNormals_;
//- When is medial axis
scalar minMedianAxisAngleCos_;
//- When to slip along wall
scalar slipFeatureAngle_;
//- Number of attempts shrinking the mesh
label nSnap_;
// Pre-calculated medial axis information // Pre-calculated medial axis information
//- normal of nearest wall. Where this normal changes direction //- normal of nearest wall. Where this normal changes direction
@ -125,13 +92,6 @@ class medialAxisMeshMover
// Private Member Functions // Private Member Functions
const dictionary& coeffDict() const;
const dictionary& meshQualityDict() const;
//- Read model coefficients
void read();
//- Extract fixed-value patchfields //- Extract fixed-value patchfields
static labelList getFixedValueBCs(const pointVectorField&); static labelList getFixedValueBCs(const pointVectorField&);
@ -201,8 +161,9 @@ class medialAxisMeshMover
const scalar minCos const scalar minCos
) const; ) const;
//- Initialise medial axis //- Initialise medial axis. Uses pointDisplacement field only
void update(); // for boundary types, not values.
void update(const dictionary&);
// Calculation of mesh movement // Calculation of mesh movement
@ -272,6 +233,7 @@ class medialAxisMeshMover
// (in pointDisplacement) // (in pointDisplacement)
void calculateDisplacement void calculateDisplacement
( (
const dictionary&,
const scalarField& minThickness, const scalarField& minThickness,
List<autoLayerDriver::extrudeMode>& extrudeStatus, List<autoLayerDriver::extrudeMode>& extrudeStatus,
pointField& patchDisp pointField& patchDisp
@ -281,10 +243,10 @@ class medialAxisMeshMover
bool shrinkMesh bool shrinkMesh
( (
const dictionary& meshQualityDict, const dictionary& meshQualityDict,
const label nAllowableErrors const label nAllowableErrors,
labelList& checkFaces
); );
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct
medialAxisMeshMover medialAxisMeshMover
( (
@ -319,8 +281,15 @@ public:
// Member Functions // Member Functions
//- Move mesh. Return true if succesful //- Move mesh using current pointDisplacement boundary values.
virtual bool move(const label nAllowableErrors); // Return true if succesful (errors on checkFaces less than
// allowable). Updates pointDisplacement.
virtual bool move
(
const dictionary&,
const label nAllowableErrors,
labelList& checkFaces
);
//- Update local data for geometry changes //- Update local data for geometry changes
virtual void movePoints(const pointField&); virtual void movePoints(const pointField&);