diff --git a/applications/utilities/mesh/generation/CV3DMesher/CV3D.C b/applications/utilities/mesh/generation/CV3DMesher/CV3D.C
index 267f88cf402bc9fb1066a0ae91a650addfbe9bd4..2263d46fce91b87e6d3fce59c9df05d7666c129e 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 4c908476ab0c62156e3f466a195674f7061916ca..041e1e37145e40f4fb06fb637d1489866964bbf0 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 4104d280cfc8aaf3c281a2fdf9370be8e0c0c7d3..cefaa6104e4247e0d29c7ea178496afaede7a196 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 0000000000000000000000000000000000000000..b495bc770eb1b4f18ee81d0c959d4a06c304b6d4
--- /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 a58f715d76086bb9df93cb0ce67e81a2f84c90dc..0f694eddaa76519f5ac6a922ab58cf83ba358adc 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 eb680b312d4630028fcb57b92637fb41ba2638d9..8df1a62b3f5a1ec94bd33b9a4d140a146ffe1d4f 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 767a9beb9553f24994928232fd0185ab41e88978..68eae2e9e2e17fe6e6284711b0816b7a1094759d 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;
+ }
+}
+
// ************************************************************************* //