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 (