From 188e43299dd174aaab1326a8d1c2ad0c947353e4 Mon Sep 17 00:00:00 2001
From: graham <graham.macpherson@strath.ac.uk>
Date: Thu, 2 Oct 2008 18:57:03 +0100
Subject: [PATCH] Basic isotropic cell relaxation calculated using the centres
of the dual faces. Feature control points stored after construction. Point
motion achieved by storing the location of the new points, clearing the whole
triangulation, reinserting feature points, then reinserting the new
locations. Much faster than removing than moving the points and avoids the
unexpectedly costly remove operation on the surface points.
---
.../mesh/generation/CV3DMesher/CV3D.C | 189 ++++++++----
.../mesh/generation/CV3DMesher/CV3D.H | 10 +
.../mesh/generation/CV3DMesher/CV3DMesher.C | 2 +-
.../indexedVertex_with_displacementSum.H | 289 ++++++++++++++++++
.../mesh/generation/CV3DMesher/calcDualMesh.C | 11 +-
.../generation/CV3DMesher/indexedVertex.H | 17 +-
.../CV3DMesher/insertFeaturePoints.C | 171 ++++++++---
7 files changed, 579 insertions(+), 110 deletions(-)
create mode 100644 applications/utilities/mesh/generation/CV3DMesher/backup/indexedVertex_with_displacementSum.H
diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3D.C b/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
index 267f88cf402..2263d46fce9 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
@@ -47,6 +47,25 @@ void Foam::CV3D::insertBoundingBox()
}
+void Foam::CV3D::reinsertPoints(const pointField& points)
+{
+ Info<< "Reinserting points after motion. ";
+
+ startOfInternalPoints_ = number_of_vertices();
+ label nVert = startOfInternalPoints_;
+
+ // Add the points and index them
+ forAll(points, i)
+ {
+ const point& p = points[i];
+
+ insert(toPoint(p))->index() = nVert++;
+ }
+
+ Info<< nVert << " vertices reinserted" << endl;
+}
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::CV3D::CV3D
@@ -60,7 +79,8 @@ Foam::CV3D::CV3D
controls_(controlDict),
tols_(controlDict, controls_.minCellSize, qSurf.bb()),
startOfInternalPoints_(0),
- startOfSurfacePointPairs_(0)
+ startOfSurfacePointPairs_(0),
+ featureConstrainingVertices_(0)
{
// insertBoundingBox();
insertFeaturePoints();
@@ -218,94 +238,143 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
{
Info<< "Calculating new points: " << endl;
- vector totalDisp = vector::zero;
- scalar totalDist = 0;
+ pointField dualVertices(number_of_cells());
+
+ pointField displacementAccumulator(startOfSurfacePointPairs_, vector::zero);
+
+ label dualVerti = 0;
+
+ // Find the dual point of each tetrahedron and assign it an index.
for
(
- Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
- vit != finite_vertices_end();
- ++vit
+ Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+ cit != finite_cells_end();
+ ++cit
)
{
- if (vit->internalPoint())
+ if
+ (
+ cit->vertex(0)->internalOrBoundaryPoint()
+ || cit->vertex(1)->internalOrBoundaryPoint()
+ || cit->vertex(2)->internalOrBoundaryPoint()
+ || cit->vertex(3)->internalOrBoundaryPoint()
+ )
{
- std::list<Facet> facets;
- incident_facets(vit, std::back_inserter(facets));
+ cit->cellIndex() = dualVerti;
+ dualVertices[dualVerti] = topoint(dual(cit));
+ dualVerti++;
+ }
+ else
+ {
+ cit->cellIndex() = -1;
+ }
+ }
- label maxIncidentFacets = 120;
- List<point> vertices(maxIncidentFacets);
- List<vector> edges(maxIncidentFacets);
+ dualVertices.setSize(dualVerti);
- point vd(topoint(vit->point()));
+ // loop around the Delaunay edges to construct the dual faces.
+ // Find the face-centre and use it to calculate the displacement vector
+ // contribution to the Delaunay vertices (Dv) attached to the edge. Add the
+ // contribution to the running displacement vector of each Dv.
- point vi0 = topoint(dual(facets.begin()->first));
+ for
+ (
+ Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+ eit != finite_edges_end();
+ ++eit
+ )
+ {
+ if
+ (
+ eit->first->vertex(eit->second)->internalOrBoundaryPoint()
+ && eit->first->vertex(eit->third)->internalOrBoundaryPoint()
+ )
+ {
+ Cell_circulator ccStart = incident_cells(*eit);
+ Cell_circulator cc = ccStart;
- label edgei = 0;
+ DynamicList<label> verticesOnFace;
- for
- (
- std::list<Facet>::iterator fit=facets.begin();
- fit != facets.end();
- ++fit
- )
+ do
{
- if
- (
- is_infinite(fit->first)
- || is_infinite(fit->first->neighbor(fit->second))
- )
+ if (!is_infinite(cc))
{
- FatalErrorIn("relaxPoints")
- << "Finite cell attached to facet incident on vertex"
- << exit(FatalError);
+ if (cc->cellIndex() < 0)
+ {
+ FatalErrorIn("Foam::CV3D::relaxPoints")
+ << "Dual face uses circumcenter defined by a "
+ << " Delaunay tetrahedron with no internal "
+ << "or boundary points."
+ << exit(FatalError);
+ }
+
+ verticesOnFace.append(cc->cellIndex());
}
+ } while (++cc != ccStart);
- point vi1 = topoint(dual(fit->first->neighbor(fit->second)));
+ verticesOnFace.shrink();
- edges[edgei] = vi1 - vi0;
+ face dualFace(verticesOnFace);
- vertices[edgei] = 0.5*(vi1 + vi0);
+ Cell_handle c = eit->first;
+ Vertex_handle vA = c->vertex(eit->second);
+ Vertex_handle vB = c->vertex(eit->third);
- vi0 = vi1;
+ point dVA = topoint(vA->point());
+ point dVB = topoint(vB->point());
- edgei++;
- }
-
- edgei = 0;
+ point dualFaceCentre(dualFace.centre(dualVertices));
- // Initialise the displacement for the centre and sum-weights
- vector disp = vector::zero;
- scalar sumw = 0;
+ scalar weight = 1.0;
- for
- (
- std::list<Facet>::iterator fit=facets.begin();
- fit != facets.end();
- ++fit
- )
+ if (vA->internalPoint())
+ {
+ //displacementAccumulator[vA->index()] = vA->index()*vector::one;
+ displacementAccumulator[vA->index()] +=
+ relaxation*weight*(dualFaceCentre - dVA);
+ }
+ if (vB->internalPoint())
{
- vector deltai = vertices[edgei] - vd;
+ //displacementAccumulator[vB->index()] = vB->index()*vector::one;
+ displacementAccumulator[vB->index()] +=
+ relaxation*weight*(dualFaceCentre - dVB);
+ }
+ }
+ }
- scalar w = 1;
+ vector totalDisp = sum(displacementAccumulator);
+ scalar totalDist = sum(mag(displacementAccumulator));
- disp += w*deltai;
+ Info<< "Total displacement = " << totalDisp
+ << " total distance = " << totalDist << endl;
- sumw += w;
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->internalPoint())
+ {
+ displacementAccumulator[vit->index()] += topoint(vit->point());
+ }
+ }
- edgei++;
- }
+ pointField internalDelaunayVertices = SubField<point>
+ (
+ displacementAccumulator,
+ displacementAccumulator.size() - startOfInternalPoints_,
+ startOfInternalPoints_
+ );
- disp /= sumw;
- totalDisp += disp;
- totalDist += mag(disp);
+ // Remove the entire triangulation
+ this->clear();
- movePoint(vit, vd + relaxation*disp);
- }
- }
+ reinsertFeaturePoints();
- Info<< "Total displacement = " << totalDisp
- << " total distance = " << totalDist << endl;
+ reinsertPoints(internalDelaunayVertices);
}
diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3D.H b/applications/utilities/mesh/generation/CV3DMesher/CV3D.H
index 4c908476ab0..041e1e37145 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/CV3D.H
+++ b/applications/utilities/mesh/generation/CV3DMesher/CV3D.H
@@ -209,6 +209,10 @@ private:
// removing and insertin the surface point-pairs
label startOfBoundaryConformPointPairs_;
+ //- Store the feature constraining points to be reinserted after a
+ //- triangulation clear
+ List<Vb> featureConstrainingVertices_;
+
// Private Member Functions
@@ -248,6 +252,12 @@ private:
//- Insert point groups at the feature points.
void insertFeaturePoints();
+ //- Reinsert stored feature point defining points.
+ void reinsertFeaturePoints();
+
+ //- Reinsert points that have been saved to clear the mesh.
+ void reinsertPoints (const pointField& points);
+
//- Insert point-pairs at the given set of points using the surface
// normals corresponding to the given set of surface triangles
// and write the inserted point locations to the given file.
diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C b/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C
index 4104d280cfc..cefaa6104e4 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C
@@ -110,7 +110,7 @@ int main(int argc, char *argv[])
mesh.relaxPoints(relaxation);
- mesh.removeSurfacePointPairs();
+ //mesh.removeSurfacePointPairs();
mesh.insertSurfacePointPairs();
mesh.boundaryConform();
diff --git a/applications/utilities/mesh/generation/CV3DMesher/backup/indexedVertex_with_displacementSum.H b/applications/utilities/mesh/generation/CV3DMesher/backup/indexedVertex_with_displacementSum.H
new file mode 100644
index 00000000000..b495bc770eb
--- /dev/null
+++ b/applications/utilities/mesh/generation/CV3DMesher/backup/indexedVertex_with_displacementSum.H
@@ -0,0 +1,289 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
+ \\/ 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 2 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, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+ indexedVertex
+
+Description
+ An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep
+ track of the vertices in the triangulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef indexedVertex_H
+#define indexedVertex_H
+
+#include <CGAL/Triangulation_3.h>
+#include "vector.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace CGAL
+{
+
+/*---------------------------------------------------------------------------*\
+ Class indexedVertex Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Gt, class Vb=CGAL::Triangulation_vertex_base_3<Gt> >
+class indexedVertex
+:
+ public Vb
+{
+ // Private data
+
+ //- The index for this triangle vertex
+ int index_;
+
+ //- Index of pair-point :
+ // NEAR_BOUNDARY_POINT : internal near boundary point.
+ // INTERNAL_POINT : internal point.
+ // FAR_POINT : far-point.
+ // >= 0 : part of point-pair. Index of other point.
+ // Lowest numbered is inside one (master).
+ int type_;
+
+ Foam::vector displacementSum_;
+
+
+public:
+
+ enum pointTypes
+ {
+ NEAR_BOUNDARY_POINT = -4,
+ INTERNAL_POINT = -3,
+ MIRROR_POINT = -2,
+ FAR_POINT = -1
+ };
+
+ typedef typename Vb::Vertex_handle Vertex_handle;
+ typedef typename Vb::Cell_handle Cell_handle;
+ typedef typename Vb::Point Point;
+
+ template<typename TDS2>
+ struct Rebind_TDS
+ {
+ typedef typename Vb::template Rebind_TDS<TDS2>::Other Vb2;
+ typedef indexedVertex<Gt,Vb2> Other;
+ };
+
+ indexedVertex()
+ :
+ Vb(),
+ index_(INTERNAL_POINT),
+ type_(INTERNAL_POINT),
+ displacementSum_(Foam::vector::zero)
+ {}
+
+ indexedVertex(const Point& p)
+ :
+ Vb(p),
+ index_(INTERNAL_POINT),
+ type_(INTERNAL_POINT),
+ displacementSum_(Foam::vector::zero)
+ {}
+
+ indexedVertex(const Point& p, Cell_handle f)
+ :
+ Vb(f, p),
+ index_(INTERNAL_POINT),
+ type_(INTERNAL_POINT),
+ displacementSum_(Foam::vector::zero)
+ {}
+
+ indexedVertex(Cell_handle f)
+ :
+ Vb(f),
+ index_(INTERNAL_POINT),
+ type_(INTERNAL_POINT),
+ displacementSum_(Foam::vector::zero)
+ {}
+
+
+ int& index()
+ {
+ return index_;
+ }
+
+ int index() const
+ {
+ return index_;
+ }
+
+
+ int& type()
+ {
+ return type_;
+ }
+
+ int type() const
+ {
+ return type_;
+ }
+
+
+ Foam::vector& displacementSum()
+ {
+ return displacementSum_;
+ }
+
+ const Foam::vector& displacementSum() const
+ {
+ return displacementSum_;
+ }
+
+ //- Is point a far-point
+ inline bool farPoint() const
+ {
+ return type_ == FAR_POINT;
+ }
+
+ //- Is point internal, i.e. not on boundary
+ inline bool internalPoint() const
+ {
+ return type_ <= INTERNAL_POINT;
+ }
+
+ //- Is point internal and near the boundary
+ inline bool nearBoundary() const
+ {
+ return type_ == NEAR_BOUNDARY_POINT;
+ }
+
+ //- Set the point to be near the boundary
+ inline void setNearBoundary()
+ {
+ type_ = NEAR_BOUNDARY_POINT;
+ }
+
+ //- Is point a mirror point
+ inline bool mirrorPoint() const
+ {
+ return type_ == MIRROR_POINT;
+ }
+
+ //- Either master or slave of pointPair.
+ inline bool pairPoint() const
+ {
+ return type_ >= 0;
+ }
+
+ //- Master of a pointPair is the lowest numbered one.
+ inline bool ppMaster() const
+ {
+ if (type_ > index_)
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ }
+
+ //- Slave of a pointPair is the highest numbered one.
+ inline bool ppSlave() const
+ {
+ if (type_ >= 0 && type_ < index_)
+ {
+ return true;
+ }
+ else
+ {
+ return false;
+ }
+ }
+
+ //- Either original internal point or master of pointPair.
+ inline bool internalOrBoundaryPoint() const
+ {
+ return internalPoint() || ppMaster();
+ }
+
+ //- Is point near the boundary or part of the boundary definition
+ inline bool nearOrOnBoundary() const
+ {
+ return pairPoint() || mirrorPoint() || nearBoundary();
+ }
+
+ //- Do the two given vertices consitute a boundary point-pair
+ inline friend bool pointPair
+ (
+ const indexedVertex& v0,
+ const indexedVertex& v1
+ )
+ {
+ return v0.index_ == v1.type_ || v1.index_ == v0.type_;
+ }
+
+ //- Do the three given vertices consitute a boundary triangle
+ inline friend bool boundaryTriangle
+ (
+ const indexedVertex& v0,
+ const indexedVertex& v1,
+ const indexedVertex& v2
+ )
+ {
+ return (v0.pairPoint() && pointPair(v1, v2))
+ || (v1.pairPoint() && pointPair(v2, v0))
+ || (v2.pairPoint() && pointPair(v0, v1));
+ }
+
+ //- Do the three given vertices consitute an outside triangle
+ inline friend bool outsideTriangle
+ (
+ const indexedVertex& v0,
+ const indexedVertex& v1,
+ const indexedVertex& v2
+ )
+ {
+ return (v0.farPoint() || v0.ppSlave())
+ || (v1.farPoint() || v1.ppSlave())
+ || (v2.farPoint() || v2.ppSlave());
+ }
+
+
+ // inline void operator=(const Triangulation::Finite_vertices_iterator vit)
+ // {
+ // Vb::operator=indexedVertex(vit->point());
+
+ // this->index_ = vit->index();
+
+ // this->type_ = vit->type();
+
+ // this->displacementSum_ = vit->displacementSum();
+ // }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace CGAL
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
index a58f715d760..0f694eddaa7 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C
@@ -149,15 +149,14 @@ void Foam::CV3D::calcDualMesh
{
if (cc->cellIndex() < 0)
{
- FatalErrorIn
- (
- "Foam::CV3D::calcDualMesh"
- )<< "Dual face uses circumcenter defined by a Delaunay"
- " tetrahedron with no internal or boundary points."
+ FatalErrorIn("Foam::CV3D::calcDualMesh")
+ << "Dual face uses circumcenter defined by a "
+ << " Delaunay tetrahedron with no internal "
+ << "or boundary points."
<< exit(FatalError);
}
- verticesOnFace.append(cc->cellIndex());
+ verticesOnFace.append(cc->cellIndex());
}
} while (++cc != ccStart);
diff --git a/applications/utilities/mesh/generation/CV3DMesher/indexedVertex.H b/applications/utilities/mesh/generation/CV3DMesher/indexedVertex.H
index eb680b312d4..8df1a62b3f5 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/indexedVertex.H
+++ b/applications/utilities/mesh/generation/CV3DMesher/indexedVertex.H
@@ -46,7 +46,7 @@ namespace CGAL
\*---------------------------------------------------------------------------*/
template<class Gt, class Vb=CGAL::Triangulation_vertex_base_3<Gt> >
-class indexedVertex
+class indexedVertex
:
public Vb
{
@@ -81,11 +81,10 @@ public:
template<typename TDS2>
struct Rebind_TDS
{
- typedef typename Vb::template Rebind_TDS<TDS2>::Other Vb2;
+ typedef typename Vb::template Rebind_TDS<TDS2>::Other Vb2;
typedef indexedVertex<Gt,Vb2> Other;
};
-
indexedVertex()
:
Vb(),
@@ -246,6 +245,18 @@ public:
|| (v1.farPoint() || v1.ppSlave())
|| (v2.farPoint() || v2.ppSlave());
}
+
+
+ // inline void operator=(const Triangulation::Finite_vertices_iterator vit)
+ // {
+ // Vb::operator=indexedVertex(vit->point());
+
+ // this->index_ = vit->index();
+
+ // this->type_ = vit->type();
+
+ // this->displacementSum_ = vit->displacementSum();
+ // }
};
diff --git a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
index 767a9beb955..68eae2e9e2e 100644
--- a/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
+++ b/applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C
@@ -40,7 +40,8 @@ void Foam::CV3D::insertFeaturePoints()
scalar planeErrorAngle = 0.1*(180.0 - controls_.includedAngle);
- scalar planeErrorAngleCos = cos(mathematicalConstant::pi*planeErrorAngle/180.0);
+ scalar planeErrorAngleCos =
+ cos(mathematicalConstant::pi*planeErrorAngle/180.0);
forAll(featPoints, i)
{
@@ -169,7 +170,8 @@ void Foam::CV3D::insertFeaturePoints()
cornerNormal /= mag(cornerNormal);
point internalPt = featPt - tols_.ppDist*cornerNormal;
- label internalPtIndex = insertPoint(internalPt, number_of_vertices() + 1);
+ label internalPtIndex =
+ insertPoint(internalPt, number_of_vertices() + 1);
forAll (uniquePlaneNormals, uPN)
{
@@ -177,7 +179,8 @@ void Foam::CV3D::insertFeaturePoints()
plane planeN = plane(featPt, n);
- point externalPt = internalPt + 2.0 * planeN.distance(internalPt) * n;
+ point externalPt =
+ internalPt + 2.0 * planeN.distance(internalPt) * n;
insertPoint(externalPt, internalPtIndex);
}
@@ -192,6 +195,7 @@ void Foam::CV3D::insertFeaturePoints()
cornerNormal /= mag(cornerNormal);
point externalPt = featPt + tols_.ppDist*cornerNormal;
+
label externalPtIndex = number_of_vertices() + concaveEdges.size();
label internalPtIndex = -1;
@@ -202,7 +206,8 @@ void Foam::CV3D::insertFeaturePoints()
plane planeN = plane(featPt, n);
- point internalPt = externalPt - 2.0 * planeN.distance(externalPt) * n;
+ point internalPt = externalPt
+ - 2.0 * planeN.distance(externalPt) * n;
internalPtIndex = insertPoint(internalPt, externalPtIndex);
}
@@ -216,7 +221,8 @@ void Foam::CV3D::insertFeaturePoints()
if (convexEdges.size() + concaveEdges.size() > 3)
{
- Info<< concaveEdges.size() + convexEdges.size() << " mixed edge feature."
+ Info<< concaveEdges.size() + convexEdges.size()
+ << " mixed edge feature."
<< " NOT IMPLEMENTED." << endl;
}
else if (convexEdges.size() > concaveEdges.size())
@@ -231,6 +237,7 @@ void Foam::CV3D::insertFeaturePoints()
const labelList& eFaces = qSurf_.edgeFaces()[concaveEdgeI];
label faceA = eFaces[0];
+
vector nA = qSurf_.faceNormals()[faceA];
scalar maxNormalDotProduct = -SMALL;
@@ -265,9 +272,10 @@ void Foam::CV3D::insertFeaturePoints()
}
const vector& concaveEdgePlaneANormal =
- uniquePlaneNormals[concaveEdgePlanes[0]];
+ uniquePlaneNormals[concaveEdgePlanes[0]];
+
const vector& concaveEdgePlaneBNormal =
- uniquePlaneNormals[concaveEdgePlanes[1]];
+ uniquePlaneNormals[concaveEdgePlanes[1]];
label convexEdgesPlaneI;
@@ -299,8 +307,9 @@ void Foam::CV3D::insertFeaturePoints()
edgeDirection *= -1.0;
}
- // Intersect planes parallel to the concave edge planes offset by ppDist
- // and the plane defined by featPt and the edge vector.
+ // Intersect planes parallel to the concave edge planes offset
+ // by ppDist and the plane defined by featPt and the edge
+ // vector.
plane planeA
(
featPt + tols_.ppDist*concaveEdgePlaneANormal,
@@ -317,35 +326,48 @@ void Foam::CV3D::insertFeaturePoints()
+ tols_.ppDist*edgeDirection
* concaveEdge.vec(localPts)/concaveEdge.mag(localPts);
- // Finding the nearest point on the intersecting line to the edge point.
- // Floating point errors often encountered using planePlaneIntersect
+ // Finding the nearest point on the intersecting line to the
+ // edge point. Floating point errors often encountered using
+ // planePlaneIntersect
plane planeF(concaveEdgeLocalFeatPt, concaveEdge.vec(localPts));
- point concaveEdgeExternalPt = planeF.planePlaneIntersect(planeA,planeB);
+
+ point concaveEdgeExternalPt =
+ planeF.planePlaneIntersect(planeA,planeB);
label concaveEdgeExternalPtI = number_of_vertices() + 4;
- // Redefine planes to be on the feature surfaces to project through
+ // Redefine planes to be on the feature surfaces to project
+ // through
planeA = plane(featPt, concaveEdgePlaneANormal);
+
planeB = plane(featPt, concaveEdgePlaneBNormal);
- point internalPtA = concaveEdgeExternalPt -
- 2*planeA.distance(concaveEdgeExternalPt) * concaveEdgePlaneANormal;
- label internalPtAI = insertPoint(internalPtA, concaveEdgeExternalPtI);
+ point internalPtA = concaveEdgeExternalPt
+ - 2*planeA.distance(concaveEdgeExternalPt)
+ *concaveEdgePlaneANormal;
- point internalPtB = concaveEdgeExternalPt -
- 2*planeB.distance(concaveEdgeExternalPt) * concaveEdgePlaneBNormal;
- label internalPtBI = insertPoint(internalPtB, concaveEdgeExternalPtI);
+ label internalPtAI =
+ insertPoint(internalPtA, concaveEdgeExternalPtI);
+
+ point internalPtB = concaveEdgeExternalPt
+ - 2*planeB.distance(concaveEdgeExternalPt)
+ *concaveEdgePlaneBNormal;
+
+ label internalPtBI =
+ insertPoint(internalPtB, concaveEdgeExternalPtI);
plane planeC(featPt, convexEdgesPlaneNormal);
point externalPtD = internalPtA +
- 2*planeC.distance(internalPtA) * convexEdgesPlaneNormal;
+ 2*planeC.distance(internalPtA) * convexEdgesPlaneNormal;
+
insertPoint(externalPtD, internalPtAI);
point externalPtE = internalPtB +
- 2*planeC.distance(internalPtB) * convexEdgesPlaneNormal;
+ 2*planeC.distance(internalPtB) * convexEdgesPlaneNormal;
+
insertPoint(externalPtE, internalPtBI);
insertPoint(concaveEdgeExternalPt, internalPtAI);
@@ -361,16 +383,24 @@ void Foam::CV3D::insertFeaturePoints()
// Add additional mitering points
vector concaveEdgeNormal =
- edgeDirection*concaveEdge.vec(localPts)/concaveEdge.mag(localPts);
+ edgeDirection*concaveEdge.vec(localPts)
+ /concaveEdge.mag(localPts);
scalar angleSign = 1.0;
- if (qSurf_.outside(featPt - convexEdgesPlaneNormal*tols_.ppDist))
+ if
+ (
+ qSurf_.outside
+ (
+ featPt - convexEdgesPlaneNormal*tols_.ppDist
+ )
+ )
{
angleSign = -1.0;
}
- scalar phi = angleSign*acos(concaveEdgeNormal & -convexEdgesPlaneNormal);
+ scalar phi = angleSign
+ *acos(concaveEdgeNormal & -convexEdgesPlaneNormal);
scalar guard =
(
@@ -380,14 +410,15 @@ void Foam::CV3D::insertFeaturePoints()
)
)/cos(phi) - 1;
- point internalPtF = concaveEdgeExternalPt + (2 + guard) *
- (concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
+ point internalPtF = concaveEdgeExternalPt + (2 + guard)
+ *(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
- label internalPtFI = insertPoint(internalPtF, number_of_vertices() + 1);
+ label internalPtFI =
+ insertPoint(internalPtF, number_of_vertices() + 1);
point externalPtG = internalPtF +
- 2*planeC.distance(internalPtF) * convexEdgesPlaneNormal;
- insertPoint(externalPtG, internalPtFI);
+ 2*planeC.distance(internalPtF) * convexEdgesPlaneNormal;
+ insertPoint(externalPtG, internalPtFI);
}
}
else
@@ -470,8 +501,8 @@ void Foam::CV3D::insertFeaturePoints()
edgeDirection *= -1.0;
}
- // Intersect planes parallel to the convex edge planes offset by ppDist
- // and the plane defined by featPt and the edge vector.
+ // Intersect planes parallel to the convex edge planes offset by
+ // ppDist and the plane defined by featPt and the edge vector.
plane planeA
(
featPt - tols_.ppDist*convexEdgePlaneANormal,
@@ -489,21 +520,28 @@ void Foam::CV3D::insertFeaturePoints()
* convexEdge.vec(localPts)/convexEdge.mag(localPts);
- // Finding the nearest point on the intersecting line to the edge point.
- // Floating point errors often encountered using planePlaneIntersect
+ // Finding the nearest point on the intersecting line to the
+ // edge point. Floating point errors often encountered using
+ // planePlaneIntersect
plane planeF(convexEdgeLocalFeatPt, convexEdge.vec(localPts));
- point convexEdgeInternalPt = planeF.planePlaneIntersect(planeA,planeB);
+
+ point convexEdgeInternalPt =
+ planeF.planePlaneIntersect(planeA,planeB);
planeA = plane(featPt, convexEdgePlaneANormal);
+
planeB = plane(featPt, convexEdgePlaneBNormal);
- point externalPtA = convexEdgeInternalPt +
- 2*planeA.distance(convexEdgeInternalPt) * convexEdgePlaneANormal;
+ point externalPtA = convexEdgeInternalPt
+ + 2*planeA.distance(convexEdgeInternalPt)
+ * convexEdgePlaneANormal;
label externalPtAI = number_of_vertices() + 3;
- point externalPtB = convexEdgeInternalPt +
- 2*planeB.distance(convexEdgeInternalPt) * convexEdgePlaneBNormal;
+ point externalPtB = convexEdgeInternalPt
+ + 2*planeB.distance(convexEdgeInternalPt)
+ * convexEdgePlaneBNormal;
+
label externalPtBI = number_of_vertices() + 4;
label convexEdgeInternalPtI = insertPoint
@@ -522,13 +560,35 @@ void Foam::CV3D::insertFeaturePoints()
2*planeC.distance(externalPtB) * concaveEdgesPlaneNormal;
insertPoint(internalPtE, externalPtBI);
- insertPoint(externalPtA,convexEdgeInternalPtI);
+ insertPoint(externalPtA, convexEdgeInternalPtI);
- insertPoint(externalPtB,convexEdgeInternalPtI);
+ insertPoint(externalPtB, convexEdgeInternalPtI);
}
}
}
+ featureConstrainingVertices_.setSize(number_of_vertices());
+
+ label featPtI = 0;
+
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ vit++
+ )
+ {
+ // featureConstrainingVertices_[featPtI] = vit;
+
+ featureConstrainingVertices_[featPtI] = Vb(vit->point());
+
+ featureConstrainingVertices_[featPtI].index() = vit->index();
+
+ featureConstrainingVertices_[featPtI].type() = vit->type();
+
+ featPtI++;
+ }
+
if (controls_.writeFeatureTriangulation)
{
writePoints("feat_allPoints.obj", false);
@@ -538,4 +598,35 @@ void Foam::CV3D::insertFeaturePoints()
}
+void Foam::CV3D::reinsertFeaturePoints()
+{
+ if (featureConstrainingVertices_.size())
+ {
+
+ forAll(featureConstrainingVertices_, f)
+ {
+ const Point& fPt(featureConstrainingVertices_[f].point());
+
+ uint nVert = number_of_vertices();
+
+ Vertex_handle vh = insert(fPt);
+
+ if (nVert == number_of_vertices())
+ {
+ FatalErrorIn("Foam::CV3D::reinsertFeaturePoints")
+ << "Failed to reinsert feature point " << topoint(fPt)
+ << endl;
+ }
+
+ vh->index() = featureConstrainingVertices_[f].index();
+ vh->type() = featureConstrainingVertices_[f].type();
+ }
+ }
+ else
+ {
+ WarningIn("Foam::CV3D::reinsertFeaturePoints")
+ << "No stored feature points to reinsert." << endl;
+ }
+}
+
// ************************************************************************* //
--
GitLab