Basic isotropic cell relaxation calculated using the incident_facets of a vertex. Not the way to do it, use the incident edges of the vertex to calculate the poly (dual) faces, use their centres and areas for forcing points and weights.

This commit is contained in:
graham
2008-10-01 19:10:55 +01:00
parent b992b6e1cc
commit 33afefccde
7 changed files with 133 additions and 18 deletions

View File

@ -216,7 +216,10 @@ void Foam::CV3D::insertGrid()
void Foam::CV3D::relaxPoints(const scalar relaxation) void Foam::CV3D::relaxPoints(const scalar relaxation)
{ {
Info<< "Calculating new points: " << nl << endl; Info<< "Calculating new points: " << endl;
vector totalDisp = vector::zero;
scalar totalDist = 0;
for for
( (
@ -227,10 +230,82 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
{ {
if (vit->internalPoint()) if (vit->internalPoint())
{ {
// movePoint(vit, newPoint); std::list<Facet> facets;
incident_facets(vit, std::back_inserter(facets));
label maxIncidentFacets = 20;
List<point> vertices(maxIncidentFacets);
List<vector> edges(maxIncidentFacets);
point vd(topoint(vit->point()));
point vi0 = topoint(dual(facets.begin()->first));
label edgei = 0;
for
(
std::list<Facet>::iterator fit=facets.begin();
fit != facets.end();
++fit
)
{
if
(
is_infinite(fit->first)
|| is_infinite(fit->first->neighbor(fit->second))
)
{
FatalErrorIn("relaxPoints")
<< "Finite cell attached to facet incident on vertex"
<< exit(FatalError);
}
point vi1 = topoint(dual(fit->first->neighbor(fit->second)));
edges[edgei] = vi1 - vi0;
vertices[edgei] = 0.5*(vi1 + vi0);
vi0 = vi1;
edgei++;
}
edgei = 0;
// Initialise the displacement for the centre and sum-weights
vector disp = vector::zero;
scalar sumw = 0;
for
(
std::list<Facet>::iterator fit=facets.begin();
fit != facets.end();
++fit
)
{
vector deltai = vertices[edgei] - vd;
scalar w = 1;
disp += w*deltai;
sumw += w;
edgei++;
}
disp /= sumw;
totalDisp += disp;
totalDist += mag(disp);
movePoint(vit, vd + relaxation*disp);
} }
} }
Info<< "Total displacement = " << totalDisp
<< " total distance = " << totalDist << endl;
} }

View File

@ -79,6 +79,12 @@ public:
{ {
public: public:
//- Relaxation factor at the start of the iteration
scalar relaxationFactorStart;
//- Relaxation factor at the end of the iteration
scalar relaxationFactorEnd;
//- Minimum cell size below which protusions through the surface are //- Minimum cell size below which protusions through the surface are
// not split // not split
scalar minCellSize; scalar minCellSize;
@ -95,6 +101,19 @@ public:
// additional "mitering" lines are added // additional "mitering" lines are added
scalar maxQuadAngle; scalar maxQuadAngle;
//- Should the mesh be square-dominated or of unbiased hexagons
Switch squares;
//- Near-wall region where cells are aligned with the wall
scalar nearWallAlignedDist;
//- Square of nearWallAlignedDist
scalar nearWallAlignedDist2;
//- Chose if the cell orientation should relax during the iterations
// or remain fixed to the x-y directions
Switch relaxOrientation;
//- Insert near-boundary point mirror or point-pairs //- Insert near-boundary point mirror or point-pairs
Switch insertSurfaceNearestPointPairs; Switch insertSurfaceNearestPointPairs;
@ -190,6 +209,7 @@ private:
// removing and insertin the surface point-pairs // removing and insertin the surface point-pairs
label startOfBoundaryConformPointPairs_; label startOfBoundaryConformPointPairs_;
// Private Member Functions // Private Member Functions
//- Disallow default bitwise copy construct //- Disallow default bitwise copy construct

View File

@ -90,16 +90,22 @@ int main(int argc, char *argv[])
mesh.boundaryConform(); mesh.boundaryConform();
} }
scalar relaxation = 1.0; scalar relaxation =
mesh.meshingControls().relaxationFactorStart;
for (int iter=1; iter<=nIterations; iter++) scalar relaxationDelta =
(
mesh.meshingControls().relaxationFactorStart
- mesh.meshingControls().relaxationFactorEnd
)
/max(nIterations, 1);
for (label iter = 0; iter < nIterations; iter++)
{ {
Info<< nl Info<< nl
<< "Relaxation iteration " << iter << nl << "Relaxation iteration " << iter << nl
<< "~~~~~~~~~~~~~~~~~~~~~~~~" << endl; << "~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
relaxation -= 0.02;
Info<< "relaxation = " << relaxation << endl; Info<< "relaxation = " << relaxation << endl;
mesh.relaxPoints(relaxation); mesh.relaxPoints(relaxation);
@ -107,6 +113,8 @@ int main(int argc, char *argv[])
mesh.removeSurfacePointPairs(); mesh.removeSurfacePointPairs();
mesh.insertSurfacePointPairs(); mesh.insertSurfacePointPairs();
mesh.boundaryConform(); mesh.boundaryConform();
relaxation -= relaxationDelta;
} }
mesh.write(); mesh.write();

View File

@ -30,10 +30,22 @@ License
Foam::CV3D::controls::controls(const dictionary& controlDict) Foam::CV3D::controls::controls(const dictionary& controlDict)
: :
relaxationFactorStart
(
readScalar(controlDict.lookup("relaxationFactorStart"))
),
relaxationFactorEnd(readScalar(controlDict.lookup("relaxationFactorEnd"))),
minCellSize(readScalar(controlDict.lookup("minCellSize"))), minCellSize(readScalar(controlDict.lookup("minCellSize"))),
minCellSize2(Foam::sqr(minCellSize)), minCellSize2(Foam::sqr(minCellSize)),
includedAngle(readScalar(controlDict.lookup("includedAngle"))), includedAngle(readScalar(controlDict.lookup("includedAngle"))),
maxQuadAngle(readScalar(controlDict.lookup("maxQuadAngle"))), maxQuadAngle(readScalar(controlDict.lookup("maxQuadAngle"))),
squares(controlDict.lookup("squares")),
nearWallAlignedDist
(
readScalar(controlDict.lookup("nearWallAlignedDist"))*minCellSize
),
nearWallAlignedDist2(Foam::sqr(nearWallAlignedDist)),
relaxOrientation(controlDict.lookup("relaxOrientation")),
insertSurfaceNearestPointPairs insertSurfaceNearestPointPairs
( (
controlDict.lookup("insertSurfaceNearestPointPairs") controlDict.lookup("insertSurfaceNearestPointPairs")

View File

@ -442,8 +442,8 @@ void Foam::CV3D::smoothEdge
label nInsertions = label(gap/tols_.maxEdgeSpacing); label nInsertions = label(gap/tols_.maxEdgeSpacing);
Info<< "Gap at start of edge of " << gap // Info<< "Gap at start of edge of " << gap
<< ". Inserting " << nInsertions << " points" << endl; // << ". Inserting " << nInsertions << " points" << endl;
scalar spacing = gap / (nInsertions + 1); scalar spacing = gap / (nInsertions + 1);
@ -473,8 +473,8 @@ void Foam::CV3D::smoothEdge
label nInsertions = label(gap/tols_.maxEdgeSpacing); label nInsertions = label(gap/tols_.maxEdgeSpacing);
Info<< "Gap in edge of " << gap // Info<< "Gap in edge of " << gap
<< ". Inserting " << nInsertions << " points" << endl; // << ". Inserting " << nInsertions << " points" << endl;
scalar spacing = gap / (nInsertions + 1); scalar spacing = gap / (nInsertions + 1);
@ -493,7 +493,7 @@ void Foam::CV3D::smoothEdge
tempEdgePoints.append(edgePoints[eP]); tempEdgePoints.append(edgePoints[eP]);
} }
// Special treatment for gaps between closest point to start // Special treatment for gaps between closest point to end
if if
( (
@ -510,8 +510,8 @@ void Foam::CV3D::smoothEdge
label nInsertions = label(gap/tols_.maxEdgeSpacing); label nInsertions = label(gap/tols_.maxEdgeSpacing);
Info<< "Gap at end of edge of " << gap // Info<< "Gap at end of edge of " << gap
<< ". Inserting " << nInsertions << " points" << endl; // << ". Inserting " << nInsertions << " points" << endl;
scalar spacing = gap / (nInsertions + 1); scalar spacing = gap / (nInsertions + 1);
@ -579,7 +579,7 @@ void Foam::CV3D::smoothEdge
} }
} }
Info<< edgeI << tab << nPointsRemoved << " points removed." << endl; // Info<< edgeI << tab << nPointsRemoved << " points removed." << endl;
edgePoints.transfer(tempEdgePoints.shrink()); edgePoints.transfer(tempEdgePoints.shrink());
} }