From 5879a730fc59c1e90de7f4fa35d2257255648887 Mon Sep 17 00:00:00 2001 From: laurence Date: Thu, 3 Jan 2013 17:07:25 +0000 Subject: [PATCH 01/16] cv2DMesh: Initial Commit with tutorials --- .../utilities/mesh/generation/Allwmake | 1 + .../mesh/generation/cv2DMesh/Allwclean | 8 + .../mesh/generation/cv2DMesh/Allwmake | 8 + .../CGALTriangulation2DKernel.H} | 40 +- .../cv2DMesh/CGALTriangulation2Ddefs.H | 71 + .../utilities/mesh/generation/cv2DMesh/CV2D.C | 997 + .../utilities/mesh/generation/cv2DMesh/CV2D.H | 474 + .../mesh/generation/cv2DMesh/CV2DI.H | 227 + .../mesh/generation/cv2DMesh/CV2DIO.C | 386 + .../mesh/generation/cv2DMesh/Make/files | 12 + .../mesh/generation/cv2DMesh/Make/options | 43 + .../conformalVoronoi2DMesh/Make/files | 3 + .../conformalVoronoi2DMesh/Make/options | 3 + .../cv2DControls/cv2DControls.C | 154 + .../cv2DControls/cv2DControls.H | 260 + .../cv2DControls/cv2DControlsI.H | 158 + .../mesh/generation/cv2DMesh/cv2DMesh.C | 224 + .../mesh/generation/cv2DMesh/cv2DMeshDict | 208 + .../mesh/generation/cv2DMesh/indexedFace.H | 135 + .../mesh/generation/cv2DMesh/indexedFaceI.H | 110 + .../mesh/generation/cv2DMesh/indexedVertex.H | 210 + .../mesh/generation/cv2DMesh/indexedVertexI.H | 233 + .../insertBoundaryConformPointPairs.C | 323 + .../generation/cv2DMesh/insertFeaturePoints.C | 394 + .../cv2DMesh/insertSurfaceNearPointPairs.C | 114 + .../cv2DMesh/insertSurfaceNearestPointPairs.C | 244 + .../generation/cv2DMesh/shortEdgeFilter2D.C | 533 + .../generation/cv2DMesh/shortEdgeFilter2D.H | 133 + .../cellSizeControlSurfaces.C | 804 - .../cellSizeControlSurfaces.H | 266 - tutorials/mesh/cv2DMesh/OpenCFD/0.org/T | 50 + tutorials/mesh/cv2DMesh/OpenCFD/0.org/U | 59 + tutorials/mesh/cv2DMesh/OpenCFD/0.org/p | 51 + tutorials/mesh/cv2DMesh/OpenCFD/Allclean | 15 + tutorials/mesh/cv2DMesh/OpenCFD/Allrun | 15 + .../cv2DMesh/OpenCFD/Allrun-rhoCentralFoam | 18 + .../OpenCFD/constant/thermophysicalProperties | 43 + .../constant/triSurface/opencfd_box.stl | 96 + .../constant/triSurface/opencfd_text.stl | 16466 ++++++++++++++++ .../OpenCFD/constant/turbulenceProperties | 21 + .../mesh/cv2DMesh/OpenCFD/system/controlDict | 55 + .../OpenCFD/system/controlDict.mesher | 55 + .../OpenCFD/system/controlDict.rhoCentralFoam | 55 + .../mesh/cv2DMesh/OpenCFD/system/cv2DMeshDict | 158 + .../cv2DMesh/OpenCFD/system/decomposeParDict | 45 + .../cv2DMesh/OpenCFD/system/extrude2DMeshDict | 42 + .../mesh/cv2DMesh/OpenCFD/system/fvSchemes | 61 + .../mesh/cv2DMesh/OpenCFD/system/fvSolution | 42 + .../OpenCFD/system/surfaceFeatureExtractDict | 87 + .../mesh/cv2DMesh/jaggedBoundary/Allclean | 16 + tutorials/mesh/cv2DMesh/jaggedBoundary/Allrun | 13 + .../constant/triSurface/jaggedBoundary.stl | 954 + .../jaggedBoundary/system/controlDict | 47 + .../jaggedBoundary/system/cv2DMeshDict | 147 + .../jaggedBoundary/system/extrude2DMeshDict | 42 + .../cv2DMesh/jaggedBoundary/system/fvSchemes | 54 + .../cv2DMesh/jaggedBoundary/system/fvSolution | 22 + .../system/surfaceFeatureExtractDict | 52 + tutorials/mesh/cv2DMesh/square/Allclean | 16 + tutorials/mesh/cv2DMesh/square/Allrun | 13 + .../square/constant/triSurface/unit_cube.stl | 88 + .../mesh/cv2DMesh/square/system/controlDict | 52 + .../mesh/cv2DMesh/square/system/cv2DMeshDict | 165 + .../cv2DMesh/square/system/extrude2DMeshDict | 42 + .../mesh/cv2DMesh/square/system/fvSchemes | 54 + .../mesh/cv2DMesh/square/system/fvSolution | 22 + .../square/system/surfaceFeatureExtractDict | 52 + 67 files changed, 24975 insertions(+), 1086 deletions(-) create mode 100755 applications/utilities/mesh/generation/cv2DMesh/Allwclean create mode 100755 applications/utilities/mesh/generation/cv2DMesh/Allwmake rename applications/utilities/mesh/generation/{cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfacesI.H => cv2DMesh/CGALTriangulation2DKernel.H} (59%) 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/cv2DControls/cv2DControls.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControls.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/conformalVoronoi2DMesh/cv2DControls/cv2DControlsI.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/cv2DMeshDict create mode 100644 applications/utilities/mesh/generation/cv2DMesh/indexedFace.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/indexedFaceI.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/indexedVertex.H create mode 100644 applications/utilities/mesh/generation/cv2DMesh/indexedVertexI.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/shortEdgeFilter2D.C create mode 100644 applications/utilities/mesh/generation/cv2DMesh/shortEdgeFilter2D.H delete mode 100644 applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C delete mode 100644 applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.H create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/0.org/T create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/0.org/U create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/0.org/p create mode 100755 tutorials/mesh/cv2DMesh/OpenCFD/Allclean create mode 100755 tutorials/mesh/cv2DMesh/OpenCFD/Allrun create mode 100755 tutorials/mesh/cv2DMesh/OpenCFD/Allrun-rhoCentralFoam create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/constant/thermophysicalProperties create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/constant/triSurface/opencfd_box.stl create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/constant/triSurface/opencfd_text.stl create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/constant/turbulenceProperties create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/system/controlDict create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/system/controlDict.mesher create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/system/controlDict.rhoCentralFoam create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/system/cv2DMeshDict create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/system/decomposeParDict create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/system/extrude2DMeshDict create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/system/fvSchemes create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/system/fvSolution create mode 100644 tutorials/mesh/cv2DMesh/OpenCFD/system/surfaceFeatureExtractDict create mode 100755 tutorials/mesh/cv2DMesh/jaggedBoundary/Allclean create mode 100755 tutorials/mesh/cv2DMesh/jaggedBoundary/Allrun create mode 100644 tutorials/mesh/cv2DMesh/jaggedBoundary/constant/triSurface/jaggedBoundary.stl create mode 100644 tutorials/mesh/cv2DMesh/jaggedBoundary/system/controlDict create mode 100644 tutorials/mesh/cv2DMesh/jaggedBoundary/system/cv2DMeshDict create mode 100644 tutorials/mesh/cv2DMesh/jaggedBoundary/system/extrude2DMeshDict create mode 100644 tutorials/mesh/cv2DMesh/jaggedBoundary/system/fvSchemes create mode 100644 tutorials/mesh/cv2DMesh/jaggedBoundary/system/fvSolution create mode 100644 tutorials/mesh/cv2DMesh/jaggedBoundary/system/surfaceFeatureExtractDict create mode 100755 tutorials/mesh/cv2DMesh/square/Allclean create mode 100755 tutorials/mesh/cv2DMesh/square/Allrun create mode 100644 tutorials/mesh/cv2DMesh/square/constant/triSurface/unit_cube.stl create mode 100644 tutorials/mesh/cv2DMesh/square/system/controlDict create mode 100644 tutorials/mesh/cv2DMesh/square/system/cv2DMeshDict create mode 100644 tutorials/mesh/cv2DMesh/square/system/extrude2DMeshDict create mode 100644 tutorials/mesh/cv2DMesh/square/system/fvSchemes create mode 100644 tutorials/mesh/cv2DMesh/square/system/fvSolution create mode 100644 tutorials/mesh/cv2DMesh/square/system/surfaceFeatureExtractDict diff --git a/applications/utilities/mesh/generation/Allwmake b/applications/utilities/mesh/generation/Allwmake index a1141cfb1b..48ca530c96 100755 --- a/applications/utilities/mesh/generation/Allwmake +++ b/applications/utilities/mesh/generation/Allwmake @@ -12,6 +12,7 @@ wmake snappyHexMesh if [ -d "$CGAL_ARCH_PATH" ] then cvMesh/Allwmake + cv2DMesh/Allwmake fi # ----------------------------------------------------------------- end-of-file diff --git a/applications/utilities/mesh/generation/cv2DMesh/Allwclean b/applications/utilities/mesh/generation/cv2DMesh/Allwclean new file mode 100755 index 0000000000..d0ae53e415 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/Allwclean @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wclean libso conformalVoronoi2DMesh +wclean + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/utilities/mesh/generation/cv2DMesh/Allwmake b/applications/utilities/mesh/generation/cv2DMesh/Allwmake new file mode 100755 index 0000000000..5486849957 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/Allwmake @@ -0,0 +1,8 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake libso conformalVoronoi2DMesh +wmake + +# ----------------------------------------------------------------- end-of-file diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfacesI.H b/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2DKernel.H similarity index 59% rename from applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfacesI.H rename to applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2DKernel.H index f59a9cbdc0..8bbbe59920 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfacesI.H +++ b/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2DKernel.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -21,26 +21,34 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . +Typedefs + CGALTriangulation2DKernel + +Description + \*---------------------------------------------------------------------------*/ -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +#ifndef CGALTriangulation2DKernel_H +#define CGALTriangulation2DKernel_H -const Foam::searchableSurfaces& Foam::cellSizeControlSurfaces::geometry() const -{ - return allGeometry_; -} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "CGAL/Delaunay_triangulation_2.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 -const Foam::labelList& Foam::cellSizeControlSurfaces::surfaces() const -{ - return surfaces_; -} - - -Foam::scalar Foam::cellSizeControlSurfaces::defaultCellSize() const -{ - return defaultCellSize_; -} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#endif // ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2Ddefs.H b/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2Ddefs.H new file mode 100644 index 0000000000..6a17732356 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/CGALTriangulation2Ddefs.H @@ -0,0 +1,71 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Typedefs + CGALTriangulation2Ddefs + +Description + CGAL data structures used for 2D Delaunay meshing. + + Define CGAL_INEXACT to use Exact_predicates_inexact_constructions kernel + otherwise the more robust but much less efficient + Exact_predicates_exact_constructions will be used. + + Define CGAL_HIERARCHY to use hierarchical Delaunay triangulation which is + faster but uses more memory than the standard Delaunay triangulation. + +\*---------------------------------------------------------------------------*/ + +#ifndef CGALTriangulation2Ddefs_H +#define CGALTriangulation2Ddefs_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "CGALTriangulation2DKernel.H" + +#include "indexedVertex.H" +#include "indexedFace.H" + +typedef CGAL::indexedVertex Vb; +typedef CGAL::indexedFace Fb; + +#ifdef CGAL_HIERARCHY + // Data structures for hierarchical Delaunay triangulation which is more + // efficient but also uses more storage + #include "CGAL/Triangulation_hierarchy_2.h" + typedef CGAL::Triangulation_hierarchy_vertex_base_2 Vbh; + typedef CGAL::Triangulation_data_structure_2 Tds; + typedef CGAL::Delaunay_triangulation_2 DT; + typedef CGAL::Triangulation_hierarchy_2
Delaunay; +#else + // Data structures for standard Delaunay triangulation + typedef CGAL::Triangulation_data_structure_2 Tds; + typedef CGAL::Delaunay_triangulation_2 Delaunay; +#endif + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.C b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C new file mode 100644 index 0000000000..8a5931df61 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C @@ -0,0 +1,997 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*----------------------------------------------------------------------------*/ + +#include "CV2D.H" +#include "Random.H" +#include "transform.H" +#include "IFstream.H" +#include "uint.H" +#include "ulong.H" + +namespace Foam +{ + defineTypeNameAndDebug(CV2D, 0); +} + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::CV2D::insertBoundingBox() +{ + Info<< "insertBoundingBox: creating bounding mesh" << endl; + scalar bigSpan = 10*meshControls().span(); + insertPoint(point2D(-bigSpan, -bigSpan), Vb::FAR_POINT); + insertPoint(point2D(-bigSpan, bigSpan), Vb::FAR_POINT); + insertPoint(point2D(bigSpan, -bigSpan), Vb::FAR_POINT); + insertPoint(point2D(bigSpan, bigSpan), Vb::FAR_POINT); +} + + +void Foam::CV2D::fast_restore_Delaunay(Vertex_handle vh) +{ + int i; + Face_handle f = vh->face(), next, start(f); + + do + { + i=f->index(vh); + if (!is_infinite(f)) + { + if (!internal_flip(f, cw(i))) external_flip(f, i); + if (f->neighbor(i) == start) start = f; + } + f = f->neighbor(cw(i)); + } while (f != start); +} + + +void Foam::CV2D::external_flip(Face_handle& f, int i) +{ + Face_handle n = f->neighbor(i); + + if + ( + CGAL::ON_POSITIVE_SIDE + != side_of_oriented_circle(n, f->vertex(i)->point()) + ) return; + + flip(f, i); + i = n->index(f->vertex(i)); + external_flip(n, i); +} + + +bool Foam::CV2D::internal_flip(Face_handle& f, int i) +{ + Face_handle n = f->neighbor(i); + + if + ( + CGAL::ON_POSITIVE_SIDE + != side_of_oriented_circle(n, f->vertex(i)->point()) + ) + { + return false; + } + + flip(f, i); + + return true; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::CV2D::CV2D +( + const Time& runTime, + const dictionary& cvMeshDict +) +: + Delaunay(), + runTime_(runTime), + rndGen_(64293*Pstream::myProcNo()), + allGeometry_ + ( + IOobject + ( + "cvSearchableSurfaces", + runTime_.constant(), + "triSurface", + runTime_, + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + cvMeshDict.subDict("geometry") + ), + qSurf_ + ( + runTime_, + rndGen_, + allGeometry_, + cvMeshDict.subDict("surfaceConformation") + ), + controls_(cvMeshDict, qSurf_.globalBounds()), + cellSizeControl_ + ( + runTime, + cvMeshDict.subDict("motionControl").subDict("shapeControlFunctions"), + qSurf_, + controls_.minCellSize() + ), + relaxationModel_ + ( + relaxationModel::New + ( + cvMeshDict.subDict("motionControl"), + runTime + ) + ), + z_ + ( + point + ( + cvMeshDict.subDict("surfaceConformation").lookup("locationInMesh") + ).z() + ), + startOfInternalPoints_(0), + startOfSurfacePointPairs_(0), + startOfBoundaryConformPointPairs_(0), + featurePoints_() +{ + Info<< meshControls() << endl; + + insertBoundingBox(); + insertFeaturePoints(); +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::CV2D::~CV2D() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::CV2D::insertPoints +( + const point2DField& points, + const scalar nearness +) +{ + Info<< "insertInitialPoints(const point2DField& points): "; + + startOfInternalPoints_ = number_of_vertices(); + label nVert = startOfInternalPoints_; + + // Add the points and index them + forAll(points, i) + { + const point2D& p = points[i]; + + if (qSurf_.wellInside(toPoint3D(p), nearness)) + { + insert(toPoint(p))->index() = nVert++; + } + else + { + Warning + << "Rejecting point " << p << " outside surface" << endl; + } + } + + Info<< nVert << " vertices inserted" << endl; + + if (meshControls().objOutput()) + { + // Checking validity of triangulation + assert(is_valid()); + + writeTriangles("initial_triangles.obj", true); + writeFaces("initial_faces.obj", true); + } +} + + +void Foam::CV2D::insertPoints(const fileName& pointFileName) +{ + IFstream pointsFile(pointFileName); + + if (pointsFile.good()) + { + insertPoints + ( + point2DField(pointsFile), + 0.5*meshControls().minCellSize2() + ); + } + else + { + FatalErrorIn("insertInitialPoints") + << "Could not open pointsFile " << pointFileName + << exit(FatalError); + } +} + + +void Foam::CV2D::insertGrid() +{ + Info<< "insertInitialGrid: "; + + startOfInternalPoints_ = number_of_vertices(); + label nVert = startOfInternalPoints_; + + scalar x0 = qSurf_.globalBounds().min().x(); + scalar xR = qSurf_.globalBounds().max().x() - x0; + int ni = int(xR/meshControls().minCellSize()) + 1; + scalar deltax = xR/ni; + + scalar y0 = qSurf_.globalBounds().min().y(); + scalar yR = qSurf_.globalBounds().max().y() - y0; + int nj = int(yR/meshControls().minCellSize()) + 1; + scalar deltay = yR/nj; + + Random rndGen(1321); + scalar pert = meshControls().randomPerturbation()*min(deltax, deltay); + + for (int i=0; iindex() = nVert++; + } + } + } + + Info<< nVert << " vertices inserted" << endl; + + if (meshControls().objOutput()) + { + // Checking validity of triangulation + assert(is_valid()); + + writeTriangles("initial_triangles.obj", true); + writeFaces("initial_faces.obj", true); + } +} + + +void Foam::CV2D::insertSurfacePointPairs() +{ + startOfSurfacePointPairs_ = number_of_vertices(); + + if (meshControls().insertSurfaceNearestPointPairs()) + { + insertSurfaceNearestPointPairs(); + } + + write("nearest"); + + // Insertion of point-pairs for near-points may cause protrusions + // so insertBoundaryConformPointPairs must be executed last + if (meshControls().insertSurfaceNearPointPairs()) + { + insertSurfaceNearPointPairs(); + } + + startOfBoundaryConformPointPairs_ = number_of_vertices(); +} + + +void Foam::CV2D::boundaryConform() +{ + if (!meshControls().insertSurfaceNearestPointPairs()) + { + markNearBoundaryPoints(); + } + + // Mark all the faces as SAVE_CHANGED + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + fit++ + ) + { + fit->faceIndex() = Fb::SAVE_CHANGED; + } + + for (label iter=1; iter<=meshControls().maxBoundaryConformingIter(); iter++) + { + label nIntersections = insertBoundaryConformPointPairs + ( + "surfaceIntersections_" + Foam::name(iter) + ".obj" + ); + + if (nIntersections == 0) + { + break; + } + else + { + Info<< "BC iteration " << iter << ": " + << nIntersections << " point-pairs inserted" << endl; + } + + // Any faces changed by insertBoundaryConformPointPairs will now + // be marked CHANGED, mark those as SAVE_CHANGED and those that + // remained SAVE_CHANGED as UNCHANGED + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + fit++ + ) + { + if (fit->faceIndex() == Fb::SAVE_CHANGED) + { + fit->faceIndex() = Fb::UNCHANGED; + } + else if (fit->faceIndex() == Fb::CHANGED) + { + fit->faceIndex() = Fb::SAVE_CHANGED; + } + } + } + + Info<< nl; + + write("boundary"); +} + + +void Foam::CV2D::removeSurfacePointPairs() +{ + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->index() >= startOfSurfacePointPairs_) + { + remove(vit); + } + } +} + + +void Foam::CV2D::newPoints() +{ + const scalar relaxation = relaxationModel_->relaxation(); + + Info<< "Relaxation = " << relaxation << endl; + + Field dualVertices(number_of_faces()); + + label dualVerti = 0; + + // Find the dual point of each tetrahedron and assign it an index. + for + ( + Triangulation::Finite_faces_iterator fit = finite_faces_begin(); + fit != finite_faces_end(); + ++fit + ) + { + fit->faceIndex() = -1; + + if + ( + fit->vertex(0)->internalOrBoundaryPoint() + || fit->vertex(1)->internalOrBoundaryPoint() + || fit->vertex(2)->internalOrBoundaryPoint() + ) + { + fit->faceIndex() = dualVerti; + + dualVertices[dualVerti] = toPoint2D(circumcenter(fit)); + + dualVerti++; + } + } + + dualVertices.setSize(dualVerti); + + Field displacementAccumulator + ( + startOfSurfacePointPairs_, + vector2D::zero + ); + + // Calculate target size and alignment for vertices + scalarField sizes + ( + number_of_vertices(), + meshControls().minCellSize() + ); + + Field alignments + ( + number_of_vertices(), + vector2D(1, 0) + ); + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->internalOrBoundaryPoint()) + { + point2D vert = toPoint2D(vit->point()); + + // alignment and size determination + pointIndexHit pHit; + label hitSurface = -1; + + qSurf_.findSurfaceNearest + ( + toPoint3D(vert), + meshControls().span2(), + pHit, + hitSurface + ); + + if (pHit.hit()) + { + vectorField norm(1); + allGeometry_[hitSurface].getNormal + ( + List(1, pHit), + norm + ); + + alignments[vit->index()] = toPoint2D(norm[0]); + + sizes[vit->index()] = + cellSizeControl_.cellSize(toPoint3D(vit->point())); + } + } + } + + // Info<< "Calculated alignments" << endl; + + scalar cosAlignmentAcceptanceAngle = 0.68; + + // Upper and lower edge length ratios for weight + scalar u = 1.0; + scalar l = 0.7; + + PackedBoolList pointToBeRetained(startOfSurfacePointPairs_, true); + + std::list pointsToInsert; + + for + ( + Triangulation::Finite_edges_iterator eit = finite_edges_begin(); + eit != finite_edges_end(); + eit++ + ) + { + Vertex_handle vA = eit->first->vertex(cw(eit->second)); + Vertex_handle vB = eit->first->vertex(ccw(eit->second)); + + if (!vA->internalOrBoundaryPoint() || !vB->internalOrBoundaryPoint()) + { + continue; + } + + const point2D& dualV1 = dualVertices[eit->first->faceIndex()]; + const point2D& dualV2 = + dualVertices[eit->first->neighbor(eit->second)->faceIndex()]; + + scalar dualEdgeLength = mag(dualV1 - dualV2); + + point2D dVA = toPoint2D(vA->point()); + point2D dVB = toPoint2D(vB->point()); + + Field alignmentDirsA(2); + + alignmentDirsA[0] = alignments[vA->index()]; + alignmentDirsA[1] = vector2D + ( + -alignmentDirsA[0].y(), + alignmentDirsA[0].x() + ); + + Field alignmentDirsB(2); + + alignmentDirsB[0] = alignments[vB->index()]; + alignmentDirsB[1] = vector2D + ( + -alignmentDirsB[0].y(), + alignmentDirsB[0].x() + ); + + Field alignmentDirs(2); + + forAll(alignmentDirsA, aA) + { + const vector2D& a(alignmentDirsA[aA]); + + scalar maxDotProduct = 0.0; + + forAll(alignmentDirsB, aB) + { + const vector2D& b(alignmentDirsB[aB]); + + scalar dotProduct = a & b; + + if (mag(dotProduct) > maxDotProduct) + { + maxDotProduct = mag(dotProduct); + + alignmentDirs[aA] = a + sign(dotProduct)*b; + + alignmentDirs[aA] /= mag(alignmentDirs[aA]); + } + } + } + + vector2D rAB = dVA - dVB; + + scalar rABMag = mag(rAB); + + forAll(alignmentDirs, aD) + { + vector2D& alignmentDir = alignmentDirs[aD]; + + if ((rAB & alignmentDir) < 0) + { + // swap the direction of the alignment so that has the + // same sense as rAB + alignmentDir *= -1; + } + + scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir); + + if (alignmentDotProd > cosAlignmentAcceptanceAngle) + { + scalar targetFaceSize = + 0.5*(sizes[vA->index()] + sizes[vB->index()]); + + // Test for changing aspect ratio on second alignment (first + // alignment is neartest surface normal) + // if (aD == 1) + // { + // targetFaceSize *= 2.0; + // } + + alignmentDir *= 0.5*targetFaceSize; + + vector2D delta = alignmentDir - 0.5*rAB; + + if (dualEdgeLength < 0.7*targetFaceSize) + { + delta *= 0; + } + else if (dualEdgeLength < targetFaceSize) + { + delta *= + ( + dualEdgeLength + /(targetFaceSize*(u - l)) + - 1/((u/l) - 1) + ); + } + + if + ( + vA->internalPoint() + && vB->internalPoint() + && rABMag > 1.75*targetFaceSize + && dualEdgeLength > 0.05*targetFaceSize + && alignmentDotProd > 0.93 + ) + { + // Point insertion + pointsToInsert.push_back(toPoint(0.5*(dVA + dVB))); + } + else if + ( + (vA->internalPoint() || vB->internalPoint()) + && rABMag < 0.65*targetFaceSize + ) + { + // Point removal + + // Only insert a point at the midpoint of the short edge + // if neither attached point has already been identified + // to be removed. + if + ( + pointToBeRetained[vA->index()] == true + && pointToBeRetained[vB->index()] == true + ) + { + pointsToInsert.push_back(toPoint(0.5*(dVA + dVB))); + } + + if (vA->internalPoint()) + { + pointToBeRetained[vA->index()] = false; + } + + if (vB->internalPoint()) + { + pointToBeRetained[vB->index()] = false; + } + } + else + { + if (vA->internalPoint()) + { + displacementAccumulator[vA->index()] += delta; + } + + if (vB->internalPoint()) + { + displacementAccumulator[vB->index()] += -delta; + } + } + } + } + } + + vector2D totalDisp = sum(displacementAccumulator); + scalar totalDist = sum(mag(displacementAccumulator)); + + // Relax the calculated displacement + displacementAccumulator *= relaxation; + + label numberOfNewPoints = pointsToInsert.size(); + + for + ( + Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if (vit->internalPoint()) + { + if (pointToBeRetained[vit->index()]) + { + pointsToInsert.push_front + ( + toPoint + ( + toPoint2D(vit->point()) + + displacementAccumulator[vit->index()] + ) + ); + } + } + } + + // Clear the triangulation and reinsert the bounding box and feature points. + // This is faster than removing and moving points. + this->clear(); + + insertBoundingBox(); + + reinsertFeaturePoints(); + + startOfInternalPoints_ = number_of_vertices(); + + label nVert = startOfInternalPoints_; + + Info<< "Inserting " << numberOfNewPoints << " new points" << endl; + + // Use the range insert as it is faster than individually inserting points. + insert(pointsToInsert.begin(), pointsToInsert.end()); + + for + ( + Delaunay::Finite_vertices_iterator vit = finite_vertices_begin(); + vit != finite_vertices_end(); + ++vit + ) + { + if + ( + vit->type() == Vb::INTERNAL_POINT + && vit->index() == Vb::INTERNAL_POINT + ) + { + vit->index() = nVert++; + } + } + + Info<< " Total displacement = " << totalDisp << nl + << " Total distance = " << totalDist << nl + << " Points added = " << pointsToInsert.size() + << endl; + + write("internal"); + + insertSurfacePointPairs(); + + boundaryConform(); + + +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// Old Method +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +// for +// ( +// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); +// vit != finite_vertices_end(); +// ++vit +// ) +// { +// if (vit->internalPoint()) +// { +// // Current dual-cell defining vertex ("centre") +// point2DFromPoint defVert0 = toPoint2D(vit->point()); + +// Triangulation::Edge_circulator ec = incident_edges(vit); +// Triangulation::Edge_circulator ecStart = ec; + +// // Circulate around the edges to find the first which is not +// // infinite +// do +// { +// if (!is_infinite(ec)) break; +// } while (++ec != ecStart); + +// // Store the start-end of the first non-infinte edge +// point2D de0 = toPoint2D(circumcenter(ec->first)); + +// // Keep track of the maximum edge length^2 +// scalar maxEdgeLen2 = 0.0; + +// // Keep track of the index of the longest edge +// label edgecd0i = -1; + +// // Edge counter +// label edgei = 0; + +// do +// { +// if (!is_infinite(ec)) +// { +// // Get the end of the current edge +// point2D de1 = toPoint2D +// ( +// circumcenter(ec->first->neighbor(ec->second)) +// ); + +// // Store the current edge vector +// edges[edgei] = de1 - de0; + +// // Store the edge mid-point in the vertices array +// vertices[edgei] = 0.5*(de1 + de0); + +// // Move the current edge end into the edge start for the +// // next iteration +// de0 = de1; + +// // Keep track of the longest edge + +// scalar edgeLen2 = magSqr(edges[edgei]); + +// if (edgeLen2 > maxEdgeLen2) +// { +// maxEdgeLen2 = edgeLen2; +// edgecd0i = edgei; +// } + +// edgei++; +// } +// } while (++ec != ecStart); + +// // Initialise cd0 such that the mesh will align +// // in in the x-y directions +// vector2D cd0(1, 0); + +// if (meshControls().relaxOrientation()) +// { +// // Get the longest edge from the array and use as the primary +// // direction of the coordinate system of the "square" cell +// cd0 = edges[edgecd0i]; +// } + +// if (meshControls().nearWallAlignedDist() > 0) +// { +// pointIndexHit pHit = qSurf_.tree().findNearest +// ( +// toPoint3D(defVert0), +// meshControls().nearWallAlignedDist2() +// ); + +// if (pHit.hit()) +// { +// cd0 = toPoint2D(faceNormals[pHit.index()]); +// } +// } + +// // Rotate by 45deg needed to create an averaging procedure which +// // encourages the cells to be square +// cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x()); + +// // Normalise the primary coordinate direction +// cd0 /= mag(cd0); + +// // Calculate the orthogonal coordinate direction +// vector2D cd1(-cd0.y(), cd0.x()); + + +// // Restart the circulator +// ec = ecStart; + +// // ... and the counter +// edgei = 0; + +// // Initialise the displacement for the centre and sum-weights +// vector2D disp = vector2D::zero; +// scalar sumw = 0; + +// do +// { +// if (!is_infinite(ec)) +// { +// // Pick up the current edge +// const vector2D& ei = edges[edgei]; + +// // Calculate the centre to edge-centre vector +// vector2D deltai = vertices[edgei] - defVert0; + +// // Set the weight for this edge contribution +// scalar w = 1; + +// if (meshControls().squares()) +// { +// w = magSqr(deltai.x()*ei.y() - deltai.y()*ei.x()); +// // alternative weights +// //w = mag(deltai.x()*ei.y() - deltai.y()*ei.x()); +// //w = magSqr(ei)*mag(deltai); + +// // Use the following for an ~square mesh +// // Find the coordinate contributions for this edge delta +// scalar cd0deltai = cd0 & deltai; +// scalar cd1deltai = cd1 & deltai; + +// // Create a "square" displacement +// if (mag(cd0deltai) > mag(cd1deltai)) +// { +// disp += (w*cd0deltai)*cd0; +// } +// else +// { +// disp += (w*cd1deltai)*cd1; +// } +// } +// else +// { +// // Use this for a hexagon/pentagon mesh +// disp += w*deltai; +// } + +// // Sum the weights +// sumw += w; +// } +// else +// { +// FatalErrorIn("CV2D::newPoints() const") +// << "Infinite triangle found in internal mesh" +// << exit(FatalError); +// } + +// edgei++; + +// } while (++ec != ecStart); + +// // Calculate the average displacement +// disp /= sumw; +// totalDisp += disp; +// totalDist += mag(disp); + +// // Move the point by a fraction of the average displacement +// movePoint(vit, defVert0 + relaxation*disp); +// } +// } + +// Info << "\nTotal displacement = " << totalDisp +// << " total distance = " << totalDist << endl; +} + + +//void Foam::CV2D::moveInternalPoints(const point2DField& newPoints) +//{ +// label pointI = 0; + +// for +// ( +// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); +// vit != finite_vertices_end(); +// ++vit +// ) +// { +// if (vit->internalPoint()) +// { +// movePoint(vit, newPoints[pointI++]); +// } +// } +//} + + +void Foam::CV2D::write() const +{ + if (meshControls().objOutput()) + { + writeFaces("allFaces.obj", false); + writeFaces("faces.obj", true); + writeTriangles("allTriangles.obj", false); + writeTriangles("triangles.obj", true); + writePatch("patch.pch"); + } +} + + +void Foam::CV2D::write(const word& stage) const +{ + if (meshControls().objOutput()) + { + Foam::mkDir(stage + "Faces"); + Foam::mkDir(stage + "Triangles"); + + writeFaces + ( + stage + + "Faces/allFaces_" + + runTime_.timeName() + + ".obj", + false + ); + + writeTriangles + ( + stage + + "Triangles/allTriangles_" + + runTime_.timeName() + + ".obj", + false + ); + } +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.H b/applications/utilities/mesh/generation/cv2DMesh/CV2D.H new file mode 100644 index 0000000000..4948196fc2 --- /dev/null +++ b/applications/utilities/mesh/generation/cv2DMesh/CV2D.H @@ -0,0 +1,474 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 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 . + +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 "cv2DControls.H" +#include "tolerances.H" +#include "meshTools.H" +#include "triSurface.H" +#include "searchableSurfaces.H" +#include "conformationSurfaces.H" +#include "relaxationModel.H" +#include "cellSizeAndAlignmentControls.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 + cv2DControls controls_; + + //- The cell size control object + cellSizeAndAlignmentControls cellSizeControl_; + + //- Relaxation coefficient model. Runtime selectable. + autoPtr relaxationModel_; + + //- 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_; + + //- Store the feature points + std::list featurePoints_; + + //- 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 + ); + + //- Insert point and return it's index + inline label insertPoint + ( + const point2D& pt, + const label index, + const label type + ); + + inline label insertPoint + ( + const Point& p, + const label index, + 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(); + + //- Check if a point is within a line. + bool on2DLine(const point2D& p, const linePointRef& line); + + //- Insert point groups at the feature points. + void insertFeaturePoints(); + + //- Re-insert point groups at the feature points. + void reinsertFeaturePoints(); + + //- 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& nearSurfacePoints, + const DynamicList& surfacePoints, + const DynamicList
Delaunay; + #else + // Data structures for standard Delaunay triangulation typedef CGAL::Triangulation_data_structure_2 Tds; typedef CGAL::Delaunay_triangulation_2 Delaunay; + #endif diff --git a/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C b/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C index 98ec79e102..e44a249f28 100644 --- a/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C +++ b/applications/utilities/mesh/generation/cv2DMesh/cv2DMesh.C @@ -163,20 +163,26 @@ int main(int argc, char *argv[]) Info<< "Constructing patches." << endl; List patches(poly2DMesh.patchNames().size()); + label countPatches = 0; forAll(patches, patchI) { - patches[patchI] = new polyPatch - ( - poly2DMesh.patchNames()[patchI], - poly2DMesh.patchSizes()[patchI], - poly2DMesh.patchStarts()[patchI], - patchI, - pMesh.boundaryMesh(), - word::null - ); - } + if (poly2DMesh.patchSizes()[patchI] != 0) + { + patches[countPatches] = new polyPatch + ( + poly2DMesh.patchNames()[patchI], + poly2DMesh.patchSizes()[patchI], + poly2DMesh.patchStarts()[patchI], + countPatches, + pMesh.boundaryMesh(), + word::null + ); + countPatches++; + } + } + patches.setSize(countPatches); pMesh.addPatches(patches); if (extrude) From 45ce8ab68dade441aaa3f14479c72c81a3d982b9 Mon Sep 17 00:00:00 2001 From: laurence Date: Fri, 4 Jan 2013 14:30:07 +0000 Subject: [PATCH 04/16] ENH: Update cellShapeControl to work with 2D mesher and remove dependency on the Triangulation --- .../cellSizeAndAlignmentGrid.C | 183 +++++++++--------- .../cellShapeControl/cellShapeControl.C | 166 +++++++++------- .../cellShapeControl/cellShapeControl.H | 5 +- .../cellSizeAndAlignmentControl.C | 5 +- .../cellSizeAndAlignmentControl.H | 14 +- .../cellSizeAndAlignmentControls.C | 105 +++++++++- .../cellSizeAndAlignmentControls.H | 17 +- .../fileControl/fileControl.C | 8 +- .../fileControl/fileControl.H | 10 +- .../searchableSurfaceControl.C | 11 +- .../searchableSurfaceControl.H | 22 ++- .../cellSizeFunction/cellSizeFunction.C | 5 +- .../cellSizeFunction/cellSizeFunction.H | 9 +- .../conformalVoronoiMesh.C | 29 ++- .../conformalVoronoiMesh.H | 8 +- .../conformalVoronoiMeshCalcDualMesh.C | 4 +- .../conformalVoronoiMeshI.H | 92 ++++++++- .../cvControls/cvControls.C | 6 +- .../cvControls/cvControls.H | 18 +- .../cvControls/cvControlsI.H | 14 +- .../mesh/generation/cvMesh/cvMeshDict | 3 + .../extrude2DMesh/extrude2DMeshApp.C | 14 +- .../mesh/cvMesh/blob/system/collapseDict | 2 +- tutorials/mesh/cvMesh/blob/system/controlDict | 2 +- .../mesh/cvMesh/flange/system/collapseDict | 2 +- .../mesh/cvMesh/flange/system/controlDict | 2 +- .../mesh/cvMesh/flange/system/cvMeshDict | 4 +- 27 files changed, 524 insertions(+), 236 deletions(-) diff --git a/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/cellSizeAndAlignmentGrid.C b/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/cellSizeAndAlignmentGrid.C index 18109761a5..454416299d 100644 --- a/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/cellSizeAndAlignmentGrid.C +++ b/applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/cellSizeAndAlignmentGrid.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -174,12 +174,7 @@ Foam::tmp buildAlignmentField(const T& mesh) continue; } - alignments[vit->index()] = triad - ( - vit->alignment().x(), - vit->alignment().y(), - vit->alignment().z() - ); + alignments[vit->index()] = vit->alignment(); } return tAlignments; @@ -214,6 +209,61 @@ Foam::tmp buildPointField(const T& mesh) } +void refine +( + cellShapeControlMesh& mesh, + const conformationSurfaces& geometryToConformTo, + const label maxRefinementIterations, + const scalar defaultCellSize +) +{ + for (label iter = 0; iter < maxRefinementIterations; ++iter) + { + DynamicList ptsToInsert; + + for + ( + CellSizeDelaunay::Finite_cells_iterator cit = + mesh.finite_cells_begin(); + cit != mesh.finite_cells_end(); + ++cit + ) + { + const point newPoint = + topoint + ( + CGAL::centroid + ( + cit->vertex(0)->point(), + cit->vertex(1)->point(), + cit->vertex(2)->point(), + cit->vertex(3)->point() + ) + ); + + if (geometryToConformTo.inside(newPoint)) + { + ptsToInsert.append(newPoint); + } + } + + Info<< " Adding " << returnReduce(ptsToInsert.size(), sumOp
Delaunay; - -#else - - // Data structures for standard Delaunay triangulation - typedef CGAL::Triangulation_data_structure_2 Tds; - typedef CGAL::Delaunay_triangulation_2 Delaunay; - -#endif - - -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // - -#endif - -// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.C b/applications/utilities/mesh/generation/cv2DMesh/CV2D.C deleted file mode 100644 index 8a5931df61..0000000000 --- a/applications/utilities/mesh/generation/cv2DMesh/CV2D.C +++ /dev/null @@ -1,997 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation - \\/ M anipulation | -------------------------------------------------------------------------------- -License - This file is part of OpenFOAM. - - OpenFOAM is free software: you can redistribute it and/or modify it - under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OpenFOAM is distributed in the hope that it will be useful, but WITHOUT - ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or - FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License - for more details. - - You should have received a copy of the GNU General Public License - along with OpenFOAM. If not, see . - -\*----------------------------------------------------------------------------*/ - -#include "CV2D.H" -#include "Random.H" -#include "transform.H" -#include "IFstream.H" -#include "uint.H" -#include "ulong.H" - -namespace Foam -{ - defineTypeNameAndDebug(CV2D, 0); -} - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -void Foam::CV2D::insertBoundingBox() -{ - Info<< "insertBoundingBox: creating bounding mesh" << endl; - scalar bigSpan = 10*meshControls().span(); - insertPoint(point2D(-bigSpan, -bigSpan), Vb::FAR_POINT); - insertPoint(point2D(-bigSpan, bigSpan), Vb::FAR_POINT); - insertPoint(point2D(bigSpan, -bigSpan), Vb::FAR_POINT); - insertPoint(point2D(bigSpan, bigSpan), Vb::FAR_POINT); -} - - -void Foam::CV2D::fast_restore_Delaunay(Vertex_handle vh) -{ - int i; - Face_handle f = vh->face(), next, start(f); - - do - { - i=f->index(vh); - if (!is_infinite(f)) - { - if (!internal_flip(f, cw(i))) external_flip(f, i); - if (f->neighbor(i) == start) start = f; - } - f = f->neighbor(cw(i)); - } while (f != start); -} - - -void Foam::CV2D::external_flip(Face_handle& f, int i) -{ - Face_handle n = f->neighbor(i); - - if - ( - CGAL::ON_POSITIVE_SIDE - != side_of_oriented_circle(n, f->vertex(i)->point()) - ) return; - - flip(f, i); - i = n->index(f->vertex(i)); - external_flip(n, i); -} - - -bool Foam::CV2D::internal_flip(Face_handle& f, int i) -{ - Face_handle n = f->neighbor(i); - - if - ( - CGAL::ON_POSITIVE_SIDE - != side_of_oriented_circle(n, f->vertex(i)->point()) - ) - { - return false; - } - - flip(f, i); - - return true; -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::CV2D::CV2D -( - const Time& runTime, - const dictionary& cvMeshDict -) -: - Delaunay(), - runTime_(runTime), - rndGen_(64293*Pstream::myProcNo()), - allGeometry_ - ( - IOobject - ( - "cvSearchableSurfaces", - runTime_.constant(), - "triSurface", - runTime_, - IOobject::MUST_READ, - IOobject::NO_WRITE - ), - cvMeshDict.subDict("geometry") - ), - qSurf_ - ( - runTime_, - rndGen_, - allGeometry_, - cvMeshDict.subDict("surfaceConformation") - ), - controls_(cvMeshDict, qSurf_.globalBounds()), - cellSizeControl_ - ( - runTime, - cvMeshDict.subDict("motionControl").subDict("shapeControlFunctions"), - qSurf_, - controls_.minCellSize() - ), - relaxationModel_ - ( - relaxationModel::New - ( - cvMeshDict.subDict("motionControl"), - runTime - ) - ), - z_ - ( - point - ( - cvMeshDict.subDict("surfaceConformation").lookup("locationInMesh") - ).z() - ), - startOfInternalPoints_(0), - startOfSurfacePointPairs_(0), - startOfBoundaryConformPointPairs_(0), - featurePoints_() -{ - Info<< meshControls() << endl; - - insertBoundingBox(); - insertFeaturePoints(); -} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::CV2D::~CV2D() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::CV2D::insertPoints -( - const point2DField& points, - const scalar nearness -) -{ - Info<< "insertInitialPoints(const point2DField& points): "; - - startOfInternalPoints_ = number_of_vertices(); - label nVert = startOfInternalPoints_; - - // Add the points and index them - forAll(points, i) - { - const point2D& p = points[i]; - - if (qSurf_.wellInside(toPoint3D(p), nearness)) - { - insert(toPoint(p))->index() = nVert++; - } - else - { - Warning - << "Rejecting point " << p << " outside surface" << endl; - } - } - - Info<< nVert << " vertices inserted" << endl; - - if (meshControls().objOutput()) - { - // Checking validity of triangulation - assert(is_valid()); - - writeTriangles("initial_triangles.obj", true); - writeFaces("initial_faces.obj", true); - } -} - - -void Foam::CV2D::insertPoints(const fileName& pointFileName) -{ - IFstream pointsFile(pointFileName); - - if (pointsFile.good()) - { - insertPoints - ( - point2DField(pointsFile), - 0.5*meshControls().minCellSize2() - ); - } - else - { - FatalErrorIn("insertInitialPoints") - << "Could not open pointsFile " << pointFileName - << exit(FatalError); - } -} - - -void Foam::CV2D::insertGrid() -{ - Info<< "insertInitialGrid: "; - - startOfInternalPoints_ = number_of_vertices(); - label nVert = startOfInternalPoints_; - - scalar x0 = qSurf_.globalBounds().min().x(); - scalar xR = qSurf_.globalBounds().max().x() - x0; - int ni = int(xR/meshControls().minCellSize()) + 1; - scalar deltax = xR/ni; - - scalar y0 = qSurf_.globalBounds().min().y(); - scalar yR = qSurf_.globalBounds().max().y() - y0; - int nj = int(yR/meshControls().minCellSize()) + 1; - scalar deltay = yR/nj; - - Random rndGen(1321); - scalar pert = meshControls().randomPerturbation()*min(deltax, deltay); - - for (int i=0; iindex() = nVert++; - } - } - } - - Info<< nVert << " vertices inserted" << endl; - - if (meshControls().objOutput()) - { - // Checking validity of triangulation - assert(is_valid()); - - writeTriangles("initial_triangles.obj", true); - writeFaces("initial_faces.obj", true); - } -} - - -void Foam::CV2D::insertSurfacePointPairs() -{ - startOfSurfacePointPairs_ = number_of_vertices(); - - if (meshControls().insertSurfaceNearestPointPairs()) - { - insertSurfaceNearestPointPairs(); - } - - write("nearest"); - - // Insertion of point-pairs for near-points may cause protrusions - // so insertBoundaryConformPointPairs must be executed last - if (meshControls().insertSurfaceNearPointPairs()) - { - insertSurfaceNearPointPairs(); - } - - startOfBoundaryConformPointPairs_ = number_of_vertices(); -} - - -void Foam::CV2D::boundaryConform() -{ - if (!meshControls().insertSurfaceNearestPointPairs()) - { - markNearBoundaryPoints(); - } - - // Mark all the faces as SAVE_CHANGED - for - ( - Triangulation::Finite_faces_iterator fit = finite_faces_begin(); - fit != finite_faces_end(); - fit++ - ) - { - fit->faceIndex() = Fb::SAVE_CHANGED; - } - - for (label iter=1; iter<=meshControls().maxBoundaryConformingIter(); iter++) - { - label nIntersections = insertBoundaryConformPointPairs - ( - "surfaceIntersections_" + Foam::name(iter) + ".obj" - ); - - if (nIntersections == 0) - { - break; - } - else - { - Info<< "BC iteration " << iter << ": " - << nIntersections << " point-pairs inserted" << endl; - } - - // Any faces changed by insertBoundaryConformPointPairs will now - // be marked CHANGED, mark those as SAVE_CHANGED and those that - // remained SAVE_CHANGED as UNCHANGED - for - ( - Triangulation::Finite_faces_iterator fit = finite_faces_begin(); - fit != finite_faces_end(); - fit++ - ) - { - if (fit->faceIndex() == Fb::SAVE_CHANGED) - { - fit->faceIndex() = Fb::UNCHANGED; - } - else if (fit->faceIndex() == Fb::CHANGED) - { - fit->faceIndex() = Fb::SAVE_CHANGED; - } - } - } - - Info<< nl; - - write("boundary"); -} - - -void Foam::CV2D::removeSurfacePointPairs() -{ - for - ( - Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); - vit != finite_vertices_end(); - ++vit - ) - { - if (vit->index() >= startOfSurfacePointPairs_) - { - remove(vit); - } - } -} - - -void Foam::CV2D::newPoints() -{ - const scalar relaxation = relaxationModel_->relaxation(); - - Info<< "Relaxation = " << relaxation << endl; - - Field dualVertices(number_of_faces()); - - label dualVerti = 0; - - // Find the dual point of each tetrahedron and assign it an index. - for - ( - Triangulation::Finite_faces_iterator fit = finite_faces_begin(); - fit != finite_faces_end(); - ++fit - ) - { - fit->faceIndex() = -1; - - if - ( - fit->vertex(0)->internalOrBoundaryPoint() - || fit->vertex(1)->internalOrBoundaryPoint() - || fit->vertex(2)->internalOrBoundaryPoint() - ) - { - fit->faceIndex() = dualVerti; - - dualVertices[dualVerti] = toPoint2D(circumcenter(fit)); - - dualVerti++; - } - } - - dualVertices.setSize(dualVerti); - - Field displacementAccumulator - ( - startOfSurfacePointPairs_, - vector2D::zero - ); - - // Calculate target size and alignment for vertices - scalarField sizes - ( - number_of_vertices(), - meshControls().minCellSize() - ); - - Field alignments - ( - number_of_vertices(), - vector2D(1, 0) - ); - - for - ( - Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); - vit != finite_vertices_end(); - ++vit - ) - { - if (vit->internalOrBoundaryPoint()) - { - point2D vert = toPoint2D(vit->point()); - - // alignment and size determination - pointIndexHit pHit; - label hitSurface = -1; - - qSurf_.findSurfaceNearest - ( - toPoint3D(vert), - meshControls().span2(), - pHit, - hitSurface - ); - - if (pHit.hit()) - { - vectorField norm(1); - allGeometry_[hitSurface].getNormal - ( - List(1, pHit), - norm - ); - - alignments[vit->index()] = toPoint2D(norm[0]); - - sizes[vit->index()] = - cellSizeControl_.cellSize(toPoint3D(vit->point())); - } - } - } - - // Info<< "Calculated alignments" << endl; - - scalar cosAlignmentAcceptanceAngle = 0.68; - - // Upper and lower edge length ratios for weight - scalar u = 1.0; - scalar l = 0.7; - - PackedBoolList pointToBeRetained(startOfSurfacePointPairs_, true); - - std::list pointsToInsert; - - for - ( - Triangulation::Finite_edges_iterator eit = finite_edges_begin(); - eit != finite_edges_end(); - eit++ - ) - { - Vertex_handle vA = eit->first->vertex(cw(eit->second)); - Vertex_handle vB = eit->first->vertex(ccw(eit->second)); - - if (!vA->internalOrBoundaryPoint() || !vB->internalOrBoundaryPoint()) - { - continue; - } - - const point2D& dualV1 = dualVertices[eit->first->faceIndex()]; - const point2D& dualV2 = - dualVertices[eit->first->neighbor(eit->second)->faceIndex()]; - - scalar dualEdgeLength = mag(dualV1 - dualV2); - - point2D dVA = toPoint2D(vA->point()); - point2D dVB = toPoint2D(vB->point()); - - Field alignmentDirsA(2); - - alignmentDirsA[0] = alignments[vA->index()]; - alignmentDirsA[1] = vector2D - ( - -alignmentDirsA[0].y(), - alignmentDirsA[0].x() - ); - - Field alignmentDirsB(2); - - alignmentDirsB[0] = alignments[vB->index()]; - alignmentDirsB[1] = vector2D - ( - -alignmentDirsB[0].y(), - alignmentDirsB[0].x() - ); - - Field alignmentDirs(2); - - forAll(alignmentDirsA, aA) - { - const vector2D& a(alignmentDirsA[aA]); - - scalar maxDotProduct = 0.0; - - forAll(alignmentDirsB, aB) - { - const vector2D& b(alignmentDirsB[aB]); - - scalar dotProduct = a & b; - - if (mag(dotProduct) > maxDotProduct) - { - maxDotProduct = mag(dotProduct); - - alignmentDirs[aA] = a + sign(dotProduct)*b; - - alignmentDirs[aA] /= mag(alignmentDirs[aA]); - } - } - } - - vector2D rAB = dVA - dVB; - - scalar rABMag = mag(rAB); - - forAll(alignmentDirs, aD) - { - vector2D& alignmentDir = alignmentDirs[aD]; - - if ((rAB & alignmentDir) < 0) - { - // swap the direction of the alignment so that has the - // same sense as rAB - alignmentDir *= -1; - } - - scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir); - - if (alignmentDotProd > cosAlignmentAcceptanceAngle) - { - scalar targetFaceSize = - 0.5*(sizes[vA->index()] + sizes[vB->index()]); - - // Test for changing aspect ratio on second alignment (first - // alignment is neartest surface normal) - // if (aD == 1) - // { - // targetFaceSize *= 2.0; - // } - - alignmentDir *= 0.5*targetFaceSize; - - vector2D delta = alignmentDir - 0.5*rAB; - - if (dualEdgeLength < 0.7*targetFaceSize) - { - delta *= 0; - } - else if (dualEdgeLength < targetFaceSize) - { - delta *= - ( - dualEdgeLength - /(targetFaceSize*(u - l)) - - 1/((u/l) - 1) - ); - } - - if - ( - vA->internalPoint() - && vB->internalPoint() - && rABMag > 1.75*targetFaceSize - && dualEdgeLength > 0.05*targetFaceSize - && alignmentDotProd > 0.93 - ) - { - // Point insertion - pointsToInsert.push_back(toPoint(0.5*(dVA + dVB))); - } - else if - ( - (vA->internalPoint() || vB->internalPoint()) - && rABMag < 0.65*targetFaceSize - ) - { - // Point removal - - // Only insert a point at the midpoint of the short edge - // if neither attached point has already been identified - // to be removed. - if - ( - pointToBeRetained[vA->index()] == true - && pointToBeRetained[vB->index()] == true - ) - { - pointsToInsert.push_back(toPoint(0.5*(dVA + dVB))); - } - - if (vA->internalPoint()) - { - pointToBeRetained[vA->index()] = false; - } - - if (vB->internalPoint()) - { - pointToBeRetained[vB->index()] = false; - } - } - else - { - if (vA->internalPoint()) - { - displacementAccumulator[vA->index()] += delta; - } - - if (vB->internalPoint()) - { - displacementAccumulator[vB->index()] += -delta; - } - } - } - } - } - - vector2D totalDisp = sum(displacementAccumulator); - scalar totalDist = sum(mag(displacementAccumulator)); - - // Relax the calculated displacement - displacementAccumulator *= relaxation; - - label numberOfNewPoints = pointsToInsert.size(); - - for - ( - Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); - vit != finite_vertices_end(); - ++vit - ) - { - if (vit->internalPoint()) - { - if (pointToBeRetained[vit->index()]) - { - pointsToInsert.push_front - ( - toPoint - ( - toPoint2D(vit->point()) - + displacementAccumulator[vit->index()] - ) - ); - } - } - } - - // Clear the triangulation and reinsert the bounding box and feature points. - // This is faster than removing and moving points. - this->clear(); - - insertBoundingBox(); - - reinsertFeaturePoints(); - - startOfInternalPoints_ = number_of_vertices(); - - label nVert = startOfInternalPoints_; - - Info<< "Inserting " << numberOfNewPoints << " new points" << endl; - - // Use the range insert as it is faster than individually inserting points. - insert(pointsToInsert.begin(), pointsToInsert.end()); - - for - ( - Delaunay::Finite_vertices_iterator vit = finite_vertices_begin(); - vit != finite_vertices_end(); - ++vit - ) - { - if - ( - vit->type() == Vb::INTERNAL_POINT - && vit->index() == Vb::INTERNAL_POINT - ) - { - vit->index() = nVert++; - } - } - - Info<< " Total displacement = " << totalDisp << nl - << " Total distance = " << totalDist << nl - << " Points added = " << pointsToInsert.size() - << endl; - - write("internal"); - - insertSurfacePointPairs(); - - boundaryConform(); - - -// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -// Old Method -// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -// for -// ( -// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); -// vit != finite_vertices_end(); -// ++vit -// ) -// { -// if (vit->internalPoint()) -// { -// // Current dual-cell defining vertex ("centre") -// point2DFromPoint defVert0 = toPoint2D(vit->point()); - -// Triangulation::Edge_circulator ec = incident_edges(vit); -// Triangulation::Edge_circulator ecStart = ec; - -// // Circulate around the edges to find the first which is not -// // infinite -// do -// { -// if (!is_infinite(ec)) break; -// } while (++ec != ecStart); - -// // Store the start-end of the first non-infinte edge -// point2D de0 = toPoint2D(circumcenter(ec->first)); - -// // Keep track of the maximum edge length^2 -// scalar maxEdgeLen2 = 0.0; - -// // Keep track of the index of the longest edge -// label edgecd0i = -1; - -// // Edge counter -// label edgei = 0; - -// do -// { -// if (!is_infinite(ec)) -// { -// // Get the end of the current edge -// point2D de1 = toPoint2D -// ( -// circumcenter(ec->first->neighbor(ec->second)) -// ); - -// // Store the current edge vector -// edges[edgei] = de1 - de0; - -// // Store the edge mid-point in the vertices array -// vertices[edgei] = 0.5*(de1 + de0); - -// // Move the current edge end into the edge start for the -// // next iteration -// de0 = de1; - -// // Keep track of the longest edge - -// scalar edgeLen2 = magSqr(edges[edgei]); - -// if (edgeLen2 > maxEdgeLen2) -// { -// maxEdgeLen2 = edgeLen2; -// edgecd0i = edgei; -// } - -// edgei++; -// } -// } while (++ec != ecStart); - -// // Initialise cd0 such that the mesh will align -// // in in the x-y directions -// vector2D cd0(1, 0); - -// if (meshControls().relaxOrientation()) -// { -// // Get the longest edge from the array and use as the primary -// // direction of the coordinate system of the "square" cell -// cd0 = edges[edgecd0i]; -// } - -// if (meshControls().nearWallAlignedDist() > 0) -// { -// pointIndexHit pHit = qSurf_.tree().findNearest -// ( -// toPoint3D(defVert0), -// meshControls().nearWallAlignedDist2() -// ); - -// if (pHit.hit()) -// { -// cd0 = toPoint2D(faceNormals[pHit.index()]); -// } -// } - -// // Rotate by 45deg needed to create an averaging procedure which -// // encourages the cells to be square -// cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x()); - -// // Normalise the primary coordinate direction -// cd0 /= mag(cd0); - -// // Calculate the orthogonal coordinate direction -// vector2D cd1(-cd0.y(), cd0.x()); - - -// // Restart the circulator -// ec = ecStart; - -// // ... and the counter -// edgei = 0; - -// // Initialise the displacement for the centre and sum-weights -// vector2D disp = vector2D::zero; -// scalar sumw = 0; - -// do -// { -// if (!is_infinite(ec)) -// { -// // Pick up the current edge -// const vector2D& ei = edges[edgei]; - -// // Calculate the centre to edge-centre vector -// vector2D deltai = vertices[edgei] - defVert0; - -// // Set the weight for this edge contribution -// scalar w = 1; - -// if (meshControls().squares()) -// { -// w = magSqr(deltai.x()*ei.y() - deltai.y()*ei.x()); -// // alternative weights -// //w = mag(deltai.x()*ei.y() - deltai.y()*ei.x()); -// //w = magSqr(ei)*mag(deltai); - -// // Use the following for an ~square mesh -// // Find the coordinate contributions for this edge delta -// scalar cd0deltai = cd0 & deltai; -// scalar cd1deltai = cd1 & deltai; - -// // Create a "square" displacement -// if (mag(cd0deltai) > mag(cd1deltai)) -// { -// disp += (w*cd0deltai)*cd0; -// } -// else -// { -// disp += (w*cd1deltai)*cd1; -// } -// } -// else -// { -// // Use this for a hexagon/pentagon mesh -// disp += w*deltai; -// } - -// // Sum the weights -// sumw += w; -// } -// else -// { -// FatalErrorIn("CV2D::newPoints() const") -// << "Infinite triangle found in internal mesh" -// << exit(FatalError); -// } - -// edgei++; - -// } while (++ec != ecStart); - -// // Calculate the average displacement -// disp /= sumw; -// totalDisp += disp; -// totalDist += mag(disp); - -// // Move the point by a fraction of the average displacement -// movePoint(vit, defVert0 + relaxation*disp); -// } -// } - -// Info << "\nTotal displacement = " << totalDisp -// << " total distance = " << totalDist << endl; -} - - -//void Foam::CV2D::moveInternalPoints(const point2DField& newPoints) -//{ -// label pointI = 0; - -// for -// ( -// Triangulation::Finite_vertices_iterator vit = finite_vertices_begin(); -// vit != finite_vertices_end(); -// ++vit -// ) -// { -// if (vit->internalPoint()) -// { -// movePoint(vit, newPoints[pointI++]); -// } -// } -//} - - -void Foam::CV2D::write() const -{ - if (meshControls().objOutput()) - { - writeFaces("allFaces.obj", false); - writeFaces("faces.obj", true); - writeTriangles("allTriangles.obj", false); - writeTriangles("triangles.obj", true); - writePatch("patch.pch"); - } -} - - -void Foam::CV2D::write(const word& stage) const -{ - if (meshControls().objOutput()) - { - Foam::mkDir(stage + "Faces"); - Foam::mkDir(stage + "Triangles"); - - writeFaces - ( - stage - + "Faces/allFaces_" - + runTime_.timeName() - + ".obj", - false - ); - - writeTriangles - ( - stage - + "Triangles/allTriangles_" - + runTime_.timeName() - + ".obj", - false - ); - } -} - - -// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cv2DMesh/CV2D.H b/applications/utilities/mesh/generation/cv2DMesh/CV2D.H deleted file mode 100644 index 4948196fc2..0000000000 --- a/applications/utilities/mesh/generation/cv2DMesh/CV2D.H +++ /dev/null @@ -1,474 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2013 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 . - -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 "cv2DControls.H" -#include "tolerances.H" -#include "meshTools.H" -#include "triSurface.H" -#include "searchableSurfaces.H" -#include "conformationSurfaces.H" -#include "relaxationModel.H" -#include "cellSizeAndAlignmentControls.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 - cv2DControls controls_; - - //- The cell size control object - cellSizeAndAlignmentControls cellSizeControl_; - - //- Relaxation coefficient model. Runtime selectable. - autoPtr relaxationModel_; - - //- 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_; - - //- Store the feature points - std::list featurePoints_; - - //- 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 - ); - - //- Insert point and return it's index - inline label insertPoint - ( - const point2D& pt, - const label index, - const label type - ); - - inline label insertPoint - ( - const Point& p, - const label index, - 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(); - - //- Check if a point is within a line. - bool on2DLine(const point2D& p, const linePointRef& line); - - //- Insert point groups at the feature points. - void insertFeaturePoints(); - - //- Re-insert point groups at the feature points. - void reinsertFeaturePoints(); - - //- 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& nearSurfacePoints, - const DynamicList& surfacePoints, - const DynamicList