From ce651255742b9f4d2fb9305cad2126a44dd7ecd5 Mon Sep 17 00:00:00 2001
From: graham <g.macpherson@opencfd.co.uk>
Date: Fri, 10 Apr 2009 17:31:59 +0100
Subject: [PATCH] Surface conformation point pair insertions and polyMesh
output (without patching) included.
---
.../utilities/mesh/generation/cvMesh/cvMesh.C | 1 -
.../conformalVoronoiMesh.C | 817 +++++++++++++++++-
.../conformalVoronoiMesh.H | 42 +
.../conformalVoronoiMeshIO.C | 119 +++
.../conformationSurfaces.C | 38 +
.../conformationSurfaces.H | 11 +
6 files changed, 993 insertions(+), 35 deletions(-)
diff --git a/applications/utilities/mesh/generation/cvMesh/cvMesh.C b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
index 4fc80a526d5..8cf41aed2c8 100644
--- a/applications/utilities/mesh/generation/cvMesh/cvMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
@@ -64,7 +64,6 @@ int main(int argc, char *argv[])
Info<< nl << "Time = " << runTime.timeName()
<< endl;
-
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
diff --git a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 94054bf4d79..f03249f78ae 100644
--- a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -66,7 +66,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceIntersection
point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
// If other edge end is outside bounding box then edge cuts boundary
- if (!!geometryToConformTo_.bounds().contains(dE1))
+ if (!geometryToConformTo_.bounds().contains(dE1))
{
return true;
}
@@ -82,6 +82,655 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceIntersection
}
+void Foam::conformalVoronoiMesh::calcDualMesh
+(
+ pointField& points,
+ faceList& faces,
+ labelList& owner,
+ labelList& neighbour,
+ wordList& patchNames,
+ labelList& patchSizes,
+ labelList& patchStarts
+)
+{
+ // ~~~~~~~~~~~ removing short edges by indexing dual vertices ~~~~~~~~~~~~~~
+
+ for
+ (
+ Triangulation::Finite_cells_iterator cit = finite_cells_begin();
+ cit != finite_cells_end();
+ ++cit
+ )
+ {
+ cit->cellIndex() = -1;
+ }
+
+ points.setSize(number_of_cells());
+
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Looking up details from a dictionary, in future the will be available
+ // from a controls class.
+
+ const dictionary& cvMeshDict( cvMeshControls_.cvMeshDict());
+
+ scalar defaultCellSize
+ (
+ readScalar
+ (
+ cvMeshDict.subDict("motionControl").lookup("defaultCellSize")
+ )
+ );
+
+ scalar minimumEdgeLengthCoeff
+ (
+ readScalar
+ (
+ cvMeshDict.subDict("polyMeshFiltering").lookup
+ (
+ "minimumEdgeLengthCoeff"
+ )
+ )
+ );
+
+ scalar minEdgeLenSqr = Foam::sqr(defaultCellSize*minimumEdgeLengthCoeff);
+
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ 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
+ (
+ 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.
+
+ 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);
+
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~ dual cell indexing ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // assigns an index to the Delaunay vertices which will be the dual cell
+ // index used for owner neighbour assignment.
+
+ // The indices of the points are reset which destroys the point-pair
+ // matching, so the type of each vertex are reset to avoid any ambiguity.
+
+ label dualCelli = 0;
+
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->internalOrBoundaryPoint())
+ {
+ vit->type() = Vb::INTERNAL_POINT;
+ vit->index() = dualCelli;
+ dualCelli++;
+ }
+ else
+ {
+ vit->type() = Vb::FAR_POINT;
+ vit->index() = -1;
+ }
+ }
+
+ // ~~~~~~~~~~~~ dual face and owner neighbour construction ~~~~~~~~~~~~~~~~~
+
+ //label nPatches = qSurf_.patches().size() + 1;
+
+ //label defaultPatchIndex = qSurf_.patches().size();
+
+ label nPatches = 1;
+
+ label defaultPatchIndex = 0;
+
+ patchNames.setSize(nPatches);
+
+ //const geometricSurfacePatchList& surfacePatches = qSurf_.patches();
+
+ // forAll(surfacePatches, sP)
+ // {
+ // patchNames[sP] = surfacePatches[sP].name();
+ // }
+
+ patchNames[defaultPatchIndex] = "cvMesh_defaultPatch";
+
+ patchSizes.setSize(nPatches);
+
+ patchStarts.setSize(nPatches);
+
+ List<DynamicList<face> > patchFaces(nPatches, DynamicList<face>(0));
+
+ List<DynamicList<label> > patchOwners(nPatches, DynamicList<label>(0));
+
+ faces.setSize(number_of_edges());
+
+ owner.setSize(number_of_edges());
+
+ neighbour.setSize(number_of_edges());
+
+ label dualFacei = 0;
+
+ for
+ (
+ Triangulation::Finite_edges_iterator eit = finite_edges_begin();
+ eit != finite_edges_end();
+ ++eit
+ )
+ {
+ Cell_handle c = eit->first;
+ Vertex_handle vA = c->vertex(eit->second);
+ Vertex_handle vB = c->vertex(eit->third);
+
+ if
+ (
+ vA->internalOrBoundaryPoint()
+ || vB->internalOrBoundaryPoint()
+ )
+ {
+ Cell_circulator ccStart = incident_cells(*eit);
+ Cell_circulator cc1 = ccStart;
+ Cell_circulator cc2 = cc1;
+
+ // Advance the second circulator so that it always stays on the next
+ // cell around the edge;
+ cc2++;
+
+ DynamicList<label> verticesOnFace;
+
+ do
+ {
+ label cc1I = cc1->cellIndex();
+
+ label cc2I = cc2->cellIndex();
+
+
+ if (cc1I < 0 || cc2I < 0)
+ {
+ FatalErrorIn("Foam::conformalVoronoiMesh::calcDualMesh")
+ << "Dual face uses circumcenter defined by a "
+ << "Delaunay tetrahedron with no internal "
+ << "or boundary points. Defining Delaunay edge ends: "
+ << topoint(vA->point()) << " "
+ << topoint(vB->point()) << nl
+ << exit(FatalError);
+ }
+
+ if (cc1I != cc2I)
+ {
+ verticesOnFace.append(cc1I);
+ }
+
+ cc1++;
+
+ cc2++;
+ } while (cc1 != ccStart);
+
+ verticesOnFace.shrink();
+
+ if (verticesOnFace.size() >= 3)
+ {
+ face newDualFace(verticesOnFace);
+
+ label dcA = vA->index();
+
+ if (!vA->internalOrBoundaryPoint())
+ {
+ dcA = -1;
+ }
+
+ label dcB = vB->index();
+
+ if (!vB->internalOrBoundaryPoint())
+ {
+ dcB = -1;
+ }
+
+ label dcOwn = -1;
+ label dcNei = -1;
+
+ if (dcA == -1 && dcB == -1)
+ {
+ FatalErrorIn("calcDualMesh")
+ << "Attempting to create a face joining "
+ << "two external dual cells "
+ << exit(FatalError);
+ }
+ else if (dcA == -1 || dcB == -1)
+ {
+ // boundary face, find which is the owner
+
+ if (dcA == -1)
+ {
+ dcOwn = dcB;
+
+ // reverse face order to correctly orientate normal
+ reverse(newDualFace);
+ }
+ else
+ {
+ dcOwn = dcA;
+ }
+
+ // Find which patch this face is on by finding the
+ // intersection with the surface of the Delaunay edge
+ // generating the face and identify the region of the
+ // intersection.
+
+ point ptA = topoint(vA->point());
+
+ point ptB = topoint(vB->point());
+
+ //pointIndexHit pHit = qSurf_.tree().findLineAny(ptA, ptB);
+
+ //label patchIndex = qSurf_[pHit.index()].region();
+
+ label patchIndex = defaultPatchIndex;
+
+ if (patchIndex == -1)
+ {
+ patchIndex = defaultPatchIndex;
+
+ WarningIn("Foam::conformalVoronoiMesh::calcDualMesh")
+ << "Dual face found that is not on a surface "
+ << "patch. Adding to "
+ << patchNames[defaultPatchIndex]
+ << endl;
+ }
+
+ patchFaces[patchIndex].append(newDualFace);
+ patchOwners[patchIndex].append(dcOwn);
+ }
+ else
+ {
+ // internal face, find the lower cell to be the owner
+
+ if (dcB > dcA)
+ {
+ dcOwn = dcA;
+ dcNei = dcB;
+ }
+ else
+ {
+ dcOwn = dcB;
+ dcNei = dcA;
+
+ // reverse face order to correctly orientate normal
+ reverse(newDualFace);
+ }
+
+ faces[dualFacei] = newDualFace;
+
+ owner[dualFacei] = dcOwn;
+
+ neighbour[dualFacei] = dcNei;
+
+ dualFacei++;
+ }
+ }
+ // else
+ // {
+ // Info<< verticesOnFace.size()
+ // << " size face not created." << endl;
+ // }
+ }
+ }
+
+ label nInternalFaces = dualFacei;
+
+ faces.setSize(nInternalFaces);
+
+ owner.setSize(nInternalFaces);
+
+ neighbour.setSize(nInternalFaces);
+
+ // ~~~~~~~~ sort owner, reordinging neighbour and faces to match ~~~~~~~~~~~
+ // two stage sort for upper triangular order: sort by owner first, then for
+ // each block of owners sort by neighbour
+
+ labelList sortingIndices;
+
+ // Stage 1
+
+ {
+ SortableList<label> sortedOwner(owner);
+
+ sortingIndices = sortedOwner.indices();
+ }
+
+ {
+ labelList copyOwner(owner.size());
+
+ forAll(sortingIndices, sI)
+ {
+ copyOwner[sI] = owner[sortingIndices[sI]];
+ }
+
+ owner = copyOwner;
+ }
+
+ {
+ labelList copyNeighbour(neighbour.size());
+
+ forAll(sortingIndices, sI)
+ {
+ copyNeighbour[sI] = neighbour[sortingIndices[sI]];
+ }
+
+ neighbour = copyNeighbour;
+ }
+
+ {
+ faceList copyFaces(faces.size());
+
+ forAll(sortingIndices, sI)
+ {
+ copyFaces[sI] = faces[sortingIndices[sI]];
+ }
+
+ faces = copyFaces;
+ }
+
+ // Stage 2
+
+ sortingIndices = -1;
+
+ DynamicList<label> ownerCellJumps;
+
+ // Force first owner entry to be a jump
+ ownerCellJumps.append(0);
+
+ for (label o = 1; o < owner.size(); o++)
+ {
+ if (owner[o] > owner[o-1])
+ {
+ ownerCellJumps.append(o);
+ }
+ }
+
+ ownerCellJumps.shrink();
+
+ forAll(ownerCellJumps, oCJ)
+ {
+ label start = ownerCellJumps[oCJ];
+
+ label length;
+
+ if (oCJ == ownerCellJumps.size() - 1)
+ {
+ length = owner.size() - start;
+ }
+ else
+ {
+ length = ownerCellJumps[oCJ + 1] - start;
+ }
+
+ SubList<label> neighbourBlock(neighbour, length, start);
+
+ SortableList<label> sortedNeighbourBlock(neighbourBlock);
+
+ forAll(sortedNeighbourBlock, sNB)
+ {
+ sortingIndices[start + sNB] =
+ sortedNeighbourBlock.indices()[sNB] + start;
+ }
+ }
+
+ // Perform sort
+
+ {
+ labelList copyOwner(owner.size());
+
+ forAll(sortingIndices, sI)
+ {
+ copyOwner[sI] = owner[sortingIndices[sI]];
+ }
+
+ owner = copyOwner;
+ }
+
+ {
+ labelList copyNeighbour(neighbour.size());
+
+ forAll(sortingIndices, sI)
+ {
+ copyNeighbour[sI] = neighbour[sortingIndices[sI]];
+ }
+
+ neighbour = copyNeighbour;
+ }
+
+ {
+ faceList copyFaces(faces.size());
+
+ forAll(sortingIndices, sI)
+ {
+ copyFaces[sI] = faces[sortingIndices[sI]];
+ }
+
+ faces = copyFaces;
+ }
+
+ // ~~~~~~~~ add patch information ~~~~~~~~~~~
+
+ label nBoundaryFaces = 0;
+
+ forAll(patchFaces, p)
+ {
+ patchFaces[p].shrink();
+
+ patchOwners[p].shrink();
+
+ patchSizes[p] = patchFaces[p].size();
+
+ patchStarts[p] = nInternalFaces + nBoundaryFaces;
+
+ nBoundaryFaces += patchSizes[p];
+ }
+
+ faces.setSize(nInternalFaces + nBoundaryFaces);
+
+ owner.setSize(nInternalFaces + nBoundaryFaces);
+
+ forAll(patchFaces, p)
+ {
+ forAll(patchFaces[p], f)
+ {
+ faces[dualFacei] = patchFaces[p][f];
+
+ owner[dualFacei] = patchOwners[p][f];
+
+ dualFacei++;
+ }
+ }
+}
+
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::conformalVoronoiMesh::conformalVoronoiMesh
@@ -126,11 +775,18 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
timeCheck();
insertFeaturePoints();
-
timeCheck();
insertInitialPoints();
+ timeCheck();
+
+ conformToSurface();
+ timeCheck();
+ writePoints("allPoints.obj", false);
+ timeCheck();
+
+ writeMesh();
timeCheck();
}
@@ -145,7 +801,48 @@ Foam::conformalVoronoiMesh::~conformalVoronoiMesh()
void Foam::conformalVoronoiMesh::timeCheck() const
{
- Info<< nl << "--- [ " << runTime_.elapsedCpuTime() << "s ] --- " << endl;
+ Info<< nl << "--- [ " << runTime_.elapsedCpuTime() << "s, delta "
+ << runTime_.cpuTimeIncrement()<< "s ] --- " << endl;
+}
+
+
+void Foam::conformalVoronoiMesh::insertSurfacePointPairs
+(
+ const List<scalar>& surfacePpDist,
+ const List<point>& surfacePoints,
+ const List<vector>& surfaceNormals,
+ const fileName fName
+)
+{
+ if
+ (
+ surfacePpDist.size() != surfacePoints.size()
+ || surfacePpDist.size() != surfaceNormals.size()
+ )
+ {
+ FatalErrorIn("Foam::conformalVoronoiMesh::insertPointPairs")
+ << "surfacePpDist, surfacePoints and surfaceNormals are not "
+ << "the same size. Sizes"
+ << surfacePpDist.size() << ' '
+ << surfacePoints.size() << ' '
+ << surfaceNormals.size()
+ << exit(FatalError);
+ }
+
+ forAll(surfacePoints, p)
+ {
+ insertPointPair
+ (
+ surfacePpDist[p],
+ surfacePoints[p],
+ surfaceNormals[p]
+ );
+ }
+
+ if (fName != fileName::null)
+ {
+ writePoints(fName, surfacePoints);
+ }
}
@@ -203,60 +900,112 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
void Foam::conformalVoronoiMesh::conformToSurface()
{
+ Info<< nl << "Conforming to surfaces" << endl;
+
startOfSurfacePointPairs_ = number_of_vertices();
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Looking up details from a dictionary, in future the will be available
+ // from a controls class.
+
const dictionary& cvMeshDict( cvMeshControls_.cvMeshDict());
scalar defaultCellSize
(
- cvMeshDict.subDict("motionControl").lookup("defaultCellSize")
+ readScalar
+ (
+ cvMeshDict.subDict("motionControl").lookup("defaultCellSize")
+ )
);
scalar surfDepthCoeff
(
- cvMeshDict.subDict("surfaceConformation").lookup
+ readScalar
(
- "surfacePointSearchDepthCoeff"
+ cvMeshDict.subDict("surfaceConformation").lookup
+ (
+ "surfacePointSearchDepthCoeff"
+ )
)
);
- DynamicList<point> surfacePoints;
- DynamicList<point> surfaceNormals;
-
- for
+ scalar ppDistCoeff
(
- Triangulation::Finite_vertices_iterator vit = finite_vertices_begin();
- vit != finite_vertices_end();
- vit++
- )
- {
- if (vit->internalPoint())
- {
- point vert(topoint(vit->point()));
-
- // TODO Need to have a function to recover the local cell size, use
- // the defaultCellSize for the moment
+ readScalar
+ (
+ cvMeshDict.subDict("surfaceConformation").lookup
+ (
+ "pointPairDistanceCoeff"
+ )
+ )
+ );
- scalar searchDistanceSqr = sqr(defaultCellSize*surfDepthCoeff);
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- pointIndexHit pHit = geometryToConformTo_.findNearest
- (
- vert,
- searchDistanceSqr
- );
+ for(label iterationNo = 0; iterationNo < 1; iterationNo++)
+ {
+ DynamicList<scalar> surfacePpDist;
+ DynamicList<point> surfacePoints;
+ DynamicList<vector> surfaceNormals;
- if (pHit.hit())
+ for
+ (
+ Triangulation::Finite_vertices_iterator vit =
+ finite_vertices_begin();
+ vit != finite_vertices_end();
+ vit++
+ )
+ {
+ if (vit->internalPoint())
{
- vit->setNearBoundary();
-
- if (dualCellSurfaceIntersection(vit))
+ point vert(topoint(vit->point()));
+
+ // TODO Need to have a function to recover the local cell size,
+ // use the defaultCellSize for the moment
+
+ scalar searchDistanceSqr = Foam::sqr
+ (
+ defaultCellSize*surfDepthCoeff
+ );
+ pointIndexHit pHit;
+ vector normal;
+
+ geometryToConformTo_.findNearestAndNormal
+ (
+ vert,
+ searchDistanceSqr,
+ pHit,
+ normal
+ );
+
+ if (pHit.hit())
{
- allNearSurfacePoints.append(vert);
- allSurfacePoints.append(pHit.hitPoint());
- allSurfaceTris.append(pHit.index());
+ vit->setNearBoundary();
+
+ if (dualCellSurfaceIntersection(vit))
+ {
+ surfacePpDist.append(defaultCellSize*ppDistCoeff);
+ surfacePoints.append(pHit.hitPoint());
+ surfaceNormals.append(normal);
+ }
}
}
}
+
+ Info<< nl <<iterationNo << ": "
+ << number_of_vertices() << ": "
+ << surfacePoints.size() << endl;
+
+ insertSurfacePointPairs
+ (
+ surfacePpDist,
+ surfacePoints,
+ surfaceNormals,
+ fileName
+ (
+ "surfaceConformationLocations_" + name(iterationNo) + ".obj"
+ )
+ );
}
}
diff --git a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 249adabe288..571e1d0e619 100644
--- a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -48,6 +48,7 @@ SourceFiles
#include "DynamicList.H"
#include "Time.H"
#include "polyMesh.H"
+#include "SortableList.H"
#include "meshTools.H"
#include "triSurfaceTools.H"
#include "mathematicalConstants.H"
@@ -122,6 +123,16 @@ class conformalVoronoiMesh
//- Insert a Vb (a typedef of CGAL::indexedVertex<K>)
inline void insertVb(const Vb& v);
+ //- Insert pairs of points on the surface with the given normals, at the
+ // specified spacing
+ void insertSurfacePointPairs
+ (
+ const List<scalar>& surfacePpDist,
+ const List<point>& surfacePoints,
+ const List<vector>& surfaceNormals,
+ const fileName fName = fileName::null
+ );
+
//- Insert point groups at the feature points.
void insertFeaturePoints();
@@ -139,6 +150,18 @@ class conformalVoronoiMesh
const Triangulation::Finite_vertices_iterator& vit
) const;
+ //- Dual calculation
+ void calcDualMesh
+ (
+ pointField& points,
+ faceList& faces,
+ labelList& owner,
+ labelList& neighbour,
+ wordList& patchNames,
+ labelList& patchSizes,
+ labelList& patchStarts
+ );
+
//- Disallow default bitwise copy construct
conformalVoronoiMesh(const conformalVoronoiMesh&);
@@ -201,6 +224,25 @@ public:
//- Write Delaunay points to .obj file
void writePoints(const fileName& fName, bool internalOnly) const;
+
+ //- Write list of points to file
+ void writePoints
+ (
+ const fileName& fName,
+ const List<point>& points
+ ) const;
+
+ //- Write polyMesh
+ void writeMesh();
+
+ //- Write dual points and faces as .obj file
+ void writeDual
+ (
+ const pointField& points,
+ const faceList& faces,
+ const fileName& fName
+ ) const;
+
};
diff --git a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 3c81334eae5..ea741b1d634 100644
--- a/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/src/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -55,4 +55,123 @@ void Foam::conformalVoronoiMesh::writePoints
}
+void Foam::conformalVoronoiMesh::writePoints
+(
+ const fileName& fName,
+ const List<point>& points
+) const
+{
+ Info<< nl << "Writing " << points.size() << " points from pointList to "
+ << fName << endl;
+
+ OFstream str(fName);
+
+ forAll(points, p)
+ {
+ meshTools::writeOBJ(str, points[p]);
+ }
+}
+
+
+void Foam::conformalVoronoiMesh::writeMesh()
+{
+ pointField points(0);
+ faceList faces(0);
+ labelList owner(0);
+ labelList neighbour(0);
+ wordList patchNames(0);
+ labelList patchSizes(0);
+ labelList patchStarts(0);
+
+ calcDualMesh
+ (
+ points,
+ faces,
+ owner,
+ neighbour,
+ patchNames,
+ patchSizes,
+ patchStarts
+ );
+
+ writeDual(points, faces, "dualMesh.obj");
+
+ IOobject io
+ (
+ Foam::polyMesh::defaultRegion,
+ runTime_.constant(),
+ runTime_,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ );
+
+ Info<< nl << "Writing polyMesh to constant." << endl;
+
+
+ polyMesh pMesh
+ (
+ io,
+ xferMove(points),
+ xferMove(faces),
+ xferMove(owner),
+ xferMove(neighbour)
+ );
+
+ List<polyPatch*> patches(patchStarts.size());
+
+ forAll (patches, p)
+ {
+ patches[p] = new polyPatch
+ (
+ patchNames[p],
+ patchSizes[p],
+ patchStarts[p],
+ p,
+ pMesh.boundaryMesh()
+ );
+ }
+
+ pMesh.addPatches(patches);
+
+ if (!pMesh.write())
+ {
+ FatalErrorIn("Foam::conformalVoronoiMesh::writeMesh()")
+ << "Failed writing polyMesh."
+ << exit(FatalError);
+ }
+}
+
+
+void Foam::conformalVoronoiMesh::writeDual
+(
+ const pointField& points,
+ const faceList& faces,
+ const fileName& fName
+) const
+{
+ Info<< nl << "Writing dual points and faces to " << fName << endl;
+
+ OFstream str(fName);
+
+ forAll(points, p)
+ {
+ meshTools::writeOBJ(str, points[p]);
+ }
+
+ forAll (faces, f)
+ {
+ str<< 'f';
+
+ const face& fP = faces[f];
+
+ forAll(fP, p)
+ {
+ str<< ' ' << fP[p] + 1;
+ }
+
+ str<< nl;
+ }
+}
+
+
// ************************************************************************* //
diff --git a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 759810229e0..abde13e3f16 100644
--- a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -256,4 +256,42 @@ bool Foam::conformationSurfaces::findAnyIntersection
}
+void Foam::conformationSurfaces::findNearestAndNormal
+(
+ const point& sample,
+ scalar nearestDistSqr,
+ pointIndexHit& pHit,
+ vector& normal
+) const
+{
+ labelList hitSurfaces;
+ List<pointIndexHit> hitInfo;
+
+ searchableSurfacesQueries::findNearest
+ (
+ allGeometry_,
+ surfaces_,
+ pointField(1, sample),
+ scalarField(1, nearestDistSqr),
+ hitSurfaces,
+ hitInfo
+ );
+
+ pHit = hitInfo[0];
+
+ if (pHit.hit())
+ {
+ vectorField normals;
+
+ // hitSurfaces has returned the index of the entry in surfaces_ that was
+ // found, not the index of the surface in allGeometry_, translating this
+ // on access to allGeometry_ for the normal lookup.
+
+ allGeometry_[surfaces_[hitSurfaces[0]]].getNormal(hitInfo, normals);
+
+ normal = normals[0];
+ }
+}
+
+
// ************************************************************************* //
diff --git a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
index cd75b1f29ff..1ec6734fd9e 100644
--- a/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
+++ b/src/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
@@ -40,6 +40,7 @@ SourceFiles
#include "searchableSurfaces.H"
#include "searchableSurfacesQueries.H"
#include "featureEdgeMesh.H"
+#include "triSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -146,6 +147,16 @@ public:
// Finding if the line joining start and end intersects the surface
bool findAnyIntersection(point start, point end) const;
+ //- Find the nearest point to the sample and return it to the
+ // pointIndexHit, and the normal at the hit location, if found.
+ void findNearestAndNormal
+ (
+ const point& sample,
+ scalar nearestDistSqr,
+ pointIndexHit& pHit,
+ vector& normal
+ ) const;
+
// Member Operators
--
GitLab