Commit 59ba9207 authored by laurence's avatar laurence
Browse files

BUG: foamyHexMesh: Remove internal points that are outside the surface.

parent 7c00bbe7
......@@ -787,8 +787,8 @@ Foam::face Foam::conformalVoronoiMesh::buildDualFace
<< "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
<< vA->info() << " "
<< vB->info() << nl
<< exit(FatalError);
}
......@@ -1630,6 +1630,8 @@ void Foam::conformalVoronoiMesh::move()
);
}
DynamicList<Vertex_handle> pointsToRemove;
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
......@@ -1640,15 +1642,18 @@ void Foam::conformalVoronoiMesh::move()
if
(
(vit->internalPoint() || vit->internalBoundaryPoint())
&& !vit->referred()
//&& !vit->referred()
)
{
bool inside = geometryToConformTo_.inside
(
topoint(vit->point())
);
const Foam::point& pt = topoint(vit->point());
bool inside = geometryToConformTo_.inside(pt);
if (!inside)
if
(
!inside
|| !geometryToConformTo_.globalBounds().contains(pt)
)
{
if
(
......@@ -1658,13 +1663,16 @@ void Foam::conformalVoronoiMesh::move()
{
str().write(topoint(vit->point()));
}
remove(vit);
pointsToRemove.append(vit);
internalPtIsOutside++;
}
}
}
Info<< " " << internalPtIsOutside
remove(pointsToRemove.begin(), pointsToRemove.end());
Info<< " " << returnReduce(internalPtIsOutside, sumOp<label>())
<< " internal points were inserted outside the domain. "
<< "They have been removed." << endl;
}
......
......@@ -2458,6 +2458,54 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
|| (vB->internalOrBoundaryPoint() && !vB->referred())
)
{
if
(
(vA->internalPoint() && vB->externalBoundaryPoint())
|| (vB->internalPoint() && vA->externalBoundaryPoint())
)
{
Cell_circulator ccStart = incident_cells(*eit);
Cell_circulator cc1 = ccStart;
Cell_circulator cc2 = cc1;
cc2++;
bool skipEdge = false;
do
{
if
(
cc1->hasFarPoint() || cc2->hasFarPoint()
|| is_infinite(cc1) || is_infinite(cc2)
)
{
Pout<< "Ignoring edge between internal and external: "
<< vA->info()
<< vB->info();
skipEdge = true;
break;
}
cc1++;
cc2++;
} while (cc1 != ccStart);
// Do not create faces if the internal point is outside!
// This occurs because the internal point is not determined to
// be outside in the inside/outside test. This is most likely
// due to the triangle.nearestPointClassify test not returning
// edge/point as the nearest type.
if (skipEdge)
{
continue;
}
}
face newDualFace = buildDualFace(eit);
if (newDualFace.size() >= 3)
......
......@@ -569,7 +569,20 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
Vertex_handle vA = c->vertex(eit->second);
Vertex_handle vB = c->vertex(eit->third);
if (vA->referred() || vB->referred())
if
(
vA->referred()
|| vB->referred()
)
{
continue;
}
if
(
(vA->internalPoint() && vA->referred())
|| (vB->internalPoint() && vB->referred())
)
{
continue;
}
......@@ -2183,59 +2196,59 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
pointIndexHitAndFeature(edHit, featureHit)
);
}
else if (firstPass)
{
label hitIndex = nearestEdgeHit.index();
// Calc new edge location
Foam::point newPt =
0.5
*(
nearestEdgeHit.hitPoint()
+ edHit.hitPoint()
);
pointIndexHit pHitOld =
edgeLocationTreePtr_().findNearest
(
nearestEdgeHit.hitPoint(), GREAT
);
pointIndexHit pHitNew =
edgeLocationTreePtr_().findNearest
(
edHit.hitPoint(), GREAT
);
if
(
magSqr(pHitNew.hitPoint() - edHit.hitPoint())
< magSqr
(
pHitOld.hitPoint()
- nearestEdgeHit.hitPoint()
)
)
{
edHit.setPoint(edHit.hitPoint());
featureEdgeHits[hitIndex] =
pointIndexHitAndFeature(edHit, featureHit);
existingEdgeLocations_[hitIndex] =
edHit.hitPoint();
// Change edge location in featureEdgeHits
// remove index from edge tree
// reinsert new point into tree
edgeLocationTreePtr_().remove(hitIndex);
edgeLocationTreePtr_().insert
(
hitIndex,
hitIndex + 1
);
}
}
// else if (firstPass)
// {
// label hitIndex = nearestEdgeHit.index();
//
// // Calc new edge location
//// Foam::point newPt =
//// 0.5
//// *(
//// nearestEdgeHit.hitPoint()
//// + edHit.hitPoint()
//// );
//
// pointIndexHit pHitOld =
// edgeLocationTreePtr_().findNearest
// (
// nearestEdgeHit.hitPoint(), GREAT
// );
//
// pointIndexHit pHitNew =
// edgeLocationTreePtr_().findNearest
// (
// edHit.hitPoint(), GREAT
// );
//
// if
// (
// magSqr(pHitNew.hitPoint() - edHit.hitPoint())
// < magSqr
// (
// pHitOld.hitPoint()
// - nearestEdgeHit.hitPoint()
// )
// )
// {
// edHit.setPoint(edHit.hitPoint());
//
// featureEdgeHits[hitIndex] =
// pointIndexHitAndFeature(edHit, featureHit);
//
// existingEdgeLocations_[hitIndex] =
// edHit.hitPoint();
//
// // Change edge location in featureEdgeHits
// // remove index from edge tree
// // reinsert new point into tree
// edgeLocationTreePtr_().remove(hitIndex);
// edgeLocationTreePtr_().insert
// (
// hitIndex,
// hitIndex + 1
// );
// }
// }
}
}
}
......
......@@ -1458,8 +1458,6 @@ void Foam::conformalVoronoiMesh::createFeaturePoints(DynamicList<Vb>& pts)
forAll(feMeshes, i)
{
Info<< indent << "Edge mesh = " << feMeshes[i].name() << nl << endl;
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
for
......
......@@ -101,7 +101,9 @@ void Foam::conformalVoronoiMesh::drawDelaunayCell
// Supply offset as tet number
offset *= 4;
os << "# cell index: " << label(c->cellIndex()) << endl;
os << "# cell index: " << label(c->cellIndex())
<< " INT_MIN = " << INT_MIN
<< endl;
os << "# circumradius "
<< mag(c->dual() - topoint(c->vertex(0)->point()))
......@@ -112,7 +114,15 @@ void Foam::conformalVoronoiMesh::drawDelaunayCell
os << "# index / type / procIndex: "
<< label(c->vertex(i)->index()) << " "
<< label(c->vertex(i)->type()) << " "
<< label(c->vertex(i)->procIndex()) << endl;
<< label(c->vertex(i)->procIndex())
<< (is_infinite(c->vertex(i)) ? " # This vertex is infinite!" : "")
<<
(
c->vertex(i)->uninitialised()
? " # This vertex is uninitialised!"
: ""
)
<< endl;
meshTools::writeOBJ(os, topoint(c->vertex(i)->point()));
}
......
......@@ -706,8 +706,8 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInside
// Info<< surface.name() << " = "
// << volumeType::names[surfaceVolumeTests[s][i]] << endl;
//if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
if (surfaceVolumeTests[s][i] != volumeType::INSIDE)
if (surfaceVolumeTests[s][i] == volumeType::OUTSIDE)
// if (surfaceVolumeTests[s][i] != volumeType::INSIDE)
{
insidePoint[i] = false;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment