diff --git a/applications/utilities/mesh/generation/Allwmake b/applications/utilities/mesh/generation/Allwmake index ef28f9aaf5..48ca530c96 100755 --- a/applications/utilities/mesh/generation/Allwmake +++ b/applications/utilities/mesh/generation/Allwmake @@ -4,12 +4,15 @@ set -x wmake blockMesh wmake all extrude -wmake extrude2DMesh + +extrude2DMesh/Allwmake + wmake snappyHexMesh if [ -d "$CGAL_ARCH_PATH" ] then - cd cvMesh && ./Allwmake + cvMesh/Allwmake + cv2DMesh/Allwmake fi # ----------------------------------------------------------------- end-of-file diff --git a/applications/utilities/mesh/generation/cv2DMesh/Allwclean b/applications/utilities/mesh/generation/cv2DMesh/Allwclean new file mode 100755 index 0000000000..d0ae53e415 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/Allwclean @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wclean libso conformalVoronoi2DMesh +wclean + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/utilities/mesh/generation/cv2DMesh/Allwmake b/applications/utilities/mesh/generation/cv2DMesh/Allwmake new file mode 100755 index 0000000000..5486849957 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/Allwmake @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake libso conformalVoronoi2DMesh +wmake + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2Ddefs.H b/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2Ddefs.H new file mode 100644 index 0000000000..1485226d98 --- /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-2011 OpenFOAM Foundation + \\/ 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..aa45a88b3e --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C @@ -0,0 +1,1050 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2007-2011 OpenFOAM Foundation + \\/ 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" + +namespace Foam +{ + +defineTypeNameAndDebug(CV2D, 0); + +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::CV2D::insertBoundingBox() +{ + Info<< "insertBoundingBox: creating bounding mesh" << endl; + scalar bigSpan = 10*meshControls().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, qSurf_.globalBounds()), + cellSizeControl_ + ( + allGeometry_, + cvMeshDict.subDict("motionControl") + ), + z_ + ( + point + ( + cvMeshDict.subDict("surfaceConformation").lookup("locationInMesh") + ).z() + ), + startOfInternalPoints_(0), + startOfSurfacePointPairs_(0), + startOfBoundaryConformPointPairs_(0), + featurePoints_() +{ + Info<< meshControls() << endl; + + 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 (meshControls().objOutput()) + { + // 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*meshControls().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/meshControls().minCellSize()) + 1; + scalar deltax = xR/ni; + + scalar y0 = qSurf_.globalBounds().min().y(); + scalar yR = qSurf_.globalBounds().max().y() - y0; + int nj = int(yR/meshControls().minCellSize()) + 1; + scalar deltay = yR/nj; + + Random rndGen(1321); + scalar pert = meshControls().randomPerturbation()*min(deltax, deltay); + + for (int i=0; iindex() = nVert++; + } + } + } + + Info<< nVert << " vertices inserted" << endl; + + if (meshControls().objOutput()) + { + // 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 (meshControls().insertSurfaceNearestPointPairs()) + { + insertSurfaceNearestPointPairs(); + } + + write("nearest"); + + // Insertion of point-pairs for near-points may cause protrusions + // so insertBoundaryConformPointPairs must be executed last + if (meshControls().insertSurfaceNearPointPairs()) + { + insertSurfaceNearPointPairs(); + } + + startOfBoundaryConformPointPairs_ = number_of_vertices(); +} + + +void Foam::CV2D::boundaryConform() +{ + if (!meshControls().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<=meshControls().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; + + write("boundary"); +} + + +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(), + meshControls().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), + meshControls().span2(), + pHit, + hitSurface + ); + + if (pHit.hit()) + { + vectorField norm(1); + allGeometry_[hitSurface].getNormal + ( + List(1, pHit), + norm + ); + + alignments[vit->index()] = toPoint2D(norm[0]); + + sizes[vit->index()] = + cellSizeControl_.cellSize(toPoint3D(vit->point())); + } + } + } + + // 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); + + std::list 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() + ); + + Field alignmentDirsB(2); + + alignmentDirsB[0] = alignments[vB->index()]; + alignmentDirsB[1] = vector2D + ( + -alignmentDirsB[0].y(), + alignmentDirsB[0].x() + ); + + 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; + + if (mag(dotProduct) > maxDotProduct) + { + maxDotProduct = mag(dotProduct); + + alignmentDirs[aA] = a + sign(dotProduct)*b; + + alignmentDirs[aA] /= mag(alignmentDirs[aA]); + } + } + } + + 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.push_back(toPoint(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.push_back(toPoint(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; + + label numberOfNewPoints = pointsToInsert.size(); + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->internalPoint()) + { + if (pointToBeRetained[vit->index()]) + { + pointsToInsert.push_front + ( + toPoint + ( + toPoint2D(vit->point()) + + displacementAccumulator[vit->index()] + ) + ); + } + } + } + + // Clear the triangulation and reinsert the bounding box and feature points. + // This is faster than removing and moving points. + this->clear(); + + insertBoundingBox(); + + reinsertFeaturePoints(); + + startOfInternalPoints_ = number_of_vertices(); + + label nVert = startOfInternalPoints_; + + Info<< "Inserting " << numberOfNewPoints << " new points" << endl; + + // Use the range insert as it is faster than individually inserting points. + insert(pointsToInsert.begin(), pointsToInsert.end()); + + for + ( + Delaunay::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if + ( + vit->type() == Vb::INTERNAL_POINT + && vit->index() == Vb::INTERNAL_POINT + ) + { + vit->index() = nVert++; + } + } + + Info<< " Total displacement = " << totalDisp << nl + << " Total distance = " << totalDist << nl + << " Points added = " << pointsToInsert.size() + << endl; + + write("internal"); + + 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 (meshControls().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 (meshControls().nearWallAlignedDist() > 0) +// { +// pointIndexHit pHit = qSurf_.tree().findNearest +// ( +// toPoint3D(defVert0), +// meshControls().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 (meshControls().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