Commit 5154e387 authored by mattijs's avatar mattijs
Browse files

ENH: autoLayerShrink: added medial axis smoothing

parent 35b9d636
......@@ -1473,6 +1473,11 @@ void Foam::autoLayerDriver::getPatchDisplacement
// Start off from same thickness everywhere (except where no extrusion)
patchDisp = thickness*pointNormals;
label nNoVisNormal = 0;
label nExtrudeRemove = 0;
// Check if no extrude possible.
forAll(pointNormals, patchPointI)
{
......@@ -1491,12 +1496,17 @@ void Foam::autoLayerDriver::getPatchDisplacement
if (!meshTools::visNormal(n, faceNormals, pointFaces[patchPointI]))
{
Pout<< "No valid normal for point " << meshPointI
<< ' ' << pp.points()[meshPointI]
<< "; setting displacement to " << patchDisp[patchPointI]
<< endl;
if (debug&meshRefinement::OBJINTERSECTIONS)
{
Pout<< "No valid normal for point " << meshPointI
<< ' ' << pp.points()[meshPointI]
<< "; setting displacement to "
<< patchDisp[patchPointI]
<< endl;
}
extrudeStatus[patchPointI] = EXTRUDEREMOVE;
nNoVisNormal++;
}
}
}
......@@ -1526,18 +1536,30 @@ void Foam::autoLayerDriver::getPatchDisplacement
if (nPoints > 0)
{
Pout<< "Displacement at illegal point "
<< localPoints[patchPointI]
<< " set to " << (avg / nPoints - localPoints[patchPointI])
<< endl;
if (debug&meshRefinement::OBJINTERSECTIONS)
{
Pout<< "Displacement at illegal point "
<< localPoints[patchPointI]
<< " set to "
<< (avg / nPoints - localPoints[patchPointI])
<< endl;
}
patchDisp[patchPointI] =
avg / nPoints
- localPoints[patchPointI];
nExtrudeRemove++;
}
}
}
Info<< "Detected " << returnReduce(nNoVisNormal, sumOp<label>())
<< " points with point normal pointing through faces." << nl
<< "Reset displacement at "
<< returnReduce(nExtrudeRemove, sumOp<label>())
<< " points to average of surrounding points." << endl;
// Make sure displacement is equal on both sides of coupled patches.
syncPatchDisplacement
(
......@@ -2849,6 +2871,7 @@ void Foam::autoLayerDriver::addLayers
baffles,
layerParams.nSmoothThickness(),
layerParams.nSmoothDisplacement(),
layerParams.maxThicknessToMedialRatio(),
nAllowableErrors,
layerParams.nSnap(),
......
......@@ -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
......@@ -29,6 +29,7 @@ Description
SourceFiles
autoLayerDriver.C
autoLayerDriverShrink.C
\*---------------------------------------------------------------------------*/
......@@ -482,7 +483,8 @@ class autoLayerDriver
motionSmoother& meshMover,
const dictionary& meshQualityDict,
const List<labelPair>& baffles,
const label nSmoothThickness,
const label nSmoothPatchThickness,
const label nSmoothDisplacement,
const scalar maxThicknessToMedialRatio,
const label nAllowableErrors,
const label nSnap,
......
......@@ -1210,7 +1210,8 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
motionSmoother& meshMover,
const dictionary& meshQualityDict,
const List<labelPair>& baffles,
const label nSmoothThickness,
const label nSmoothPatchThickness,
const label nSmoothDisplacement,
const scalar maxThicknessToMedialRatio,
const label nAllowableErrors,
const label nSnap,
......@@ -1412,7 +1413,7 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
isMasterEdge,
meshEdges,
minThickness,
nSmoothThickness,
nSmoothPatchThickness,
thickness
);
......@@ -1473,6 +1474,90 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
// Smear displacement away from fixed values (medialRatio=0 or 1)
if (nSmoothDisplacement > 0)
{
scalarField invSumWeight(mesh.nPoints());
sumWeights
(
isMasterEdge,
identity(mesh.nEdges()),
identity(mesh.nPoints()),
mesh.edges(),
invSumWeight
);
// Get smoothly varying patch field.
Info<< "shrinkMeshDistance : Smoothing displacement ..." << endl;
const scalar lambda = 0.33;
const scalar mu = -0.34;
pointField average(mesh.nPoints());
for (label iter = 0; iter < nSmoothDisplacement; iter++)
{
// Calculate average of field
averageNeighbours
(
mesh,
isMasterEdge,
identity(mesh.nEdges()), //meshEdges,
identity(mesh.nPoints()), //meshPoints,
mesh.edges(), //edges,
invSumWeight,
displacement,
average
);
forAll(displacement, i)
{
if (medialRatio[i] > SMALL && medialRatio[i] < 1-SMALL)
{
displacement[i] =
(1-lambda)*displacement[i]
+lambda*average[i];
}
}
// Calculate average of field
averageNeighbours
(
mesh,
isMasterEdge,
identity(mesh.nEdges()), //meshEdges,
identity(mesh.nPoints()), //meshPoints,
mesh.edges(), //edges,
invSumWeight,
displacement,
average
);
forAll(displacement, i)
{
if (medialRatio[i] > SMALL && medialRatio[i] < 1-SMALL)
{
displacement[i] = (1-mu)*displacement[i]+mu*average[i];
}
}
// Do residual calculation every so often.
if ((iter % 10) == 0)
{
Info<< " Iteration " << iter << " residual "
<< gSum(mag(displacement-average))
/returnReduce(average.size(), sumOp<label>())
<< endl;
}
}
}
// Make sure displacement boundary conditions is uptodate with
// internal field
meshMover.setDisplacementPatchFields();
// Check a bit the sync of displacements
if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
......@@ -1509,92 +1594,6 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance
}
//XXXXX
// // Smear displacement away from fixed values (medialRatio=0 or 1)
// {
// const edgeList& edges = mesh.edges();
// scalarField edgeWeight(edges.size(), 0.0);
// forAll(edges, edgeI)
// {
// if (isMasterEdge[edgeI])
// {
// scalar eMag = edges[edgeI].mag(mesh.points());
// if (eMag > VSMALL)
// {
// edgeWeight[edgeI] = 1.0/eMag;
// }
// else
// {
// edgeWeight[edgeI] = GREAT;
// }
// }
// }
// scalarField invSumWeight(mesh.nPoints());
// sumWeights(isMasterEdge, edgeWeight, invSumWeight);
//
//
// // Get smoothly varying patch field.
// Info<< "shrinkMeshDistance : Smoothing displacement ..." << endl;
//
// const scalar lambda = 0.33;
// const scalar mu = -0.34;
//
// pointField average(mesh.nPoints());
// for (label iter = 0; iter < 90; iter++)
// {
// // Calculate average of field
// averageNeighbours
// (
// mesh,
// edgeWeight,
// invSumWeight,
// displacement,
// average
// );
//
// forAll(displacement, i)
// {
// if (medialRatio[i] > SMALL && medialRatio[i] < 1-SMALL)
// {
// displacement[i] =
// (1-lambda)*displacement[i]
// +lambda*average[i];
// }
// }
//
//
// // Calculate average of field
// averageNeighbours
// (
// mesh,
// edgeWeight,
// invSumWeight,
// displacement,
// average
// );
//
// forAll(displacement, i)
// {
// if (medialRatio[i] > SMALL && medialRatio[i] < 1-SMALL)
// {
// displacement[i] = (1-mu)*displacement[i]+mu*average[i];
// }
// }
//
//
// // Do residual calculation every so often.
// if ((iter % 10) == 0)
// {
// Info<< " Iteration " << iter << " residual "
// << gSum(mag(displacement-average))
// /returnReduce(average.size(), sumOp<label>())
// << endl;
// }
// }
// }
//XXXXX
if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
{
const_cast<Time&>(mesh.time())++;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment