1051 lines
27 KiB
C
1051 lines
27 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
\*----------------------------------------------------------------------------*/
|
|
|
|
#include "CV2D.H"
|
|
#include "Random.H"
|
|
#include "transform.H"
|
|
#include "IFstream.H"
|
|
#include "uint.H"
|
|
#include "ulong.H"
|
|
|
|
namespace Foam
|
|
{
|
|
defineTypeNameAndDebug(CV2D, 0);
|
|
}
|
|
|
|
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
|
|
|
void Foam::CV2D::insertBoundingBox()
|
|
{
|
|
Info<< "insertBoundingBox: creating bounding mesh" << endl;
|
|
scalar bigSpan = 10*meshControls().span();
|
|
insertPoint(point2D(-bigSpan, -bigSpan), Vb::FAR_POINT);
|
|
insertPoint(point2D(-bigSpan, bigSpan), Vb::FAR_POINT);
|
|
insertPoint(point2D(bigSpan, -bigSpan), Vb::FAR_POINT);
|
|
insertPoint(point2D(bigSpan, bigSpan), Vb::FAR_POINT);
|
|
}
|
|
|
|
|
|
void Foam::CV2D::fast_restore_Delaunay(Vertex_handle vh)
|
|
{
|
|
int i;
|
|
Face_handle f = vh->face(), next, start(f);
|
|
|
|
do
|
|
{
|
|
i=f->index(vh);
|
|
if (!is_infinite(f))
|
|
{
|
|
if (!internal_flip(f, cw(i))) external_flip(f, i);
|
|
if (f->neighbor(i) == start) start = f;
|
|
}
|
|
f = f->neighbor(cw(i));
|
|
} while (f != start);
|
|
}
|
|
|
|
|
|
void Foam::CV2D::external_flip(Face_handle& f, int i)
|
|
{
|
|
Face_handle n = f->neighbor(i);
|
|
|
|
if
|
|
(
|
|
CGAL::ON_POSITIVE_SIDE
|
|
!= side_of_oriented_circle(n, f->vertex(i)->point())
|
|
) return;
|
|
|
|
flip(f, i);
|
|
i = n->index(f->vertex(i));
|
|
external_flip(n, i);
|
|
}
|
|
|
|
|
|
bool Foam::CV2D::internal_flip(Face_handle& f, int i)
|
|
{
|
|
Face_handle n = f->neighbor(i);
|
|
|
|
if
|
|
(
|
|
CGAL::ON_POSITIVE_SIDE
|
|
!= side_of_oriented_circle(n, f->vertex(i)->point())
|
|
)
|
|
{
|
|
return false;
|
|
}
|
|
|
|
flip(f, i);
|
|
|
|
return true;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
|
|
|
|
Foam::CV2D::CV2D
|
|
(
|
|
const Time& runTime,
|
|
const dictionary& cvMeshDict
|
|
)
|
|
:
|
|
Delaunay(),
|
|
runTime_(runTime),
|
|
rndGen_(64293*Pstream::myProcNo()),
|
|
allGeometry_
|
|
(
|
|
IOobject
|
|
(
|
|
"cvSearchableSurfaces",
|
|
runTime_.constant(),
|
|
"triSurface",
|
|
runTime_,
|
|
IOobject::MUST_READ,
|
|
IOobject::NO_WRITE
|
|
),
|
|
cvMeshDict.subDict("geometry")
|
|
),
|
|
qSurf_
|
|
(
|
|
runTime_,
|
|
rndGen_,
|
|
allGeometry_,
|
|
cvMeshDict.subDict("surfaceConformation")
|
|
),
|
|
controls_(cvMeshDict, qSurf_.globalBounds()),
|
|
cellSizeControl_
|
|
(
|
|
allGeometry_,
|
|
cvMeshDict.subDict("motionControl")
|
|
),
|
|
z_
|
|
(
|
|
point
|
|
(
|
|
cvMeshDict.subDict("surfaceConformation").lookup("locationInMesh")
|
|
).z()
|
|
),
|
|
startOfInternalPoints_(0),
|
|
startOfSurfacePointPairs_(0),
|
|
startOfBoundaryConformPointPairs_(0),
|
|
featurePoints_()
|
|
{
|
|
Info<< meshControls() << endl;
|
|
|
|
insertBoundingBox();
|
|
insertFeaturePoints();
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
|
|
|
|
Foam::CV2D::~CV2D()
|
|
{}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
|
|
|
|
void Foam::CV2D::insertPoints
|
|
(
|
|
const point2DField& points,
|
|
const scalar nearness
|
|
)
|
|
{
|
|
Info<< "insertInitialPoints(const point2DField& points): ";
|
|
|
|
startOfInternalPoints_ = number_of_vertices();
|
|
label nVert = startOfInternalPoints_;
|
|
|
|
// Add the points and index them
|
|
forAll(points, i)
|
|
{
|
|
const point2D& p = points[i];
|
|
|
|
if (qSurf_.wellInside(toPoint3D(p), nearness))
|
|
{
|
|
insert(toPoint(p))->index() = nVert++;
|
|
}
|
|
else
|
|
{
|
|
Warning
|
|
<< "Rejecting point " << p << " outside surface" << endl;
|
|
}
|
|
}
|
|
|
|
Info<< nVert << " vertices inserted" << endl;
|
|
|
|
if (meshControls().objOutput())
|
|
{
|
|
// Checking validity of triangulation
|
|
assert(is_valid());
|
|
|
|
writeTriangles("initial_triangles.obj", true);
|
|
writeFaces("initial_faces.obj", true);
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::CV2D::insertPoints(const fileName& pointFileName)
|
|
{
|
|
IFstream pointsFile(pointFileName);
|
|
|
|
if (pointsFile.good())
|
|
{
|
|
insertPoints
|
|
(
|
|
point2DField(pointsFile),
|
|
0.5*meshControls().minCellSize2()
|
|
);
|
|
}
|
|
else
|
|
{
|
|
FatalErrorIn("insertInitialPoints")
|
|
<< "Could not open pointsFile " << pointFileName
|
|
<< exit(FatalError);
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::CV2D::insertGrid()
|
|
{
|
|
Info<< "insertInitialGrid: ";
|
|
|
|
startOfInternalPoints_ = number_of_vertices();
|
|
label nVert = startOfInternalPoints_;
|
|
|
|
scalar x0 = qSurf_.globalBounds().min().x();
|
|
scalar xR = qSurf_.globalBounds().max().x() - x0;
|
|
int ni = int(xR/meshControls().minCellSize()) + 1;
|
|
scalar deltax = xR/ni;
|
|
|
|
scalar y0 = qSurf_.globalBounds().min().y();
|
|
scalar yR = qSurf_.globalBounds().max().y() - y0;
|
|
int nj = int(yR/meshControls().minCellSize()) + 1;
|
|
scalar deltay = yR/nj;
|
|
|
|
Random rndGen(1321);
|
|
scalar pert = meshControls().randomPerturbation()*min(deltax, deltay);
|
|
|
|
for (int i=0; i<ni; i++)
|
|
{
|
|
for (int j=0; j<nj; j++)
|
|
{
|
|
point p(x0 + i*deltax, y0 + j*deltay, 0);
|
|
|
|
if (meshControls().randomiseInitialGrid())
|
|
{
|
|
p.x() += pert*(rndGen.scalar01() - 0.5);
|
|
p.y() += pert*(rndGen.scalar01() - 0.5);
|
|
}
|
|
|
|
if (qSurf_.wellInside(p, 0.5*meshControls().minCellSize2()))
|
|
{
|
|
insert(Point(p.x(), p.y()))->index() = nVert++;
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< nVert << " vertices inserted" << endl;
|
|
|
|
if (meshControls().objOutput())
|
|
{
|
|
// Checking validity of triangulation
|
|
assert(is_valid());
|
|
|
|
writeTriangles("initial_triangles.obj", true);
|
|
writeFaces("initial_faces.obj", true);
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::CV2D::insertSurfacePointPairs()
|
|
{
|
|
startOfSurfacePointPairs_ = number_of_vertices();
|
|
|
|
if (meshControls().insertSurfaceNearestPointPairs())
|
|
{
|
|
insertSurfaceNearestPointPairs();
|
|
}
|
|
|
|
write("nearest");
|
|
|
|
// Insertion of point-pairs for near-points may cause protrusions
|
|
// so insertBoundaryConformPointPairs must be executed last
|
|
if (meshControls().insertSurfaceNearPointPairs())
|
|
{
|
|
insertSurfaceNearPointPairs();
|
|
}
|
|
|
|
startOfBoundaryConformPointPairs_ = number_of_vertices();
|
|
}
|
|
|
|
|
|
void Foam::CV2D::boundaryConform()
|
|
{
|
|
if (!meshControls().insertSurfaceNearestPointPairs())
|
|
{
|
|
markNearBoundaryPoints();
|
|
}
|
|
|
|
// Mark all the faces as SAVE_CHANGED
|
|
for
|
|
(
|
|
Triangulation::Finite_faces_iterator fit = finite_faces_begin();
|
|
fit != finite_faces_end();
|
|
fit++
|
|
)
|
|
{
|
|
fit->faceIndex() = Fb::SAVE_CHANGED;
|
|
}
|
|
|
|
for (label iter=1; iter<=meshControls().maxBoundaryConformingIter(); iter++)
|
|
{
|
|
label nIntersections = insertBoundaryConformPointPairs
|
|
(
|
|
"surfaceIntersections_" + Foam::name(iter) + ".obj"
|
|
);
|
|
|
|
if (nIntersections == 0)
|
|
{
|
|
break;
|
|
}
|
|
else
|
|
{
|
|
Info<< "BC iteration " << iter << ": "
|
|
<< nIntersections << " point-pairs inserted" << endl;
|
|
}
|
|
|
|
// Any faces changed by insertBoundaryConformPointPairs will now
|
|
// be marked CHANGED, mark those as SAVE_CHANGED and those that
|
|
// remained SAVE_CHANGED as UNCHANGED
|
|
for
|
|
(
|
|
Triangulation::Finite_faces_iterator fit = finite_faces_begin();
|
|
fit != finite_faces_end();
|
|
fit++
|
|
)
|
|
{
|
|
if (fit->faceIndex() == Fb::SAVE_CHANGED)
|
|
{
|
|
fit->faceIndex() = Fb::UNCHANGED;
|
|
}
|
|
else if (fit->faceIndex() == Fb::CHANGED)
|
|
{
|
|
fit->faceIndex() = Fb::SAVE_CHANGED;
|
|
}
|
|
}
|
|
}
|
|
|
|
Info<< nl;
|
|
|
|
write("boundary");
|
|
}
|
|
|
|
|
|
void Foam::CV2D::removeSurfacePointPairs()
|
|
{
|
|
for
|
|
(
|
|
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
|
|
vit != finite_vertices_end();
|
|
++vit
|
|
)
|
|
{
|
|
if (vit->index() >= startOfSurfacePointPairs_)
|
|
{
|
|
remove(vit);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::CV2D::newPoints(const scalar relaxation)
|
|
{
|
|
Info<< "newPointsFromVertices: ";
|
|
|
|
Field<point2D> dualVertices(number_of_faces());
|
|
|
|
label dualVerti = 0;
|
|
|
|
// Find the dual point of each tetrahedron and assign it an index.
|
|
for
|
|
(
|
|
Triangulation::Finite_faces_iterator fit = finite_faces_begin();
|
|
fit != finite_faces_end();
|
|
++fit
|
|
)
|
|
{
|
|
fit->faceIndex() = -1;
|
|
|
|
if
|
|
(
|
|
fit->vertex(0)->internalOrBoundaryPoint()
|
|
|| fit->vertex(1)->internalOrBoundaryPoint()
|
|
|| fit->vertex(2)->internalOrBoundaryPoint()
|
|
)
|
|
{
|
|
fit->faceIndex() = dualVerti;
|
|
|
|
dualVertices[dualVerti] = toPoint2D(circumcenter(fit));
|
|
|
|
dualVerti++;
|
|
}
|
|
}
|
|
|
|
dualVertices.setSize(dualVerti);
|
|
|
|
Field<vector2D> displacementAccumulator
|
|
(
|
|
startOfSurfacePointPairs_,
|
|
vector2D::zero
|
|
);
|
|
|
|
// Calculate target size and alignment for vertices
|
|
scalarField sizes
|
|
(
|
|
number_of_vertices(),
|
|
meshControls().minCellSize()
|
|
);
|
|
|
|
Field<vector2D> alignments
|
|
(
|
|
number_of_vertices(),
|
|
vector2D(1, 0)
|
|
);
|
|
|
|
for
|
|
(
|
|
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
|
|
vit != finite_vertices_end();
|
|
++vit
|
|
)
|
|
{
|
|
if (vit->internalOrBoundaryPoint())
|
|
{
|
|
point2D vert = toPoint2D(vit->point());
|
|
|
|
// alignment and size determination
|
|
pointIndexHit pHit;
|
|
label hitSurface = -1;
|
|
|
|
qSurf_.findSurfaceNearest
|
|
(
|
|
toPoint3D(vert),
|
|
meshControls().span2(),
|
|
pHit,
|
|
hitSurface
|
|
);
|
|
|
|
if (pHit.hit())
|
|
{
|
|
vectorField norm(1);
|
|
allGeometry_[hitSurface].getNormal
|
|
(
|
|
List<pointIndexHit>(1, pHit),
|
|
norm
|
|
);
|
|
|
|
alignments[vit->index()] = toPoint2D(norm[0]);
|
|
|
|
sizes[vit->index()] =
|
|
cellSizeControl_.cellSize(toPoint3D(vit->point()));
|
|
}
|
|
}
|
|
}
|
|
|
|
// Info<< "Calculated alignments" << endl;
|
|
|
|
scalar cosAlignmentAcceptanceAngle = 0.68;
|
|
|
|
// Upper and lower edge length ratios for weight
|
|
scalar u = 1.0;
|
|
scalar l = 0.7;
|
|
|
|
PackedBoolList pointToBeRetained(startOfSurfacePointPairs_, true);
|
|
|
|
std::list<Point> pointsToInsert;
|
|
|
|
for
|
|
(
|
|
Triangulation::Finite_edges_iterator eit = finite_edges_begin();
|
|
eit != finite_edges_end();
|
|
eit++
|
|
)
|
|
{
|
|
Vertex_handle vA = eit->first->vertex(cw(eit->second));
|
|
Vertex_handle vB = eit->first->vertex(ccw(eit->second));
|
|
|
|
if (!vA->internalOrBoundaryPoint() || !vB->internalOrBoundaryPoint())
|
|
{
|
|
continue;
|
|
}
|
|
|
|
const point2D& dualV1 = dualVertices[eit->first->faceIndex()];
|
|
const point2D& dualV2 =
|
|
dualVertices[eit->first->neighbor(eit->second)->faceIndex()];
|
|
|
|
scalar dualEdgeLength = mag(dualV1 - dualV2);
|
|
|
|
point2D dVA = toPoint2D(vA->point());
|
|
point2D dVB = toPoint2D(vB->point());
|
|
|
|
Field<vector2D> alignmentDirsA(2);
|
|
|
|
alignmentDirsA[0] = alignments[vA->index()];
|
|
alignmentDirsA[1] = vector2D
|
|
(
|
|
-alignmentDirsA[0].y(),
|
|
alignmentDirsA[0].x()
|
|
);
|
|
|
|
Field<vector2D> alignmentDirsB(2);
|
|
|
|
alignmentDirsB[0] = alignments[vB->index()];
|
|
alignmentDirsB[1] = vector2D
|
|
(
|
|
-alignmentDirsB[0].y(),
|
|
alignmentDirsB[0].x()
|
|
);
|
|
|
|
Field<vector2D> alignmentDirs(2);
|
|
|
|
forAll(alignmentDirsA, aA)
|
|
{
|
|
const vector2D& a(alignmentDirsA[aA]);
|
|
|
|
scalar maxDotProduct = 0.0;
|
|
|
|
forAll(alignmentDirsB, aB)
|
|
{
|
|
const vector2D& b(alignmentDirsB[aB]);
|
|
|
|
scalar dotProduct = a & b;
|
|
|
|
if (mag(dotProduct) > maxDotProduct)
|
|
{
|
|
maxDotProduct = mag(dotProduct);
|
|
|
|
alignmentDirs[aA] = a + sign(dotProduct)*b;
|
|
|
|
alignmentDirs[aA] /= mag(alignmentDirs[aA]);
|
|
}
|
|
}
|
|
}
|
|
|
|
vector2D rAB = dVA - dVB;
|
|
|
|
scalar rABMag = mag(rAB);
|
|
|
|
forAll(alignmentDirs, aD)
|
|
{
|
|
vector2D& alignmentDir = alignmentDirs[aD];
|
|
|
|
if ((rAB & alignmentDir) < 0)
|
|
{
|
|
// swap the direction of the alignment so that has the
|
|
// same sense as rAB
|
|
alignmentDir *= -1;
|
|
}
|
|
|
|
scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir);
|
|
|
|
if (alignmentDotProd > cosAlignmentAcceptanceAngle)
|
|
{
|
|
scalar targetFaceSize =
|
|
0.5*(sizes[vA->index()] + sizes[vB->index()]);
|
|
|
|
// Test for changing aspect ratio on second alignment (first
|
|
// alignment is neartest surface normal)
|
|
// if (aD == 1)
|
|
// {
|
|
// targetFaceSize *= 2.0;
|
|
// }
|
|
|
|
alignmentDir *= 0.5*targetFaceSize;
|
|
|
|
vector2D delta = alignmentDir - 0.5*rAB;
|
|
|
|
if (dualEdgeLength < 0.7*targetFaceSize)
|
|
{
|
|
delta *= 0;
|
|
}
|
|
else if (dualEdgeLength < targetFaceSize)
|
|
{
|
|
delta *=
|
|
(
|
|
dualEdgeLength
|
|
/(targetFaceSize*(u - l))
|
|
- 1/((u/l) - 1)
|
|
);
|
|
}
|
|
|
|
if
|
|
(
|
|
vA->internalPoint()
|
|
&& vB->internalPoint()
|
|
&& rABMag > 1.75*targetFaceSize
|
|
&& dualEdgeLength > 0.05*targetFaceSize
|
|
&& alignmentDotProd > 0.93
|
|
)
|
|
{
|
|
// Point insertion
|
|
pointsToInsert.push_back(toPoint(0.5*(dVA + dVB)));
|
|
}
|
|
else if
|
|
(
|
|
(vA->internalPoint() || vB->internalPoint())
|
|
&& rABMag < 0.65*targetFaceSize
|
|
)
|
|
{
|
|
// Point removal
|
|
|
|
// Only insert a point at the midpoint of the short edge
|
|
// if neither attached point has already been identified
|
|
// to be removed.
|
|
if
|
|
(
|
|
pointToBeRetained[vA->index()] == true
|
|
&& pointToBeRetained[vB->index()] == true
|
|
)
|
|
{
|
|
pointsToInsert.push_back(toPoint(0.5*(dVA + dVB)));
|
|
}
|
|
|
|
if (vA->internalPoint())
|
|
{
|
|
pointToBeRetained[vA->index()] = false;
|
|
}
|
|
|
|
if (vB->internalPoint())
|
|
{
|
|
pointToBeRetained[vB->index()] = false;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (vA->internalPoint())
|
|
{
|
|
displacementAccumulator[vA->index()] += delta;
|
|
}
|
|
|
|
if (vB->internalPoint())
|
|
{
|
|
displacementAccumulator[vB->index()] += -delta;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
vector2D totalDisp = sum(displacementAccumulator);
|
|
scalar totalDist = sum(mag(displacementAccumulator));
|
|
|
|
// Relax the calculated displacement
|
|
displacementAccumulator *= relaxation;
|
|
|
|
label numberOfNewPoints = pointsToInsert.size();
|
|
|
|
for
|
|
(
|
|
Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
|
|
vit != finite_vertices_end();
|
|
++vit
|
|
)
|
|
{
|
|
if (vit->internalPoint())
|
|
{
|
|
if (pointToBeRetained[vit->index()])
|
|
{
|
|
pointsToInsert.push_front
|
|
(
|
|
toPoint
|
|
(
|
|
toPoint2D(vit->point())
|
|
+ displacementAccumulator[vit->index()]
|
|
)
|
|
);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Clear the triangulation and reinsert the bounding box and feature points.
|
|
// This is faster than removing and moving points.
|
|
this->clear();
|
|
|
|
insertBoundingBox();
|
|
|
|
reinsertFeaturePoints();
|
|
|
|
startOfInternalPoints_ = number_of_vertices();
|
|
|
|
label nVert = startOfInternalPoints_;
|
|
|
|
Info<< "Inserting " << numberOfNewPoints << " new points" << endl;
|
|
|
|
// Use the range insert as it is faster than individually inserting points.
|
|
insert(pointsToInsert.begin(), pointsToInsert.end());
|
|
|
|
for
|
|
(
|
|
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
|
|
vit != finite_vertices_end();
|
|
++vit
|
|
)
|
|
{
|
|
if
|
|
(
|
|
vit->type() == Vb::INTERNAL_POINT
|
|
&& vit->index() == Vb::INTERNAL_POINT
|
|
)
|
|
{
|
|
vit->index() = nVert++;
|
|
}
|
|
}
|
|
|
|
Info<< " Total displacement = " << totalDisp << nl
|
|
<< " Total distance = " << totalDist << nl
|
|
<< " Points added = " << pointsToInsert.size()
|
|
<< endl;
|
|
|
|
write("internal");
|
|
|
|
insertSurfacePointPairs();
|
|
|
|
boundaryConform();
|
|
|
|
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
// Old Method
|
|
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
|
|
|
// for
|
|
// (
|
|
// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
|
|
// vit != finite_vertices_end();
|
|
// ++vit
|
|
// )
|
|
// {
|
|
// if (vit->internalPoint())
|
|
// {
|
|
// // Current dual-cell defining vertex ("centre")
|
|
// point2DFromPoint defVert0 = toPoint2D(vit->point());
|
|
|
|
// Triangulation::Edge_circulator ec = incident_edges(vit);
|
|
// Triangulation::Edge_circulator ecStart = ec;
|
|
|
|
// // Circulate around the edges to find the first which is not
|
|
// // infinite
|
|
// do
|
|
// {
|
|
// if (!is_infinite(ec)) break;
|
|
// } while (++ec != ecStart);
|
|
|
|
// // Store the start-end of the first non-infinte edge
|
|
// point2D de0 = toPoint2D(circumcenter(ec->first));
|
|
|
|
// // Keep track of the maximum edge length^2
|
|
// scalar maxEdgeLen2 = 0.0;
|
|
|
|
// // Keep track of the index of the longest edge
|
|
// label edgecd0i = -1;
|
|
|
|
// // Edge counter
|
|
// label edgei = 0;
|
|
|
|
// do
|
|
// {
|
|
// if (!is_infinite(ec))
|
|
// {
|
|
// // Get the end of the current edge
|
|
// point2D de1 = toPoint2D
|
|
// (
|
|
// circumcenter(ec->first->neighbor(ec->second))
|
|
// );
|
|
|
|
// // Store the current edge vector
|
|
// edges[edgei] = de1 - de0;
|
|
|
|
// // Store the edge mid-point in the vertices array
|
|
// vertices[edgei] = 0.5*(de1 + de0);
|
|
|
|
// // Move the current edge end into the edge start for the
|
|
// // next iteration
|
|
// de0 = de1;
|
|
|
|
// // Keep track of the longest edge
|
|
|
|
// scalar edgeLen2 = magSqr(edges[edgei]);
|
|
|
|
// if (edgeLen2 > maxEdgeLen2)
|
|
// {
|
|
// maxEdgeLen2 = edgeLen2;
|
|
// edgecd0i = edgei;
|
|
// }
|
|
|
|
// edgei++;
|
|
// }
|
|
// } while (++ec != ecStart);
|
|
|
|
// // Initialise cd0 such that the mesh will align
|
|
// // in in the x-y directions
|
|
// vector2D cd0(1, 0);
|
|
|
|
// if (meshControls().relaxOrientation())
|
|
// {
|
|
// // Get the longest edge from the array and use as the primary
|
|
// // direction of the coordinate system of the "square" cell
|
|
// cd0 = edges[edgecd0i];
|
|
// }
|
|
|
|
// if (meshControls().nearWallAlignedDist() > 0)
|
|
// {
|
|
// pointIndexHit pHit = qSurf_.tree().findNearest
|
|
// (
|
|
// toPoint3D(defVert0),
|
|
// meshControls().nearWallAlignedDist2()
|
|
// );
|
|
|
|
// if (pHit.hit())
|
|
// {
|
|
// cd0 = toPoint2D(faceNormals[pHit.index()]);
|
|
// }
|
|
// }
|
|
|
|
// // Rotate by 45deg needed to create an averaging procedure which
|
|
// // encourages the cells to be square
|
|
// cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x());
|
|
|
|
// // Normalise the primary coordinate direction
|
|
// cd0 /= mag(cd0);
|
|
|
|
// // Calculate the orthogonal coordinate direction
|
|
// vector2D cd1(-cd0.y(), cd0.x());
|
|
|
|
|
|
// // Restart the circulator
|
|
// ec = ecStart;
|
|
|
|
// // ... and the counter
|
|
// edgei = 0;
|
|
|
|
// // Initialise the displacement for the centre and sum-weights
|
|
// vector2D disp = vector2D::zero;
|
|
// scalar sumw = 0;
|
|
|
|
// do
|
|
// {
|
|
// if (!is_infinite(ec))
|
|
// {
|
|
// // Pick up the current edge
|
|
// const vector2D& ei = edges[edgei];
|
|
|
|
// // Calculate the centre to edge-centre vector
|
|
// vector2D deltai = vertices[edgei] - defVert0;
|
|
|
|
// // Set the weight for this edge contribution
|
|
// scalar w = 1;
|
|
|
|
// if (meshControls().squares())
|
|
// {
|
|
// w = magSqr(deltai.x()*ei.y() - deltai.y()*ei.x());
|
|
// // alternative weights
|
|
// //w = mag(deltai.x()*ei.y() - deltai.y()*ei.x());
|
|
// //w = magSqr(ei)*mag(deltai);
|
|
|
|
// // Use the following for an ~square mesh
|
|
// // Find the coordinate contributions for this edge delta
|
|
// scalar cd0deltai = cd0 & deltai;
|
|
// scalar cd1deltai = cd1 & deltai;
|
|
|
|
// // Create a "square" displacement
|
|
// if (mag(cd0deltai) > mag(cd1deltai))
|
|
// {
|
|
// disp += (w*cd0deltai)*cd0;
|
|
// }
|
|
// else
|
|
// {
|
|
// disp += (w*cd1deltai)*cd1;
|
|
// }
|
|
// }
|
|
// else
|
|
// {
|
|
// // Use this for a hexagon/pentagon mesh
|
|
// disp += w*deltai;
|
|
// }
|
|
|
|
// // Sum the weights
|
|
// sumw += w;
|
|
// }
|
|
// else
|
|
// {
|
|
// FatalErrorIn("CV2D::newPoints() const")
|
|
// << "Infinite triangle found in internal mesh"
|
|
// << exit(FatalError);
|
|
// }
|
|
|
|
// edgei++;
|
|
|
|
// } while (++ec != ecStart);
|
|
|
|
// // Calculate the average displacement
|
|
// disp /= sumw;
|
|
// totalDisp += disp;
|
|
// totalDist += mag(disp);
|
|
|
|
// // Move the point by a fraction of the average displacement
|
|
// movePoint(vit, defVert0 + relaxation*disp);
|
|
// }
|
|
// }
|
|
|
|
// Info << "\nTotal displacement = " << totalDisp
|
|
// << " total distance = " << totalDist << endl;
|
|
}
|
|
|
|
|
|
//void Foam::CV2D::moveInternalPoints(const point2DField& newPoints)
|
|
//{
|
|
// label pointI = 0;
|
|
|
|
// for
|
|
// (
|
|
// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
|
|
// vit != finite_vertices_end();
|
|
// ++vit
|
|
// )
|
|
// {
|
|
// if (vit->internalPoint())
|
|
// {
|
|
// movePoint(vit, newPoints[pointI++]);
|
|
// }
|
|
// }
|
|
//}
|
|
|
|
|
|
void Foam::CV2D::extractPatches
|
|
(
|
|
wordList& patchNames,
|
|
labelList& patchSizes,
|
|
EdgeMap<label>& mapEdgesRegion
|
|
) const
|
|
{
|
|
label nPatches = qSurf_.patchNames().size() + 1;
|
|
label defaultPatchIndex = qSurf_.patchNames().size();
|
|
|
|
patchNames.setSize(nPatches);
|
|
patchSizes.setSize(nPatches, 0);
|
|
mapEdgesRegion.clear();
|
|
|
|
const wordList& existingPatches = qSurf_.patchNames();
|
|
|
|
forAll(existingPatches, sP)
|
|
{
|
|
patchNames[sP] = existingPatches[sP];
|
|
}
|
|
|
|
patchNames[defaultPatchIndex] = "CV2D_default_patch";
|
|
|
|
for
|
|
(
|
|
Triangulation::Finite_edges_iterator eit = finite_edges_begin();
|
|
eit != finite_edges_end();
|
|
++eit
|
|
)
|
|
{
|
|
Face_handle fOwner = eit->first;
|
|
Face_handle fNeighbor = fOwner->neighbor(eit->second);
|
|
|
|
Vertex_handle vA = fOwner->vertex(cw(eit->second));
|
|
Vertex_handle vB = fOwner->vertex(ccw(eit->second));
|
|
|
|
if
|
|
(
|
|
(vA->internalOrBoundaryPoint() && !vB->internalOrBoundaryPoint())
|
|
|| (vB->internalOrBoundaryPoint() && !vA->internalOrBoundaryPoint())
|
|
)
|
|
{
|
|
point ptA = toPoint3D(vA->point());
|
|
point ptB = toPoint3D(vB->point());
|
|
|
|
label patchIndex = qSurf_.findPatch(ptA, ptB);
|
|
|
|
if (patchIndex == -1)
|
|
{
|
|
patchIndex = defaultPatchIndex;
|
|
|
|
WarningIn("Foam::CV2D::extractPatches")
|
|
<< "Dual face found that is not on a surface "
|
|
<< "patch. Adding to CV2D_default_patch."
|
|
<< endl;
|
|
}
|
|
|
|
edge e(fOwner->faceIndex(), fNeighbor->faceIndex());
|
|
patchSizes[patchIndex]++;
|
|
mapEdgesRegion.insert(e, patchIndex);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::CV2D::write() const
|
|
{
|
|
if (meshControls().objOutput())
|
|
{
|
|
writeFaces("allFaces.obj", false);
|
|
writeFaces("faces.obj", true);
|
|
writeTriangles("allTriangles.obj", false);
|
|
writeTriangles("triangles.obj", true);
|
|
writePatch("patch.pch");
|
|
}
|
|
}
|
|
|
|
|
|
void Foam::CV2D::write(const word& stage) const
|
|
{
|
|
if (meshControls().objOutput())
|
|
{
|
|
Foam::mkDir(stage + "Faces");
|
|
Foam::mkDir(stage + "Triangles");
|
|
|
|
writeFaces
|
|
(
|
|
stage
|
|
+ "Faces/allFaces_"
|
|
+ runTime_.timeName()
|
|
+ ".obj",
|
|
false
|
|
);
|
|
|
|
writeTriangles
|
|
(
|
|
stage
|
|
+ "Triangles/allTriangles_"
|
|
+ runTime_.timeName()
|
|
+ ".obj",
|
|
false
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|