From e1f00987f04d009c0772af4872d223783edc1dfb Mon Sep 17 00:00:00 2001 From: graham Date: Wed, 19 Jan 2011 10:04:10 +0000 Subject: [PATCH] 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). --- .../cellSizeControlSurfaces.C | 230 +++++++++++------- .../cellSizeControlSurfaces.H | 11 +- .../conformalVoronoiMesh.C | 28 +-- .../conformalVoronoiMesh.H | 12 +- .../conformalVoronoiMeshConformToSurface.C | 190 ++++++++++++++- .../conformalVoronoiMeshI.H | 20 +- 6 files changed, 375 insertions(+), 116 deletions(-) diff --git a/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C b/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C index 6ad77155d6..19169d04f5 100644 --- a/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C +++ b/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C @@ -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) + cvMesh_.geometryToConformTo().findSurfaceNearest + ( + pt, + cvMesh_.geometryToConformTo().spanMagSqr(), + surfHit, + hitSurface + ); + + if (!surfHit.hit()) { - Info<< "size function " - << allGeometry_.names()[surfaces_[i]] - << " priority " << cSF.priority() - << endl; + 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; } - - if (cSF.priority() < previousPriority) + else { - return minSize; - } - - scalar sizeI; - - if (cSF.cellSize(pt, sizeI, isSurfacePoint)) - { - if (cSF.priority() == previousPriority) - { - if (sizeI < minSize) - { - minSize = sizeI; - } - } - 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; } diff --git a/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H b/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H index a2ddfb8c9e..f12547e7e2 100644 --- a/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H +++ b/src/mesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H @@ -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&); diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C index 692b5baef6..2829cf62f6 100644 --- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C +++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C @@ -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()] + ) ); } } diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H index 0c3a124f2f..c8647b4bd3 100644 --- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H +++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H @@ -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 diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C index 438d861453..8c160a031a 100644 --- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C +++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C @@ -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 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::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 facets; + incident_facets(vit, std::back_inserter(facets)); + + Foam::point vert(topoint(vit->point())); + + scalar minIncursionDistance = -maxSurfaceProtrusion(vert); + + for + ( + std::list::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(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, diff --git a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H index e3d92c06fc..ca0c67a21c 100644 --- a/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H +++ b/src/mesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H @@ -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(p); } - #else + inline Foam::conformalVoronoiMesh::pointFromPoint Foam::conformalVoronoiMesh::topoint (