mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
Surface conformation point pair insertions and polyMesh output (without patching) included.
This commit is contained in:
graham
@ -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
|
||||
|
||||
scalar searchDistanceSqr = sqr(defaultCellSize*surfDepthCoeff);
|
||||
|
||||
pointIndexHit pHit = geometryToConformTo_.findNearest
|
||||
readScalar
|
||||
(
|
||||
cvMeshDict.subDict("surfaceConformation").lookup
|
||||
(
|
||||
vert,
|
||||
searchDistanceSqr
|
||||
);
|
||||
"pointPairDistanceCoeff"
|
||||
)
|
||||
)
|
||||
);
|
||||
|
||||
if (pHit.hit())
|
||||
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
|
||||
for(label iterationNo = 0; iterationNo < 1; iterationNo++)
|
||||
{
|
||||
DynamicList<scalar> surfacePpDist;
|
||||
DynamicList<point> surfacePoints;
|
||||
DynamicList<vector> surfaceNormals;
|
||||
|
||||
for
|
||||
(
|
||||
Triangulation::Finite_vertices_iterator vit =
|
||||
finite_vertices_begin();
|
||||
vit != finite_vertices_end();
|
||||
vit++
|
||||
)
|
||||
{
|
||||
if (vit->internalPoint())
|
||||
{
|
||||
vit->setNearBoundary();
|
||||
point vert(topoint(vit->point()));
|
||||
|
||||
if (dualCellSurfaceIntersection(vit))
|
||||
// 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"
|
||||
)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -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;
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -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];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user