From ba05ddc58c2cc4a1168f8dac9eeb0f8e00859fd9 Mon Sep 17 00:00:00 2001 From: Graham Date: Tue, 17 Mar 2009 10:13:46 +0000 Subject: [PATCH 001/484] Reintroducing CV meshers after removal by a commit in master, originating in the dsmc branch. --- .../CV2DMesher/CGALTriangulation2Ddefs.H | 84 + .../mesh/generation/CV2DMesher/CV2D.C | 567 ++++++ .../mesh/generation/CV2DMesher/CV2D.H | 514 +++++ .../mesh/generation/CV2DMesher/CV2DI.H | 169 ++ .../mesh/generation/CV2DMesher/CV2DIO.C | 285 +++ .../mesh/generation/CV2DMesher/CV2DMesher.C | 116 ++ .../mesh/generation/CV2DMesher/Make/files | 14 + .../mesh/generation/CV2DMesher/Make/options | 16 + .../mesh/generation/CV2DMesher/controls.C | 76 + .../mesh/generation/CV2DMesher/indexedFace.H | 148 ++ .../generation/CV2DMesher/indexedVertex.H | 261 +++ .../insertBoundaryConformPointPairs.C | 290 +++ .../CV2DMesher/insertFeaturePoints.C | 162 ++ .../CV2DMesher/insertSurfaceNearPointPairs.C | 103 + .../insertSurfaceNearestPointPairs.C | 223 +++ .../mesh/generation/CV2DMesher/querySurface.C | 229 +++ .../mesh/generation/CV2DMesher/querySurface.H | 149 ++ .../mesh/generation/CV2DMesher/tolerances.C | 62 + .../CV3DMesher/CGALTriangulation3Ddefs.H | 84 + .../mesh/generation/CV3DMesher/CV3D.C | 1714 +++++++++++++++++ .../mesh/generation/CV3DMesher/CV3D.H | 460 +++++ .../mesh/generation/CV3DMesher/CV3DI.H | 177 ++ .../mesh/generation/CV3DMesher/CV3DIO.C | 244 +++ .../mesh/generation/CV3DMesher/CV3DMesher.C | 130 ++ .../mesh/generation/CV3DMesher/Make/files | 15 + .../mesh/generation/CV3DMesher/Make/options | 21 + ...0080723_edgeAimingFeatureReconstuction.tar | Bin 0 -> 3441 bytes .../20080808_3EdgeSpecificImplementation.tar | Bin 0 -> 4814 bytes .../backup/20080820_featurePointsDone.tar | Bin 0 -> 16761 bytes .../hardCodedSimpleCubeForPolyTopoChange.H | 96 + .../indexedVertex_with_displacementSum.H | 289 +++ .../mesh/generation/CV3DMesher/calcDualMesh.C | 642 ++++++ .../mesh/generation/CV3DMesher/controls.C | 68 + .../mesh/generation/CV3DMesher/indexedCell.H | 154 ++ .../generation/CV3DMesher/indexedVertex.H | 315 +++ .../insertBoundaryConformPointPairs.C | 56 + .../CV3DMesher/insertFeaturePoints.C | 615 ++++++ .../CV3DMesher/insertSurfaceNearPointPairs.C | 45 + .../insertSurfaceNearestPointPairs.C | 988 ++++++++++ .../mesh/generation/CV3DMesher/querySurface.C | 248 +++ .../mesh/generation/CV3DMesher/querySurface.H | 162 ++ .../mesh/generation/CV3DMesher/tolerances.C | 77 + 42 files changed, 10068 insertions(+) create mode 100644 applications/utilities/mesh/generation/CV2DMesher/CGALTriangulation2Ddefs.H create mode 100644 applications/utilities/mesh/generation/CV2DMesher/CV2D.C create mode 100644 applications/utilities/mesh/generation/CV2DMesher/CV2D.H create mode 100644 applications/utilities/mesh/generation/CV2DMesher/CV2DI.H create mode 100644 applications/utilities/mesh/generation/CV2DMesher/CV2DIO.C create mode 100644 applications/utilities/mesh/generation/CV2DMesher/CV2DMesher.C create mode 100755 applications/utilities/mesh/generation/CV2DMesher/Make/files create mode 100755 applications/utilities/mesh/generation/CV2DMesher/Make/options create mode 100644 applications/utilities/mesh/generation/CV2DMesher/controls.C create mode 100644 applications/utilities/mesh/generation/CV2DMesher/indexedFace.H create mode 100644 applications/utilities/mesh/generation/CV2DMesher/indexedVertex.H create mode 100644 applications/utilities/mesh/generation/CV2DMesher/insertBoundaryConformPointPairs.C create mode 100644 applications/utilities/mesh/generation/CV2DMesher/insertFeaturePoints.C create mode 100644 applications/utilities/mesh/generation/CV2DMesher/insertSurfaceNearPointPairs.C create mode 100644 applications/utilities/mesh/generation/CV2DMesher/insertSurfaceNearestPointPairs.C create mode 100644 applications/utilities/mesh/generation/CV2DMesher/querySurface.C create mode 100644 applications/utilities/mesh/generation/CV2DMesher/querySurface.H create mode 100644 applications/utilities/mesh/generation/CV2DMesher/tolerances.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/CGALTriangulation3Ddefs.H create mode 100644 applications/utilities/mesh/generation/CV3DMesher/CV3D.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/CV3D.H create mode 100644 applications/utilities/mesh/generation/CV3DMesher/CV3DI.H create mode 100644 applications/utilities/mesh/generation/CV3DMesher/CV3DIO.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/CV3DMesher.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/Make/files create mode 100644 applications/utilities/mesh/generation/CV3DMesher/Make/options create mode 100644 applications/utilities/mesh/generation/CV3DMesher/backup/20080723_edgeAimingFeatureReconstuction.tar create mode 100644 applications/utilities/mesh/generation/CV3DMesher/backup/20080808_3EdgeSpecificImplementation.tar create mode 100644 applications/utilities/mesh/generation/CV3DMesher/backup/20080820_featurePointsDone.tar create mode 100644 applications/utilities/mesh/generation/CV3DMesher/backup/hardCodedSimpleCubeForPolyTopoChange.H create mode 100644 applications/utilities/mesh/generation/CV3DMesher/backup/indexedVertex_with_displacementSum.H create mode 100644 applications/utilities/mesh/generation/CV3DMesher/calcDualMesh.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/controls.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/indexedCell.H create mode 100644 applications/utilities/mesh/generation/CV3DMesher/indexedVertex.H create mode 100644 applications/utilities/mesh/generation/CV3DMesher/insertBoundaryConformPointPairs.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/insertFeaturePoints.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearPointPairs.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/insertSurfaceNearestPointPairs.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/querySurface.C create mode 100644 applications/utilities/mesh/generation/CV3DMesher/querySurface.H create mode 100644 applications/utilities/mesh/generation/CV3DMesher/tolerances.C diff --git a/applications/utilities/mesh/generation/CV2DMesher/CGALTriangulation2Ddefs.H b/applications/utilities/mesh/generation/CV2DMesher/CGALTriangulation2Ddefs.H new file mode 100644 index 0000000000..0ecbed8082 --- /dev/null +++ b/applications/utilities/mesh/generation/CV2DMesher/CGALTriangulation2Ddefs.H @@ -0,0 +1,84 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Typedefs + CGALTriangulation2Ddefs + +Description + CGAL data structures used for 2D Delaunay meshing. + + Define CGAL_INEXACT to use Exact_predicates_inexact_constructions kernel + otherwise the more robust but much less efficient + Exact_predicates_exact_constructions will be used. + + Define CGAL_HIERARCHY to use hierarchical Delaunay triangulation which is + faster but uses more memory than the standard Delaunay triangulation. + +\*---------------------------------------------------------------------------*/ + +#ifndef CGALTriangulation2Ddefs_H +#define CGALTriangulation2Ddefs_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "CGAL/Delaunay_triangulation_2.h" + +#include "indexedVertex.H" +#include "indexedFace.H" + +#ifdef CGAL_INEXACT + // Fast kernel using a double as the storage type but the triangulation + // may fail + #include "CGAL/Exact_predicates_inexact_constructions_kernel.h" + typedef CGAL::Exact_predicates_inexact_constructions_kernel K; +#else + // Very robust but expensive kernel + #include "CGAL/Exact_predicates_exact_constructions_kernel.h" + typedef CGAL::Exact_predicates_exact_constructions_kernel K; +#endif + +typedef CGAL::indexedVertex Vb; +typedef CGAL::indexedFace Fb; + +#ifdef CGAL_HIERARCHY + // Data structures for hierarchical Delaunay triangulation which is more + // efficient but also uses more storage + #include "CGAL/Triangulation_hierarchy_2.h" + typedef CGAL::Triangulation_hierarchy_vertex_base_2 Vbh; + typedef CGAL::Triangulation_data_structure_2 Tds; + typedef CGAL::Delaunay_triangulation_2 Triangulation; + typedef CGAL::Triangulation_hierarchy_2 HTriangulation; +#else + // Data structures for standard Delaunay triangulation + typedef CGAL::Triangulation_data_structure_2 Tds; + typedef CGAL::Delaunay_triangulation_2 Triangulation; + typedef Triangulation HTriangulation; +#endif + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/CV2DMesher/CV2D.C b/applications/utilities/mesh/generation/CV2DMesher/CV2D.C new file mode 100644 index 0000000000..e03fccbf47 --- /dev/null +++ b/applications/utilities/mesh/generation/CV2DMesher/CV2D.C @@ -0,0 +1,567 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*----------------------------------------------------------------------------*/ + +#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 dictionary& controlDict, + const querySurface& qSurf +) +: + HTriangulation(), + qSurf_(qSurf), + controls_(controlDict), + tols_(controlDict, controls_.minCellSize, qSurf.bb()), + z_((1.0/3.0)*(qSurf_.bb().min().z() + qSurf_.bb().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_.bb().min().x(); + scalar xR = qSurf_.bb().max().x() - x0; + int ni = int(xR/controls_.minCellSize) + 1; + scalar deltax = xR/ni; + + scalar y0 = qSurf_.bb().min().y(); + scalar yR = qSurf_.bb().max().y() - y0; + int nj = int(yR/controls_.minCellSize) + 1; + scalar deltay = yR/nj; + + Random rndGen(1321); + scalar pert = controls_.randomPurturbation*min(deltax, deltay); + + for (int i=0; iindex() = nVert++; + } + } + } + + Info<< nVert << " vertices inserted" << endl; + + if (controls_.writeInitialTriangulation) + { + // Checking validity of triangulation + assert(is_valid()); + + writeTriangles("initial_triangles.obj", true); + writeFaces("initial_faces.obj", true); + } +} + + +void Foam::CV2D::insertSurfacePointPairs() +{ + startOfSurfacePointPairs_ = number_of_vertices(); + + if (controls_.insertSurfaceNearestPointPairs) + { + insertSurfaceNearestPointPairs(); + } + + if (controls_.writeNearestTriangulation) + { + writeFaces("near_allFaces.obj", false); + writeFaces("near_faces.obj", true); + writeTriangles("near_triangles.obj", true); + } + + // Insertion of point-pais 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: "; + + const vectorField& faceNormals = qSurf_.faceNormals(); + + // Initialise the total displacement and its distance for writing out + vector2D totalDisp = vector2D::zero; + scalar totalDist = 0; + + 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::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/CV2DMesher/CV2D.H b/applications/utilities/mesh/generation/CV2DMesher/CV2D.H new file mode 100644 index 0000000000..171623fb04 --- /dev/null +++ b/applications/utilities/mesh/generation/CV2DMesher/CV2D.H @@ -0,0 +1,514 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +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 "querySurface.H" +#include "point2DFieldFwd.H" +#include "dictionary.H" +#include "Switch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class CV2D Declaration +\*---------------------------------------------------------------------------*/ + +class CV2D +: + public HTriangulation +{ +public: + + class 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; + + controls(const dictionary& controlDict); + }; + + + class 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; + + tolerances + ( + const dictionary& controlDict, + scalar minCellSize, + const boundBox& + ); + }; + + +private: + + // Private data + + //- The surface to mesh + const querySurface& qSurf_; + + //- Meshing controls + controls controls_; + + //- 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& nearSurfacePoints, + const DynamicList& surfacePoints, + const DynamicList