diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C index f9516c09a7615646cfc4b6a6ba4074108e9b3579..4e97af7a7a7583dea1535bae83092571ea12bbb0 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C @@ -178,8 +178,7 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment "(" "const Foam::point& pt" ") const" - ) - << "No secondary surface hit found in spoke search " + ) << "No secondary surface hit found in spoke search " << "using " << s << " spokes, try increasing alignmentSearchSpokes." << endl; diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H index a51e5a595f7451c0cdafb9dcb7022ed58453ea86..4ee0b2a42d6ab720f3a3b81d46fecf4119cb6140 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H @@ -124,8 +124,15 @@ public: private: + // Static data + static const scalar tolParallel; + static const scalar searchConeAngle; + + static const scalar searchAngleOppositeSurface; + + // Private data //- The time registry of the application @@ -420,6 +427,8 @@ private: //- Determine and insert point groups at the feature points void insertFeaturePoints(); + bool edgesShareNormal(const label e1, const label e2) const; + //- Create point groups at convex feature points void createConvexFeaturePoints ( @@ -470,6 +479,11 @@ private: //- Store the locations of all of the features to be conformed to void constructFeaturePointLocations(); + List<pointIndexHit> findSurfacePtLocationsNearFeaturePoint + ( + const Foam::point& featurePoint + ) const; + //- Reinsert stored feature point defining points void reinsertFeaturePoints(bool distribute = false); @@ -582,6 +596,7 @@ private: // intersect. bool clipLineToProc ( + const Foam::point& pt, Foam::point& a, Foam::point& b ) const; @@ -699,9 +714,9 @@ private: ) const; //- Check if a point is near any feature edge points. - bool pointIsNearFeatureEdge(const Foam::point& pt) const; + bool pointIsNearFeatureEdgeLocation(const Foam::point& pt) const; - bool pointIsNearFeatureEdge + bool pointIsNearFeatureEdgeLocation ( const Foam::point& pt, pointIndexHit& info diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C index 20ae246d279ebeccecc7c8f0ce77dfd3f1ec3b9f..925a19c099c057fa710b77bca13264838251740f 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C @@ -29,6 +29,12 @@ License using namespace Foam::vectorTools; +const Foam::scalar Foam::conformalVoronoiMesh::searchConeAngle + = Foam::cos(degToRad(30)); + +const Foam::scalar Foam::conformalVoronoiMesh::searchAngleOppositeSurface + = Foam::cos(degToRad(150)); + // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // void Foam::conformalVoronoiMesh::conformToSurface() @@ -277,79 +283,6 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation timeCheck("After initial conformation"); initialTotalHits = nSurfHits + nFeatEdHits; - - // Filter small edges at the boundary by inserting surface point pairs -// for -// ( -// Delaunay::Finite_cells_iterator cit = finite_cells_begin(); -// cit != finite_cells_end(); -// cit++ -// ) -// { -// const Foam::point dualVertex = topoint(dual(cit)); -// -// -// } -// for -// ( -// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin(); -// vit != finite_vertices_end(); -// vit++ -// ) -// { -// if (vit->nearBoundary()) -// { -// std::list<Facet> faces; -// incident_facets(vit, std::back_inserter(faces)); -// -// List<scalar> edgeLengthList(faces.size()); -// scalar totalLength = 0; -// label count = 0; -// -// for -// ( -// std::list<Facet>::iterator fit=faces.begin(); -// fit != faces.end(); -// ++fit -// ) -// { -// if -// ( -// is_infinite(fit->first) -// || is_infinite(fit->first->neighbor(fit->second)) -// ) -// { -// continue; -// } -// -// const Point dE0 = dual(fit->first); -// const Point dE1 = dual(fit->first->neighbor(fit->second)); -// -// const Segment s(dE0, dE1); -// -// const scalar length = Foam::sqrt(s.squared_length()); -// -// edgeLengthList[count++] = length; -// -// totalLength += length; -// -// //Info<< length << " / " << totalLength << endl; -// } -// -// const scalar averageLength = totalLength/edgeLengthList.size(); -// -// forAll(edgeLengthList, eI) -// { -// const scalar& el = edgeLengthList[eI]; -// -// if (el < 0.1*averageLength) -// { -// //Info<< "SMALL" << endl; -// } -// } -// } -// } - } // Remember which vertices were referred to each processor so only updates @@ -598,7 +531,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection Foam::point& a = dE0; Foam::point& b = dE1; - bool inProc = clipLineToProc(a, b); + bool inProc = clipLineToProc(topoint(vit->point()), a, b); // Check for the edge passing through a surface if @@ -642,7 +575,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections for ( - std::list<Facet>::iterator fit=facets.begin(); + std::list<Facet>::iterator fit = facets.begin(); fit != facets.end(); ++fit ) @@ -666,7 +599,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections if (Pstream::parRun()) { - bool inProc = clipLineToProc(dE0, dE1); + bool inProc = clipLineToProc(topoint(vit->point()), dE0, dE1); if (!inProc) { @@ -684,8 +617,6 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections if (infoIntersection.hit()) { - flagIntersection = true; - vectorField norm(1); allGeometry_[hitSurfaceIntersection].getNormal @@ -709,10 +640,10 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections pointIndexHit info; label hitSurface = -1; - geometryToConformTo_.findSurfaceNearestIntersection + geometryToConformTo_.findSurfaceNearest ( - vertex, - newPoint + SMALL*n, + newPoint, + surfaceSearchDistanceSqr(newPoint), info, hitSurface ); @@ -738,12 +669,12 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections = mag(p - info.hitPoint()); const scalar minSepDist - = cvMeshControls().removalDistCoeff() - *targetCellSize(p); + = sqr(cvMeshControls().removalDistCoeff() + *targetCellSize(p)); // Reject the point if it is too close to another // surface point. - // Could merge? + // Could merge the points? if (separationDistance < minSepDist) { rejectPoint = true; @@ -757,6 +688,7 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections // because another surface may be in the way. if (!rejectPoint && info.hit()) { + flagIntersection = true; infoList.append(info); hitSurfaceList.append(hitSurface); } @@ -769,72 +701,45 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections bool Foam::conformalVoronoiMesh::clipLineToProc ( + const Foam::point& pt, Foam::point& a, Foam::point& b ) const { - bool intersects = false; + bool inProc = false; + + pointIndexHit findAnyIntersection = decomposition_().findLine(a, b); - if (decomposition_().positionOnThisProcessor(a)) + if (!findAnyIntersection.hit()) { - // a is inside this processor + pointIndexHit info = decomposition_().findLine(a, pt); - if (decomposition_().positionOnThisProcessor(b)) + if (!info.hit()) { - // both a and b inside, no clip required - intersects = true; + inProc = true; } else { - // b is outside, clip b to bounding box. - - pointIndexHit info = decomposition_().findLine(b, a); - - if (info.hit()) - { - intersects = true; - b = info.hitPoint(); - } + inProc = false; } } else { - // a is outside this processor + pointIndexHit info = decomposition_().findLine(a, pt); - if (decomposition_().positionOnThisProcessor(b)) + if (!info.hit()) { - // b is inside this processor, clip a to processor - - pointIndexHit info = decomposition_().findLine(a, b); - - if (info.hit()) - { - intersects = true; - a = info.hitPoint(); - } + inProc = true; + b = findAnyIntersection.hitPoint(); } else { - // both a and b outside, but they can still intersect the processor - - pointIndexHit info = decomposition_().findLine(a, b); - - if (info.hit()) - { - intersects = true; - a = info.hitPoint(); - - info = decomposition_().findLine(b, a); - - if (info.hit()) - { - b = info.hitPoint(); - } - } + inProc = true; + a = findAnyIntersection.hitPoint(); } } - return intersects; + return inProc; } @@ -1877,7 +1782,7 @@ Foam::scalar Foam::conformalVoronoiMesh::angleBetweenSurfacePoints const vector nB = norm[0]; - return vectorTools::degAngleBetween(nA, nB); + return vectorTools::cosPhi(nA, nB); } @@ -1896,11 +1801,11 @@ bool Foam::conformalVoronoiMesh::nearSurfacePoint if (closeToSurfacePt) { - const scalar angle + const scalar cosAngle = angleBetweenSurfacePoints(pt, closePoint.hitPoint()); // @todo make this tolerance run-time selectable? - if (angle > 170.0) + if (cosAngle < searchAngleOppositeSurface) { pointIndexHit pCloseHit; label pCloseSurfaceHit = -1; @@ -1999,7 +1904,7 @@ bool Foam::conformalVoronoiMesh::appendToEdgeLocationTree } -bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdge +bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation ( const Foam::point& pt ) const @@ -2013,7 +1918,7 @@ bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdge } -bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdge +bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation ( const Foam::point& pt, pointIndexHit& info @@ -2066,7 +1971,7 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation pointIndexHit info; - bool closeToFeatureEdge = pointIsNearFeatureEdge(pt, info); + bool closeToFeatureEdge = pointIsNearFeatureEdgeLocation(pt, info); if (!closeToFeatureEdge) { @@ -2075,7 +1980,7 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation else { // Check if the edge location that the new edge location is near to - // might be on a different edge. If so, add it anyway. + // "might" be on a different edge. If so, add it anyway. pointIndexHit edgeHit; label featureHit = -1; @@ -2094,14 +1999,17 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation const vector lineBetweenPoints = pt - info.hitPoint(); - const scalar angle = degAngleBetween(edgeDir, lineBetweenPoints); + const scalar cosAngle = vectorTools::cosPhi(edgeDir, lineBetweenPoints); // Allow the point to be added if it is almost at right angles to the // other point. Also check it is not the same point. +// Info<< cosAngle<< " " +// << radToDeg(acos(cosAngle)) << " " +// << searchConeAngle << " " +// << radToDeg(acos(searchConeAngle)) << endl; if ( - angle < 100.0 - && angle > 80.0 + mag(cosAngle) < searchConeAngle && mag(lineBetweenPoints) > SMALL ) { @@ -2274,7 +2182,16 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits continue; } - if (nearFeaturePt(surfHitI.hitPoint())) + bool isNearFeaturePt = nearFeaturePt(surfHitI.hitPoint()); + + bool isNearSurfacePt = nearSurfacePoint + ( + surfHitI, + hitSurfaceI, + existingSurfacePtLocations + ); + + if (isNearFeaturePt || isNearSurfacePt) { keepSurfacePoint = false; @@ -2289,19 +2206,6 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits } } - if - ( - nearSurfacePoint - ( - surfHitI, - hitSurfaceI, - existingSurfacePtLocations - ) - ) - { - keepSurfacePoint = false; - } - List<List<pointIndexHit> > edHitsByFeature; labelList featuresHit; @@ -2332,7 +2236,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits { pointIndexHit& edHit = edHits[eHitI]; - if (!nearFeaturePt(edHit.hitPoint())) + if (!nearFeaturePt(edHit.hitPoint()) && keepSurfacePoint) { if ( @@ -2345,6 +2249,8 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits // this will prevent "pits" forming. keepSurfacePoint = false; + + // NEED TO REMOVE FROM THE SURFACE TREE... } if diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C index aa8e6722fe0bbc4f064d81fd00c0426a43bcc5e3..cb2a0c9340c3bd5c88c13c794c28f2c459400ff2 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C @@ -84,12 +84,12 @@ void Foam::conformalVoronoiMesh::createEdgePointGroup } case extendedFeatureEdgeMesh::OPEN: { - createOpenEdgePointGroup(feMesh, edHit, pts, indices, types); + //createOpenEdgePointGroup(feMesh, edHit, pts, indices, types); break; } case extendedFeatureEdgeMesh::MULTIPLE: { - createMultipleEdgePointGroup(feMesh, edHit, pts, indices, types); + //createMultipleEdgePointGroup(feMesh, edHit, pts, indices, types); break; } case extendedFeatureEdgeMesh::NONE: @@ -124,7 +124,6 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup { // The normals are nearly parallel, so this is too sharp a feature to // conform to. - return; } @@ -448,6 +447,40 @@ void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute) } +bool Foam::conformalVoronoiMesh::edgesShareNormal +( + const label e1, + const label e2 +) const +{ + const PtrList<extendedFeatureEdgeMesh>& feMeshes + ( + geometryToConformTo_.features() + ); + + forAll(feMeshes, i) + { + const extendedFeatureEdgeMesh& feMesh(feMeshes[i]); + + const vectorField& e1normals = feMesh.edgeNormals(e1); + const vectorField& e2normals = feMesh.edgeNormals(e2); + + forAll(e1normals, nI1) + { + forAll(e2normals, nI2) + { + if (degAngleBetween(e1normals[nI1], e2normals[nI2]) < 1) + { + return true; + } + } + } + } + + return false; +} + + void Foam::conformalVoronoiMesh::createConvexFeaturePoints ( DynamicList<Foam::point>& pts, @@ -479,6 +512,7 @@ void Foam::conformalVoronoiMesh::createConvexFeaturePoints } const vectorField& featPtNormals = feMesh.featurePointNormals(ptI); + const scalar ppDist = - pointPairDistance(featPt); vector cornerNormal = sum(featPtNormals); @@ -706,6 +740,61 @@ void Foam::conformalVoronoiMesh::createMixedFeaturePoints } } +// +//void Foam::conformalVoronoiMesh::createFeaturePoints +//( +// DynamicList<Foam::point>& pts, +// DynamicList<label>& indices, +// DynamicList<label>& types +//) +//{ +// const PtrList<extendedFeatureEdgeMesh>& feMeshes +// ( +// geometryToConformTo_.features() +// ); +// +// forAll(feMeshes, i) +// { +// const extendedFeatureEdgeMesh& feMesh(feMeshes[i]); +// +// for +// ( +// label ptI = feMesh.convexStart(); +// ptI < feMesh.mixedStart(); +// ++ptI +// ) +// { +// const Foam::point& featPt = feMesh.points()[ptI]; +// +// if (!positionOnThisProc(featPt)) +// { +// continue; +// } +// +// const scalar searchRadiusSqr = 5.0*targetCellSize(featPt); +// +// labelList indices = +// surfacePtLocationTreePtr_().findSphere(featPt, searchRadiusSqr); +// +// pointField nearestSurfacePoints(indices.size()); +// +// forAll(indices, pI) +// { +// nearestSurfacePoints[pI] = +// surfaceConformationVertices_[indices[pI]]; +// } +// +// forAll() +// +// // Now find the nearest points within the edge cones. +// +// // Calculate preliminary surface point locations +// +// +// } +// } +//} + void Foam::conformalVoronoiMesh::insertFeaturePoints() { @@ -721,6 +810,8 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints() createConcaveFeaturePoints(pts, indices, types); +// createFeaturePoints(pts, indices, types); + createMixedFeaturePoints(pts, indices, types); // Insert the created points, distributing to the appropriate processor @@ -798,3 +889,35 @@ void Foam::conformalVoronoiMesh::constructFeaturePointLocations() featurePointLocations_.transfer(ftPtLocs); } + + +Foam::List<Foam::pointIndexHit> +Foam::conformalVoronoiMesh::findSurfacePtLocationsNearFeaturePoint +( + const Foam::point& featurePoint +) const +{ + DynamicList<pointIndexHit> dynPointList; + + const scalar searchRadiusSqr = 3*targetCellSize(featurePoint); + + labelList surfacePtList = surfacePtLocationTreePtr_().findSphere + ( + featurePoint, + searchRadiusSqr + ); + + forAll(surfacePtList, elemI) + { + label index = surfacePtList[elemI]; + + const Foam::point& p + = surfacePtLocationTreePtr_().shapes().shapePoints()[index]; + + pointIndexHit nearHit(true, p, index); + + dynPointList.append(nearHit); + } + + return dynPointList.shrink(); +}