mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
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.
This commit is contained in:
graham
@ -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,8 +238,116 @@ 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_cells_iterator cit = finite_cells_begin();
|
||||
cit != finite_cells_end();
|
||||
++cit
|
||||
)
|
||||
{
|
||||
if
|
||||
(
|
||||
cit->vertex(0)->internalOrBoundaryPoint()
|
||||
|| cit->vertex(1)->internalOrBoundaryPoint()
|
||||
|| cit->vertex(2)->internalOrBoundaryPoint()
|
||||
|| cit->vertex(3)->internalOrBoundaryPoint()
|
||||
)
|
||||
{
|
||||
cit->cellIndex() = dualVerti;
|
||||
dualVertices[dualVerti] = topoint(dual(cit));
|
||||
dualVerti++;
|
||||
}
|
||||
else
|
||||
{
|
||||
cit->cellIndex() = -1;
|
||||
}
|
||||
}
|
||||
|
||||
dualVertices.setSize(dualVerti);
|
||||
|
||||
// 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.
|
||||
|
||||
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;
|
||||
|
||||
DynamicList<label> verticesOnFace;
|
||||
|
||||
do
|
||||
{
|
||||
if (!is_infinite(cc))
|
||||
{
|
||||
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);
|
||||
|
||||
verticesOnFace.shrink();
|
||||
|
||||
face dualFace(verticesOnFace);
|
||||
|
||||
Cell_handle c = eit->first;
|
||||
Vertex_handle vA = c->vertex(eit->second);
|
||||
Vertex_handle vB = c->vertex(eit->third);
|
||||
|
||||
point dVA = topoint(vA->point());
|
||||
point dVB = topoint(vB->point());
|
||||
|
||||
point dualFaceCentre(dualFace.centre(dualVertices));
|
||||
|
||||
scalar weight = 1.0;
|
||||
|
||||
if (vA->internalPoint())
|
||||
{
|
||||
//displacementAccumulator[vA->index()] = vA->index()*vector::one;
|
||||
displacementAccumulator[vA->index()] +=
|
||||
relaxation*weight*(dualFaceCentre - dVA);
|
||||
}
|
||||
if (vB->internalPoint())
|
||||
{
|
||||
//displacementAccumulator[vB->index()] = vB->index()*vector::one;
|
||||
displacementAccumulator[vB->index()] +=
|
||||
relaxation*weight*(dualFaceCentre - dVB);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
vector totalDisp = sum(displacementAccumulator);
|
||||
scalar totalDist = sum(mag(displacementAccumulator));
|
||||
|
||||
Info<< "Total displacement = " << totalDisp
|
||||
<< " total distance = " << totalDist << endl;
|
||||
|
||||
for
|
||||
(
|
||||
@ -230,82 +358,23 @@ void Foam::CV3D::relaxPoints(const scalar relaxation)
|
||||
{
|
||||
if (vit->internalPoint())
|
||||
{
|
||||
std::list<Facet> facets;
|
||||
incident_facets(vit, std::back_inserter(facets));
|
||||
|
||||
label maxIncidentFacets = 120;
|
||||
List<point> vertices(maxIncidentFacets);
|
||||
List<vector> edges(maxIncidentFacets);
|
||||
|
||||
point vd(topoint(vit->point()));
|
||||
|
||||
point vi0 = topoint(dual(facets.begin()->first));
|
||||
|
||||
label edgei = 0;
|
||||
|
||||
for
|
||||
(
|
||||
std::list<Facet>::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<Facet>::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);
|
||||
displacementAccumulator[vit->index()] += topoint(vit->point());
|
||||
}
|
||||
}
|
||||
|
||||
Info<< "Total displacement = " << totalDisp
|
||||
<< " total distance = " << totalDist << endl;
|
||||
pointField internalDelaunayVertices = SubField<point>
|
||||
(
|
||||
displacementAccumulator,
|
||||
displacementAccumulator.size() - startOfInternalPoints_,
|
||||
startOfInternalPoints_
|
||||
);
|
||||
|
||||
// Remove the entire triangulation
|
||||
this->clear();
|
||||
|
||||
reinsertFeaturePoints();
|
||||
|
||||
reinsertPoints(internalDelaunayVertices);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -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.
|
||||
|
||||
@ -110,7 +110,7 @@ int main(int argc, char *argv[])
|
||||
|
||||
mesh.relaxPoints(relaxation);
|
||||
|
||||
mesh.removeSurfacePointPairs();
|
||||
//mesh.removeSurfacePointPairs();
|
||||
mesh.insertSurfacePointPairs();
|
||||
mesh.boundaryConform();
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -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);
|
||||
|
||||
|
||||
@ -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();
|
||||
// }
|
||||
};
|
||||
|
||||
|
||||
|
||||
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
Reference in New Issue
Block a user