From 7b88ff0824bec6e0d999d52608f9624b61084c38 Mon Sep 17 00:00:00 2001 From: laurence <laurence> Date: Fri, 11 Nov 2011 09:55:46 +0000 Subject: [PATCH] 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 | 67 +- .../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, 5631 insertions(+), 78 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 00000000000..9528f060d57 --- /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 <http://www.gnu.org/licenses/>. + +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<K> Vb; +typedef CGAL::indexedFace<K> 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<Vb> Vbh; + typedef CGAL::Triangulation_data_structure_2<Vbh, Fb> Tds; + typedef CGAL::Delaunay_triangulation_2<K, Tds> DT; + typedef CGAL::Triangulation_hierarchy_2<DT> Delaunay; +#else + // Data structures for standard Delaunay triangulation + typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds; + typedef CGAL::Delaunay_triangulation_2<K, Tds> 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 00000000000..f6ab73cc3f8 --- /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 <http://www.gnu.org/licenses/>. + +\*----------------------------------------------------------------------------*/ + +#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; i<ni; i++) + { + for (int j=0; j<nj; j++) + { + point p(x0 + i*deltax, y0 + j*deltay, 0); + + if (controls_.randomiseInitialGrid) + { + p.x() += pert*(rndGen.scalar01() - 0.5); + p.y() += pert*(rndGen.scalar01() - 0.5); + } + + if (qSurf_.wellInside(p, 0.5*controls_.minCellSize2)) + { + insert(Point(p.x(), p.y()))->index() = 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<point2D> 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<vector2D> displacementAccumulator + ( + startOfSurfacePointPairs_, + vector2D::zero + ); + + // Calculate target size and alignment for vertices + scalarField sizes + ( + number_of_vertices(), + controls_.minCellSize + ); + + Field<vector2D> 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<pointIndexHit>(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<point2D> 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<vector2D> alignmentDirsA(2); + + alignmentDirsA[0] = alignments[vA->index()]; + alignmentDirsA[1] = vector2D + ( + -alignmentDirsA[0].y(), + alignmentDirsA[0].x() + ); + + //Info<< "DirsA " << alignmentDirsA[0] << " " << alignmentDirsA[1] << endl; + + Field<vector2D> alignmentDirsB(2); + + alignmentDirsB[0] = alignments[vB->index()]; + alignmentDirsB[1] = vector2D + ( + -alignmentDirsB[0].y(), + alignmentDirsB[0].x() + ); + + //Info<< "DirsB " << alignmentDirsB[0] << " " << alignmentDirsB[1] << endl; + + Field<vector2D> 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<label>& mapEdgesRegion +) const +{ + label nPatches = qSurf_.patchNames().size() + 1; + label defaultPatchIndex = qSurf_.patchNames().size(); + + patchNames.setSize(nPatches); + patchSizes.setSize(nPatches, 0); + mapEdgesRegion.clear(); + + const wordList& existingPatches = qSurf_.patchNames(); + + forAll(existingPatches, sP) + { + patchNames[sP] = existingPatches[sP]; + } + + patchNames[defaultPatchIndex] = "CV2D_default_patch"; + + for + ( + Triangulation::Finite_edges_iterator eit = finite_edges_begin(); + eit != finite_edges_end(); + ++eit + ) + { + Face_handle fOwner = eit->first; + Face_handle fNeighbor = fOwner->neighbor(eit->second); + + Vertex_handle vA = fOwner->vertex(cw(eit->second)); + Vertex_handle vB = fOwner->vertex(ccw(eit->second)); + + if + ( + (vA->internalOrBoundaryPoint() && !vB->internalOrBoundaryPoint()) + || (vB->internalOrBoundaryPoint() && !vA->internalOrBoundaryPoint()) + ) + { + point ptA = toPoint3D(vA->point()); + point ptB = toPoint3D(vB->point()); + + label patchIndex = qSurf_.findPatch(ptA, ptB); + + if (patchIndex == -1) + { + patchIndex = defaultPatchIndex; + + WarningIn("Foam::CV2D::extractPatches") + << "Dual face found that is not on a surface " + << "patch. Adding to CV2D_default_patch." + << endl; + } + + edge e(fOwner->faceIndex(), fNeighbor->faceIndex()); + patchSizes[patchIndex]++; + mapEdgesRegion.insert(e, patchIndex); + } + } +} + + +void Foam::CV2D::write() const +{ + if (controls_.writeFinalTriangulation) + { + writeFaces("allFaces.obj", false); + writeFaces("faces.obj", true); + writeTriangles("allTriangles.obj", false); + writeTriangles("triangles.obj", true); + writePatch("patch.pch"); + } +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.H b/applications/utilities/mesh/generation/cv2DMesh/CV2D.H new file mode 100644 index 00000000000..e295300b128 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/CV2D.H @@ -0,0 +1,443 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +Class + CV2D + +Description + Conformal-Voronoi 2D automatic mesher with grid or read initial points + and point position relaxation with optional "squarification". + + There are a substantial number of options to this mesher read from + CV2DMesherDict file e.g.: + + // Min cell size used in tolerances when inserting points for + // boundary conforming. + // Also used to as the grid spacing usind in insertGrid. + minCellSize 0.05; + + // Feature angle used to inser feature points + // 0 = all features, 180 = no features + featureAngle 45; + + // Maximum quadrant angle allowed at a concave corner before + // additional "mitering" lines are added + maxQuadAngle 110; + + // Should the mesh be square-dominated or of unbiased hexagons + squares yes; + + // Near-wall region where cells are aligned with the wall specified as a + // number of cell layers + nearWallAlignedDist 3; + + // Chose if the cell orientation should relax during the iterations + // or remain fixed to the x-y directions + relaxOrientation no; + + // Insert near-boundary point mirror or point-pairs + insertSurfaceNearestPointPairs yes; + + // Mirror near-boundary points rather than insert point-pairs + mirrorPoints no; + + // Insert point-pairs vor dual-cell vertices very near the surface + insertSurfaceNearPointPairs yes; + + // Choose if to randomise the initial grid created by insertGrid. + randomiseInitialGrid yes; + + // Perturbation fraction, 1 = cell-size. + randomPurturbation 0.1; + + // Number of relaxation iterations. + nIterations 5; + + // Relaxation factor at the start of the iteration sequence. + // 0.5 is a sensible maximum and < 0.2 converges better. + relaxationFactorStart 0.8; + + // Relaxation factor at the end of the iteration sequence. + // Should be <= relaxationFactorStart + relaxationFactorEnd 0; + + writeInitialTriangulation no; + writeFeatureTriangulation no; + writeNearestTriangulation no; + writeInsertedPointPairs no; + writeFinalTriangulation yes; + + // Maximum number of iterations used in boundaryConform. + maxBoundaryConformingIter 5; + + minEdgeLenCoeff 0.5; + maxNotchLenCoeff 0.3; + minNearPointDistCoeff 0.25; + ppDistCoeff 0.05; + +SourceFiles + CGALTriangulation2Ddefs.H + indexedVertex.H + indexedFace.H + CV2DI.H + CV2D.C + CV2DIO.C + tolerances.C + controls.C + insertFeaturePoints.C + insertSurfaceNearestPointPairs.C + insertSurfaceNearPointPairs.C + insertBoundaryConformPointPairs.C + +\*---------------------------------------------------------------------------*/ + +#ifndef CV2D_H +#define CV2D_H + +#define CGAL_INEXACT +#define CGAL_HIERARCHY + +#include "CGALTriangulation2Ddefs.H" + +#include "Time.H" +#include "point2DFieldFwd.H" +#include "dictionary.H" +#include "Switch.H" +#include "PackedBoolList.H" +#include "EdgeMap.H" +#include "controls.H" +#include "tolerances.H" +#include "meshTools.H" +#include "triSurface.H" +#include "searchableSurfaces.H" +#include "conformationSurfaces.H" +#include "cellSizeControlSurfaces.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class CV2D Declaration +\*---------------------------------------------------------------------------*/ + +class CV2D +: + public Delaunay +{ + +private: + + // Private data + + //- The time registry of the application + const Time& runTime_; + + mutable Random rndGen_; + + //- The surface to mesh + //const querySurface& qSurf_; + //- All geometry of the meshing process, including surfaces to be + // conformed to and those to be used for refinement + searchableSurfaces allGeometry_; + + conformationSurfaces qSurf_; + + //- Meshing controls + controls controls_; + + //- The cell size control object + cellSizeControlSurfaces cellSizeControl_; + + //- Meshing tolerances + tolerances tols_; + + //- z-level + scalar z_; + + //- Keep track of the start of the internal points + label startOfInternalPoints_; + + //- Keep track of the start of the surface point-pairs + label startOfSurfacePointPairs_; + + //- Keep track of the boundary conform point-pairs + // stored after the insertion of the surface point-pairs in case + // the boundary conform function is called more than once without + // removing and insertin the surface point-pairs + label startOfBoundaryConformPointPairs_; + + //- Temporary storage for a dual-cell + static const label maxNvert = 20; + mutable point2D vertices[maxNvert+1]; + mutable vector2D edges[maxNvert+1]; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + CV2D(const CV2D&); + + //- Disallow default bitwise assignment + void operator=(const CV2D&); + + + //- Insert point and return it's index + inline label insertPoint + ( + const point2D& pt, + const label type + ); + + inline bool insertMirrorPoint + ( + const point2D& nearSurfPt, + const point2D& surfPt + ); + + //- Insert a point-pair at a distance ppDist either side of + // surface point point surfPt in the direction n + inline void insertPointPair + ( + const scalar mirrorDist, + const point2D& surfPt, + const vector2D& n + ); + + //- Create the initial mesh from the bounding-box + void insertBoundingBox(); + + //- Insert point groups at the feature points. + void insertFeaturePoints(); + + //- Insert point-pairs at the given set of points using the surface + // normals corresponding to the given set of surface triangles + // and write the inserted point locations to the given file. + void insertPointPairs + ( + const DynamicList<point2D>& nearSurfacePoints, + const DynamicList<point2D>& surfacePoints, + const DynamicList<label>& surfaceTris, + const fileName fName + ); + + //- Check to see if dual cell specified by given vertex iterator + // intersects the boundary and hence reqires a point-pair. + bool dualCellSurfaceIntersection + ( + const Triangulation::Finite_vertices_iterator& vit + ) const; + + //- Insert point-pairs at the nearest points on the surface to the + // control vertex of dual-cells which intersect the boundary in order + // to provide a boundary-layer mesh. + // NB: This is not guaranteed to close the boundary + void insertSurfaceNearestPointPairs(); + + //- Insert point-pairs at small dual-cell edges on the surface in order + // to improve the boundary-layer mesh generated by + // insertSurfaceNearestPointPairs. + void insertSurfaceNearPointPairs(); + + //- Insert point-pair and correcting the Finite_vertices_iterator + // to account for the additional vertices + void insertPointPair + ( + Triangulation::Finite_vertices_iterator& vit, + const point2D& p, + const label trii + ); + + //- Insert point-pair at the best intersection point between the lines + // from the dual-cell real centroid and it's vertices and the surface. + bool insertPointPairAtIntersection + ( + Triangulation::Finite_vertices_iterator& vit, + const point2D& defVert, + const point2D vertices[], + const scalar maxProtSize + ); + + //- Insert point-pairs corresponding to dual-cells which intersect + // the boundary surface + label insertBoundaryConformPointPairs(const fileName& fName); + + void markNearBoundaryPoints(); + + //- Restore the Delaunay contraint + void fast_restore_Delaunay(Vertex_handle vh); + + // Flip operations used by fast_restore_Delaunay + void external_flip(Face_handle& f, int i); + bool internal_flip(Face_handle& f, int i); + + +public: + + // Constructors + + //- Construct for given surface + CV2D + ( + const Time& runTime, + const dictionary& controlDict + ); + + + //- Destructor + ~CV2D(); + + + // Member Functions + + // Access + + const controls& meshingControls() const + { + return controls_; + } + + + // Conversion functions between point2D, point and Point + + inline const point2D& toPoint2D(const point&) const; + inline point toPoint3D(const point2D&) const; + + #ifdef CGAL_INEXACT + typedef const point2D& point2DFromPoint; + typedef const Point& PointFromPoint2D; + #else + typedef point2D point2DFromPoint; + typedef Point PointFromPoint2D; + #endif + + inline point2DFromPoint toPoint2D(const Point&) const; + inline PointFromPoint2D toPoint(const point2D&) const; + inline point toPoint3D(const Point&) const; + + + // Point insertion + + //- Create the initial mesh from the given internal points. + // Points must be inside the boundary by at least nearness + // otherwise they are ignored. + void insertPoints + ( + const point2DField& points, + const scalar nearness + ); + + //- Create the initial mesh from the internal points in the given + // file. Points outside the geometry are ignored. + void insertPoints(const fileName& pointFileName); + + //- Create the initial mesh as a regular grid of points. + // Points outside the geometry are ignored. + void insertGrid(); + + //- Insert all surface point-pairs from + // insertSurfaceNearestPointPairs and + // findIntersectionForOutsideCentroid + void insertSurfacePointPairs(); + + //- Insert point-pairs where there are protrusions into + // or out of the surface + void boundaryConform(); + + + // Point removal + + //- Remove the point-pairs introduced by insertSurfacePointPairs + // and boundaryConform + void removeSurfacePointPairs(); + + + // Point motion + + inline void movePoint(const Vertex_handle& vh, const Point& P); + + //- Move the internal points to the given new locations and update + // the triangulation to ensure it is Delaunay + // void moveInternalPoints(const point2DField& newPoints); + + //- Calculate the displacements to create the new points + void newPoints(const scalar relaxation); + + //- Extract patch names and sizes. + void extractPatches + ( + wordList& patchNames, + labelList& patchSizes, + EdgeMap<label>& mapEdgesRegion + ) const; + + // Write + + //- Write internal points to .obj file + void writePoints(const fileName& fName, bool internalOnly) const; + + //- Write triangles as .obj file + void writeTriangles(const fileName& fName, bool internalOnly) const; + + //- Write dual faces as .obj file + void writeFaces(const fileName& fName, bool internalOnly) const; + + //- Calculates dual points (circumcentres of tets) and faces + // (point-cell walk of tets). + // Returns: + // - dualPoints (in triangle ordering) + // - dualFaces (compacted) + void calcDual + ( + point2DField& dualPoints, + faceList& dualFaces, + wordList& patchNames, + labelList& patchSizes, + EdgeMap<label>& mapEdgesRegion + ) const; + + //- Write patch + void writePatch(const fileName& fName) const; + + void write() const; +}; + + +inline bool boundaryTriangle(const CV2D::Face_handle fc); +inline bool outsideTriangle(const CV2D::Face_handle fc); + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "CV2DI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H b/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H new file mode 100644 index 00000000000..a389c32891d --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/CV2DI.H @@ -0,0 +1,182 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +inline Foam::label Foam::CV2D::insertPoint +( + const point2D& p, + const label type +) +{ + uint nVert = number_of_vertices(); + + Vertex_handle vh = insert(toPoint(p)); + + if (nVert == number_of_vertices()) + { + WarningIn("Foam::CV2D::insertPoint") + << "Failed to insert point " << p << endl; + } + else + { + vh->index() = nVert; + vh->type() = type; + } + + return vh->index(); +} + + +inline bool Foam::CV2D::insertMirrorPoint +( + const point2D& nearSurfPt, + const point2D& surfPt +) +{ + point2D mirrorPoint(2*surfPt - nearSurfPt); + + if (qSurf_.outside(toPoint3D(mirrorPoint))) + { + insertPoint(mirrorPoint, Vb::MIRROR_POINT); + return true; + } + else + { + return false; + } +} + + +inline void Foam::CV2D::insertPointPair +( + const scalar ppDist, + const point2D& surfPt, + const vector2D& n +) +{ + vector2D ppDistn = ppDist*n; + + label master = insertPoint + ( + surfPt - ppDistn, + number_of_vertices() + 1 + ); + + insertPoint(surfPt + ppDistn, master); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline const Foam::point2D& Foam::CV2D::toPoint2D(const point& p) const +{ + return reinterpret_cast<const point2D&>(p); +} + +inline Foam::point Foam::CV2D::toPoint3D(const point2D& p) const +{ + return point(p.x(), p.y(), z_); +} + + +#ifdef CGAL_INEXACT + +inline Foam::CV2D::point2DFromPoint Foam::CV2D::toPoint2D(const Point& P) const +{ + return reinterpret_cast<point2DFromPoint>(P); +} + +inline Foam::CV2D::PointFromPoint2D Foam::CV2D::toPoint(const point2D& p) const +{ + return reinterpret_cast<PointFromPoint2D>(p); +} + +#else + +inline Foam::CV2D::point2DFromPoint Foam::CV2D::toPoint2D(const Point& P) const +{ + return point2D(CGAL::to_double(P.x()), CGAL::to_double(P.y())); +} + +inline Foam::CV2D::PointFromPoint2D Foam::CV2D::toPoint(const point2D& p) const +{ + return Point(p.x(), p.y()); +} + +#endif + + +inline Foam::point Foam::CV2D::toPoint3D(const Point& P) const +{ + return point(CGAL::to_double(P.x()), CGAL::to_double(P.y()), z_); +} + + +inline void Foam::CV2D::movePoint(const Vertex_handle& vh, const Point& P) +{ + int i = vh->index(); + int t = vh->type(); + + remove(vh); + + Vertex_handle newVh = insert(P); + + newVh->index() = i; + newVh->type() = t; + + // label i = vh->index(); + // move(vh, P); + // vh->index() = i; + + // vh->set_point(P); + // fast_restore_Delaunay(vh); +} + + +// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // + +inline bool Foam::boundaryTriangle(const CV2D::Face_handle fc) +{ + return boundaryTriangle + ( + *fc->vertex(0), + *fc->vertex(1), + *fc->vertex(2) + ); +} + +inline bool Foam::outsideTriangle(const CV2D::Face_handle fc) +{ + return outsideTriangle + ( + *fc->vertex(0), + *fc->vertex(1), + *fc->vertex(2) + ); +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2DIO.C b/applications/utilities/mesh/generation/cv2DMesh/CV2DIO.C new file mode 100644 index 00000000000..6038bbefe0d --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/CV2DIO.C @@ -0,0 +1,304 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "CV2D.H" +#include "OFstream.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::CV2D::writePoints(const fileName& fName, bool internalOnly) const +{ + Info<< "Writing points to " << fName << nl << endl; + OFstream str(fName); + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (!internalOnly || vit->internalOrBoundaryPoint()) + { + meshTools::writeOBJ(str, toPoint3D(vit->point())); + } + } +} + + +void Foam::CV2D::writeTriangles(const fileName& fName, bool internalOnly) const +{ + Info<< "Writing triangles to " << fName << nl << endl; + OFstream str(fName); + + labelList vertexMap(number_of_vertices(), -2); + label verti = 0; + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (!internalOnly || !vit->farPoint()) + { + vertexMap[vit->index()] = verti++; + meshTools::writeOBJ(str, toPoint3D(vit->point())); + } + } + + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + ++fit + ) + { + if + ( + !internalOnly + || ( + fit->vertex(0)->internalOrBoundaryPoint() + || fit->vertex(1)->internalOrBoundaryPoint() + || fit->vertex(2)->internalOrBoundaryPoint() + ) + ) + { + str << "f"; + for (label i = 0; i < 3; ++i) + { + str << " " << vertexMap[fit->vertex(i)->index()] + 1; + } + str << nl; + } + } +} + + +void Foam::CV2D::writeFaces(const fileName& fName, bool internalOnly) const +{ + Info<< "Writing dual faces to " << fName << nl << endl; + OFstream str(fName); + + label dualVerti = 0; + + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + ++fit + ) + { + if + ( + !internalOnly + || ( + fit->vertex(0)->internalOrBoundaryPoint() + || fit->vertex(1)->internalOrBoundaryPoint() + || fit->vertex(2)->internalOrBoundaryPoint() + ) + ) + { + fit->faceIndex() = dualVerti++; + meshTools::writeOBJ(str, toPoint3D(circumcenter(fit))); + } + else + { + fit->faceIndex() = -1; + } + } + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (!internalOnly || vit->internalOrBoundaryPoint()) + { + Face_circulator fcStart = incident_faces(vit); + Face_circulator fc = fcStart; + + str<< 'f'; + + do + { + if (!is_infinite(fc)) + { + if (fc->faceIndex() < 0) + { + FatalErrorIn + ( + "Foam::CV2D::writeFaces" + "(const fileName& fName, bool internalOnly)" + )<< "Dual face uses vertex defined by a triangle" + " defined by an external point" + << exit(FatalError); + } + + str<< ' ' << fc->faceIndex() + 1; + } + } while (++fc != fcStart); + + str<< nl; + } + } +} + + +void Foam::CV2D::calcDual +( + point2DField& dualPoints, + faceList& dualFaces, + wordList& patchNames, + labelList& patchSizes, + EdgeMap<label>& mapEdgesRegion +) const +{ + // Dual points stored in triangle order. + dualPoints.setSize(number_of_faces()); + label dualVerti = 0; + + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + ++fit + ) + { + if + ( + fit->vertex(0)->internalOrBoundaryPoint() + || fit->vertex(1)->internalOrBoundaryPoint() + || fit->vertex(2)->internalOrBoundaryPoint() + ) + { + fit->faceIndex() = dualVerti; + + dualPoints[dualVerti++] = toPoint2D(circumcenter(fit)); + } + else + { + fit->faceIndex() = -1; + } + } + + dualPoints.setSize(dualVerti); + + extractPatches(patchNames, patchSizes, mapEdgesRegion); + + forAll(patchNames, patchI) + { + Info<< "Patch " << patchNames[patchI] + << " has size " << patchSizes[patchI] << endl; + } + + // Create dual faces + // ~~~~~~~~~~~~~~~~~ + + dualFaces.setSize(number_of_vertices()); + label dualFacei = 0; + labelList faceVerts(maxNvert); + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->internalOrBoundaryPoint()) + { + Face_circulator fcStart = incident_faces(vit); + Face_circulator fc = fcStart; + label verti = 0; + + do + { + if (!is_infinite(fc)) + { + if (fc->faceIndex() < 0) + { + FatalErrorIn + ( + "Foam::CV2D::calcDual" + "(point2DField& dualPoints, faceList& dualFaces)" + )<< "Dual face uses vertex defined by a triangle" + " defined by an external point" + << exit(FatalError); + } + + // Look up the index of the triangle + faceVerts[verti++] = fc->faceIndex(); + } + } while (++fc != fcStart); + + if (faceVerts.size() > 2) + { + dualFaces[dualFacei++] = + face(labelList::subList(faceVerts, verti)); + } + else + { + Info<< "From triangle point:" << vit->index() + << " coord:" << toPoint2D(vit->point()) + << " generated illegal dualFace:" << faceVerts + << endl; + } + } + } + + dualFaces.setSize(dualFacei); +} + + +void Foam::CV2D::writePatch(const fileName& fName) const +{ + point2DField dual2DPoints; + faceList dualFaces; + wordList patchNames; + labelList patchSizes; + EdgeMap<label> mapEdgesRegion; + + calcDual(dual2DPoints, dualFaces, patchNames, patchSizes, mapEdgesRegion); + pointField dualPoints(dual2DPoints.size()); + forAll(dualPoints, ip) + { + dualPoints[ip] = toPoint3D(dual2DPoints[ip]); + } + + // Dump as primitive patch to be read by extrudeMesh. + OFstream str(fName); + + Info<< "Writing patch to be used with extrudeMesh to file " << fName + << endl; + + str << dualPoints << nl << dualFaces << nl; +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/Make/files b/applications/utilities/mesh/generation/cv2DMesh/Make/files new file mode 100755 index 00000000000..75e58f25a15 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/Make/files @@ -0,0 +1,15 @@ +#include CGAL_FILES + +CV2D.C +controls.C +tolerances.C +insertFeaturePoints.C +insertSurfaceNearestPointPairs.C +insertSurfaceNearPointPairs.C +insertBoundaryConformPointPairs.C +CV2DIO.C +shortEdgeFilter2D.C +cv2DMesh.C + +EXE = $(FOAM_USER_APPBIN)/cv2DMesh + diff --git a/applications/utilities/mesh/generation/cv2DMesh/Make/options b/applications/utilities/mesh/generation/cv2DMesh/Make/options new file mode 100755 index 00000000000..2344ae6a691 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/Make/options @@ -0,0 +1,35 @@ +EXE_DEBUG = -DFULLDEBUG -g -O0 +EXE_FROUNDING_MATH = -frounding-math +EXE_NDEBUG = -DNDEBUG + +include $(GENERAL_RULES)/CGAL +FFLAGS = -DCGAL_FILES='"${CGAL_ARCH_PATH}/share/files"' + +EXE_INC = \ + ${EXE_FROUNDING_MATH} \ + ${EXE_NDEBUG} \ + ${CGAL_INC} \ + -I$(FOAM_APP)/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/lnInclude \ + -IconformalVoronoi2DMesh/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/surfMesh/lnInclude \ + -I$(LIB_SRC)/edgeMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/lnInclude \ + -I$(FOAM_APP)/utilities/mesh/generation/extrude/extrudeModel/lnInclude \ + -I$(LIB_SRC)/triSurface/lnInclude -DFULLDEBUG -g -O0 + +EXE_LIBS = \ + -L$(FOAM_USER_LIBBIN) \ + $(CGAL_LIBS) \ + -lmeshTools \ + -lsurfMesh \ + -ledgeMesh \ + -ltriSurface \ + -ldynamicMesh \ + -lextrude2DMesh \ + -lextrudeModel \ + -lconformalVoronoiMesh \ + -lcv2DMesh \ + -ldecompositionMethods diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files new file mode 100755 index 00000000000..df6cc213776 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/files @@ -0,0 +1,4 @@ +patchTo2DpolyMesh/patchTo2DpolyMesh.C + +LIB = $(FOAM_LIBBIN)/libcv2DMesh + diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options new file mode 100755 index 00000000000..298044fab80 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/Make/options @@ -0,0 +1,14 @@ +EXE_INC = \ + -I$(FOAM_APP)/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/surfMesh/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/triSurface/lnInclude -DFULLDEBUG -g -O0 + +LIB_LIBS = \ + -lmeshTools \ + -lsurfMesh \ + -ltriSurface \ + -ldynamicMesh \ + -lextrude2DMesh diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.C b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.C new file mode 100644 index 00000000000..c0121533bb7 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.C @@ -0,0 +1,334 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "patchTo2DpolyMesh.H" +#include "PatchTools.H" +#include "polyMesh.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::patchTo2DpolyMesh::flipFaceOrder() +{ + const edgeList& edges = patch_.edges(); + const faceList& localFaces = patch_.localFaces(); + const labelList& meshPoints = patch_.meshPoints(); + + Info<< "Flipping face order if necessary." << endl; + forAll(edges, edgeI) + { + const edge& e = edges[edgeI]; + + faces_[edgeI].setSize(2); + + label edgeOwner = owner_[edgeI]; + + const face& f = localFaces[edgeOwner]; + + label fp = findIndex(f, e[0]); + + if (f.nextLabel(fp) != e[1]) + { + Info<< "Flipping face " << faces_[edgeI] << endl; + faces_[edgeI][0] = meshPoints[e[1]]; + faces_[edgeI][1] = meshPoints[e[0]]; + } + else + { + faces_[edgeI][0] = meshPoints[e[0]]; + faces_[edgeI][1] = meshPoints[e[1]]; + } + } +} + + +void Foam::patchTo2DpolyMesh::createNeighbours() +{ + const edgeList& edges = patch_.edges(); + const labelListList& edgeFaces = patch_.edgeFaces(); + + Info<< "Calculating neighbours." << endl; + forAll(edges, edgeI) + { + const labelList& eFaces = edgeFaces[edgeI]; + if (eFaces.size() == 2) + { + if (owner_[edgeI] == eFaces[0]) + { + neighbour_[edgeI] = eFaces[1]; + } + else + { + neighbour_[edgeI] = eFaces[0]; + } + } + else if (eFaces.size() == 1) + { + continue; + } + else + { + FatalErrorIn("polyMesh neighbour construction") + << abort(FatalError); + } + } +} + + +Foam::labelList Foam::patchTo2DpolyMesh::internalFaceOrder() +{ + const labelListList& cellFaces = patch_.faceEdges(); + + labelList oldToNew(owner_.size(), -1); + + label newFaceI = 0; + + forAll(cellFaces, cellI) + { + const labelList& cFaces = cellFaces[cellI]; + // Neighbouring cells + SortableList<label> nbr(cFaces.size(), -1); + + forAll(cFaces, cfI) + { + if (cFaces[cfI] < neighbour_.size()) + { + label nbrFaceI = neighbour_[cFaces[cfI]]; + + // Internal face. Get cell on other side. + if (nbrFaceI == cellI) + { + nbrFaceI = owner_[cellI]; + } + + if (cellI < nbrFaceI) + { + // CellI is master + nbr[cfI] = nbrFaceI; + } + } + } + + nbr.sort(); + + forAll(nbr, i) + { + if (nbr[i] != -1) + { + oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++; + } + } + } + + return oldToNew; +} + + +void Foam::patchTo2DpolyMesh::addPatchFacesToFaces() +{ + const labelList& meshPoints = patch_.meshPoints(); + + label offset = patch_.nInternalEdges(); + face f(2); + + forAll(patchNames_, patchI) + { + forAllConstIter(EdgeMap<label>, mapEdgesRegion_, eIter) + { + if (eIter() == patchI) + { + f[0] = meshPoints[eIter.key().start()]; + f[1] = meshPoints[eIter.key().end()]; + faces_[offset++] = f; + } + } + } + + f.clear(); +} + + +void Foam::patchTo2DpolyMesh::addPatchFacesToOwner() +{ + const label nInternalEdges = patch_.nInternalEdges(); + const faceList& faces = patch_.faces(); + const label nExternalEdges = patch_.edges().size() - nInternalEdges; + const labelList& meshPoints = patch_.meshPoints(); + + // Reorder patch faces on owner list. + labelList newOwner = owner_; + + label nMatched = 0; + + for + ( + label bFaceI = nInternalEdges; + bFaceI < faces_.size(); + ++bFaceI + ) + { + const face& e = faces_[bFaceI]; + + bool matched = false; + + for + ( + label bEdgeI = nInternalEdges; + bEdgeI < faces_.size(); + ++bEdgeI + ) + { + if + ( + e[0] == meshPoints[patch_.edges()[bEdgeI][0]] + && e[1] == meshPoints[patch_.edges()[bEdgeI][1]] + ) + { + const face& f = faces[owner_[bEdgeI]]; + + label fp = findIndex(f, e[0]); + + newOwner[bFaceI] = owner_[bEdgeI]; + + if (f.nextLabel(fp) != e[1]) + { + Info<< "Flipping" << endl; + + faces_[bFaceI][0] = e[1]; + faces_[bFaceI][1] = e[0]; + } + + nMatched++; + + matched = true; + } + else if + ( + e[0] == patch_.edges()[bEdgeI][1] + && e[1] == patch_.edges()[bEdgeI][0] + ) + { + Info<< "Warning: Wrong orientation." << endl; + nMatched++; + matched = true; + } + } + if (!matched) + { + Info<< "No match for edge." << endl; + } + } + + if (nMatched != nExternalEdges) + { + Info<< "Number of matched edges, " << nMatched + << ", does not match number of external edges, " + << nExternalEdges << endl; + } + + owner_ = newOwner.xfer(); +} + + +void Foam::patchTo2DpolyMesh::createPolyMeshComponents() +{ + flipFaceOrder(); + + createNeighbours(); + + // New function for returning a map of old faces to new faces. + labelList oldToNew = internalFaceOrder(); + + inplaceReorder(oldToNew, faces_); + inplaceReorder(oldToNew, owner_); + inplaceReorder(oldToNew, neighbour_); + + // Add patches. + addPatchFacesToFaces(); + + addPatchFacesToOwner(); +} + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +//- Construct from a primitivePatch +Foam::patchTo2DpolyMesh::patchTo2DpolyMesh +( + const MeshedSurface<face>& patch, + const wordList& patchNames, + const labelList& patchSizes, + const EdgeMap<label>& mapEdgesRegion +) +: + patch_(patch), + patchNames_(patchNames), + patchSizes_(patchSizes), + patchStarts_(patchNames.size(), 0), + mapEdgesRegion_(mapEdgesRegion), + points_(patch.points()), + faces_(patch.nEdges()), + owner_(PatchTools::edgeOwner(patch)), + neighbour_(patch.nInternalEdges()) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::patchTo2DpolyMesh::~patchTo2DpolyMesh() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void Foam::patchTo2DpolyMesh::createMesh() +{ + createPolyMeshComponents(); + + label startFace = patch_.nInternalEdges(); + forAll(patchNames_, patchI) + { + patchStarts_[patchI] = startFace; + startFace += patchSizes_[patchI]; + } +} + +// * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * // + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.H b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.H new file mode 100644 index 00000000000..4b831b574e7 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/patchTo2DpolyMesh/patchTo2DpolyMesh.H @@ -0,0 +1,169 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-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 <http://www.gnu.org/licenses/>. + +Class + Foam::patchTo2DpolyMesh + +Description + Convert a primitivePatch into a 2D polyMesh. + +SourceFiles + patchTo2DpolyMeshI.H + patchTo2DpolyMesh.C + patchTo2DpolyMeshIO.C + +\*---------------------------------------------------------------------------*/ + +#ifndef patchTo2DpolyMesh_H +#define patchTo2DpolyMesh_H + +#include "EdgeMap.H" +#include "polyMesh.H" +#include "MeshedSurface.H" +#include "Time.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class patchTo2DpolyMesh Declaration +\*---------------------------------------------------------------------------*/ + +class patchTo2DpolyMesh +{ + // Private data + + const MeshedSurface<face>& patch_; + + const wordList& patchNames_; + + const labelList& patchSizes_; + + labelList patchStarts_; + + const EdgeMap<label>& mapEdgesRegion_; + + pointField points_; + faceList faces_; + labelList owner_; + labelList neighbour_; + + //- Description of data_ + void flipFaceOrder(); + + void createNeighbours(); + + labelList internalFaceOrder(); + + void addPatchFacesToFaces(); + + void addPatchFacesToOwner(); + + void createPolyMeshComponents(); + + + // Private Member Functions + + //- Disallow default bitwise copy construct + patchTo2DpolyMesh(const patchTo2DpolyMesh&); + + //- Disallow default bitwise assignment + void operator=(const patchTo2DpolyMesh&); + + +public: + + // Constructors + + //- Construct from a primitivePatch + patchTo2DpolyMesh + ( + const MeshedSurface<face>& patch, + const wordList& patchNames, + const labelList& patchSizes, + const EdgeMap<label>& mapEdgesRegion + ); + + + //- Destructor + ~patchTo2DpolyMesh(); + + + // Member Functions + + // Access + pointField& points() + { + return points_; + } + + faceList& faces() + { + return faces_; + } + + labelList& owner() + { + return owner_; + } + + labelList& neighbour() + { + return neighbour_; + } + + const wordList& patchNames() const + { + return patchNames_; + } + + const labelList& patchSizes() const + { + return patchSizes_; + } + + const labelList& patchStarts() const + { + return patchStarts_; + } + // Check + + // Edit + void createMesh(); + + // Write + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/controls.C b/applications/utilities/mesh/generation/cv2DMesh/controls.C new file mode 100644 index 00000000000..fd90a2d2ddd --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/controls.C @@ -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 <http://www.gnu.org/licenses/>. + +\*----------------------------------------------------------------------------*/ + +#include "controls.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::controls::controls(const dictionary& controlDict) +: + dict_(controlDict), + minCellSize(readScalar(controlDict.lookup("minCellSize"))), + minCellSize2(Foam::sqr(minCellSize)), + + featAngle(readScalar(controlDict.lookup("featureAngle"))), + maxQuadAngle(readScalar(controlDict.lookup("maxQuadAngle"))), + squares(controlDict.lookup("squares")), + + nearWallAlignedDist + ( + readScalar(controlDict.lookup("nearWallAlignedDist"))*minCellSize + ), + nearWallAlignedDist2(Foam::sqr(nearWallAlignedDist)), + + relaxOrientation(controlDict.lookup("relaxOrientation")), + + insertSurfaceNearestPointPairs + ( + controlDict.lookup("insertSurfaceNearestPointPairs") + ), + mirrorPoints(controlDict.lookup("mirrorPoints")), + insertSurfaceNearPointPairs + ( + controlDict.lookup("insertSurfaceNearPointPairs") + ), + writeInitialTriangulation(controlDict.lookup("writeInitialTriangulation")), + writeFeatureTriangulation(controlDict.lookup("writeFeatureTriangulation")), + writeNearestTriangulation(controlDict.lookup("writeNearestTriangulation")), + writeInsertedPointPairs(controlDict.lookup("writeInsertedPointPairs")), + writeFinalTriangulation(controlDict.lookup("writeFinalTriangulation")), + randomiseInitialGrid(controlDict.lookup("randomiseInitialGrid")), + randomPurturbation(readScalar(controlDict.lookup("randomPurturbation"))), + maxBoundaryConformingIter + ( + readLabel(controlDict.lookup("maxBoundaryConformingIter")) + ), + relaxationFactorStart + ( + readScalar(controlDict.lookup("relaxationFactorStart")) + ), + relaxationFactorEnd(readScalar(controlDict.lookup("relaxationFactorEnd"))) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::controls::~controls() +{} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/controls.H b/applications/utilities/mesh/generation/cv2DMesh/controls.H new file mode 100644 index 00000000000..1e6c7d79cd3 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/controls.H @@ -0,0 +1,164 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-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 <http://www.gnu.org/licenses/>. + +Class + Foam::controls + +Description + Controls for the 2D cv mesh generator. + +SourceFiles + controls.C + +\*---------------------------------------------------------------------------*/ + +#ifndef controls_H +#define controls_H + +#include "Switch.H" +#include "dictionary.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class controls Declaration +\*---------------------------------------------------------------------------*/ + +class controls +{ + // Private data + + //- Description of data_ + const dictionary& dict_; + + // Private Member Functions + + //- Disallow default bitwise copy construct + controls(const controls&); + + //- Disallow default bitwise assignment + void operator=(const controls&); + + +public: + + //- Minimum cell size below which protusions through the surface are + // not split + scalar minCellSize; + + //- Square of minCellSize + scalar minCellSize2; + + //- The feature angle used to select corners to be + // explicitly represented in the mesh. + // 0 = all features, 180 = no features + scalar featAngle; + + //- Maximum quadrant angle allowed at a concave corner before + // additional "mitering" lines are added + scalar maxQuadAngle; + + //- Should the mesh be square-dominated or of unbiased hexagons + Switch squares; + + //- Near-wall region where cells are aligned with the wall + scalar nearWallAlignedDist; + + //- Square of nearWallAlignedDist + scalar nearWallAlignedDist2; + + //- Chose if the cell orientation should relax during the iterations + // or remain fixed to the x-y directions + Switch relaxOrientation; + + //- Insert near-boundary point mirror or point-pairs + Switch insertSurfaceNearestPointPairs; + + //- Mirror near-boundary points rather than insert point-pairs + Switch mirrorPoints; + + //- Insert point-pairs vor dual-cell vertices very near the surface + Switch insertSurfaceNearPointPairs; + + Switch writeInitialTriangulation; + Switch writeFeatureTriangulation; + Switch writeNearestTriangulation; + Switch writeInsertedPointPairs; + Switch writeFinalTriangulation; + + Switch randomiseInitialGrid; + scalar randomPurturbation; + + label maxBoundaryConformingIter; + + //- Relaxation factor at the start of the iteration + scalar relaxationFactorStart; + + //- Relaxation factor at the end of the iteration + scalar relaxationFactorEnd; + + + // Constructors + + controls(const dictionary& controlDict); + + + //- Destructor + ~controls(); + + + // Member Functions + + // Access + + // Check + + // Edit + + // Write + + + // Friend Functions + + // Friend Operators + + // IOstream Operators +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//#include "controlsI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C b/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C new file mode 100644 index 00000000000..4ed77e79333 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C @@ -0,0 +1,214 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +Application + cv2DMesh + +Description + Conformal-Voronoi 2D extruding automatic mesher with grid or read + initial points and point position relaxation with optional + "squarification". + +\*---------------------------------------------------------------------------*/ + +#include "CV2D.H" +#include "argList.H" +#include "IFstream.H" + +#include "MeshedSurfaces.H" +#include "shortEdgeFilter2D.H" +#include "extrude2DMesh.H" +#include "polyMesh.H" +#include "PatchTools.H" +#include "patchTo2DpolyMesh.H" +#include "extrudeModel.H" +#include "polyTopoChange.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// Main program: + +int main(int argc, char *argv[]) +{ + argList::noParallel(); + argList::validArgs.clear(); + argList::validArgs.append("surface"); + argList::validOptions.insert("pointsFile", "<filename>"); + + #include "setRootCase.H" + #include "createTime.H" + + // Read control dictionary + // ~~~~~~~~~~~~~~~~~~~~~~~ + dictionary controlDict(IFstream("system/" + args.executable() + "Dict")()); + dictionary shortEdgeFilterDict(controlDict.subDict("shortEdgeFilter")); + dictionary extrusionDict(controlDict.subDict("extrusion")); + + Switch extrude(extrusionDict.lookup("extrude")); + label nIterations(readLabel(controlDict.lookup("nIterations"))); + label sefDebug(shortEdgeFilterDict.lookupOrDefault<label>("debug", 0)); + + // Read the surface to conform to + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// querySurface surf(args.args()[1]); +// surf.writeTreeOBJ(); + +// Info<< nl +// << "Read surface with " << surf.size() << " triangles from file " +// << args.args()[1] << nl << endl; + + // surf.write("surface.obj"); + + // Read and triangulation + // ~~~~~~~~~~~~~~~~~~~~~~ + CV2D mesh(runTime, controlDict); + if (args.options().found("pointsFile")) + { + fileName pointFileName(IStringStream(args.options()["pointsFile"])()); + mesh.insertPoints(pointFileName); + mesh.insertSurfacePointPairs(); + mesh.boundaryConform(); + } + else + { + mesh.insertGrid(); + mesh.insertSurfacePointPairs(); + mesh.boundaryConform(); + } + + for (int iter=1; iter<=nIterations; iter++) + { + Info<< nl + << "Relaxation iteration " << iter << nl + << "~~~~~~~~~~~~~~~~~~~~~~~~" << endl; + + scalar relax = + mesh.meshingControls().relaxationFactorStart + + + ( + mesh.meshingControls().relaxationFactorEnd + - mesh.meshingControls().relaxationFactorStart + ) + *scalar(iter)/scalar(nIterations); + + Info<< "Relaxation = " << relax << endl; + + mesh.newPoints(relax); + } + + mesh.write(); + + Info<< "Finished Delaunay in = " + << runTime.cpuTimeIncrement() << " s." << endl; + + shortEdgeFilter2D sef(mesh, shortEdgeFilterDict); + + shortEdgeFilter2D::debug = sefDebug; + + sef.filter(); + + Info<< "Meshed surface after edge filtering :" << endl; + sef.fMesh().writeStats(Info); + + Info<< "Write .obj file : MeshedSurface.obj" << endl; + sef.fMesh().write("MeshedSurface.obj"); + + Info<< "Finished filtering in = " + << runTime.cpuTimeIncrement() << " s." << endl; + + patchTo2DpolyMesh poly2DMesh + ( + sef.fMesh(), + sef.patchNames(), + sef.patchSizes(), + sef.mapEdgesRegion() + ); + + poly2DMesh.createMesh(); + + polyMesh pMesh + ( + IOobject + ( + polyMesh::defaultRegion, + runTime.constant(), + runTime, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + xferMove(poly2DMesh.points()), + xferMove(poly2DMesh.faces()), + xferMove(poly2DMesh.owner()), + xferMove(poly2DMesh.neighbour()) + ); + + Info<< "Constructing patches." << endl; + List<polyPatch*> patches(poly2DMesh.patchNames().size()); + + forAll(patches, patchI) + { + patches[patchI] = new polyPatch + ( + poly2DMesh.patchNames()[patchI], + poly2DMesh.patchSizes()[patchI], + poly2DMesh.patchStarts()[patchI], + patchI, + pMesh.boundaryMesh() + ); + } + + pMesh.addPatches(patches); + + if (extrude) + { + // Point generator + autoPtr<extrudeModel> model(extrudeModel::New(extrusionDict)); + + extrude2DMesh extruder(pMesh, extrusionDict, model()); + + polyTopoChange meshMod(pMesh.boundaryMesh().size()); + + extruder.addFrontBackPatches(); + extruder.setRefinement(meshMod); + + autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(pMesh, false); + + pMesh.updateMesh(morphMap); + + pMesh.setInstance("constant"); + } + + pMesh.write(); + + Info<< "Finished extruding in = " + << runTime.cpuTimeIncrement() << " s." << endl; + + Info<< nl << "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/indexedFace.H b/applications/utilities/mesh/generation/cv2DMesh/indexedFace.H new file mode 100644 index 00000000000..4d1dcce5747 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/indexedFace.H @@ -0,0 +1,147 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +Class + indexedFace + +Description + An indexed form of CGAL::Triangulation_face_base_2<K> used to keep + track of the vertices in the triangulation. + +\*---------------------------------------------------------------------------*/ + +#ifndef indexedFace_H +#define indexedFace_H + +#include <CGAL/Triangulation_2.h> +#include "indexedVertex.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace CGAL +{ + +/*---------------------------------------------------------------------------*\ + Class indexedFace Declaration +\*---------------------------------------------------------------------------*/ + +template<class Gt, class Fb=CGAL::Triangulation_face_base_2<Gt> > +class indexedFace +: + public Fb +{ + // Private data + + //- The index for this triangle face + // -1: triangle and changed and associated data needs updating + // >=0: index of triangles, set by external update algorithm + int index_; + + +public: + + enum faceTypes + { + UNCHANGED = 0, + CHANGED = -1, + SAVE_CHANGED = -2 + }; + + typedef typename Fb::Vertex_handle Vertex_handle; + typedef typename Fb::Face_handle Face_handle; + + template < typename TDS2 > + struct Rebind_TDS + { + typedef typename Fb::template Rebind_TDS<TDS2>::Other Fb2; + typedef indexedFace<Gt, Fb2> Other; + }; + + + indexedFace() + : + Fb(), + index_(CHANGED) + {} + + indexedFace(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2) + : + Fb(v0, v1, v2), + index_(CHANGED) + {} + + indexedFace + ( + Vertex_handle v0, + Vertex_handle v1, + Vertex_handle v2, + Face_handle n0, + Face_handle n1, + Face_handle n2 + ) + : + Fb(v0, v1, v2, n0, n1, n2), + index_(CHANGED) + {} + + + void set_vertex(int i, Vertex_handle v) + { + index_ = CHANGED; + Fb::set_vertex(i, v); + } + + void set_vertices() + { + index_ = CHANGED; + Fb::set_vertices(); + } + + void set_vertices(Vertex_handle v0, Vertex_handle v1, Vertex_handle v2) + { + index_ = CHANGED; + Fb::set_vertices(v0, v1, v2); + } + + + int& faceIndex() + { + return index_; + } + + int faceIndex() const + { + return index_; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace CGAL + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H b/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H new file mode 100644 index 00000000000..edb8d4c9541 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H @@ -0,0 +1,260 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +Class + indexedVertex + +Description + An indexed form of CGAL::Triangulation_vertex_base_2<K> used to keep + track of the vertices in the triangulation. + +\*---------------------------------------------------------------------------*/ + +#ifndef indexedVertex_H +#define indexedVertex_H + +#include <CGAL/Triangulation_2.h> + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace CGAL +{ + +/*---------------------------------------------------------------------------*\ + Class indexedVertex Declaration +\*---------------------------------------------------------------------------*/ + +template<class Gt, class Vb=CGAL::Triangulation_vertex_base_2<Gt> > +class indexedVertex +: + public Vb +{ + // Private data + + //- The index for this triangle vertex + int index_; + + //- Index of pair-point : + // NEAR_BOUNDARY_POINT : internal near boundary point. + // INTERNAL_POINT : internal point. + // FAR_POINT : far-point. + // >= 0 : part of point-pair. Index of other point. + // Lowest numbered is inside one (master). + int type_; + + +public: + + enum pointTypes + { + NEAR_BOUNDARY_POINT = -4, + INTERNAL_POINT = -3, + MIRROR_POINT = -2, + FAR_POINT = -1 + }; + + typedef typename Vb::Vertex_handle Vertex_handle; + typedef typename Vb::Face_handle Face_handle; + typedef typename Vb::Point Point; + + template<typename TDS2> + struct Rebind_TDS + { + typedef typename Vb::template Rebind_TDS<TDS2>::Other Vb2; + typedef indexedVertex<Gt,Vb2> Other; + }; + + + indexedVertex() + : + Vb(), + index_(INTERNAL_POINT), + type_(INTERNAL_POINT) + {} + + indexedVertex(const Point& p) + : + Vb(p), + index_(INTERNAL_POINT), + type_(INTERNAL_POINT) + {} + + indexedVertex(const Point& p, Face_handle f) + : + Vb(f, p), + index_(INTERNAL_POINT), + type_(INTERNAL_POINT) + {} + + indexedVertex(Face_handle f) + : + Vb(f), + index_(INTERNAL_POINT), + type_(INTERNAL_POINT) + {} + + + int& index() + { + return index_; + } + + int index() const + { + return index_; + } + + + int& type() + { + return type_; + } + + int type() const + { + return type_; + } + + + //- Is point a far-point + inline bool farPoint() const + { + return type_ == FAR_POINT; + } + + //- Is point internal, i.e. not on boundary + inline bool internalPoint() const + { + return type_ <= INTERNAL_POINT; + } + + //- Is point internal and near the boundary + inline bool nearBoundary() const + { + return type_ == NEAR_BOUNDARY_POINT; + } + + //- Set the point to be near the boundary + inline void setNearBoundary() + { + type_ = NEAR_BOUNDARY_POINT; + } + + //- Is point a mirror point + inline bool mirrorPoint() const + { + return type_ == MIRROR_POINT; + } + + //- Either master or slave of pointPair. + inline bool pairPoint() const + { + return type_ >= 0; + } + + //- Master of a pointPair is the lowest numbered one. + inline bool ppMaster() const + { + if (type_ > index_) + { + return true; + } + else + { + return false; + } + } + + //- Slave of a pointPair is the highest numbered one. + inline bool ppSlave() const + { + if (type_ >= 0 && type_ < index_) + { + return true; + } + else + { + return false; + } + } + + //- Either original internal point or master of pointPair. + inline bool internalOrBoundaryPoint() const + { + return internalPoint() || ppMaster(); + } + + //- Is point near the boundary or part of the boundary definition + inline bool nearOrOnBoundary() const + { + return pairPoint() || mirrorPoint() || nearBoundary(); + } + + //- Do the two given vertices consitute a boundary point-pair + inline friend bool pointPair + ( + const indexedVertex& v0, + const indexedVertex& v1 + ) + { + return v0.index_ == v1.type_ || v1.index_ == v0.type_; + } + + //- Do the three given vertices consitute a boundary triangle + inline friend bool boundaryTriangle + ( + const indexedVertex& v0, + const indexedVertex& v1, + const indexedVertex& v2 + ) + { + return (v0.pairPoint() && pointPair(v1, v2)) + || (v1.pairPoint() && pointPair(v2, v0)) + || (v2.pairPoint() && pointPair(v0, v1)); + } + + //- Do the three given vertices consitute an outside triangle + inline friend bool outsideTriangle + ( + const indexedVertex& v0, + const indexedVertex& v1, + const indexedVertex& v2 + ) + { + return (v0.farPoint() || v0.ppSlave()) + || (v1.farPoint() || v1.ppSlave()) + || (v2.farPoint() || v2.ppSlave()); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace CGAL + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C new file mode 100644 index 00000000000..dec7a5d38c7 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/insertBoundaryConformPointPairs.C @@ -0,0 +1,316 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +\*----------------------------------------------------------------------------*/ + +#include "CV2D.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::CV2D::insertPointPair +( + Triangulation::Finite_vertices_iterator& vit, + const point2D& p, + const label trii +) +{ + if + ( + !controls_.mirrorPoints + || !insertMirrorPoint(toPoint2D(vit->point()), p) + ) + { + pointIndexHit pHit + ( + true, + toPoint3D(p), + trii + ); + + vectorField norm(1); + qSurf_.geometry()[trii].getNormal + ( + List<pointIndexHit>(1, pHit), + norm + ); + + insertPointPair + ( + tols_.ppDist, + p, + toPoint2D(norm[0]) + ); + } + + vit = Triangulation::Finite_vertices_iterator + ( + CGAL::Filter_iterator + < + Triangulation::All_vertices_iterator, + Triangulation::Infinite_tester + >(finite_vertices_end(), vit.predicate(), vit.base()) + ); +} + + +bool Foam::CV2D::insertPointPairAtIntersection +( + Triangulation::Finite_vertices_iterator& vit, + const point2D& defVert, + const point2D vertices[], + const scalar maxProtSize2 +) +{ + bool found = false; + point2D interPoint; + label interTri = -1; + scalar interDist2 = 0; + + Face_circulator fcStart = incident_faces(vit); + Face_circulator fc = fcStart; + label vi = 0; + + do + { + if (!is_infinite(fc)) + { + pointIndexHit pHit; + label hitSurface = -1; + + qSurf_.findSurfaceAnyIntersection + ( + toPoint3D(defVert), + toPoint3D(vertices[vi]), + pHit, + hitSurface + ); + + if (pHit.hit()) + { + scalar dist2 = + magSqr(toPoint2D(pHit.hitPoint()) - vertices[vi]); + + // Check the point is further away than the furthest so far + if (dist2 > interDist2) + { + scalar mps2 = maxProtSize2; + + // If this is a boundary triangle reset the tolerance + // to avoid finding a hit point very close to a boundary + // vertex + if (boundaryTriangle(fc)) + { + mps2 = tols_.maxNotchLen2; + } + + if (dist2 > mps2) + { + found = true; + interPoint = toPoint2D(pHit.hitPoint()); + interTri = hitSurface; + interDist2 = dist2; + } + } + } + + vi++; + } + } while (++fc != fcStart); + + if (found) + { + insertPointPair(vit, interPoint, interTri); + return true; + } + else + { + return false; + } +} + + +Foam::label Foam::CV2D::insertBoundaryConformPointPairs +( + const fileName& fName +) +{ + label nIntersections = 0; + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + vit++ + ) + { + // Consider only those points part of point-pairs or near boundary + if (!vit->nearOrOnBoundary()) + { + continue; + } + + // Counter-clockwise circulator + Face_circulator fcStart = incident_faces(vit); + Face_circulator fc = fcStart; + + bool infinite = false; + bool changed = false; + + do + { + if (is_infinite(fc)) + { + infinite = true; + break; + } + else if (fc->faceIndex() < Fb::UNCHANGED) + { + changed = true; + break; + } + } while (++fc != fcStart); + + // If the dual-cell is connected to the infinite point or none of the + // faces whose circumcentres it uses have changed ignore + if (infinite || !changed) continue; + + fc = fcStart; + label nVerts = 0; + + do + { + vertices[nVerts++] = toPoint2D(circumcenter(fc)); + + if (nVerts == maxNvert) + { + break; + } + } while (++fc != fcStart); + + // Check if dual-cell has a large number of faces in which case + // assumed to be in the far-field and reject + if (nVerts == maxNvert) continue; + + // Set n+1 vertex to the first vertex for easy circulating + vertices[nVerts] = vertices[0]; + + // Convert triangle vertex to OpenFOAM point + point2DFromPoint defVert = toPoint2D(vit->point()); + + scalar maxProtSize2 = tols_.maxNotchLen2; + + if (vit->internalOrBoundaryPoint()) + { + // Calculate metrics of the dual-cell + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + // The perimeter of the dual-cell + scalar perimeter = 0; + + // Twice the area of the dual-cell + scalar areaT2 = 0; + + for (int vi=0; vi<nVerts; vi++) + { + vector2D edge(vertices[vi+1] - vertices[vi]); + perimeter += mag(edge); + vector2D otherEdge = defVert - vertices[vi]; + areaT2 += mag(edge.x()*otherEdge.y() - edge.y()*otherEdge.x()); + } + + // If the dual-cell is very small reject refinement + if (areaT2 < tols_.minEdgeLen2) continue; + + // Estimate the cell width + scalar cellWidth = areaT2/perimeter; + + + // Check dimensions of dual-cell + /* + // Quick rejection of dual-cell refinement based on it's perimeter + if (perimeter < 2*tols_.minCellSize) continue; + + // Also check the area of the cell and reject refinement + // if it is less than that allowed + if (areaT2 < tols_.minCellSize2) continue; + + // Estimate the cell width and reject refinement if it is less than + // that allowed + if (cellWidth < 0.5*tols_.minEdgeLen) continue; + */ + + if + ( + perimeter > 2*controls_.minCellSize + && areaT2 > controls_.minCellSize2 + && cellWidth > 0.5*tols_.minEdgeLen + ) + { + maxProtSize2 = 0.25*tols_.maxNotchLen2; + } + } + + if (insertPointPairAtIntersection(vit, defVert, vertices, maxProtSize2)) + { + nIntersections++; + } + } + + return nIntersections; +} + + +void Foam::CV2D::markNearBoundaryPoints() +{ + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + vit++ + ) + { + if (vit->internalPoint()) + { + point vert(toPoint3D(vit->point())); + + pointIndexHit pHit; + label hitSurface = -1; + + qSurf_.findSurfaceNearest + ( + vert, + 4*controls_.minCellSize2, + pHit, + hitSurface + ); + + if (pHit.hit()) + { + vit->setNearBoundary(); + } + } + } +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C b/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C new file mode 100644 index 00000000000..a1512f66f20 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/insertFeaturePoints.C @@ -0,0 +1,197 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +\*----------------------------------------------------------------------------*/ + +#include "CV2D.H" +#include "plane.H" +#include "triSurfaceTools.H" +#include "unitConversion.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +// Create feature points/edges by creating a triplet in the corner. +// (this triplet will have as its circumcentre the feature) +void Foam::CV2D::insertFeaturePoints() +{ + //labelList featEdges(qSurf_.extractFeatures2D(controls_.featAngle)); + + /* const PtrList<extendedFeatureEdgeMesh>& feMeshes + ( + qSurf_.features() + ); + + // const pointField& localPts = qSurf_.localPoints(); + + forAll(feMeshes, i) + { + const extendedFeatureEdgeMesh& feMesh(feMeshes[i]); + + // Loop over convex points + for + ( + label ptI = feMesh.convexStart(); + ptI < feMesh.concaveStart(); + ptI++ + ) + { + + const Foam::point& featPt = feMesh.points()[ptI]; + } + + // Loop over concave points + for + ( + label ptI = feMesh.concaveStart(); + ptI < feMesh.mixedStart(); + ptI++ + ) + { + + } + + // Loop over mixed points + for + ( + label ptI = feMesh.mixedStart(); + ptI < feMesh.nonFeatureStart(); + ptI++ + ) + { + + } + + //label edgeI = featEdges[i]; + //const edge& featEdge = qSurf_.edges()[edgeI]; + + // Get the feature point as the mid-point of the edge and convert to 2D + //point2D featPt = toPoint2D(featEdge.centre(qSurf_.localPoints())); + + // Pick up the two faces adjacent to the feature edge + const labelList& eFaces = qSurf_.edgeFaces()[edgeI]; + + label faceA = eFaces[0]; + vector2D nA = toPoint2D(qSurf_.faceNormals()[faceA]); + + label faceB = eFaces[1]; + vector2D nB = toPoint2D(qSurf_.faceNormals()[faceB]); + + // Intersect planes parallel to faceA and faceB offset by ppDist. + plane planeA(toPoint3D(featPt - tols_.ppDist*nA), toPoint3D(nA)); + plane planeB(toPoint3D(featPt - tols_.ppDist*nB), toPoint3D(nB)); + plane::ray interLine(planeA.planeIntersect(planeB)); + + // The reference point is where this line intersects the z_ plane + point2D refPt = toPoint2D + ( + interLine.refPoint() + + ((z_ - interLine.refPoint().z())/interLine.dir().z()) + *interLine.dir() + ); + + point2D faceAVert = toPoint2D + ( + localPts[triSurfaceTools::oppositeVertex(qSurf_, faceA, edgeI)] + ); + + // Determine convex or concave angle + if (((faceAVert - featPt) & nB) < 0) + { + // Convex. So refPt will be inside domain and hence a master point + + // Insert the master point refering the the first slave + label masterPtIndex = insertPoint(refPt, number_of_vertices() + 1); + + // Insert the slave points by reflecting refPt in both faces. + // with each slave refering to the master + + point2D reflectedA = refPt + 2*((featPt - refPt) & nA)*nA; + insertPoint(reflectedA, masterPtIndex); + + point2D reflectedB = refPt + 2*((featPt - refPt) & nB)*nB; + insertPoint(reflectedB, masterPtIndex); + } + else + { + // Concave. master and reflected points inside the domain. + // Generate reflected master to be outside. + point2D reflMasterPt = refPt + 2*(featPt - refPt); + + // Reflect refPt in both faces. + point2D reflectedA = + reflMasterPt + 2*((featPt - reflMasterPt) & nA)*nA; + + point2D reflectedB = + reflMasterPt + 2*((featPt - reflMasterPt) & nB)*nB; + + // Total angle around the concave feature + scalar totalAngle = + radToDeg(constant::mathematical::pi + acos(mag(nA & nB))); + + // Number of quadrants the angle should be split into + int nQuads = int(totalAngle/controls_.maxQuadAngle) + 1; + + // The number of additional master points needed to obtain the + // required number of quadrants. + int nAddPoints = min(max(nQuads - 2, 0), 2); + + // index of reflMaster + label reflectedMaster = number_of_vertices() + 2 + nAddPoints; + + // Master A is inside. + label reflectedAI = insertPoint(reflectedA, reflectedMaster); + + // Master B is inside. + insertPoint(reflectedB, reflectedMaster); + + if (nAddPoints == 1) + { + // One additinal point is the reflection of the slave point, + // i.e. the original reference point + insertPoint(refPt, reflectedMaster); + } + else if (nAddPoints == 2) + { + point2D reflectedAa = refPt - ((featPt - reflMasterPt) & nB)*nB; + insertPoint(reflectedAa, reflectedMaster); + + point2D reflectedBb = refPt - ((featPt - reflMasterPt) & nA)*nA; + insertPoint(reflectedBb, reflectedMaster); + } + + // Slave is outside. + insertPoint(reflMasterPt, reflectedAI); + } + } + + if (controls_.writeFeatureTriangulation) + { + writePoints("feat_allPoints.obj", false); + writeFaces("feat_allFaces.obj", false); + writeFaces("feat_faces.obj", true); + writeTriangles("feat_triangles.obj", true); + }*/ +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C new file mode 100644 index 00000000000..81349b8df55 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearPointPairs.C @@ -0,0 +1,115 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +\*----------------------------------------------------------------------------*/ + +#include "CV2D.H" +#include "treeDataTriSurface.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::CV2D::insertSurfaceNearPointPairs() +{ + Info<< "insertSurfaceNearPointPairs: "; + + label nNearPoints = 0; + + for + ( + Triangulation::Finite_edges_iterator eit = finite_edges_begin(); + eit != finite_edges_end(); + eit++ + ) + { + Vertex_handle v0h = eit->first->vertex(cw(eit->second)); + Vertex_handle v1h = eit->first->vertex(ccw(eit->second)); + + if (v0h->ppMaster() && v1h->ppMaster()) + { + point2DFromPoint v0(toPoint2D(v0h->point())); + point2DFromPoint v1(toPoint2D(v1h->point())); + + // Check that the two triangle vertices are further apart than the + // minimum cell size + if (magSqr(v1 - v0) > controls_.minCellSize2) + { + point2D e0(toPoint2D(circumcenter(eit->first))); + + point2D e1 + ( + toPoint2D(circumcenter(eit->first->neighbor(eit->second))) + ); + + // Calculate the length^2 of the edge normal to the surface + scalar edgeLen2 = magSqr(e0 - e1); + + if (edgeLen2 < tols_.minNearPointDist2) + { + pointIndexHit pHit; + label hitSurface = -1; + + qSurf_.findSurfaceNearest + ( + toPoint3D(e0), + tols_.minEdgeLen2, + pHit, + hitSurface + ); + + if (pHit.hit()) + { + vectorField norm(1); + qSurf_.geometry()[hitSurface].getNormal + ( + List<pointIndexHit>(1, pHit), + norm + ); + + insertPointPair + ( + tols_.ppDist, + toPoint2D(pHit.hitPoint()), + toPoint2D(norm[0]) + ); + + nNearPoints++; + + // Correct the edge iterator for the change in the + // number od edges following the point-pair insertion + eit = Finite_edges_iterator + ( + finite_edges_end().base(), + eit.predicate(), + eit.base() + ); + } + } + } + } + } + + Info<< nNearPoints << " point-pairs inserted" << endl; +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C new file mode 100644 index 00000000000..63ba5061d98 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/insertSurfaceNearestPointPairs.C @@ -0,0 +1,239 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +\*----------------------------------------------------------------------------*/ + +#include "CV2D.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +bool Foam::CV2D::dualCellSurfaceIntersection +( + const Triangulation::Finite_vertices_iterator& vit +) const +{ + Triangulation::Edge_circulator ecStart = incident_edges(vit); + Triangulation::Edge_circulator ec = ecStart; + + do + { + if (!is_infinite(ec)) + { + point e0 = toPoint3D(circumcenter(ec->first)); + + // If edge end is outside bounding box then edge cuts boundary + if (!qSurf_.globalBounds().contains(e0)) + { + return true; + } + + point e1 = toPoint3D(circumcenter(ec->first->neighbor(ec->second))); + + // If other edge end is ouside bounding box then edge cuts boundary + if (!qSurf_.globalBounds().contains(e1)) + { + return true; + } + + if (magSqr(e1 - e0) > tols_.minEdgeLen2) + { + if (qSurf_.findSurfaceAnyIntersection(e0, e1)) + { + return true; + } + } + } + + } while (++ec != ecStart); + + return false; +} + + +void Foam::CV2D::insertPointPairs +( + const DynamicList<point2D>& nearSurfacePoints, + const DynamicList<point2D>& surfacePoints, + const DynamicList<label>& surfaceTris, + const fileName fName +) +{ + if (controls_.mirrorPoints) + { + forAll(surfacePoints, ppi) + { + insertMirrorPoint + ( + nearSurfacePoints[ppi], + surfacePoints[ppi] + ); + } + } + else + { + forAll(surfacePoints, ppi) + { + pointIndexHit pHit + ( + true, + toPoint3D(surfacePoints[ppi]), + surfaceTris[ppi] + ); + + vectorField norm(1); + qSurf_.geometry()[surfaceTris[ppi]].getNormal + ( + List<pointIndexHit>(1, pHit), + norm + ); + + insertPointPair + ( + tols_.ppDist, + surfacePoints[ppi], + toPoint2D(norm[0]) + ); + } + } + + Info<< surfacePoints.size() << " point-pairs inserted" << endl; + + if (controls_.writeInsertedPointPairs) + { + OFstream str(fName); + label vertI = 0; + + forAll(surfacePoints, ppi) + { + meshTools::writeOBJ(str, toPoint3D(surfacePoints[ppi])); + vertI++; + } + + Info<< "insertPointPairs: Written " << surfacePoints.size() + << " inserted point-pair locations to file " + << str.name() << endl; + } +} + + +void Foam::CV2D::insertSurfaceNearestPointPairs() +{ + Info<< "insertSurfaceNearestPointPairs: "; + + label nSurfacePointsEst = min + ( + number_of_vertices(), + size_t(10*sqrt(scalar(number_of_vertices()))) + ); + + DynamicList<point2D> nearSurfacePoints(nSurfacePointsEst); + DynamicList<point2D> surfacePoints(nSurfacePointsEst); + DynamicList<label> surfaceTris(nSurfacePointsEst); + + // Local references to surface mesh addressing +// const pointField& localPoints = qSurf_.localPoints(); +// const labelListList& edgeFaces = qSurf_.edgeFaces(); +// const vectorField& faceNormals = qSurf_.faceNormals(); +// const labelListList& faceEdges = qSurf_.faceEdges(); + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + vit++ + ) + { + if (vit->internalPoint()) + { + point2DFromPoint vert(toPoint2D(vit->point())); + + pointIndexHit pHit; + label hitSurface = -1; + + qSurf_.findSurfaceNearest + ( + toPoint3D(vert), + 4*controls_.minCellSize2, + pHit, + hitSurface + ); + + if (pHit.hit()) + { + vit->setNearBoundary(); + + // Reference to the nearest triangle +// const labelledTri& f = qSurf_[hitSurface]; + +// // Find where point is on triangle. +// // Note tolerance needed is relative one +// // (used in comparing normalized [0..1] triangle coordinates). +// label nearType, nearLabel; +// triPointRef +// ( +// localPoints[f[0]], +// localPoints[f[1]], +// localPoints[f[2]] +// ).classify(pHit.hitPoint(), nearType, nearLabel); + +// // If point is on a edge check if it is an internal feature + +// bool internalFeatureEdge = false; + +// if (nearType == triPointRef::EDGE) +// { +// label edgeI = faceEdges[hitSurface][nearLabel]; +// const labelList& eFaces = edgeFaces[edgeI]; + +// if +// ( +// eFaces.size() == 2 +// && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) +// < -0.2 +// ) +// { +// internalFeatureEdge = true; +// } +// } + +// if (!internalFeatureEdge && dualCellSurfaceIntersection(vit)) + { + nearSurfacePoints.append(vert); + surfacePoints.append(toPoint2D(pHit.hitPoint())); + surfaceTris.append(hitSurface); + } + } + } + } + + insertPointPairs + ( + nearSurfacePoints, + surfacePoints, + surfaceTris, + "surfaceNearestIntersections.obj" + ); +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/querySurface.C b/applications/utilities/mesh/generation/cv2DMesh/querySurface.C new file mode 100644 index 00000000000..6a780c7a0d4 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/querySurface.C @@ -0,0 +1,237 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +\*----------------------------------------------------------------------------*/ + +#include "querySurface.H" +#include "unitConversion.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::querySurface::querySurface(const fileName& surfaceFileName) +: + triSurface(surfaceFileName), + rndGen_(12345), + bb_(localPoints()), + tree_ + ( + treeDataTriSurface + ( + *this, + indexedOctree<treeDataTriSurface>::perturbTol() + ), + bb_.extend(rndGen_, 1e-3), // slightly randomize bb + 8, // maxLevel + 4, //10, // leafsize + 10.0 //3.0 // duplicity + ) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::querySurface::~querySurface() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::labelList Foam::querySurface::extractFeatures2D +( + const scalar featAngle +) const +{ + scalar featCos = cos(degToRad(featAngle)); + + const labelListList& edgeFaces = this->edgeFaces(); + const pointField& points = this->points(); + const edgeList& edges = this->edges(); + const vectorField& faceNormals = this->faceNormals(); + const labelList& meshPoints = this->meshPoints(); + + DynamicList<label> featEdges(edges.size()); + + forAll(edgeFaces, edgeI) + { + const edge& e = edges[edgeI]; + + point p = + points[meshPoints[e.end()]] + - points[meshPoints[e.start()]]; + + if (magSqr(p & vector(1,1,0)) < SMALL) + { + const labelList& eFaces = edgeFaces[edgeI]; + + if + ( + eFaces.size() == 2 + && (faceNormals[eFaces[0]] & faceNormals[eFaces[1]]) < featCos + ) + { + featEdges.append(edgeI); + } + } + } + + return featEdges.shrink(); +} + +Foam::indexedOctree<Foam::treeDataTriSurface>::volumeType +Foam::querySurface::insideOutside +( + const scalar searchSpan2, + const point& pt +) const +{ + if (!bb_.contains(pt)) + { + return indexedOctree<treeDataTriSurface>::OUTSIDE; + } + + pointIndexHit pHit = tree_.findNearest(pt, searchSpan2); + + if (!pHit.hit()) + { + return tree_.getVolumeType(pt); + } + else + { + return indexedOctree<treeDataTriSurface>::MIXED; + } +} + + +// Check if point is inside surface +bool Foam::querySurface::inside(const point& pt) const +{ + if (!bb_.contains(pt)) + { + return false; + } + + return + ( + tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::INSIDE + ); +} + + +// Check if point is outside surface +bool Foam::querySurface::outside(const point& pt) const +{ + if (!bb_.contains(pt)) + { + return true; + } + + return + ( + tree_.getVolumeType(pt) == indexedOctree<treeDataTriSurface>::OUTSIDE + ); +} + + +// Check if point is inside surface by at least dist2 +bool Foam::querySurface::wellInside(const point& pt, const scalar dist2) const +{ + if (!bb_.contains(pt)) + { + return false; + } + + pointIndexHit pHit = tree_.findNearest(pt, dist2); + + if (pHit.hit()) + { + return false; + } + else + { + return + tree_.getVolumeType(pt) + == indexedOctree<treeDataTriSurface>::INSIDE; + } +} + + +// Check if point is outside surface by at least dist2 +bool Foam::querySurface::wellOutside(const point& pt, const scalar dist2) const +{ + if (!bb_.contains(pt)) + { + return true; + } + + pointIndexHit pHit = tree_.findNearest(pt, dist2); + + if (pHit.hit()) + { + return false; + } + else + { + return + tree_.getVolumeType(pt) + == indexedOctree<treeDataTriSurface>::OUTSIDE; + } +} + + +void Foam::querySurface::writeTreeOBJ() const +{ + OFstream str("tree.obj"); + label vertI = 0; + + const List<indexedOctree<treeDataTriSurface>::node>& nodes = tree_.nodes(); + + forAll(nodes, nodeI) + { + const indexedOctree<treeDataTriSurface>::node& nod = nodes[nodeI]; + + const treeBoundBox& bb = nod.bb_; + + const pointField points(bb.points()); + + label startVertI = vertI; + + forAll(points, i) + { + meshTools::writeOBJ(str, points[i]); + vertI++; + } + + const edgeList edges(treeBoundBox::edges); + + forAll(edges, i) + { + const edge& e = edges[i]; + + str << "l " << e[0]+startVertI+1 << ' ' << e[1]+startVertI+1 + << nl; + } + } +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/querySurface.H b/applications/utilities/mesh/generation/cv2DMesh/querySurface.H new file mode 100644 index 00000000000..b5d9ad0185c --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/querySurface.H @@ -0,0 +1,149 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +Class + querySurface + +Description + Searchable triSurface using an octree to speed-up queries. + +SourceFiles + querySurface.C + +\*---------------------------------------------------------------------------*/ + +#ifndef querySurface_H +#define querySurface_H + +#include "triSurface.H" +#include "treeDataTriSurface.H" +#include "indexedOctree.H" +#include "Random.H" +#include "meshTools.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class querySurface Declaration +\*---------------------------------------------------------------------------*/ + +class querySurface +: + public triSurface +{ + // Private data + + Random rndGen_; + + // Bounding box of surface. Used for relative tolerances. + treeBoundBox bb_; + + // Search engine on surface + indexedOctree<treeDataTriSurface> tree_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + querySurface(const querySurface&); + + //- Disallow default bitwise assignment + void operator=(const querySurface&); + + +public: + + // Constructors + + //- Construct given file name of the surface + querySurface(const fileName& surfaceFileName); + + + // Destructor + + ~querySurface(); + + + // Member Functions + + // Access + + const treeBoundBox& bb() const + { + return bb_; + } + + const indexedOctree<treeDataTriSurface>& tree() const + { + return tree_; + } + + + // Query + + //- Extract feature edges/points(2D) + // using the given feature angle in deg + labelList extractFeatures2D(const scalar featAngle) const; + + //- Returns inside, outside or mixed (= close to surface) + indexedOctree<Foam::treeDataTriSurface>::volumeType insideOutside + ( + const scalar searchSpan2, + const point& pt + ) const; + + //- Check if point is inside surface + bool inside(const point& pt) const; + + //- Check if point is outside surface + bool outside(const point& pt) const; + + //- Check if point is inside surface by at least dist2 + bool wellInside(const point& pt, const scalar dist2) const; + + //- Check if point is outside surface by at least dist2 + bool wellOutside(const point& pt, const scalar dist2) const; + + + // Write + + void writeTreeOBJ() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//#include "querySurfaceI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C new file mode 100644 index 00000000000..1ef698f3af6 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.C @@ -0,0 +1,502 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "shortEdgeFilter2D.H" + +namespace Foam +{ + +defineTypeNameAndDebug(shortEdgeFilter2D, 0); + +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::shortEdgeFilter2D::shortEdgeFilter2D +( + const Foam::CV2D& cv2Dmesh, + const dictionary& dict +) +: + cv2Dmesh_(cv2Dmesh), + shortEdgeFilterFactor_(readScalar(dict.lookup("shortEdgeFilterFactor"))), + edgeAttachedToBoundaryFactor_ + ( + dict.lookupOrDefault<scalar>("edgeAttachedToBoundaryFactor", 2.0) + ), + patchNames_(wordList()), + patchSizes_(labelList()), + mapEdgesRegion_() +{ + point2DField points2D; + faceList faces; + + cv2Dmesh.calcDual + ( + points2D, + faces, + patchNames_, + patchSizes_, + mapEdgesRegion_ + ); + + pointField points(points2D.size()); + forAll(points, ip) + { + points[ip] = cv2Dmesh.toPoint3D(points2D[ip]); + } + + points2D.clear(); + + ms_ = MeshedSurface<face>(xferMove(points), xferMove(faces)); + + Info<< "Meshed surface stats before edge filtering :" << endl; + ms_.writeStats(Info); + + ms_.write("MeshedSurface_preFilter.obj"); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::shortEdgeFilter2D::~shortEdgeFilter2D() +{} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +void +Foam::shortEdgeFilter2D::filter() +{ + if (debug) + { + writeInfo(Info); + } + + // These are global indices. + const pointField& points = ms_.points(); + const edgeList& edges = ms_.edges(); + const faceList& faces = ms_.faces(); + const labelList& meshPoints = ms_.meshPoints(); + const labelList& boundaryPoints = ms_.boundaryPoints(); + + label maxChain = 0; + label nPointsToRemove = 0; + + labelList pointsToRemove(ms_.points().size(), -1); + + // List of number of vertices in a face. + labelList newFaceVertexCount(faces.size(), -1); + forAll(faces, faceI) + { + newFaceVertexCount[faceI] = faces[faceI].size(); + } + + // Check if the point is a boundary point. Flag if it is so that + // it will not be deleted. + boolList boundaryPointFlags(points.size(), false); + // This has been removed, otherwise small edges on the boundary are not + // removed. + /* forAll(boundaryPointFlags, pointI) + { + forAll(boundaryPoints, bPoint) + { + if (meshPoints[boundaryPoints[bPoint]] == pointI) + { + boundaryPointFlags[pointI] = true; + } + } + }*/ + + // Check if an edge has a boundary point. It it does the edge length + // will be doubled when working out its length. + Info<< " Marking edges attached to boundaries." << endl; + boolList edgeAttachedToBoundary(edges.size(), false); + forAll(edges, edgeI) + { + const edge& e = edges[edgeI]; + const label startVertex = e.start(); + const label endVertex = e.end(); + + forAll(boundaryPoints, bPoint) + { + if + ( + boundaryPoints[bPoint] == startVertex + || boundaryPoints[bPoint] == endVertex + ) + { + edgeAttachedToBoundary[edgeI] = true; + } + } + } + + forAll(edges, edgeI) + { + const edge& e = edges[edgeI]; + + // get the vertices of that edge. + const label startVertex = e.start(); + const label endVertex = e.end(); + + scalar edgeLength = mag + ( + points[meshPoints[e.start()]] + -points[meshPoints[e.end()]] + ); + + if (edgeAttachedToBoundary[edgeI]) + { + edgeLength *= edgeAttachedToBoundaryFactor_; + } + + scalar shortEdgeFilterValue = 0.0; + + const labelList& psEdges = ms_.pointEdges()[startVertex]; + const labelList& peEdges = ms_.pointEdges()[endVertex]; + + forAll(psEdges, psEdgeI) + { + const edge& psE = edges[psEdges[psEdgeI]]; + if (edgeI != psEdges[psEdgeI]) + { + shortEdgeFilterValue += mag + ( + points[meshPoints[psE.start()]] + -points[meshPoints[psE.end()]] + ); + } + } + + forAll(peEdges, peEdgeI) + { + const edge& peE = edges[peEdges[peEdgeI]]; + if (edgeI != peEdges[peEdgeI]) + { + shortEdgeFilterValue += mag + ( + points[meshPoints[peE.start()]] + -points[meshPoints[peE.end()]] + ); + } + } + + shortEdgeFilterValue *= + shortEdgeFilterFactor_ + /(psEdges.size() + peEdges.size() - 2); + + if (edgeLength < shortEdgeFilterValue) + { + bool flagDegenerateFace = false; + const labelList& pFaces = ms_.pointFaces()[startVertex]; + + forAll(pFaces, pFaceI) + { + const face& f = ms_.localFaces()[pFaces[pFaceI]]; + forAll(f, fp) + { + // If the edge is part of this face... + if (f[fp] == endVertex) + { + // If deleting vertex would create a triangle, don't! + if (newFaceVertexCount[pFaces[pFaceI]] < 4) + { + flagDegenerateFace = true; + } + else + { + newFaceVertexCount[pFaces[pFaceI]]--; + } + } + // If the edge is not part of this face... + else + { + // Deleting vertex of a triangle... + if (newFaceVertexCount[pFaces[pFaceI]] < 3) + { + flagDegenerateFace = true; + } + } + } + } + + // This if statement determines whether a point should be deleted. + if + ( + pointsToRemove[meshPoints[startVertex]] == -1 + && pointsToRemove[meshPoints[endVertex]] == -1 + && !boundaryPointFlags[meshPoints[startVertex]] + && !flagDegenerateFace + ) + { + pointsToRemove[meshPoints[startVertex]] = + meshPoints[endVertex]; + ++nPointsToRemove; + } + } + } + + label totalNewPoints = points.size() - nPointsToRemove; + + pointField newPoints(totalNewPoints, vector(0, 0, 0)); + labelList newPointNumbers(points.size(), -1); + label numberRemoved=0; + + forAll(points, pointI) + { + // If the point is NOT going to be removed. + if (pointsToRemove[pointI] == -1) + { + newPoints[pointI-numberRemoved] = points[pointI]; + newPointNumbers[pointI] = pointI-numberRemoved; + } + else + { + numberRemoved++; + } + } + + // Need a new faceList + faceList newFaces(faces.size()); + label newFaceI = 0; + + labelList newFace; + label newFaceSize = 0; + + // Now need to iterate over the faces and remove points. Global index. + forAll(faces, faceI) + { + const face& f = faces[faceI]; + + newFace.clear(); + newFace.setSize(f.size()); + newFaceSize = 0; + + forAll(f, fp) + { + label pointI = f[fp]; + // If not removing the point, then add it to the new face. + if (pointsToRemove[pointI] == -1) + { + newFace[newFaceSize++] = newPointNumbers[pointI]; + } + else + { + label newPointI = pointsToRemove[pointI]; + // Replace deleted point with point that it is being + // collapsed to. + if + ( + f.nextLabel(fp) != newPointI + && f.prevLabel(fp) != newPointI + ) + { + label pChain = newPointI; + label totalChain = 0; + for (label nChain = 0; nChain <= totalChain; ++nChain) + { + if (newPointNumbers[pChain] != -1) + { + newFace[newFaceSize++] = newPointNumbers[pChain]; + newPointNumbers[pointI] + = newPointNumbers[pChain]; + maxChain = max(totalChain, maxChain); + } + else + { + WarningIn("shortEdgeFilter") + << "Point " << pChain + << " marked for deletion as well as point " + << pointI << nl + << " Incrementing maxChain by 1 from " + << totalChain << " to " << totalChain + 1 + << endl; + totalChain++; + } + pChain = pointsToRemove[pChain]; + } + } + else + { + if (newPointNumbers[newPointI] != -1) + { + newPointNumbers[pointI] = newPointNumbers[newPointI]; + } + } + } + } + + newFace.setSize(newFaceSize); + + if (newFace.size() > 2) + { + newFaces[newFaceI++] = face(newFace); + } + else + { + FatalErrorIn("shortEdgeFilter") + << "Only " << newFace.size() << " in face " << faceI + << exit(FatalError); + } + } + + newFaces.setSize(newFaceI); + + MeshedSurface<face> fMesh + ( + xferMove(newPoints), + xferMove(newFaces), + xferCopy(List<surfZone>()) + ); + + const Map<int>& fMeshPointMap = fMesh.meshPointMap(); + + // Reset patchSizes_ + patchSizes_.clear(); + patchSizes_.setSize(patchNames_.size(), 0); + + label equalEdges = 0; + label notFound = 0; + label matches = 0; + label negativeLabels = 0; + + forAll(newPointNumbers, pointI) + { + if (newPointNumbers[pointI] == -1) + { + WarningIn("shortEdgeFilter") + << pointI << " will be deleted and " << newPointNumbers[pointI] + << ", so it will not be replaced. " + << "This will cause edges to be deleted." << endl; + } + } + + // Create new EdgeMap. + Info<< "Creating new EdgeMap." << endl; + EdgeMap<label> newMapEdgesRegion(mapEdgesRegion_.size()); + + for + ( + label bEdgeI = ms_.nInternalEdges(); + bEdgeI < edges.size(); + ++bEdgeI + ) + { + label p1 = meshPoints[edges[bEdgeI][0]]; + label p2 = meshPoints[edges[bEdgeI][1]]; + + edge e(p1, p2); + + if (mapEdgesRegion_.found(e)) + { + if + ( + newPointNumbers[p1] != -1 + && newPointNumbers[p2] != -1 + ) + { + if (newPointNumbers[p1] != newPointNumbers[p2]) + { + label region = mapEdgesRegion_.find(e)(); + newMapEdgesRegion.insert + ( + edge + ( + fMeshPointMap[newPointNumbers[p1]], + fMeshPointMap[newPointNumbers[p2]] + ), + region + ); + patchSizes_[region]++; + matches++; + } + else + { + equalEdges++; + } + } + else + { + negativeLabels++; + } + } + else + { + notFound++; + } + } + + if (debug) + { + Info<< "EdgeMapping :" << nl + << " Matches : " << matches << nl + << " Equal : " << equalEdges << nl + << " Negative : " << negativeLabels << nl + << " Not Found: " << notFound << endl; + } + + mapEdgesRegion_.transfer(newMapEdgesRegion); + + ms_.transfer(fMesh); + + Info<< " Maximum number of chained collapses = " << maxChain << endl; + + if (debug) + { + writeInfo(Info); + } +} + + +void Foam::shortEdgeFilter2D::writeInfo(Ostream& os) +{ + os << "Short Edge Filtering Information:" << nl + << " shortEdgeFilterFactor : " << shortEdgeFilterFactor_ << nl + << " edgeAttachedToBoundaryFactor : " << edgeAttachedToBoundaryFactor_ + << endl; + + forAll(patchNames_, patchI) + { + os << " Patch " << patchNames_[patchI] + << ", size " << patchSizes_[patchI] << endl; + } + + os << " There are " << mapEdgesRegion_.size() + << " boundary edges." << endl; + + os << " Mesh Info:" << nl + << " Points: " << ms_.nPoints() << nl + << " Faces: " << ms_.size() << nl + << " Edges: " << ms_.nEdges() << nl + << " Internal: " << ms_.nInternalEdges() << nl + << " External: " << ms_.nEdges() - ms_.nInternalEdges() + << endl; +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H new file mode 100644 index 00000000000..0b63de73f7f --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H @@ -0,0 +1,137 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-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 <http://www.gnu.org/licenses/>. + +Class + Foam::shortEdgeFilter2D + +Description + This class filters short edges generated by the CV2D mesher. + +SourceFiles + shortEdgeFilter2D.C + +\*---------------------------------------------------------------------------*/ + +#ifndef shortEdgeFilter2D_H +#define shortEdgeFilter2D_H + +#include "MeshedSurfaces.H" +#include "CV2D.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + + +/*---------------------------------------------------------------------------*\ + Class shortEdgeFilter2D Declaration +\*---------------------------------------------------------------------------*/ + +class shortEdgeFilter2D +{ + // Private data + + //- Description of data_ + const CV2D& cv2Dmesh_; + + MeshedSurface<face> ms_; + + const scalar shortEdgeFilterFactor_; + + const scalar edgeAttachedToBoundaryFactor_; + + wordList patchNames_; + + labelList patchSizes_; + + EdgeMap<label> mapEdgesRegion_; + + + // Private Member Functions + + //- Disallow default bitwise copy construct + shortEdgeFilter2D(const shortEdgeFilter2D&); + + //- Disallow default bitwise assignment + void operator=(const shortEdgeFilter2D&); + + +public: + + //- Runtime type information + ClassName("shortEdgeFilter2D"); + + // Constructors + + shortEdgeFilter2D + ( + const CV2D& cv2Dmesh, + const dictionary& dict + ); + + + //- Destructor + ~shortEdgeFilter2D(); + + + // Access Functions + + const wordList& patchNames() const + { + return patchNames_; + } + + const labelList& patchSizes() const + { + return patchSizes_; + } + + const EdgeMap<label>& mapEdgesRegion() const + { + return mapEdgesRegion_; + } + + const MeshedSurface<face>& fMesh() const + { + return ms_; + } + + // Member Functions + + void filter(); + + void writeInfo(Ostream& os); + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/tolerances.C b/applications/utilities/mesh/generation/cv2DMesh/tolerances.C new file mode 100644 index 00000000000..ba033c907de --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/tolerances.C @@ -0,0 +1,68 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 <http://www.gnu.org/licenses/>. + +\*----------------------------------------------------------------------------*/ + +#include "tolerances.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::tolerances::tolerances +( + const dictionary& controlDict, + const scalar minCellSize, + const boundBox& bb +) +: + span + ( + max(mag(bb.max().x()), mag(bb.min().x())) + + max(mag(bb.max().y()), mag(bb.min().y())) + ), + span2(Foam::sqr(span)), + + minEdgeLen(readScalar(controlDict.lookup("minEdgeLenCoeff"))*minCellSize), + minEdgeLen2(Foam::sqr(minEdgeLen)), + + maxNotchLen(readScalar(controlDict.lookup("maxNotchLenCoeff"))*minCellSize), + maxNotchLen2(Foam::sqr(maxNotchLen)), + + minNearPointDist + ( + readScalar(controlDict.lookup("minNearPointDistCoeff"))*minCellSize + ), + minNearPointDist2(Foam::sqr(minNearPointDist)), + + ppDist(readScalar(controlDict.lookup("ppDistCoeff"))*minCellSize), + ppDist2(Foam::sqr(ppDist)) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::tolerances::~tolerances() +{} + + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/tolerances.H b/applications/utilities/mesh/generation/cv2DMesh/tolerances.H new file mode 100644 index 00000000000..7d99c3a61bc --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/tolerances.H @@ -0,0 +1,138 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-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 <http://www.gnu.org/licenses/>. + +Class + Foam::tolerances + +Description + Tolerances for the CV 2D mesher. + +SourceFiles + tolerances.C + +\*---------------------------------------------------------------------------*/ + +#ifndef tolerances_H +#define tolerances_H + +#include "dictionary.H" +#include "boundBox.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class tolerances Declaration +\*---------------------------------------------------------------------------*/ + +class tolerances +{ + // Private data + + //- Description of data_ + + + // Private Member Functions + + //- Disallow default bitwise copy construct + tolerances(const tolerances&); + + //- Disallow default bitwise assignment + void operator=(const tolerances&); + + +public: + + //- Maximum cartesian span of the geometry + scalar span; + + //- Square of span + scalar span2; + + //- Minumum edge-length of the cell size below which protusions + // through the surface are not split + scalar minEdgeLen; + + //- Square of minEdgeLen + scalar minEdgeLen2; + + //- Maximum notch size below which protusions into the surface are + // not filled + scalar maxNotchLen; + + //- Square of maxNotchLen + scalar maxNotchLen2; + + //- The minimum distance alowed between a dual-cell vertex + // and the surface before a point-pair is introduced + scalar minNearPointDist; + + //- Square of minNearPoint + scalar minNearPointDist2; + + //- Distance between boundary conforming point-pairs + scalar ppDist; + + //- Square of ppDist + scalar ppDist2; + + + // Constructors + + //- Construct null + tolerances + ( + const dictionary& controlDict, + scalar minCellSize, + const boundBox& + ); + + + //- Destructor + ~tolerances(); + + + // Member Functions + + // Access + + // Check + + // Edit + + // Write + +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C index 8f7ed466732..2edcb2410d8 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C @@ -152,12 +152,10 @@ bool Foam::cellSizeControlSurfaces::evalCellSizeFunctions Foam::cellSizeControlSurfaces::cellSizeControlSurfaces ( - const conformalVoronoiMesh& cvMesh, const searchableSurfaces& allGeometry, const dictionary& motionControlDict ) : - cvMesh_(cvMesh), allGeometry_(allGeometry), surfaces_(), cellSizeFunctions_(), @@ -223,7 +221,6 @@ Foam::cellSizeControlSurfaces::cellSizeControlSurfaces cellSizeFunction::New ( surfaceSubDict, - cvMesh, surface ) ); @@ -286,38 +283,38 @@ Foam::scalar Foam::cellSizeControlSurfaces::cellSize bool anyFunctionFound = evalCellSizeFunctions(pt, size); - if (!anyFunctionFound) - { - // Check if the point in question was actually inside the domain, if - // not, then it may be falling back to an inappropriate default size. - - if (cvMesh_.geometryToConformTo().outside(pt)) - { - pointIndexHit surfHit; - label hitSurface; - - cvMesh_.geometryToConformTo().findSurfaceNearest - ( - pt, - sqr(GREAT), - surfHit, - hitSurface - ); - - if (!surfHit.hit()) - { - FatalErrorIn - ( - "Foam::scalar Foam::cellSizeControlSurfaces::cellSize" - "(" - "const point& pt" - ") const" - ) - << "Point " << pt << " did not find a nearest surface point" - << nl << exit(FatalError) << endl; - } - } - } +// if (!anyFunctionFound) +// { +// // Check if the point in question was actually inside the domain, if +// // not, then it may be falling back to an inappropriate default size. + +// if (cvMesh_.geometryToConformTo().outside(pt)) +// { +// pointIndexHit surfHit; +// label hitSurface; + +// cvMesh_.geometryToConformTo().findSurfaceNearest +// ( +// pt, +// sqr(GREAT), +// surfHit, +// hitSurface +// ); + +// if (!surfHit.hit()) +// { +// FatalErrorIn +// ( +// "Foam::scalar Foam::cellSizeControlSurfaces::cellSize" +// "(" +// "const point& pt" +// ") const" +// ) +// << "Point " << pt << " did not find a nearest surface point" +// << nl << exit(FatalError) << endl; +// } +// } +// } return size; } diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H index 17dfd10325b..4b6c53bc2ac 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H @@ -57,9 +57,6 @@ class cellSizeControlSurfaces { // Private data - //- Reference to the conformalVoronoiMesh holding this object - const conformalVoronoiMesh& cvMesh_; - //- Reference to the searchableSurfaces object holding all geometry data const searchableSurfaces& allGeometry_; @@ -106,7 +103,6 @@ public: // searchableSurfaces cellSizeControlSurfaces ( - const conformalVoronoiMesh& cvMesh, const searchableSurfaces& allGeometry, const dictionary& motionControlDict ); diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.C index 9d7cc3084ab..37d6c7825ae 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.C @@ -45,12 +45,10 @@ cellSizeFunction::cellSizeFunction ( const word& type, const dictionary& cellSizeFunctionDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ) : dictionary(cellSizeFunctionDict), - cvMesh_(cvMesh), surface_(surface), coeffsDict_(subDict(type + "Coeffs")), sideMode_(), @@ -99,7 +97,6 @@ cellSizeFunction::cellSizeFunction autoPtr<cellSizeFunction> cellSizeFunction::New ( const dictionary& cellSizeFunctionDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ) { @@ -130,7 +127,7 @@ autoPtr<cellSizeFunction> cellSizeFunction::New return autoPtr<cellSizeFunction> ( - cstrIter()(cellSizeFunctionDict, cvMesh, surface) + cstrIter()(cellSizeFunctionDict, surface) ); } diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.H index 82535194a85..6b33ec7ca77 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/cellSizeFunction/cellSizeFunction.H @@ -81,9 +81,6 @@ protected: // Protected data - //- Reference to the conformalVoronoiMesh holding this cvs object - const conformalVoronoiMesh& cvMesh_; - //- Reference to the searchableSurface that cellSizeFunction // relates to const searchableSurface& surface_; @@ -120,10 +117,9 @@ public: dictionary, ( const dictionary& cellSizeFunctionDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ), - (cellSizeFunctionDict, cvMesh, surface) + (cellSizeFunctionDict, surface) ); @@ -134,7 +130,6 @@ public: ( const word& type, const dictionary& cellSizeFunctionDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ); @@ -145,7 +140,6 @@ public: static autoPtr<cellSizeFunction> New ( const dictionary& cellSizeFunctionDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ); diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.C index 702ac53bc68..0874d9fc823 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.C @@ -42,11 +42,10 @@ addToRunTimeSelectionTable(cellSizeFunction, linearDistance, dictionary); linearDistance::linearDistance ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ) : - cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), + cellSizeFunction(typeName, initialPointsDict, surface), surfaceCellSize_(readScalar(coeffsDict().lookup("surfaceCellSize"))), distanceCellSize_(readScalar(coeffsDict().lookup("distanceCellSize"))), distance_(readScalar(coeffsDict().lookup("distance"))), diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.H index 9799a1f7d06..0906680150b 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearDistance/linearDistance.H @@ -87,7 +87,6 @@ public: linearDistance ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ); diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearSpatial/linearSpatial.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearSpatial/linearSpatial.C index 1caeca1ffa2..d877ffee645 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearSpatial/linearSpatial.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearSpatial/linearSpatial.C @@ -42,11 +42,10 @@ addToRunTimeSelectionTable(cellSizeFunction, linearSpatial, dictionary); linearSpatial::linearSpatial ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ) : - cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), + cellSizeFunction(typeName, initialPointsDict, surface), referencePoint_(coeffsDict().lookup("referencePoint")), referenceCellSize_(readScalar(coeffsDict().lookup("referenceCellSize"))), direction_(coeffsDict().lookup("direction")), diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearSpatial/linearSpatial.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearSpatial/linearSpatial.H index 177f92fdf93..faf4d523594 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearSpatial/linearSpatial.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/linearSpatial/linearSpatial.H @@ -85,7 +85,6 @@ public: linearSpatial ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ); diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.C index 40960a81238..bb7ffa6a9a1 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.C @@ -46,11 +46,10 @@ addToRunTimeSelectionTable surfaceOffsetLinearDistance::surfaceOffsetLinearDistance ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ) : - cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), + cellSizeFunction(typeName, initialPointsDict, surface), surfaceCellSize_(readScalar(coeffsDict().lookup("surfaceCellSize"))), distanceCellSize_(readScalar(coeffsDict().lookup("distanceCellSize"))), surfaceOffset_(readScalar(coeffsDict().lookup("surfaceOffset"))), diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.H index 2b341181028..2b1d3a6f0e2 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/surfaceOffsetLinearDistance/surfaceOffsetLinearDistance.H @@ -94,7 +94,6 @@ public: surfaceOffsetLinearDistance ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ); diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.C index b8f3ffa6f46..e9a06786b48 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.C @@ -41,11 +41,10 @@ addToRunTimeSelectionTable(cellSizeFunction, uniform, dictionary); uniform::uniform ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ) : - cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), + cellSizeFunction(typeName, initialPointsDict, surface), cellSize_(readScalar(coeffsDict().lookup("cellSize"))) {} diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.H index 156e09a24fb..ea3031aad1a 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniform/uniform.H @@ -69,7 +69,6 @@ public: uniform ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ); diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.C index cba741fce2a..e67eceaa8cd 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.C @@ -41,11 +41,10 @@ addToRunTimeSelectionTable(cellSizeFunction, uniformDistance, dictionary); uniformDistance::uniformDistance ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ) : - cellSizeFunction(typeName, initialPointsDict, cvMesh, surface), + cellSizeFunction(typeName, initialPointsDict, surface), cellSize_(readScalar(coeffsDict().lookup("cellSize"))), distance_(readScalar(coeffsDict().lookup("distance"))), distanceSqr_(sqr(distance_)) diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.H index 78f06a2054f..bbd87fd548c 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeFunction/uniformDistance/uniformDistance.H @@ -75,7 +75,6 @@ public: uniformDistance ( const dictionary& initialPointsDict, - const conformalVoronoiMesh& cvMesh, const searchableSurface& surface ); diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C index 3c7720e92ac..6228f02e4cb 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C @@ -1925,13 +1925,13 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh ), geometryToConformTo_ ( - *this, + runTime_, + rndGen_, allGeometry_, cvMeshDict.subDict("surfaceConformation") ), cellSizeControl_ ( - *this, allGeometry_, cvMeshDict.subDict("motionControl") ), diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C index 9ad5ab6b3e0..5e8a8d72974 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C @@ -30,12 +30,14 @@ License Foam::conformationSurfaces::conformationSurfaces ( - const conformalVoronoiMesh& cvMesh, + const Time& runTime, + Random& rndGen, const searchableSurfaces& allGeometry, const dictionary& surfaceConformationDict ) : - cvMesh_(cvMesh), + runTime_(runTime), + rndGen_(rndGen), allGeometry_(allGeometry), features_(), locationInMesh_(surfaceConformationDict.lookup("locationInMesh")), @@ -139,9 +141,9 @@ Foam::conformationSurfaces::conformationSurfaces IOobject ( feMeshName, - cvMesh_.time().constant(), + runTime_.time().constant(), "extendedFeatureEdgeMesh", - cvMesh_.time(), + runTime_.time(), IOobject::MUST_READ, IOobject::NO_WRITE ) @@ -211,9 +213,9 @@ Foam::conformationSurfaces::conformationSurfaces IOobject ( feMeshName, - cvMesh_.time().constant(), + runTime_.time().constant(), "extendedFeatureEdgeMesh", - cvMesh_.time(), + runTime_.time(), IOobject::MUST_READ, IOobject::NO_WRITE ) @@ -232,7 +234,7 @@ Foam::conformationSurfaces::conformationSurfaces // Extend the global bounds to stop the bound box sitting on the surfaces // to be conformed to - globalBounds_ = globalBounds_.extend(cvMesh_.rndGen(), 1e-4); + globalBounds_ = globalBounds_.extend(rndGen_, 1e-4); // Look at all surfaces at determine whether the locationInMesh point is // inside or outside each, to establish a signature for the domain to be @@ -730,7 +732,7 @@ void Foam::conformationSurfaces::findEdgeNearestByType void Foam::conformationSurfaces::writeFeatureObj(const fileName& prefix) const { - OFstream ftStr(cvMesh_.time().path()/prefix + "_allFeatures.obj"); + OFstream ftStr(runTime_.time().path()/prefix + "_allFeatures.obj"); Pout<< nl << "Writing all features to " << ftStr.name() << endl; diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H index 4db3f158626..b541bed0032 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H @@ -56,8 +56,9 @@ class conformationSurfaces { // Private data - //- Reference to the conformalVoronoiMesh holding this object - const conformalVoronoiMesh& cvMesh_; + const Time& runTime_; + + Random& rndGen_; //- Reference to the searchableSurfaces object holding all geometry data const searchableSurfaces& allGeometry_; @@ -114,7 +115,8 @@ public: // searchableSurfaces conformationSurfaces ( - const conformalVoronoiMesh& cvMesh, + const Time& runTime, + Random& rndGen, const searchableSurfaces& allGeometry, const dictionary& surfaceConformationDict ); -- GitLab