ENH: Improvements to sizing and surface conformation and memory use.

Detecting surface incursions as well as protrusions.

When a cell size is requested for a point outside the surfaces, the nearest
surface size is used to avoid an inappropriate return of the default size.

Truncating point to double before adding displacement vector to avoid memory use
growth (a linear growth for gcc 4.5 - bug, see CGAL mailing list 17/1/2011, but
grows regardless).
This commit is contained in:
graham
2011-01-19 10:04:10 +00:00
parent 2398804708
commit e1f00987f0
6 changed files with 375 additions and 116 deletions

View File

@ -37,6 +37,118 @@ defineTypeNameAndDebug(cellSizeControlSurfaces, 0);
}
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
bool Foam::cellSizeControlSurfaces::evalCellSizeFunctions
(
const point& pt,
scalar& minSize,
bool isSurfacePoint
) const
{
bool anyFunctionFound = false;
// // Regions requesting with the same priority take the average
// scalar sizeAccumulator = 0;
// scalar numberOfFunctions = 0;
// label previousPriority = defaultPriority_;
// if (cellSizeFunctions_.size())
// {
// previousPriority =
// cellSizeFunctions_[cellSizeFunctions_.size() - 1].priority();
// forAll(cellSizeFunctions_, i)
// {
// const cellSizeFunction& cSF = cellSizeFunctions_[i];
// if (cSF.priority() < previousPriority && numberOfFunctions > 0)
// {
// return sizeAccumulator/numberOfFunctions;
// }
// scalar sizeI;
// if (cSF.cellSize(pt, sizeI, isSurfacePoint))
// {
// anyFunctionFound = true;
// previousPriority = cSF.priority();
// sizeAccumulator += sizeI;
// numberOfFunctions++;
// }
// }
// }
// if (previousPriority == defaultPriority_ || numberOfFunctions == 0)
// {
// sizeAccumulator += defaultCellSize_;
// numberOfFunctions++;
// }
// minSize = sizeAccumulator/numberOfFunctions;
// return anyFunctionFound;
// Regions requesting with the same priority take the smallest
if (cellSizeFunctions_.size())
{
// Initialise to the last (lowest) priority
label previousPriority = cellSizeFunctions_.last().priority();
forAll(cellSizeFunctions_, i)
{
const cellSizeFunction& cSF = cellSizeFunctions_[i];
if (debug)
{
Info<< "size function "
<< allGeometry_.names()[surfaces_[i]]
<< " priority " << cSF.priority()
<< endl;
}
if (cSF.priority() < previousPriority)
{
return minSize;
}
scalar sizeI;
if (cSF.cellSize(pt, sizeI, isSurfacePoint))
{
anyFunctionFound = true;
if (cSF.priority() == previousPriority)
{
if (sizeI < minSize)
{
minSize = sizeI;
}
}
else
{
minSize = sizeI;
}
if (debug)
{
Info<< "sizeI " << sizeI << " minSize " << minSize << endl;
}
previousPriority = cSF.priority();
}
}
}
return anyFunctionFound;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cellSizeControlSurfaces::cellSizeControlSurfaces
@ -172,102 +284,54 @@ Foam::scalar Foam::cellSizeControlSurfaces::cellSize
bool isSurfacePoint
) const
{
scalar size = defaultCellSize_;
// // Regions requesting with the same priority take the average
bool anyFunctionFound = evalCellSizeFunctions(pt, size, isSurfacePoint);
// scalar sizeAccumulator = 0;
// scalar numberOfFunctions = 0;
// label previousPriority = defaultPriority_;
// if (cellSizeFunctions_.size())
// {
// previousPriority =
// cellSizeFunctions_[cellSizeFunctions_.size() - 1].priority();
// forAll(cellSizeFunctions_, i)
// {
// const cellSizeFunction& cSF = cellSizeFunctions_[i];
// if (cSF.priority() < previousPriority && numberOfFunctions > 0)
// {
// return sizeAccumulator/numberOfFunctions;
// }
// scalar sizeI;
// if (cSF.cellSize(pt, sizeI, isSurfacePoint))
// {
// previousPriority = cSF.priority();
// sizeAccumulator += sizeI;
// numberOfFunctions++;
// }
// }
// }
// if (previousPriority == defaultPriority_ || numberOfFunctions == 0)
// {
// sizeAccumulator += defaultCellSize_;
// numberOfFunctions++;
// }
// return sizeAccumulator/numberOfFunctions;
// Regions requesting with the same priority take the smallest
scalar minSize = defaultCellSize_;
if (cellSizeFunctions_.size())
if (!anyFunctionFound)
{
// Initialise to the last (lowest) priority
label previousPriority = cellSizeFunctions_.last().priority();
// Check if the point in question was actually inside the domain, if
// not, then it may be falling back to an inappropriate default size.
forAll(cellSizeFunctions_, i)
if (cvMesh_.geometryToConformTo().outside(pt))
{
const cellSizeFunction& cSF = cellSizeFunctions_[i];
pointIndexHit surfHit;
label hitSurface;
if (debug)
{
Info<< "size function "
<< allGeometry_.names()[surfaces_[i]]
<< " priority " << cSF.priority()
<< endl;
}
cvMesh_.geometryToConformTo().findSurfaceNearest
(
pt,
cvMesh_.geometryToConformTo().spanMagSqr(),
surfHit,
hitSurface
);
if (cSF.priority() < previousPriority)
if (!surfHit.hit())
{
return minSize;
}
scalar sizeI;
if (cSF.cellSize(pt, sizeI, isSurfacePoint))
{
if (cSF.priority() == previousPriority)
{
if (sizeI < minSize)
{
minSize = sizeI;
}
FatalErrorIn
(
"Foam::scalar Foam::cellSizeControlSurfaces::cellSize"
"("
"const point& pt, "
"bool isSurfacePoint"
") const"
)
<< "Point " << pt << "was not within "
<< cvMesh_.geometryToConformTo().spanMag()
<< " of the surface." << nl
<< "findSurfaceNearest did not find a hit across the span "
<< "of the surfaces."
<< nl << exit(FatalError) << endl;
}
else
{
minSize = sizeI;
}
if (debug)
{
Info<< "sizeI " << sizeI << " minSize " << minSize << endl;
}
previousPriority = cSF.priority();
// Evaluating the cell size at the nearest surface
evalCellSizeFunctions(surfHit.hitPoint(), size, true);
}
}
}
return minSize;
return size;
}

View File

@ -74,12 +74,21 @@ class cellSizeControlSurfaces
scalar defaultCellSize_;
//- Assigning a priority to all requests for cell sizes, the highest
// overrules.
// overrules
label defaultPriority_;
// Private Member Functions
//- Evaluate the cell size functions, returning a bool stating if a
// function was found or not
bool evalCellSizeFunctions
(
const point& pt,
scalar& minSize,
bool isSurfacePoint
) const;
//- Disallow default bitwise copy construct
cellSizeControlSurfaces(const cellSizeControlSurfaces&);

View File

@ -906,8 +906,7 @@ void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
// for
// (
// Delaunay::Finite_vertices_iterator vit =
// finite_vertices_begin();
// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
// vit != finite_vertices_end();
// vit++
// )
@ -1230,10 +1229,6 @@ void Foam::conformalVoronoiMesh::move()
pointField dualVertices(number_of_cells());
Info<< dualVertices.size() << endl;
timeCheck("Start of move - after dualVertices declaration");
label dualVertI = 0;
// Find the dual point of each tetrahedron and assign it an index.
@ -1262,14 +1257,8 @@ void Foam::conformalVoronoiMesh::move()
}
}
timeCheck("Start of move - before dualVertices.setSize");
dualVertices.setSize(dualVertI);
Info<< dualVertices.size() << endl;
timeCheck("Dual vertices indexed");
setVertexSizeAndAlignment();
timeCheck("Determined sizes and alignments");
@ -1547,10 +1536,21 @@ void Foam::conformalVoronoiMesh::move()
{
if (pointToBeRetained[vit->index()] == true)
{
// Convert vit->point() to FOAM vector (double) to do addition,
// avoids memory increase because a record of the constructions
// would be kept otherwise.
// See cgal-discuss@lists-sop.inria.fr:
// "Memory issue with openSUSE 11.3, exact kernel, adding
// points/vectors"
// 14/1/2011.
pointsToInsert.push_back
(
vit->point()
+ toCGALVector(displacementAccumulator[vit->index()])
toPoint
(
topoint(vit->point())
+ displacementAccumulator[vit->index()]
)
);
}
}

View File

@ -416,7 +416,7 @@ private:
) const;
//- Find the "worst" protrusion of a dual cell through the surface,
// subject to the tolerance
// subject to the maxSurfaceProtrusion tolerance
void dualCellLargestSurfaceProtrusion
(
const Delaunay::Finite_vertices_iterator& vit,
@ -424,6 +424,16 @@ private:
label& hitSurface
) const;
//- Find the "worst" incursion of the dual cell of a non-internal or
// boundary point through the surface, subject to the
// maxSurfaceProtrusion tolerance
void dualCellLargestSurfaceIncursion
(
const Delaunay::Finite_vertices_iterator& vit,
pointIndexHit& surfHitLargest,
label& hitSurfaceLargest
) const;
//- Limit the displacement of a point so that it doesn't penetrate the
// surface to be meshed or come too close to it
void limitDisplacement

View File

@ -132,8 +132,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
for
(
Delaunay::Finite_vertices_iterator vit =
finite_vertices_begin();
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
@ -247,8 +246,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
for
(
Delaunay::Finite_vertices_iterator vit =
finite_vertices_begin();
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
@ -265,6 +263,34 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
dualCellLargestSurfaceProtrusion(vit, surfHit, hitSurface);
if (surfHit.hit())
{
addSurfaceAndEdgeHits
(
vit,
vert,
surfHit,
hitSurface,
surfacePtReplaceDistCoeffSqr,
edgeSearchDistCoeffSqr,
surfaceHits,
hitSurfaces,
featureEdgeHits,
featureEdgeFeaturesHit,
newEdgeLocations,
existingEdgeLocations,
edgeLocationTree
);
}
}
else if (vit->ppSlave())
{
Foam::point vert(topoint(vit->point()));
pointIndexHit surfHit;
label hitSurface;
dualCellLargestSurfaceIncursion(vit, surfHit, hitSurface);
if (surfHit.hit())
{
addSurfaceAndEdgeHits
@ -336,8 +362,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
// for
// (
// Delaunay::Finite_vertices_iterator vit =
// finite_vertices_begin();
// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
// vit != finite_vertices_end();
// vit++
// )
@ -365,6 +390,81 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
// }
// }
// {
// // TEST - ASSESS CLOSE SURFACE POINTS
// setVertexSizeAndAlignment();
// for
// (
// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
// vit != finite_vertices_end();
// vit++
// )
// {
// if
// (
// vit->index() >= startOfSurfacePointPairs_
// && vit->internalOrBoundaryPoint()
// )
// {
// std::list<Vertex_handle> adjacentVertices;
// adjacent_vertices(vit, std::back_inserter(adjacentVertices));
// Foam::point pt = topoint(vit->point());
// // Info<< nl << "vit: " << vit->index() << " "
// // << topoint(vit->point())
// // << endl;
// // Info<< adjacentVertices.size() << endl;
// for
// (
// std::list<Vertex_handle>::iterator
// avit = adjacentVertices.begin();
// avit != adjacentVertices.end();
// ++avit
// )
// {
// Vertex_handle avh = *avit;
// // The lower indexed vertex will perform the assessment
// if
// (
// avh->index() >= startOfSurfacePointPairs_
// && avh->internalOrBoundaryPoint()
// && vit->index() < avh->index()
// && vit->type() != avh->type()
// )
// {
// scalar targetSize = 0.2*averageAnyCellSize(vit, avh);
// // Info<< "diff " << mag(pt - topoint(avh->point()))
// // << " " << targetSize << endl;
// if
// (
// magSqr(pt - topoint(avh->point()))
// < sqr(targetSize)
// )
// {
// Info<< nl << "vit: " << vit->index() << " "
// << topoint(vit->point())
// << endl;
// Info<< " adjacent too close: "
// << avh->index() << " "
// << topoint(avh->point())
// << endl;
// }
// }
// }
// }
// }
// }
storeSurfaceConformation();
}
@ -493,6 +593,84 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
}
void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
(
const Delaunay::Finite_vertices_iterator& vit,
pointIndexHit& surfHitLargest,
label& hitSurfaceLargest
) const
{
std::list<Facet> facets;
incident_facets(vit, std::back_inserter(facets));
Foam::point vert(topoint(vit->point()));
scalar minIncursionDistance = -maxSurfaceProtrusion(vert);
for
(
std::list<Facet>::iterator fit=facets.begin();
fit != facets.end();
++fit
)
{
if
(
!is_infinite(fit->first)
&& !is_infinite(fit->first->neighbor(fit->second))
)
{
Foam::point edgeMid =
0.5
*(
topoint(dual(fit->first))
+ topoint(dual(fit->first->neighbor(fit->second)))
);
pointIndexHit surfHit;
label hitSurface;
geometryToConformTo_.findSurfaceAnyIntersection
(
vert,
edgeMid,
surfHit,
hitSurface
);
if (surfHit.hit())
{
vectorField norm(1);
allGeometry_[hitSurface].getNormal
(
List<pointIndexHit>(1, surfHit),
norm
);
const vector& n = norm[0];
scalar normalIncursionDistance =
(edgeMid - surfHit.hitPoint()) & n;
if (normalIncursionDistance < minIncursionDistance)
{
surfHitLargest = surfHit;
hitSurfaceLargest = hitSurface;
minIncursionDistance = normalIncursionDistance;
// Info<< nl << "# Incursion: " << endl;
// meshTools::writeOBJ(Info, vert);
// meshTools::writeOBJ(Info, edgeMid);
// Info<< "l Na Nb" << endl;
}
}
}
}
}
void Foam::conformalVoronoiMesh::limitDisplacement
(
const Delaunay::Finite_vertices_iterator& vit,

View File

@ -66,17 +66,15 @@ inline Foam::scalar Foam::conformalVoronoiMesh::averageAnyCellSize
&& !vB->internalOrBoundaryPoint()
)
{
FatalErrorIn
// There are no internalOrBoundaryPoints available, determine
// size from scratch
// Geometric mean
return sqrt
(
"Foam::conformalVoronoiMesh::averageAnyCellSize"
"("
"const Vertex_handle& vA,"
"const Vertex_handle& vB"
") const"
)
<< "Requested averageCellSize for two external points"
<< nl
<< exit(FatalError);
targetCellSize(topoint(vA->point()))
*targetCellSize(topoint(vB->point()))
);
}
else if (!vB->internalOrBoundaryPoint())
{
@ -391,9 +389,9 @@ Foam::conformalVoronoiMesh::toPoint
return reinterpret_cast<PointFrompoint>(p);
}
#else
inline Foam::conformalVoronoiMesh::pointFromPoint
Foam::conformalVoronoiMesh::topoint
(