From 33afefccde66d83d73ed0eae8b0d875272375ffc Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 1 Oct 2008 19:10:55 +0100 Subject: [PATCH] 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. --- .../mesh/generation/CV3DMesher/CV3D.C | 79 ++++++++++++++++++- .../mesh/generation/CV3DMesher/CV3D.H | 20 +++++ .../mesh/generation/CV3DMesher/CV3DMesher.C | 16 +++- .../mesh/generation/CV3DMesher/calcDualMesh.C | 2 +- .../mesh/generation/CV3DMesher/controls.C | 12 +++ .../mesh/generation/CV3DMesher/indexedCell.H | 6 +- .../insertSurfaceNearestPointPairs.C | 16 ++-- 7 files changed, 133 insertions(+), 18 deletions(-) diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3D.C b/applications/utilities/mesh/generation/CV3DMesher/CV3D.C index f4c4a7155c..2fab856a81 100644 --- a/applications/utilities/mesh/generation/CV3DMesher/CV3D.C +++ b/applications/utilities/mesh/generation/CV3DMesher/CV3D.C @@ -216,7 +216,10 @@ void Foam::CV3D::insertGrid() 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 ( @@ -227,10 +230,82 @@ void Foam::CV3D::relaxPoints(const scalar relaxation) { if (vit->internalPoint()) { - // movePoint(vit, newPoint); + std::list facets; + incident_facets(vit, std::back_inserter(facets)); + + label maxIncidentFacets = 20; + List vertices(maxIncidentFacets); + List edges(maxIncidentFacets); + + point vd(topoint(vit->point())); + + point vi0 = topoint(dual(facets.begin()->first)); + + label edgei = 0; + + for + ( + std::list::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::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; } diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3D.H b/applications/utilities/mesh/generation/CV3DMesher/CV3D.H index afb324c0b2..4c908476ab 100644 --- a/applications/utilities/mesh/generation/CV3DMesher/CV3D.H +++ b/applications/utilities/mesh/generation/CV3DMesher/CV3D.H @@ -79,6 +79,12 @@ 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 // not split scalar minCellSize; @@ -95,6 +101,19 @@ public: // additional "mitering" lines are added 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 Switch insertSurfaceNearestPointPairs; @@ -190,6 +209,7 @@ private: // removing and insertin the surface point-pairs label startOfBoundaryConformPointPairs_; + // Private Member Functions //- Disallow default bitwise copy construct diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C b/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C index 90e58190d0..4104d280cf 100644 --- a/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C +++ b/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C @@ -90,16 +90,22 @@ int main(int argc, char *argv[]) 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 << "Relaxation iteration " << iter << nl << "~~~~~~~~~~~~~~~~~~~~~~~~" << endl; - relaxation -= 0.02; - Info<< "relaxation = " << relaxation << endl; mesh.relaxPoints(relaxation); @@ -107,6 +113,8 @@ int main(int argc, char *argv[]) mesh.removeSurfacePointPairs(); mesh.insertSurfacePointPairs(); mesh.boundaryConform(); + + relaxation -= relaxationDelta; } mesh.write(); diff --git a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C index 8e9e48bc66..a58f715d76 100644 --- a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C +++ b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C @@ -158,7 +158,7 @@ void Foam::CV3D::calcDualMesh } verticesOnFace.append(cc->cellIndex()); - } + } } while (++cc != ccStart); verticesOnFace.shrink(); diff --git a/applications/utilities/mesh/generation/CV3DMesher/controls.C b/applications/utilities/mesh/generation/CV3DMesher/controls.C index c133415451..41bab6789c 100644 --- a/applications/utilities/mesh/generation/CV3DMesher/controls.C +++ b/applications/utilities/mesh/generation/CV3DMesher/controls.C @@ -30,10 +30,22 @@ License Foam::CV3D::controls::controls(const dictionary& controlDict) : + relaxationFactorStart + ( + readScalar(controlDict.lookup("relaxationFactorStart")) + ), + relaxationFactorEnd(readScalar(controlDict.lookup("relaxationFactorEnd"))), minCellSize(readScalar(controlDict.lookup("minCellSize"))), minCellSize2(Foam::sqr(minCellSize)), includedAngle(readScalar(controlDict.lookup("includedAngle"))), maxQuadAngle(readScalar(controlDict.lookup("maxQuadAngle"))), + squares(controlDict.lookup("squares")), + nearWallAlignedDist + ( + readScalar(controlDict.lookup("nearWallAlignedDist"))*minCellSize + ), + nearWallAlignedDist2(Foam::sqr(nearWallAlignedDist)), + relaxOrientation(controlDict.lookup("relaxOrientation")), insertSurfaceNearestPointPairs ( controlDict.lookup("insertSurfaceNearestPointPairs") diff --git a/applications/utilities/mesh/generation/CV3DMesher/indexedCell.H b/applications/utilities/mesh/generation/CV3DMesher/indexedCell.H index a8c4949719..d1810183c9 100644 --- a/applications/utilities/mesh/generation/CV3DMesher/indexedCell.H +++ b/applications/utilities/mesh/generation/CV3DMesher/indexedCell.H @@ -47,7 +47,7 @@ namespace CGAL \*---------------------------------------------------------------------------*/ template > -class indexedCell +class indexedCell : public Cb { @@ -96,11 +96,11 @@ public: indexedCell ( - Vertex_handle v0, + Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3, - Cell_handle n0, + Cell_handle n0, Cell_handle n1, Cell_handle n2, Cell_handle n3 diff --git a/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C b/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C index c4909619e7..6864833f3d 100644 --- a/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C +++ b/applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C @@ -442,8 +442,8 @@ void Foam::CV3D::smoothEdge label nInsertions = label(gap/tols_.maxEdgeSpacing); - Info<< "Gap at start of edge of " << gap - << ". Inserting " << nInsertions << " points" << endl; + // Info<< "Gap at start of edge of " << gap + // << ". Inserting " << nInsertions << " points" << endl; scalar spacing = gap / (nInsertions + 1); @@ -473,8 +473,8 @@ void Foam::CV3D::smoothEdge label nInsertions = label(gap/tols_.maxEdgeSpacing); - Info<< "Gap in edge of " << gap - << ". Inserting " << nInsertions << " points" << endl; + // Info<< "Gap in edge of " << gap + // << ". Inserting " << nInsertions << " points" << endl; scalar spacing = gap / (nInsertions + 1); @@ -493,7 +493,7 @@ void Foam::CV3D::smoothEdge tempEdgePoints.append(edgePoints[eP]); } - // Special treatment for gaps between closest point to start + // Special treatment for gaps between closest point to end if ( @@ -510,8 +510,8 @@ void Foam::CV3D::smoothEdge label nInsertions = label(gap/tols_.maxEdgeSpacing); - Info<< "Gap at end of edge of " << gap - << ". Inserting " << nInsertions << " points" << endl; + // Info<< "Gap at end of edge of " << gap + // << ". Inserting " << nInsertions << " points" << endl; 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()); }