Skip to content

GitLab

Commit e07309ee authored by mattijs's avatar mattijs
Browse files

BUG: layering: disable mesh movement if extrusion disabled

parent 96a5814b
Hide whitespace changes
Inline Side-by-side
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
View file @ e07309ee
......@@ -48,6 +48,10 @@ Description
#include "globalIndex.H"
#include "DynamicField.H"
#include "PatchTools.H"
#include "slipPointPatchFields.H"
#include "fixedValuePointPatchFields.H"
#include "calculatedPointPatchFields.H"
#include "cyclicSlipPointPatchFields.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -71,7 +75,7 @@ void Foam::autoLayerDriver::dumpDisplacement
)
{
OFstream dispStr(prefix + "_disp.obj");
Info<< "Writing all displacements to " << dispStr.name() << nl << endl;
Info<< "Writing all displacements to " << dispStr.name() << endl;
label vertI = 0;
......@@ -87,7 +91,7 @@ void Foam::autoLayerDriver::dumpDisplacement
OFstream illStr(prefix + "_illegal.obj");
Info<< "Writing invalid displacements to " << illStr.name() << nl << endl;
Info<< "Writing invalid displacements to " << illStr.name() << endl;
vertI = 0;
......@@ -801,6 +805,81 @@ void Foam::autoLayerDriver::setNumLayers
}
// Construct pointVectorField with correct boundary conditions for adding
// layers
Foam::tmp<Foam::pointVectorField>
Foam::autoLayerDriver::makeLayerDisplacementField
(
const pointMesh& pMesh,
const labelList& numLayers
)
{
// Construct displacement field.
const pointBoundaryMesh& pointPatches = pMesh.boundary();
wordList patchFieldTypes
(
pointPatches.size(),
slipPointPatchVectorField::typeName
);
forAll(numLayers, patchI)
{
// 0 layers: do not allow lslip so fixedValue 0
// >0 layers: fixedValue which gets adapted
if (numLayers[patchI] >= 0)
{
patchFieldTypes[patchI] = fixedValuePointPatchVectorField::typeName;
}
}
forAll(pointPatches, patchI)
{
if (isA<processorPointPatch>(pointPatches[patchI]))
{
patchFieldTypes[patchI] = calculatedPointPatchVectorField::typeName;
}
else if (isA<cyclicPointPatch>(pointPatches[patchI]))
{
patchFieldTypes[patchI] = cyclicSlipPointPatchVectorField::typeName;
}
}
//Pout<< "*** makeLayerDisplacementField : boundary conditions:" << endl;
//forAll(patchFieldTypes, patchI)
//{
// Pout<< "\t" << patchI << " name:" << pointPatches[patchI].name()
// << " type:" << patchFieldTypes[patchI]
// << " nLayers:" << numLayers[patchI]
// << endl;
//}
const polyMesh& mesh = pMesh();
// Note: time().timeName() instead of meshRefinement::timeName() since
// postprocessable field.
tmp<pointVectorField> tfld
(
new pointVectorField
(
IOobject
(
"pointDisplacement",
mesh.time().timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
pMesh,
dimensionedVector("displacement", dimLength, vector::zero),
patchFieldTypes
)
);
return tfld;
}
void Foam::autoLayerDriver::growNoExtrusion
(
const indirectPrimitivePatch& pp,
......@@ -2392,10 +2471,10 @@ void Foam::autoLayerDriver::addLayers
mesh,
pp(),
patchIDs,
meshRefinement::makeDisplacementField
makeLayerDisplacementField
(
pointMesh::New(mesh),
patchIDs
layerParams.numLayers()
),
motionDict
)
......@@ -3186,7 +3265,7 @@ void Foam::autoLayerDriver::doLayers
// Merge coplanar boundary faces
mergePatchFacesUndo(layerParams, motionDict);
// Per patch the number of layers (0 if no layer)
// Per patch the number of layers (-1 or 0 if no layer)
const labelList& numLayers = layerParams.numLayers();
// Patches that need to get a layer
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H
View file @ e07309ee
......@@ -195,6 +195,19 @@ class autoLayerDriver
List<extrudeMode>& extrudeStatus
) const;
//- Helper function to make a pointVectorField with correct
// bcs for layer addition:
// - numLayers > 0 : fixedValue
// - numLayers == 0 : fixedValue (always zero)
// - processor : calculated (so free to move)
// - cyclic/wedge/symmetry : slip
// - other : slip
static tmp<pointVectorField> makeLayerDisplacementField
(
const pointMesh& pMesh,
const labelList& numLayers
);
//- Grow no-extrusion layer.
void growNoExtrusion
(
......@@ -444,7 +457,6 @@ class autoLayerDriver
const PackedBoolList& isMasterEdge,
const labelList& meshEdges,
const scalar minCosLayerTermination,
scalarField& field,
List<extrudeMode>& extrudeStatus,
pointField& patchDisp,
labelList& patchNLayers
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C
View file @ e07309ee
......@@ -99,6 +99,76 @@ void Foam::autoLayerDriver::sumWeights
// Smooth field on moving patch
//void Foam::autoLayerDriver::smoothField
//(
// const motionSmoother& meshMover,
// const PackedBoolList& isMasterEdge,
// const labelList& meshEdges,
// const scalarField& fieldMin,
// const label nSmoothDisp,
// scalarField& field
//) const
//{
// const indirectPrimitivePatch& pp = meshMover.patch();
// const edgeList& edges = pp.edges();
// const labelList& meshPoints = pp.meshPoints();
//
// scalarField invSumWeight(pp.nPoints());
// sumWeights
// (
// isMasterEdge,
// meshEdges,
// meshPoints,
// edges,
// invSumWeight
// );
//
// // Get smoothly varying patch field.
// Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
//
// for (label iter = 0; iter < nSmoothDisp; iter++)
// {
// scalarField average(pp.nPoints());
// averageNeighbours
// (
// meshMover.mesh(),
// isMasterEdge,
// meshEdges,
// meshPoints,
// pp.edges(),
// invSumWeight,
// field,
// average
// );
//
// // Transfer to field
// forAll(field, pointI)
// {
// //full smoothing neighbours + point value
// average[pointI] = 0.5*(field[pointI]+average[pointI]);
//
// // perform monotonic smoothing
// if
// (
// average[pointI] < field[pointI]
// && average[pointI] >= fieldMin[pointI]
// )
// {
// field[pointI] = average[pointI];
// }
// }
//
// // Do residual calculation every so often.
// if ((iter % 10) == 0)
// {
// Info<< " Iteration " << iter << " residual "
// << gSum(mag(field-average))
// /returnReduce(average.size(), sumOp<label>())
// << endl;
// }
// }
//}
//XXXXXXXXX
void Foam::autoLayerDriver::smoothField
(
const motionSmoother& meshMover,
......@@ -126,9 +196,15 @@ void Foam::autoLayerDriver::smoothField
// Get smoothly varying patch field.
Info<< "shrinkMeshDistance : Smoothing field ..." << endl;
for (label iter = 0; iter < nSmoothDisp; iter++)
const scalar lambda = 0.33;
const scalar mu = -0.34;
for (label iter = 0; iter < 90; iter++)
{
scalarField average(pp.nPoints());
// Calculate average of field
averageNeighbours
(
meshMover.mesh(),
......@@ -141,23 +217,37 @@ void Foam::autoLayerDriver::smoothField
average
);
// Transfer to field
forAll(field, pointI)
forAll(field, i)
{
//full smoothing neighbours + point value
average[pointI] = 0.5*(field[pointI]+average[pointI]);
if (field[i] >= fieldMin[i])
{
field[i] = (1-lambda)*field[i]+lambda*average[i];
}
}
// perform monotonic smoothing
if
(
average[pointI] < field[pointI]
&& average[pointI] >= fieldMin[pointI]
)
// Calculate average of field
averageNeighbours
(
meshMover.mesh(),
isMasterEdge,
meshEdges,
meshPoints,
pp.edges(),
invSumWeight,
field,
average
);
forAll(field, i)
{
if (field[i] >= fieldMin[i])
{
field[pointI] = average[pointI];
field[i] = (1-mu)*field[i]+mu*average[i];
}
}
// Do residual calculation every so often.
if ((iter % 10) == 0)
{
......@@ -168,7 +258,7 @@ void Foam::autoLayerDriver::smoothField
}
}
}
//XXXXXXXXX
// Smooth normals on moving patch.
void Foam::autoLayerDriver::smoothPatchNormals
......@@ -480,7 +570,6 @@ void Foam::autoLayerDriver::findIsolatedRegions
const PackedBoolList& isMasterEdge,
const labelList& meshEdges,
const scalar minCosLayerTermination,
scalarField& field,
List<extrudeMode>& extrudeStatus,
pointField& patchDisp,
labelList& patchNLayers
......@@ -661,7 +750,6 @@ void Foam::autoLayerDriver::findIsolatedRegions
)
)
{
field[f[fp]] = 0.0;
nPointCounter++;
}
}
......@@ -670,7 +758,8 @@ void Foam::autoLayerDriver::findIsolatedRegions
}
reduce(nPointCounter, sumOp<label>());
Info<< "Number isolated points extrusion stopped : "<< nPointCounter<< endl;
Info<< "Number isolated points extrusion stopped : "<< nPointCounter
<< endl;
}
......@@ -875,11 +964,14 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
//const pointPatchVectorField& pvf =
// meshMover.displacement().boundaryField()[patchI];
if
(
!pp.coupled()
&& !isA<emptyPolyPatch>(pp)
//&& pvf.constraintType() != word::null //exclude fixedValue
&& !adaptPatches.found(patchI)
)
{
......@@ -1168,12 +1260,22 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
isMasterEdge,
meshEdges,
minCosLayerTermination,
thickness,
extrudeStatus,
patchDisp,
patchNLayers
);
// Update thickess for changed extrusion
forAll(thickness, patchPointI)
{
if (extrudeStatus[patchPointI] == NOEXTRUDE)
{
thickness[patchPointI] = 0.0;
}
}
// smooth layer thickness on moving patch
smoothField
(
......@@ -1182,6 +1284,7 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
meshEdges,
minThickness,
nSmoothThickness,
thickness
);
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
View file @ e07309ee
......@@ -506,7 +506,7 @@ void Foam::autoSnapDriver::dumpMove
)
{
// Dump direction of growth into file
Pout<< nl << "Dumping move direction to " << fName << endl;
Info<< "Dumping move direction to " << fName << endl;
OFstream nearestStream(fName);
......@@ -721,14 +721,14 @@ void Foam::autoSnapDriver::preSmoothPatch
if (debug)
{
const_cast<Time&>(mesh.time())++;
Pout<< "Writing patch smoothed mesh to time "
Info<< "Writing patch smoothed mesh to time "
<< meshRefiner_.timeName() << '.' << endl;
meshRefiner_.write
(
debug,
mesh.time().path()/meshRefiner_.timeName()
);
Pout<< "Dumped mesh in = "
Info<< "Dumped mesh in = "
<< mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
}
......@@ -998,7 +998,7 @@ void Foam::autoSnapDriver::smoothDisplacement
if (debug)
{
const_cast<Time&>(mesh.time())++;
Pout<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
Info<< "Writing smoothed mesh to time " << meshRefiner_.timeName()
<< endl;
// Moving mesh creates meshPhi. Can be cleared out by a mesh.clearOut
......@@ -1010,13 +1010,12 @@ void Foam::autoSnapDriver::smoothDisplacement
debug,
mesh.time().path()/meshRefiner_.timeName()
);
Pout<< "Writing displacement field ..." << endl;
Info<< "Writing displacement field ..." << endl;
disp.write();
tmp<pointScalarField> magDisp(mag(disp));
magDisp().write();
Pout<< "Writing actual patch displacement ..." << endl;
Info<< "Writing actual patch displacement ..." << endl;
vectorField actualPatchDisp(disp, pp.meshPoints());
dumpMove
(
......@@ -1068,11 +1067,11 @@ bool Foam::autoSnapDriver::scaleMesh
if (debug)
{
const_cast<Time&>(mesh.time())++;
Pout<< "Writing scaled mesh to time " << meshRefiner_.timeName()
Info<< "Writing scaled mesh to time " << meshRefiner_.timeName()
<< endl;
mesh.write();
Pout<< "Writing displacement field ..." << endl;
Info<< "Writing displacement field ..." << endl;
meshMover.displacement().write();
tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
magDisp().write();
......@@ -1521,14 +1520,14 @@ void Foam::autoSnapDriver::doSnap
if (debug)
{
const_cast<Time&>(mesh.time())++;
Pout<< "Writing scaled mesh to time "
Info<< "Writing scaled mesh to time "
<< meshRefiner_.timeName() << endl;
meshRefiner_.write
(
debug,
mesh.time().path()/meshRefiner_.timeName()
);
Pout<< "Writing displacement field ..." << endl;
Info<< "Writing displacement field ..." << endl;
meshMover.displacement().write();
tmp<pointScalarField> magDisp(mag(meshMover.displacement()));
magDisp().write();
......@@ -1564,7 +1563,7 @@ void Foam::autoSnapDriver::doSnap
if (nChanged > 0 && debug)
{
const_cast<Time&>(mesh.time())++;
Pout<< "Writing patchFace merged mesh to time "
Info<< "Writing patchFace merged mesh to time "
<< meshRefiner_.timeName() << endl;
meshRefiner_.write
(
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriverFeature.C
View file @ e07309ee
......@@ -2565,6 +2565,76 @@ Foam::vectorField Foam::autoSnapDriver::calcNearestSurfaceFeature
}
}
// Collect additionally 'normal' boundary faces for boundaryPoints of pp
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// points on the boundary of pp should pick up non-pp normals
// as well for the feature-reconstruction to behave correctly.
// (the movement is already constrained outside correctly so it
// is only that the unconstrained attraction vector is calculated
// correctly)
{
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
labelList patchID(pbm.patchID());
// Unmark all non-coupled boundary faces
forAll(pbm, patchI)
{
const polyPatch& pp = pbm[patchI];
if (pp.coupled() || isA<emptyPolyPatch>(pp))
{
forAll(pp, i)
{
label meshFaceI = pp.start()+i;
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
}
}
// Remove any meshed faces
PackedBoolList ppFaces(mesh.nFaces());
forAll(pp.addressing(), i)
{
label meshFaceI = pp.addressing()[i];
patchID[meshFaceI-mesh.nInternalFaces()] = -1;
}
// See if pp point uses any non-meshed boundary faces
const labelList& boundaryPoints = pp.boundaryPoints();
forAll(boundaryPoints, i)
{
label pointI = boundaryPoints[i];
label meshPointI = pp.meshPoints()[pointI];
const point& pt = mesh.points()[meshPointI];
const labelList& pFaces = mesh.pointFaces()[meshPointI];
List<point>& pNormals = pointFaceSurfNormals[pointI];
List<point>& pDisp = pointFaceDisp[pointI];
List<point>& pFc = pointFaceCentres[pointI];
labelList& pFid = pointFacePatchID[pointI];
forAll(pFaces, i)
{
label meshFaceI = pFaces[i];
if (!mesh.isInternalFace(meshFaceI))
{
label patchI = patchID[meshFaceI-mesh.nInternalFaces()];
if (patchI != -1)
{
vector fn = mesh.faceAreas()[meshFaceI];
pNormals.append(fn/mag(fn));
pDisp.append(mesh.faceCentres()[meshFaceI]-pt);
pFc.append(mesh.faceCentres()[meshFaceI]);
pFid.append(patchI);
}
}
}
}
}
syncTools::syncPointList
(
mesh,
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.C
View file @ e07309ee
......@@ -43,7 +43,7 @@ Foam::layerParameters::layerParameters
const polyBoundaryMesh& boundaryMesh
)
:
numLayers_(boundaryMesh.size(), 0),
numLayers_(boundaryMesh.size(), -1),
expansionRatio_
(
boundaryMesh.size(),
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/layerParameters/layerParameters.H
View file @ e07309ee
......@@ -131,6 +131,10 @@ public:
// Per patch information
//- How many layers to add.
// -1 : no specification. Assume 0 layers but allow sliding
// to make layers
// 0 : specified to have 0 layers. No sliding allowed.
// >0 : number of layers
const labelList& numLayers() const
{
return numLayers_;
......
src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
View file @ e07309ee
......@@ -571,7 +571,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
hitInfo
);
// Start of from current points
// Start off from current points
newPoints = mesh_.points();
forAll(hitInfo, i)
......@@ -581,6 +581,17 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
}
}
if (debug)
{
const_cast<Time&>(mesh_.time())++;
pointField oldPoints(mesh_.points());
mesh_.movePoints(newPoints);
Pout<< "Writing newPoints mesh to time " << timeName()
<< endl;
write(debug, mesh_.time().path()/"newPoints");
mesh_.movePoints(oldPoints);
}
}
......
Markdown is supported
0% <