diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index ac61a70f90e21fbc6a84b891b3d0ee0f95618406..9df248939aa792da5c1f884d004ae5f1c875d97a 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -1760,7 +1760,16 @@ void Foam::conformalVoronoiMesh::calcDualMesh
labelList& patchStarts
)
{
- // ~~~~~~~~~~~ removing short edges by indexing dual vertices ~~~~~~~~~~~~~~
+ timeCheck();
+
+ setVertexSizeAndAlignment();
+
+ timeCheck();
+
+ // Indexing cells, which are Dual vertices
+ label dualVertI = 0;
+
+ points.setSize(number_of_cells());
for
(
@@ -1769,211 +1778,240 @@ void Foam::conformalVoronoiMesh::calcDualMesh
++cit
)
{
- cit->cellIndex() = -1;
- }
-
- points.setSize(number_of_cells());
-
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- // Looking up minEdgeLenSqr with a dummy point, in future will be available
- // as a local value to be looked up in-place.
-
- scalar minEdgeLenSqr = sqr(minimumEdgeLength(vector::zero));
-
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- label dualVertI = 0;
-
- // Scanning by number of short (dual) edges (nSE) attached to the
- // circumcentre of each Delaunay tet. A Delaunay tet may only have four
- // dual edges emanating from its circumcentre, assigning positions and
- // indices to those with 4 short edges attached first, then >= 3, then >= 2
- // etc.
- for (label nSE = 4; nSE >= 0; nSE--)
- {
- Info<< nl << "Scanning for dual vertices with >= "
- << nSE
- << " short edges attached." << endl;
-
- for
+ if
(
- Triangulation::Finite_cells_iterator cit = finite_cells_begin();
- cit != finite_cells_end();
- ++cit
+ cit->vertex(0)->internalOrBoundaryPoint()
+ || cit->vertex(1)->internalOrBoundaryPoint()
+ || cit->vertex(2)->internalOrBoundaryPoint()
+ || cit->vertex(3)->internalOrBoundaryPoint()
)
{
- // If the Delaunay tet has an index already then it has either
- // evaluated itself and taken action or has had its index dictated
- // by a neighbouring tet with more short edges attached.
-
- if (cit->cellIndex() == -1)
- {
- point dualVertex = topoint(dual(cit));
-
- label shortEdges = 0;
-
- List<bool> edgeIsShort(4, false);
-
- List<bool> neighbourAlreadyIndexed(4, false);
-
- // Loop over the four facets of the Delaunay tet
- for (label f = 0; f < 4; f++)
- {
- // Check that at least one of the vertices of the facet is
- // an internal or boundary point
- if
- (
- cit->vertex(vertex_triple_index(f, 0))->
- internalOrBoundaryPoint()
- || cit->vertex(vertex_triple_index(f, 1))->
- internalOrBoundaryPoint()
- || cit->vertex(vertex_triple_index(f, 2))->
- internalOrBoundaryPoint()
- )
- {
- point neighDualVertex;
-
- label cNI = cit->neighbor(f)->cellIndex();
-
- if (cNI == -1)
- {
- neighDualVertex = topoint(dual(cit->neighbor(f)));
- }
- else
- {
- neighDualVertex = points[cNI];
- }
-
- if
- (
- magSqr(dualVertex - neighDualVertex) < minEdgeLenSqr
- )
- {
- edgeIsShort[f] = true;
-
- if (cNI > -1)
- {
- neighbourAlreadyIndexed[f] = true;
- }
-
- shortEdges++;
- }
- }
- }
-
- if (nSE == 0 && shortEdges == 0)
- {
- // Final iteration and no short edges are found, index
- // remaining dual vertices.
-
- if
- (
- cit->vertex(0)->internalOrBoundaryPoint()
- || cit->vertex(1)->internalOrBoundaryPoint()
- || cit->vertex(2)->internalOrBoundaryPoint()
- || cit->vertex(3)->internalOrBoundaryPoint()
- )
- {
- cit->cellIndex() = dualVertI;
- points[dualVertI] = dualVertex;
- dualVertI++;
- }
- }
- else if
- (
- shortEdges >= nSE
- )
- {
- // Info<< neighbourAlreadyIndexed << ' '
- // << edgeIsShort << endl;
-
- label numUnindexedNeighbours = 1;
-
- for (label f = 0; f < 4; f++)
- {
- if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
- {
- dualVertex += topoint(dual(cit->neighbor(f)));
-
- numUnindexedNeighbours++;
- }
- }
-
- dualVertex /= numUnindexedNeighbours;
-
- label nearestExistingIndex = -1;
-
- point nearestIndexedNeighbourPos = vector::zero;
-
- scalar minDistSqrToNearestIndexedNeighbour = VGREAT;
-
- for (label f = 0; f < 4; f++)
- {
- if (edgeIsShort[f] && neighbourAlreadyIndexed[f])
- {
- label cNI = cit->neighbor(f)->cellIndex();
+ cit->cellIndex() = dualVertI;
+ points[dualVertI] = topoint(dual(cit));
+ dualVertI++;
+ }
+ else
+ {
+ cit->cellIndex() = -1;
+ }
+ }
- point indexedNeighbourPos = points[cNI];
+ points.setSize(dualVertI);
- if
- (
- magSqr(indexedNeighbourPos - dualVertex)
- < minDistSqrToNearestIndexedNeighbour
- )
- {
- nearestExistingIndex = cNI;
+ // // ~~~~~~~~~~~ removing short edges by indexing dual vertices ~~~~~~~~~~~~~~
- nearestIndexedNeighbourPos =
- indexedNeighbourPos;
+ // for
+ // (
+ // Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+ // cit != finite_cells_end();
+ // ++cit
+ // )
+ // {
+ // cit->cellIndex() = -1;
+ // }
- minDistSqrToNearestIndexedNeighbour =
- magSqr(indexedNeighbourPos - dualVertex);
- }
- }
- }
+ // points.setSize(number_of_cells());
- if
- (
- nearestExistingIndex > -1
- && minDistSqrToNearestIndexedNeighbour < minEdgeLenSqr
- )
- {
- points[nearestExistingIndex] =
- 0.5*(dualVertex + nearestIndexedNeighbourPos);
+ // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // // Looking up minEdgeLenSqr with a dummy point, in future will be available
+ // // as a local value to be looked up in-place.
- for (label f = 0; f < 4; f++)
- {
- if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
- {
- cit->neighbor(f)->cellIndex() =
- nearestExistingIndex;
- }
- }
+ // scalar minEdgeLenSqr = sqr(minimumEdgeLength(vector::zero));
- cit->cellIndex() = nearestExistingIndex;
- }
- else
- {
- for (label f = 0; f < 4; f++)
- {
- if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
- {
- cit->neighbor(f)->cellIndex() = dualVertI;
- }
- }
+ // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- cit->cellIndex() = dualVertI;
+ // label dualVertI = 0;
- points[dualVertI] = dualVertex;
+ // // Scanning by number of short (dual) edges (nSE) attached to the
+ // // circumcentre of each Delaunay tet. A Delaunay tet may only have four
+ // // dual edges emanating from its circumcentre, assigning positions and
+ // // indices to those with 4 short edges attached first, then >= 3, then >= 2
+ // // etc.
+ // for (label nSE = 4; nSE >= 0; nSE--)
+ // {
+ // Info<< nl << "Scanning for dual vertices with >= "
+ // << nSE
+ // << " short edges attached." << endl;
+
+ // for
+ // (
+ // Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+ // cit != finite_cells_end();
+ // ++cit
+ // )
+ // {
+ // // If the Delaunay tet has an index already then it has either
+ // // evaluated itself and taken action or has had its index dictated
+ // // by a neighbouring tet with more short edges attached.
- dualVertI++;
- }
- }
- }
- }
- }
+ // if (cit->cellIndex() == -1)
+ // {
+ // point dualVertex = topoint(dual(cit));
+
+ // label shortEdges = 0;
+
+ // List<bool> edgeIsShort(4, false);
+
+ // List<bool> neighbourAlreadyIndexed(4, false);
+
+ // // Loop over the four facets of the Delaunay tet
+ // for (label f = 0; f < 4; f++)
+ // {
+ // // Check that at least one of the vertices of the facet is
+ // // an internal or boundary point
+ // if
+ // (
+ // cit->vertex(vertex_triple_index(f, 0))->
+ // internalOrBoundaryPoint()
+ // || cit->vertex(vertex_triple_index(f, 1))->
+ // internalOrBoundaryPoint()
+ // || cit->vertex(vertex_triple_index(f, 2))->
+ // internalOrBoundaryPoint()
+ // )
+ // {
+ // point neighDualVertex;
+
+ // label cNI = cit->neighbor(f)->cellIndex();
+
+ // if (cNI == -1)
+ // {
+ // neighDualVertex = topoint(dual(cit->neighbor(f)));
+ // }
+ // else
+ // {
+ // neighDualVertex = points[cNI];
+ // }
+
+ // if
+ // (
+ // magSqr(dualVertex - neighDualVertex) < minEdgeLenSqr
+ // )
+ // {
+ // edgeIsShort[f] = true;
+
+ // if (cNI > -1)
+ // {
+ // neighbourAlreadyIndexed[f] = true;
+ // }
+
+ // shortEdges++;
+ // }
+ // }
+ // }
+
+ // if (nSE == 0 && shortEdges == 0)
+ // {
+ // // Final iteration and no short edges are found, index
+ // // remaining dual vertices.
+
+ // if
+ // (
+ // cit->vertex(0)->internalOrBoundaryPoint()
+ // || cit->vertex(1)->internalOrBoundaryPoint()
+ // || cit->vertex(2)->internalOrBoundaryPoint()
+ // || cit->vertex(3)->internalOrBoundaryPoint()
+ // )
+ // {
+ // cit->cellIndex() = dualVertI;
+ // points[dualVertI] = dualVertex;
+ // dualVertI++;
+ // }
+ // }
+ // else if
+ // (
+ // shortEdges >= nSE
+ // )
+ // {
+ // // Info<< neighbourAlreadyIndexed << ' '
+ // // << edgeIsShort << endl;
+
+ // label numUnindexedNeighbours = 1;
+
+ // for (label f = 0; f < 4; f++)
+ // {
+ // if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
+ // {
+ // dualVertex += topoint(dual(cit->neighbor(f)));
+
+ // numUnindexedNeighbours++;
+ // }
+ // }
+
+ // dualVertex /= numUnindexedNeighbours;
+
+ // label nearestExistingIndex = -1;
+
+ // point nearestIndexedNeighbourPos = vector::zero;
+
+ // scalar minDistSqrToNearestIndexedNeighbour = VGREAT;
+
+ // for (label f = 0; f < 4; f++)
+ // {
+ // if (edgeIsShort[f] && neighbourAlreadyIndexed[f])
+ // {
+ // label cNI = cit->neighbor(f)->cellIndex();
+
+ // point indexedNeighbourPos = points[cNI];
+
+ // if
+ // (
+ // magSqr(indexedNeighbourPos - dualVertex)
+ // < minDistSqrToNearestIndexedNeighbour
+ // )
+ // {
+ // nearestExistingIndex = cNI;
+
+ // nearestIndexedNeighbourPos =
+ // indexedNeighbourPos;
+
+ // minDistSqrToNearestIndexedNeighbour =
+ // magSqr(indexedNeighbourPos - dualVertex);
+ // }
+ // }
+ // }
+
+ // if
+ // (
+ // nearestExistingIndex > -1
+ // && minDistSqrToNearestIndexedNeighbour < minEdgeLenSqr
+ // )
+ // {
+ // points[nearestExistingIndex] =
+ // 0.5*(dualVertex + nearestIndexedNeighbourPos);
+
+ // for (label f = 0; f < 4; f++)
+ // {
+ // if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
+ // {
+ // cit->neighbor(f)->cellIndex() =
+ // nearestExistingIndex;
+ // }
+ // }
+
+ // cit->cellIndex() = nearestExistingIndex;
+ // }
+ // else
+ // {
+ // for (label f = 0; f < 4; f++)
+ // {
+ // if (edgeIsShort[f] && !neighbourAlreadyIndexed[f])
+ // {
+ // cit->neighbor(f)->cellIndex() = dualVertI;
+ // }
+ // }
+
+ // cit->cellIndex() = dualVertI;
+
+ // points[dualVertI] = dualVertex;
+
+ // dualVertI++;
+ // }
+ // }
+ // }
+ // }
+ // }
- points.setSize(dualVertI);
+ // points.setSize(dualVertI);
// ~~~~~~~~~~~~~~~~~~~~~~~~~ dual cell indexing ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// assigns an index to the Delaunay vertices which will be the dual cell
@@ -2082,10 +2120,17 @@ void Foam::conformalVoronoiMesh::calcDualMesh
} while (cc1 != ccStart);
- if (verticesOnFace.size() >= 3)
- {
- face newDualFace(verticesOnFace);
+ face newDualFace(verticesOnFace);
+
+ bool keepFace = assessFace
+ (
+ newDualFace,
+ averageAnyCellSize(vA, vB),
+ points
+ );
+ if (verticesOnFace.size() >= 3 && keepFace)
+ {
label dcA = vA->index();
if (!vA->internalOrBoundaryPoint())
@@ -2405,6 +2450,71 @@ void Foam::conformalVoronoiMesh::calcTetMesh
}
+bool Foam::conformalVoronoiMesh::assessFace
+(
+ const face& f,
+ scalar averageCellSize,
+ const pointField& pts
+) const
+{
+ scalar smallFaceAreaCoeff = 0.1;
+ scalar highAspectRatioFaceAreaCoeff = 0.2;
+ scalar aspectRatioLimit = 1.2;
+ scalar targetArea = sqr(averageCellSize);
+
+ const edgeList& eds = f.edges();
+
+ scalar perimeter = 0.0;
+
+ forAll(eds, i)
+ {
+ perimeter += eds[i].mag(pts);
+
+ vector edVec = eds[i].vec(pts);
+ };
+
+ scalar area = f.mag(pts);
+
+ scalar equivalentSqrPerimeter = 4.0*sqrt(area);
+
+ scalar aspectRatio = perimeter/max(equivalentSqrPerimeter, VSMALL);
+
+ bool keepFace = true;
+
+ if (area < smallFaceAreaCoeff*targetArea)
+ {
+ keepFace = false;
+ }
+
+ if
+ (
+ aspectRatio > aspectRatioLimit
+ && area < highAspectRatioFaceAreaCoeff*targetArea)
+ {
+ keepFace = false;
+ }
+
+ if (!keepFace)
+ {
+ Info<< nl << "Area " << area << nl
+ << "averageCellSize " << averageCellSize << nl
+ << "Area ratio "
+ << area/max(sqr(averageCellSize), VSMALL) << nl
+ << "aspectRatio " << aspectRatio << nl
+ << endl;
+
+ forAll(f, i)
+ {
+ meshTools::writeOBJ(Info, pts[f[i]]);
+ }
+
+ Info<< nl;
+ }
+
+ return keepFace;
+}
+
+
void Foam::conformalVoronoiMesh::sortFaces
(
faceList& faces,
@@ -2753,18 +2863,7 @@ void Foam::conformalVoronoiMesh::move()
> cvMeshControls().cosAlignmentAcceptanceAngle()
)
{
- // Arithmetic mean
- // scalar targetCellSize =
- // 0.5*(vA->targetCellSize() + vB->targetCellSize());
-
- // Geometric mean
- scalar targetCellSize =
- sqrt(vA->targetCellSize()*vB->targetCellSize());
-
- // Harmonic mean
- // scalar targetCellSize =
- // 2.0*(vA->targetCellSize()*vB->targetCellSize())
- // /(vA->targetCellSize() + vB->targetCellSize());
+ scalar targetCellSize = averageCellSize(vA, vB);
scalar targetFaceArea = sqr(targetCellSize);
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 0b3c57e249bd6d30620d50c2e4eb745abb1c1ad9..dcd5eb7b9fb232fccb0f728eeec018eb8ca88ea1 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -171,6 +171,24 @@ private:
bool isSurfacePoint = false
) const;
+ //- Return the target cell size from that stored on a pair of
+ // Delaunay vertices, using a mean function.
+ inline scalar averageCellSize
+ (
+ const Vertex_handle& vA,
+ const Vertex_handle& vB
+ ) const;
+
+ //- Return the target cell size from that stored on a pair of
+ // Delaunay vertices, including the possibility that one of
+ // them is not an internalOrBoundaryPoint, and so will not
+ // have valid data.
+ inline scalar averageAnyCellSize
+ (
+ const Vertex_handle& vA,
+ const Vertex_handle& vB
+ ) const;
+
//- Return the local point pair separation at the given location
inline scalar pointPairDistance(const point& pt) const;
@@ -425,6 +443,14 @@ private:
labelList& patchStarts
);
+ //- Assess face to see if it is a candidate for removal
+ bool assessFace
+ (
+ const face& f,
+ scalar averageCellSize,
+ const pointField& pts
+ ) const;
+
//- Sort the faces, owner and neighbour lists into
// upper-triangular order. For internal faces only, use
// before adding patch faces.
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index aed10fdd25265097513566f25bc28a8ad7362f7c..adf33d418f08ae66bb098cffa77cba59945f8472 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -46,6 +46,62 @@ inline Foam::scalar Foam::conformalVoronoiMesh::targetCellSize
}
+inline Foam::scalar Foam::conformalVoronoiMesh::averageCellSize
+(
+ const Vertex_handle& vA,
+ const Vertex_handle& vB
+) const
+{
+ // Arithmetic mean
+ // return 0.5*(vA->targetCellSize() + vB->targetCellSize());
+
+ // Geometric mean
+ return sqrt(vA->targetCellSize()*vB->targetCellSize());
+
+ // Harmonic mean
+ // return
+ // 2.0*(vA->targetCellSize()*vB->targetCellSize())
+ // /(vA->targetCellSize() + vB->targetCellSize());
+}
+
+
+inline Foam::scalar Foam::conformalVoronoiMesh::averageAnyCellSize
+(
+ const Vertex_handle& vA,
+ const Vertex_handle& vB
+) const
+{
+ if
+ (
+ !vA->internalOrBoundaryPoint()
+ && !vB->internalOrBoundaryPoint()
+ )
+ {
+ FatalErrorIn
+ (
+ "Foam::conformalVoronoiMesh::averageAnyCellSize"
+ "("
+ "const Vertex_handle& vA,"
+ "const Vertex_handle& vB"
+ ") const"
+ )
+ << "Requested averageCellSize for to external points"
+ << nl
+ << exit(FatalError);
+ }
+ else if (!vB->internalOrBoundaryPoint())
+ {
+ return vA->targetCellSize();
+ }
+ else if (!vA->internalOrBoundaryPoint())
+ {
+ return vB->targetCellSize();
+ }
+
+ return averageCellSize(vA, vB);
+}
+
+
inline Foam::scalar Foam::conformalVoronoiMesh::pointPairDistance
(
const point& pt
@@ -217,7 +273,8 @@ inline void Foam::conformalVoronoiMesh::insertVb(const Vb& v, label offset)
{
FatalErrorIn("Foam::conformalVoronoiMesh::insertVb(const Vb& v")
<< "Failed to reinsert Vb at " << topoint(Pt)
- << endl;
+ << nl
+ << exit(FatalError);
}
else
{
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
index bfc0ac3298420243f6669888b28716779b6cf0df..dfd47e49a3ba951d674a695167b2fdbb5d614adc 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
@@ -27,7 +27,7 @@ Class
Description
An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep
- track of the vertices in the triangulation.
+ track of the Delaunay cells (tets) in the tessellation.
\*---------------------------------------------------------------------------*/
diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
index 32b8c7a9e264d50eff02b118f4089a97ddfc2acf..5c71317999e23f6882a65bb8f00d4fcc4acebbbd 100644
--- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
+++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
@@ -27,7 +27,7 @@ Class
Description
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep
- track of the vertices in the triangulation.
+ track of the Delaunay vertices in the tessellation.
\*---------------------------------------------------------------------------*/
@@ -97,7 +97,7 @@ public:
index_(INTERNAL_POINT),
type_(INTERNAL_POINT),
alignment_(),
- targetCellSize_()
+ targetCellSize_(0.0)
{}
indexedVertex(const Point& p)
@@ -106,7 +106,7 @@ public:
index_(INTERNAL_POINT),
type_(INTERNAL_POINT),
alignment_(),
- targetCellSize_()
+ targetCellSize_(0.0)
{}
indexedVertex(const Point& p, Cell_handle f)
@@ -115,7 +115,7 @@ public:
index_(INTERNAL_POINT),
type_(INTERNAL_POINT),
alignment_(),
- targetCellSize_()
+ targetCellSize_(0.0)
{}
indexedVertex(Cell_handle f)
@@ -124,7 +124,7 @@ public:
index_(INTERNAL_POINT),
type_(INTERNAL_POINT),
alignment_(),
- targetCellSize_()
+ targetCellSize_(0.0)
{}