zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Added to smooth surface faces function to only operate on

Browse Source 
perpendicular faces.

Added code to attempt to determine the collapse direction of a face by
detemining the shorted distance across each triangle.  Not working,
commited for reference.

Made some functions const.
This commit is contained in:
graham
2009-12-06 12:07:39 +00:00
parent e00078aa01
commit 241f6a07d4
2 changed files with 204 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files

224
src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
View File

@ -2063,21 +2063,23 @@ void Foam::conformalVoronoiMesh::calcDualMesh
// Assess close points to be merged
Info<< nl << " Merging close points" << endl;
label nPtsMerged = 0;
do
{
Map<label> dualPtIndexMap;
Info<< nl << " Merging close points" << endl;
nPtsMerged = mergeCloseDualVertices(points, dualPtIndexMap);
label nPtsMerged = 0;
Info<< " Merged " << nPtsMerged << " points" << endl;
do
{
Map<label> dualPtIndexMap;
reindexDualVertices(dualPtIndexMap);
nPtsMerged = mergeCloseDualVertices(points, dualPtIndexMap);
} while (nPtsMerged > 0);
Info<< " Merged " << nPtsMerged << " points" << endl;
reindexDualVertices(dualPtIndexMap);
} while (nPtsMerged > 0);
}
// Smooth the surface of the mesh
@ -2091,13 +2093,32 @@ void Foam::conformalVoronoiMesh::calcDualMesh
nSmoothedVertices = smoothSurfaceDualFaces(points, dualPtIndexMap);
Info<< " Smoothed " << nPtsMerged << " points (0 HARD CODED)"
Info<< " Smoothed " << nSmoothedVertices
<< " points (0 HARD CODED)"
<< endl;
reindexDualVertices(dualPtIndexMap);
} while (nSmoothedVertices > 0);
{
Info<< nl << " Merging close points" << endl;
label nPtsMerged = 0;
do
{
Map<label> dualPtIndexMap;
nPtsMerged = mergeCloseDualVertices(points, dualPtIndexMap);
Info<< " Merged " << nPtsMerged << " points" << endl;
reindexDualVertices(dualPtIndexMap);
} while (nPtsMerged > 0);
}
// Assess faces for collapse
Info<< nl << " Collapsing unnecessary faces" << endl;
@ -2387,7 +2408,7 @@ Foam::label Foam::conformalVoronoiMesh::mergeCloseDualVertices
(
const pointField& pts,
Map<label>& dualPtIndexMap
)
) const
{
label nPtsMerged = 0;
@ -2438,10 +2459,12 @@ Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
(
pointField& pts,
Map<label>& dualPtIndexMap
)
) const
{
label nSmoothedVertices = 0;
const scalar cosPerpendicularToleranceAngle = cos(degToRad(80));
for
(
Triangulation::Finite_edges_iterator eit = finite_edges_begin();
@ -2477,29 +2500,59 @@ Foam::label Foam::conformalVoronoiMesh::smoothSurfaceDualFaces
continue;
}
forAll(dualFace, fPtI)
pointIndexHit surfHit;
label hitSurface;
geometryToConformTo_.findSurfaceNearest
(
dualFace.centre(pts),
cvMeshControls().spanSqr(),
surfHit,
hitSurface
);
vectorField norm(1);
allGeometry_[hitSurface].getNormal
(
List<pointIndexHit>(1, surfHit),
norm
);
const vector& surfaceNormal = norm[0];
vector faceNormal = dualFace.normal(pts);
faceNormal /= mag(faceNormal);
if ((faceNormal & surfaceNormal) < cosPerpendicularToleranceAngle)
{
label ptI = dualFace[fPtI];
collapseFaceToEdge(dualFace, pts, dualPtIndexMap);
point& pt = pts[ptI];
// forAll(dualFace, fPtI)
// {
// label ptI = dualFace[fPtI];
pointIndexHit surfHit;
label hitSurface;
// point& pt = pts[ptI];
geometryToConformTo_.findSurfaceNearest
(
pt,
cvMeshControls().spanSqr(),
surfHit,
hitSurface
);
// pointIndexHit surfHit;
// label hitSurface;
if (surfHit.hit())
{
pt = surfHit.hitPoint();
// geometryToConformTo_.findSurfaceNearest
// (
// pt,
// cvMeshControls().spanSqr(),
// surfHit,
// hitSurface
// );
// dualPtIndexMap.insert(ptI, ptI);
}
// if (surfHit.hit())
// {
// pt = surfHit.hitPoint();
// // dualPtIndexMap.insert(ptI, ptI);
// }
// }
}
}
}
@ -2512,7 +2565,7 @@ Foam::label Foam::conformalVoronoiMesh::collapseFaces
(
pointField& pts,
Map<label>& dualPtIndexMap
)
) const
{
label nCollapsedFaces = 0;
@ -2750,6 +2803,115 @@ Foam::label Foam::conformalVoronoiMesh::collapseFaces
}
void Foam::conformalVoronoiMesh::collapseFaceToEdge
(
const face& f,
pointField& pts,
Map<label>& dualPtIndexMap
) const
{
const vector fC = f.centre(pts);
const edgeList& faceEds = f.edges();
const scalar fArea = f.mag(pts);
scalar averageQuality = 0.0;
vector faceCollapseDirection = vector::zero;
forAll(faceEds, edI)
{
const edge e = faceEds[edI];
triPointRef tri(pts[e.start()], pts[e.end()], fC);
scalar triArea = tri.mag();
if (triArea < SMALL)
{
// Sliver triangle, do not consider
continue;
}
// 0.7697997 converts from a quality being 1 for an area ratio
// of an equilateral triangle (4*pi/(3*sqrt(3))) to the
// quality being 1 for an area ratio of a quarter of a square
// (pi), i.e. (4*pi/(3*sqrt(3)))/pi = 0.7697997
scalar triQuality = tri.quality()*0.7697997;
// Weight the quality by the tri area
averageQuality += triQuality*triArea;
// find the longest edge vector, and one of the other edges
// that isn't the longest
vector longestEdge = tri.b() - tri.a();
vector otherEdge = tri.c() - tri.a();
scalar longestEdgeMagSqr = magSqr(longestEdge);
if (magSqr(tri.c() - tri.b()) > longestEdgeMagSqr)
{
longestEdge = tri.c() - tri.b();
otherEdge = tri.a() - tri.b();
longestEdgeMagSqr = magSqr(longestEdgeMagSqr);
}
if (magSqr(tri.a() - tri.c()) > longestEdgeMagSqr)
{
longestEdge = tri.a() - tri.c();
otherEdge = tri.b() - tri.c();
longestEdgeMagSqr = magSqr(longestEdgeMagSqr);
}
vector triCollapseDirection =
(otherEdge & longestEdge)*longestEdge/longestEdgeMagSqr
- otherEdge;
triCollapseDirection /= mag(triCollapseDirection);
// Add the triCollapseDirection so that it is constructive
if
(
magSqr(faceCollapseDirection + triCollapseDirection)
> magSqr(faceCollapseDirection)
)
{
faceCollapseDirection += triCollapseDirection;
// faceCollapseDirection += triCollapseDirection/triQuality;
}
else
{
faceCollapseDirection -= triCollapseDirection;
// faceCollapseDirection -= triCollapseDirection/triQuality;
}
}
faceCollapseDirection /= mag(faceCollapseDirection);
averageQuality /= fArea;
vector collapseAxis = (faceCollapseDirection ^ (f.normal(pts)/fArea));
Info<< nl << "# averageQuality " << averageQuality << endl;
forAll(f, fPtI)
{
meshTools::writeOBJ(Info, pts[f[fPtI]]);
}
collapseAxis *= 3.0*mag(pts[f[0]] - fC);
faceCollapseDirection *= 3.0*mag(pts[f[0]] - fC);
meshTools::writeOBJ(Info, fC + faceCollapseDirection);
meshTools::writeOBJ(Info, fC - faceCollapseDirection);
meshTools::writeOBJ(Info, fC + collapseAxis);
meshTools::writeOBJ(Info, fC - collapseAxis);
}
void Foam::conformalVoronoiMesh::reindexDualVertices
(
const Map<label>& dualPtIndexMap

14
src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
View File

@ -468,7 +468,7 @@ private:
(
const pointField& pts,
Map<label>& dualPtIndexMap
);
) const;
//- Smooth the dual vertices of the dual faces on the boundary
// so that they conform to the surface and remove any
@ -477,7 +477,7 @@ private:
(
pointField& pts,
Map<label>& dualPtIndexMap
);
) const;
//- Collapse a face to a point or an edge, modifying and
// mapping the points, returns the true if the face was
@ -486,7 +486,15 @@ private:
(
pointField& pts,
Map<label>& dualPtIndexMap
);
) const;
//- Collapse a face to an edge, updating the point and point map
void collapseFaceToEdge
(
const face& f,
pointField& pts,
Map<label>& dualPtIndexMap
) const;
//- Re-index all of the the Delaunay cells
void reindexDualVertices(const Map<label>& dualPtIndexMap);