From 7b88ff0824bec6e0d999d52608f9624b61084c38 Mon Sep 17 00:00:00 2001 From: laurence Date: Fri, 11 Nov 2011 09:55:46 +0000 Subject: [PATCH 1/3] ENH: cv2DMesh and cvMesh: Addition of 2D mesher. cv2DMesh: Initial commit cvMesh: Constructors take Time argument instead of conformalVoronoiMesh --- .../cv2DMesh/CGALTriangulation2Ddefs.H | 82 ++ .../utilities/mesh/generation/cv2DMesh/CV2D.C | 1028 +++++++++++++++++ .../utilities/mesh/generation/cv2DMesh/CV2D.H | 443 +++++++ .../mesh/generation/cv2DMesh/CV2DI.H | 182 +++ .../mesh/generation/cv2DMesh/CV2DIO.C | 304 +++++ .../mesh/generation/cv2DMesh/Make/files | 15 + .../mesh/generation/cv2DMesh/Make/options | 35 + .../conformalVoronoi2DMesh/Make/files | 4 + .../conformalVoronoi2DMesh/Make/options | 14 + .../patchTo2DpolyMesh/patchTo2DpolyMesh.C | 334 ++++++ .../patchTo2DpolyMesh/patchTo2DpolyMesh.H | 169 +++ .../mesh/generation/cv2DMesh/controls.C | 82 ++ .../mesh/generation/cv2DMesh/controls.H | 164 +++ .../mesh/generation/cv2DMesh/cv2DMesh.C | 214 ++++ .../mesh/generation/cv2DMesh/indexedFace.H | 147 +++ .../mesh/generation/cv2DMesh/indexedVertex.H | 260 +++++ .../insertBoundaryConformPointPairs.C | 316 +++++ .../generation/cv2DMesh/insertFeaturePoints.C | 197 ++++ .../cv2DMesh/insertSurfaceNearPointPairs.C | 115 ++ .../cv2DMesh/insertSurfaceNearestPointPairs.C | 239 ++++ .../mesh/generation/cv2DMesh/querySurface.C | 237 ++++ .../mesh/generation/cv2DMesh/querySurface.H | 149 +++ .../generation/cv2DMesh/shortEdgeFilter2D.C | 502 ++++++++ .../generation/cv2DMesh/shortEdgeFilter2D.H | 137 +++ .../mesh/generation/cv2DMesh/tolerances.C | 68 ++ .../mesh/generation/cv2DMesh/tolerances.H | 138 +++ .../cellSizeControlSurfaces.C | 61 +- .../cellSizeControlSurfaces.H | 4 - .../cellSizeFunction/cellSizeFunction.C | 5 +- .../cellSizeFunction/cellSizeFunction.H | 8 +- .../linearDistance/linearDistance.C | 3 +- .../linearDistance/linearDistance.H | 1 - .../linearSpatial/linearSpatial.C | 3 +- .../linearSpatial/linearSpatial.H | 1 - .../surfaceOffsetLinearDistance.C | 3 +- .../surfaceOffsetLinearDistance.H | 1 - .../cellSizeFunction/uniform/uniform.C | 3 +- .../cellSizeFunction/uniform/uniform.H | 1 - .../uniformDistance/uniformDistance.C | 3 +- .../uniformDistance/uniformDistance.H | 1 - .../conformalVoronoiMesh.C | 4 +- .../conformationSurfaces.C | 18 +- .../conformationSurfaces.H | 8 +- 43 files changed, 5628 insertions(+), 75 deletions(-) create mode 100644 applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2Ddefs.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/CV2D.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/CV2D.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/CV2DI.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/CV2DIO.C create mode 100755 applications/utilities/mesh/generation/cv2DMesh/Make/files create mode 100755 applications/utilities/mesh/generation/cv2DMesh/Make/options create mode 100755 applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files create mode 100755 applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options create mode 100644 applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/controls.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/controls.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/indexedFace.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/querySurface.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/querySurface.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/tolerances.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/tolerances.H diff --git a/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2Ddefs.H b/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2Ddefs.H new file mode 100644 index 0000000000..9528f060d5 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2Ddefs.H @@ -0,0 +1,82 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2007-2010 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Typedefs + CGALTriangulation2Ddefs + +Description + CGAL data structures used for 2D Delaunay meshing. + + Define CGAL_INEXACT to use Exact_predicates_inexact_constructions kernel + otherwise the more robust but much less efficient + Exact_predicates_exact_constructions will be used. + + Define CGAL_HIERARCHY to use hierarchical Delaunay triangulation which is + faster but uses more memory than the standard Delaunay triangulation. + +\*---------------------------------------------------------------------------*/ + +#ifndef CGALTriangulation2Ddefs_H +#define CGALTriangulation2Ddefs_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "CGAL/Delaunay_triangulation_2.h" + +#include "indexedVertex.H" +#include "indexedFace.H" + +#ifdef CGAL_INEXACT + // Fast kernel using a double as the storage type but the triangulation + // may fail + #include "CGAL/Exact_predicates_inexact_constructions_kernel.h" + typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +#else + // Very robust but expensive kernel + #include "CGAL/Exact_predicates_exact_constructions_kernel.h" + typedef CGAL::Exact_predicates_exact_constructions_kernel K; +#endif + +typedef CGAL::indexedVertex Vb; +typedef CGAL::indexedFace Fb; + +#ifdef CGAL_HIERARCHY + // Data structures for hierarchical Delaunay triangulation which is more + // efficient but also uses more storage + #include "CGAL/Triangulation_hierarchy_2.h" + typedef CGAL::Triangulation_hierarchy_vertex_base_2 Vbh; + typedef CGAL::Triangulation_data_structure_2 Tds; + typedef CGAL::Delaunay_triangulation_2 DT; + typedef CGAL::Triangulation_hierarchy_2
Delaunay; +#else + // Data structures for standard Delaunay triangulation + typedef CGAL::Triangulation_data_structure_2 Tds; + typedef CGAL::Delaunay_triangulation_2 Delaunay; +#endif + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.C b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C new file mode 100644 index 0000000000..f6ab73cc3f --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C @@ -0,0 +1,1028 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2007-2011 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*----------------------------------------------------------------------------*/ + +#include "CV2D.H" +#include "Random.H" +#include "transform.H" +#include "IFstream.H" +#include "uint.H" +#include "ulong.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::CV2D::insertBoundingBox() +{ + Info<< "insertBoundingBox: creating bounding mesh" << endl; + scalar bigSpan = 10*tols_.span; + insertPoint(point2D(-bigSpan, -bigSpan), Vb::FAR_POINT); + insertPoint(point2D(-bigSpan, bigSpan), Vb::FAR_POINT); + insertPoint(point2D(bigSpan, -bigSpan), Vb::FAR_POINT); + insertPoint(point2D(bigSpan, bigSpan), Vb::FAR_POINT); +} + + +void Foam::CV2D::fast_restore_Delaunay(Vertex_handle vh) +{ + int i; + Face_handle f = vh->face(), next, start(f); + + do + { + i=f->index(vh); + if (!is_infinite(f)) + { + if (!internal_flip(f, cw(i))) external_flip(f, i); + if (f->neighbor(i) == start) start = f; + } + f = f->neighbor(cw(i)); + } while (f != start); +} + +void Foam::CV2D::external_flip(Face_handle& f, int i) +{ + Face_handle n = f->neighbor(i); + + if + ( + CGAL::ON_POSITIVE_SIDE + != side_of_oriented_circle(n, f->vertex(i)->point()) + ) return; + + flip(f, i); + i = n->index(f->vertex(i)); + external_flip(n, i); +} + +bool Foam::CV2D::internal_flip(Face_handle& f, int i) +{ + Face_handle n = f->neighbor(i); + + if + ( + CGAL::ON_POSITIVE_SIDE + != side_of_oriented_circle(n, f->vertex(i)->point()) + ) + { + return false; + } + + flip(f, i); + + return true; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::CV2D::CV2D +( + const Time& runTime, + const dictionary& cvMeshDict +) +: + Delaunay(), + runTime_(runTime), + rndGen_(64293*Pstream::myProcNo()), + allGeometry_ + ( + IOobject + ( + "cvSearchableSurfaces", + runTime_.constant(), + "triSurface", + runTime_, + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + cvMeshDict.subDict("geometry") + ), + qSurf_ + ( + runTime_, + rndGen_, + allGeometry_, + cvMeshDict.subDict("surfaceConformation") + ), + controls_(cvMeshDict), + cellSizeControl_ + ( + allGeometry_, + cvMeshDict.subDict("motionControl") + ), + tols_(cvMeshDict, controls_.minCellSize, qSurf_.globalBounds()), + z_ + ( + (1.0/3.0) + *(qSurf_.globalBounds().min().z() + qSurf_.globalBounds().max().z()) + ), + startOfInternalPoints_(0), + startOfSurfacePointPairs_(0), + startOfBoundaryConformPointPairs_(0) +{ + insertBoundingBox(); + insertFeaturePoints(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::CV2D::~CV2D() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::CV2D::insertPoints +( + const point2DField& points, + const scalar nearness +) +{ + Info<< "insertInitialPoints(const point2DField& points): "; + + startOfInternalPoints_ = number_of_vertices(); + label nVert = startOfInternalPoints_; + + // Add the points and index them + forAll(points, i) + { + const point2D& p = points[i]; + + if (qSurf_.wellInside(toPoint3D(p), nearness)) + { + insert(toPoint(p))->index() = nVert++; + } + else + { + Warning + << "Rejecting point " << p << " outside surface" << endl; + } + } + + Info<< nVert << " vertices inserted" << endl; + + if (controls_.writeInitialTriangulation) + { + // Checking validity of triangulation + assert(is_valid()); + + writeTriangles("initial_triangles.obj", true); + writeFaces("initial_faces.obj", true); + } +} + + +void Foam::CV2D::insertPoints(const fileName& pointFileName) +{ + IFstream pointsFile(pointFileName); + + if (pointsFile.good()) + { + insertPoints(point2DField(pointsFile), 0.5*controls_.minCellSize2); + } + else + { + FatalErrorIn("insertInitialPoints") + << "Could not open pointsFile " << pointFileName + << exit(FatalError); + } +} + + +void Foam::CV2D::insertGrid() +{ + Info<< "insertInitialGrid: "; + + startOfInternalPoints_ = number_of_vertices(); + label nVert = startOfInternalPoints_; + + scalar x0 = qSurf_.globalBounds().min().x(); + scalar xR = qSurf_.globalBounds().max().x() - x0; + int ni = int(xR/controls_.minCellSize) + 1; + scalar deltax = xR/ni; + + scalar y0 = qSurf_.globalBounds().min().y(); + scalar yR = qSurf_.globalBounds().max().y() - y0; + int nj = int(yR/controls_.minCellSize) + 1; + scalar deltay = yR/nj; + + Random rndGen(1321); + scalar pert = controls_.randomPurturbation*min(deltax, deltay); + + for (int i=0; iindex() = nVert++; + } + } + } + + Info<< nVert << " vertices inserted" << endl; + + if (controls_.writeInitialTriangulation) + { + // Checking validity of triangulation + assert(is_valid()); + + writeTriangles("initial_triangles.obj", true); + writeFaces("initial_faces.obj", true); + } +} + + +void Foam::CV2D::insertSurfacePointPairs() +{ + startOfSurfacePointPairs_ = number_of_vertices(); + + if (controls_.insertSurfaceNearestPointPairs) + { + insertSurfaceNearestPointPairs(); + } + + if (controls_.writeNearestTriangulation) + { + writeFaces("near_allFaces.obj", false); + writeFaces("near_faces.obj", true); + writeTriangles("near_triangles.obj", true); + } + + // Insertion of point-pairs for near-points may cause protrusions + // so insertBoundaryConformPointPairs must be executed last + if (controls_.insertSurfaceNearPointPairs) + { + insertSurfaceNearPointPairs(); + } + + startOfBoundaryConformPointPairs_ = number_of_vertices(); +} + + +void Foam::CV2D::boundaryConform() +{ + if (!controls_.insertSurfaceNearestPointPairs) + { + markNearBoundaryPoints(); + } + + // Mark all the faces as SAVE_CHANGED + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + fit++ + ) + { + fit->faceIndex() = Fb::SAVE_CHANGED; + } + + for (label iter=1; iter<=controls_.maxBoundaryConformingIter; iter++) + { + label nIntersections = insertBoundaryConformPointPairs + ( + "surfaceIntersections_" + Foam::name(iter) + ".obj" + ); + + if (nIntersections == 0) + { + break; + } + else + { + Info<< "BC iteration " << iter << ": " + << nIntersections << " point-pairs inserted" << endl; + } + + // Any faces changed by insertBoundaryConformPointPairs will now + // be marked CHANGED, mark those as SAVE_CHANGED and those that + // remained SAVE_CHANGED as UNCHANGED + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + fit++ + ) + { + if (fit->faceIndex() == Fb::SAVE_CHANGED) + { + fit->faceIndex() = Fb::UNCHANGED; + } + else if (fit->faceIndex() == Fb::CHANGED) + { + fit->faceIndex() = Fb::SAVE_CHANGED; + } + } + } + + Info<< nl; +} + + +void Foam::CV2D::removeSurfacePointPairs() +{ + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->index() >= startOfSurfacePointPairs_) + { + remove(vit); + } + } +} + + +void Foam::CV2D::newPoints(const scalar relaxation) +{ + Info<< "newPointsFromVertices: "; + + Field dualVertices(number_of_faces()); + + label dualVerti = 0; + + // Find the dual point of each tetrahedron and assign it an index. + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + ++fit + ) + { + fit->faceIndex() = -1; + + if + ( + fit->vertex(0)->internalOrBoundaryPoint() + || fit->vertex(1)->internalOrBoundaryPoint() + || fit->vertex(2)->internalOrBoundaryPoint() + ) + { + fit->faceIndex() = dualVerti; + + dualVertices[dualVerti] = toPoint2D(circumcenter(fit)); + + dualVerti++; + } + } + + dualVertices.setSize(dualVerti); + + Field displacementAccumulator + ( + startOfSurfacePointPairs_, + vector2D::zero + ); + + // Calculate target size and alignment for vertices + scalarField sizes + ( + number_of_vertices(), + controls_.minCellSize + ); + + Field alignments + ( + number_of_vertices(), + vector2D(1, 0) + ); + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->internalOrBoundaryPoint()) + { + point2D vert = toPoint2D(vit->point()); + + // alignment and size determination + pointIndexHit pHit; + label hitSurface = -1; + + qSurf_.findSurfaceNearest + ( + toPoint3D(vert), + tols_.span2, + pHit, + hitSurface + ); + + if (pHit.hit()) + { + + vectorField norm(1); + allGeometry_[hitSurface].getNormal + ( + List(1, pHit), + norm + ); + + alignments[vit->index()] = toPoint2D(norm[0]); + + scalar surfDist = mag(toPoint3D(vert) - pHit.hitPoint()); + + /*if (surfDist < 0.2) + { + sizes[vit->index()] *= 0.4; + }*/ + + if (surfDist < 0.2) + { + sizes[vit->index()] *= (1 - 0.1)*surfDist/0.2 + 0.1; + } + } + + // if (vert.x() > 0) + // { + // sizes[vit->index()] *= 0.5; + // } + } + } + + // Info<< "Calculated alignments" << endl; + + scalar cosAlignmentAcceptanceAngle = 0.68; + + // Upper and lower edge length ratios for weight + scalar u = 1.0; + scalar l = 0.7; + + PackedBoolList pointToBeRetained(startOfSurfacePointPairs_, true); + + DynamicList pointsToInsert; + + for + ( + Triangulation::Finite_edges_iterator eit = finite_edges_begin(); + eit != finite_edges_end(); + eit++ + ) + { + Vertex_handle vA = eit->first->vertex(cw(eit->second)); + Vertex_handle vB = eit->first->vertex(ccw(eit->second)); + + if (!vA->internalOrBoundaryPoint() || !vB->internalOrBoundaryPoint()) + { + continue; + } + + const point2D& dualV1 = dualVertices[eit->first->faceIndex()]; + const point2D& dualV2 = + dualVertices[eit->first->neighbor(eit->second)->faceIndex()]; + + scalar dualEdgeLength = mag(dualV1 - dualV2); + + point2D dVA = toPoint2D(vA->point()); + point2D dVB = toPoint2D(vB->point()); + + Field alignmentDirsA(2); + + alignmentDirsA[0] = alignments[vA->index()]; + alignmentDirsA[1] = vector2D + ( + -alignmentDirsA[0].y(), + alignmentDirsA[0].x() + ); + + //Info<< "DirsA " << alignmentDirsA[0] << " " << alignmentDirsA[1] << endl; + + Field alignmentDirsB(2); + + alignmentDirsB[0] = alignments[vB->index()]; + alignmentDirsB[1] = vector2D + ( + -alignmentDirsB[0].y(), + alignmentDirsB[0].x() + ); + + //Info<< "DirsB " << alignmentDirsB[0] << " " << alignmentDirsB[1] << endl; + + Field alignmentDirs(2); + + forAll(alignmentDirsA, aA) + { + const vector2D& a(alignmentDirsA[aA]); + + scalar maxDotProduct = 0.0; + + forAll(alignmentDirsB, aB) + { + const vector2D& b(alignmentDirsB[aB]); + + scalar dotProduct = a & b; + + //Info<< aA << " " << aB << " " << dotProduct << endl; + + if (mag(dotProduct) > maxDotProduct) + { + maxDotProduct = mag(dotProduct); + + alignmentDirs[aA] = a + sign(dotProduct)*b; + + alignmentDirs[aA] /= mag(alignmentDirs[aA]); + } + } + } + + //Info<< "Dirs " << alignmentDirs[0] << " " << alignmentDirs[1] << endl; + + + vector2D rAB = dVA - dVB; + + scalar rABMag = mag(rAB); + + forAll(alignmentDirs, aD) + { + vector2D& alignmentDir = alignmentDirs[aD]; + + if ((rAB & alignmentDir) < 0) + { + // swap the direction of the alignment so that has the + // same sense as rAB + alignmentDir *= -1; + } + + scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir); + + if (alignmentDotProd > cosAlignmentAcceptanceAngle) + { + scalar targetFaceSize = + 0.5*(sizes[vA->index()] + sizes[vB->index()]); + + // Test for changing aspect ratio on second alignment (first + // alignment is neartest surface normal) + // if (aD == 1) + // { + // targetFaceSize *= 2.0; + // } + + alignmentDir *= 0.5*targetFaceSize; + + vector2D delta = alignmentDir - 0.5*rAB; + + if (dualEdgeLength < 0.7*targetFaceSize) + { + delta *= 0; + } + else if (dualEdgeLength < targetFaceSize) + { + delta *= + ( + dualEdgeLength + /(targetFaceSize*(u - l)) + - 1/((u/l) - 1) + ); + } + + if + ( + vA->internalPoint() + && vB->internalPoint() + && rABMag > 1.75*targetFaceSize + && dualEdgeLength > 0.05*targetFaceSize + && alignmentDotProd > 0.93 + ) + { + // Point insertion + pointsToInsert.append(0.5*(dVA + dVB)); + } + else if + ( + (vA->internalPoint() || vB->internalPoint()) + && rABMag < 0.65*targetFaceSize + ) + { + // Point removal + + // Only insert a point at the midpoint of the short edge + // if neither attached point has already been identified + // to be removed. + if + ( + pointToBeRetained[vA->index()] == true + && pointToBeRetained[vB->index()] == true + ) + { + pointsToInsert.append(0.5*(dVA + dVB)); + } + + if (vA->internalPoint()) + { + pointToBeRetained[vA->index()] = false; + } + + if (vB->internalPoint()) + { + pointToBeRetained[vB->index()] = false; + } + } + else + { + if (vA->internalPoint()) + { + displacementAccumulator[vA->index()] += delta; + } + + if (vB->internalPoint()) + { + displacementAccumulator[vB->index()] += -delta; + } + } + } + } + } + + vector2D totalDisp = sum(displacementAccumulator); + scalar totalDist = sum(mag(displacementAccumulator)); + + // Relax the calculated displacement + displacementAccumulator *= relaxation; + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->internalPoint()) + { + if (!pointToBeRetained[vit->index()]) + { + remove(vit); + } + else + { + movePoint + ( + vit, + vit->point() + + K::Vector_2 + ( + displacementAccumulator[vit->index()].x(), + displacementAccumulator[vit->index()].y() + ) + ); + } + } + } + + removeSurfacePointPairs(); + + // Re-index internal points + + label nVert = startOfInternalPoints_; + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->internalPoint()) + { + vit->index() = nVert++; + } + } + + // Insert new points + + Info<< "Inserting " << pointsToInsert.size() << " new points" << endl; + + forAll(pointsToInsert, i) + { + insertPoint(pointsToInsert[i], Vb::INTERNAL_POINT); + } + + Info<< " Total displacement = " << totalDisp << nl + << " Total distance = " << totalDist << nl + << " Points added = " << pointsToInsert.size() + << endl; + + insertSurfacePointPairs(); + boundaryConform(); + + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Old Method + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + // for + // ( + // Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + // vit != finite_vertices_end(); + // ++vit + // ) + // { + // if (vit->internalPoint()) + // { + // // Current dual-cell defining vertex ("centre") + // point2DFromPoint defVert0 = toPoint2D(vit->point()); + + // Triangulation::Edge_circulator ec = incident_edges(vit); + // Triangulation::Edge_circulator ecStart = ec; + + // // Circulate around the edges to find the first which is not + // // infinite + // do + // { + // if (!is_infinite(ec)) break; + // } while (++ec != ecStart); + + // // Store the start-end of the first non-infinte edge + // point2D de0 = toPoint2D(circumcenter(ec->first)); + + // // Keep track of the maximum edge length^2 + // scalar maxEdgeLen2 = 0.0; + + // // Keep track of the index of the longest edge + // label edgecd0i = -1; + + // // Edge counter + // label edgei = 0; + + // do + // { + // if (!is_infinite(ec)) + // { + // // Get the end of the current edge + // point2D de1 = toPoint2D + // ( + // circumcenter(ec->first->neighbor(ec->second)) + // ); + + // // Store the current edge vector + // edges[edgei] = de1 - de0; + + // // Store the edge mid-point in the vertices array + // vertices[edgei] = 0.5*(de1 + de0); + + // // Move the current edge end into the edge start for the + // // next iteration + // de0 = de1; + + // // Keep track of the longest edge + + // scalar edgeLen2 = magSqr(edges[edgei]); + + // if (edgeLen2 > maxEdgeLen2) + // { + // maxEdgeLen2 = edgeLen2; + // edgecd0i = edgei; + // } + + // edgei++; + // } + // } while (++ec != ecStart); + + // // Initialise cd0 such that the mesh will align + // // in in the x-y directions + // vector2D cd0(1, 0); + + // if (controls_.relaxOrientation) + // { + // // Get the longest edge from the array and use as the primary + // // direction of the coordinate system of the "square" cell + // cd0 = edges[edgecd0i]; + // } + + // if (controls_.nearWallAlignedDist > 0) + // { + // pointIndexHit pHit = qSurf_.tree().findNearest + // ( + // toPoint3D(defVert0), + // controls_.nearWallAlignedDist2 + // ); + + // if (pHit.hit()) + // { + // cd0 = toPoint2D(faceNormals[pHit.index()]); + // } + // } + + // // Rotate by 45deg needed to create an averaging procedure which + // // encourages the cells to be square + // cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x()); + + // // Normalise the primary coordinate direction + // cd0 /= mag(cd0); + + // // Calculate the orthogonal coordinate direction + // vector2D cd1(-cd0.y(), cd0.x()); + + + // // Restart the circulator + // ec = ecStart; + + // // ... and the counter + // edgei = 0; + + // // Initialise the displacement for the centre and sum-weights + // vector2D disp = vector2D::zero; + // scalar sumw = 0; + + // do + // { + // if (!is_infinite(ec)) + // { + // // Pick up the current edge + // const vector2D& ei = edges[edgei]; + + // // Calculate the centre to edge-centre vector + // vector2D deltai = vertices[edgei] - defVert0; + + // // Set the weight for this edge contribution + // scalar w = 1; + + // if (controls_.squares) + // { + // w = magSqr(deltai.x()*ei.y() - deltai.y()*ei.x()); + // // alternative weights + // //w = mag(deltai.x()*ei.y() - deltai.y()*ei.x()); + // //w = magSqr(ei)*mag(deltai); + + // // Use the following for an ~square mesh + // // Find the coordinate contributions for this edge delta + // scalar cd0deltai = cd0 & deltai; + // scalar cd1deltai = cd1 & deltai; + + // // Create a "square" displacement + // if (mag(cd0deltai) > mag(cd1deltai)) + // { + // disp += (w*cd0deltai)*cd0; + // } + // else + // { + // disp += (w*cd1deltai)*cd1; + // } + // } + // else + // { + // // Use this for a hexagon/pentagon mesh + // disp += w*deltai; + // } + + // // Sum the weights + // sumw += w; + // } + // else + // { + // FatalErrorIn("CV2D::newPoints() const") + // << "Infinite triangle found in internal mesh" + // << exit(FatalError); + // } + + // edgei++; + + // } while (++ec != ecStart); + + // // Calculate the average displacement + // disp /= sumw; + // totalDisp += disp; + // totalDist += mag(disp); + + // // Move the point by a fraction of the average displacement + // movePoint(vit, defVert0 + relaxation*disp); + // } + // } + + // Info << "\nTotal displacement = " << totalDisp + // << " total distance = " << totalDist << endl; +} + + +// void Foam::CV2D::moveInternalPoints(const point2DField& newPoints) +// { +// label pointI = 0; + +// for +// ( +// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); +// vit != finite_vertices_end(); +// ++vit +// ) +// { +// if (vit->internalPoint()) +// { +// movePoint(vit, newPoints[pointI++]); +// } +// } +// } + + +void Foam::CV2D::extractPatches +( + wordList& patchNames, + labelList& patchSizes, + EdgeMap