diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C index bbe3e1f3ae..9663ae8217 100644 --- a/applications/test/Matrix/Test-Matrix.C +++ b/applications/test/Matrix/Test-Matrix.C @@ -25,6 +25,7 @@ License #include "scalarMatrices.H" #include "vector.H" +#include "IFstream.H" using namespace Foam; diff --git a/applications/test/vectorTools/Make/files b/applications/test/vectorTools/Make/files new file mode 100644 index 0000000000..0b30b98f8f --- /dev/null +++ b/applications/test/vectorTools/Make/files @@ -0,0 +1,3 @@ +Test-vectorTools.C + +EXE = $(FOAM_USER_APPBIN)/Test-vectorTools diff --git a/applications/test/vectorTools/Make/options b/applications/test/vectorTools/Make/options new file mode 100644 index 0000000000..9e015e6078 --- /dev/null +++ b/applications/test/vectorTools/Make/options @@ -0,0 +1 @@ +EXE_INC = -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/vectorTools diff --git a/applications/test/vectorTools/Test-vectorTools.C b/applications/test/vectorTools/Test-vectorTools.C new file mode 100644 index 0000000000..85d18ed989 --- /dev/null +++ b/applications/test/vectorTools/Test-vectorTools.C @@ -0,0 +1,71 @@ +#include "vector.H" +#include "IOstreams.H" +#include "vectorTools.H" +#include "unitConversion.H" + +using namespace Foam; + + +void test(const vector& a, const vector& b, const scalar tolerance) +{ + Info<< "Vectors " << a << " and " << b + << " are (to tolerance of " << tolerance << "): "; + + if (vectorTools::areParallel(a, b, tolerance)) + Info<< " parallel "; + + if (vectorTools::areOrthogonal(a, b, tolerance)) + Info<< " orthogonal "; + + if (vectorTools::areAcute(a, b)) + Info<< " acute "; + + if (vectorTools::areObtuse(a, b)) + Info<< " obtuse "; + + Info<< ", angle = " << vectorTools::degAngleBetween(a, b); + + Info<< endl; +} + + +int main() +{ + vector a(1.0, 1.0, 1.0); + vector b(2.0, 2.0, 2.0); + + test(a, b, 0.0); + test(a, b, VSMALL); + test(a, b, SMALL); + test(a, b, 1e-3); + test(a, b, 1e-1); + + a = vector(1,0,0); + b = vector(0,2,0); + + test(a, b, 0.0); + test(a, b, VSMALL); + test(a, b, SMALL); + test(a, b, 1e-3); + test(a, b, 1e-1); + + a = vector(1,0,0); + b = vector(-1,0,0); + + test(a, b, 0.0); + test(a, b, VSMALL); + test(a, b, SMALL); + test(a, b, 1e-3); + test(a, b, 1e-1); + + a = vector(1,0,0); + b = vector(-1,2,0); + + test(a, b, 0.0); + test(a, b, VSMALL); + test(a, b, SMALL); + test(a, b, 1e-3); + test(a, b, 1e-1); + + return 0; +} diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseDict b/applications/utilities/mesh/advanced/collapseEdges/collapseDict index 170a0a890d..e89aa85241 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/collapseDict +++ b/applications/utilities/mesh/advanced/collapseEdges/collapseDict @@ -2,16 +2,24 @@ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: dev | -| \\ / A nd | Web: www.OpenFOAM.org | +| \\ / A nd | Web: http://www.openfoam.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ + FoamFile { - version 2.0; - format ascii; - class dictionary; - object collapseDict; + version 2.0; + format ascii; + + root ""; + case ""; + instance ""; + local ""; + + class dictionary; + object collapseDict; } + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // If on, after collapsing check the quality of the mesh. If bad faces are @@ -72,7 +80,7 @@ controlMeshQualityCoeffs // The amount that initialFaceLengthFactor will be reduced by for each // face if its collapse generates a poor quality face - faceReductionFactor $initialFaceLengthFactor; + faceReductionFactor 0.5; // Maximum number of smoothing iterations for the reductionFactors maximumSmoothingIterations 2; diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C index c2c9d0a0d3..7a07b63d8f 100644 --- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C +++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C @@ -72,10 +72,11 @@ int main(int argc, char *argv[]) "Collapse small and sliver faces as well as small edges" ); - argList::addBoolOption + argList::addOption ( - "collapseIndirectPatchFaces", - "Collapse faces that are in the face zone indirectPatchFaces" + "collapseFaceZone", + "zoneName", + "Collapse faces that are in the supplied face zone" ); # include "addOverwriteOption.H" @@ -92,8 +93,7 @@ int main(int argc, char *argv[]) const bool overwrite = args.optionFound("overwrite"); const bool collapseFaces = args.optionFound("collapseFaces"); - const bool collapseIndirectPatchFaces = - args.optionFound("collapseIndirectPatchFaces"); + const bool collapseFaceZone = args.optionFound("collapseFaceZone"); forAll(timeDirs, timeI) { @@ -115,11 +115,15 @@ int main(int argc, char *argv[]) meshMod.changeMesh(mesh, false); } - if (collapseIndirectPatchFaces) + if (collapseFaceZone) { + const word faceZoneName = args.optionRead("collapseFaceZone"); + + const faceZone& fZone = mesh.faceZones()[faceZoneName]; + // Filter faces. Pass in the number of bad faces that are present // from the previous edge filtering to use as a stopping criterion. - meshFilter.filterIndirectPatchFaces(); + meshFilter.filterFaceZone(fZone); { polyTopoChange meshMod(newMesh); diff --git a/applications/utilities/mesh/generation/Allwmake b/applications/utilities/mesh/generation/Allwmake index d557fa34f8..55623189f5 100755 --- a/applications/utilities/mesh/generation/Allwmake +++ b/applications/utilities/mesh/generation/Allwmake @@ -9,5 +9,10 @@ extrude2DMesh/Allwmake wmake snappyHexMesh +if [ -d "$CGAL_ARCH_PATH" ] +then + foamyHexMesh/Allwmake + foamyHex2DMesh/Allwmake +fi # ----------------------------------------------------------------- end-of-file diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C index ff6f636cfd..178df8fa4a 100644 --- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C +++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.C @@ -90,6 +90,7 @@ Foam::extrude2DMesh::extrude2DMesh dict_(dict), //patchDict_(dict.subDict("patchInfo")), model_(model), + modelType_(dict.lookup("extrudeModel")), patchType_(dict.lookup("patchType")), frontPatchI_(-1), backPatchI_(-1) diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H index 62a15a44b8..20b07dae87 100644 --- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H +++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMesh/extrude2DMesh/extrude2DMesh.H @@ -53,9 +53,10 @@ class polyMesh; class polyTopoChange; class mapPolyMesh; class mapDistributePolyMesh; +class polyBoundaryMesh; /*---------------------------------------------------------------------------*\ - Class extrude2DMesh Declaration + Class extrude2DMesh Declaration \*---------------------------------------------------------------------------*/ class extrude2DMesh @@ -65,24 +66,19 @@ class extrude2DMesh //- Reference to 2D mesh polyMesh& mesh_; - //- Extrusion dictionary const dictionary dict_; //const dictionary patchDict_; - //- The extrusion model const extrudeModel& model_; - //- Type of the patches that will be created by the extrusion. + const word modelType_; + const word patchType_; - //- Patch ID of the front patch label frontPatchI_; - - //- Patch ID of the back patch label backPatchI_; - // Private Member Functions //- Check the mesh is 2D @@ -97,7 +93,6 @@ class extrude2DMesh //- Disallow default bitwise assignment void operator=(const extrude2DMesh&); - public: //- Runtime type information @@ -105,8 +100,6 @@ public: // Constructors - - //- Construct from a 2D polyMesh, a dictionary and an extrusion model extrude2DMesh ( polyMesh&, @@ -121,36 +114,30 @@ public: // Member Functions - // Access + //- Add front and back patches + void addFrontBackPatches(); - //- Return the patch ID of the front patch - inline label frontPatchI() const - { - return frontPatchI_; - } + //- Play commands into polyTopoChange to extrude mesh. + void setRefinement(polyTopoChange&); - //- Return the patch ID of the back patch - inline label backPatchI() const - { - return backPatchI_; - } + //- Force recalculation of locally stored data on topological change + void updateMesh(const mapPolyMesh&) + {} + //- Force recalculation of locally stored data for mesh distribution + void distribute(const mapDistributePolyMesh&) + {} - // Edit + label frontPatchI() const + { + return frontPatchI_; + } - //- Add front and back patches - void addFrontBackPatches(); + label backPatchI() const + { + return backPatchI_; + } - //- Play commands into polyTopoChange to extrude mesh. - void setRefinement(polyTopoChange&); - - //- Force recalculation of locally stored data on topological change - void updateMesh(const mapPolyMesh&) - {} - - //- Force recalculation of locally stored data for mesh distribution - void distribute(const mapDistributePolyMesh&) - {} }; diff --git a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshDict b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshDict index fa1b37b304..5e53c29c88 100644 --- a/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshDict +++ b/applications/utilities/mesh/generation/extrude2DMesh/extrude2DMeshDict @@ -5,48 +5,42 @@ | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ + FoamFile { - version 2.0; - format ascii; - class dictionary; - object extrude2DMeshDict; + version 2.0; + format ascii; + + root ""; + case ""; + instance ""; + local ""; + + class dictionary; + object extrude2DMeshDict; } + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Type of extrusion extrudeModel linearDirection; //extrudeModel wedge; -// Patch type the extruded patches will take patchType empty; //patchType wedge; -// Number of layers to extrude nLayers 1; -// Expansion ratio. If >1 then grows the layers expansionRatio 1.0; linearDirectionCoeffs { - // Direction of extrusion direction (0 0 1); - - // Width of newly extruded cells thickness 0.1; } wedgeCoeffs { - // Point the extrusion axis goes through - axisPt (0 0 0); - - // Axis to extrude around - axis (1 0 0); - - // Total angle of the wedge in degrees - angle 10; + axisPt (0 0 0); + axis (1 0 0); + angle 10; } - -// ************************************************************************* // \ No newline at end of file diff --git a/applications/utilities/mesh/generation/foamyHex2DMesh/Allwclean b/applications/utilities/mesh/generation/foamyHex2DMesh/Allwclean new file mode 100755 index 0000000000..d0ae53e415 --- /dev/null +++ b/applications/utilities/mesh/generation/foamyHex2DMesh/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/foamyHex2DMesh/Allwmake b/applications/utilities/mesh/generation/foamyHex2DMesh/Allwmake new file mode 100755 index 0000000000..5486849957 --- /dev/null +++ b/applications/utilities/mesh/generation/foamyHex2DMesh/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/foamyHex2DMesh/CGALTriangulation2DKernel.H b/applications/utilities/mesh/generation/foamyHex2DMesh/CGALTriangulation2DKernel.H new file mode 100644 index 0000000000..28002f962f --- /dev/null +++ b/applications/utilities/mesh/generation/foamyHex2DMesh/CGALTriangulation2DKernel.H @@ -0,0 +1,60 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + CGALTriangulation2DKernel + +Description + +\*---------------------------------------------------------------------------*/ + +#ifndef CGALTriangulation2DKernel_H +#define CGALTriangulation2DKernel_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#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 + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/foamyHex2DMesh/CGALTriangulation2Ddefs.H b/applications/utilities/mesh/generation/foamyHex2DMesh/CGALTriangulation2Ddefs.H new file mode 100644 index 0000000000..8731926548 --- /dev/null +++ b/applications/utilities/mesh/generation/foamyHex2DMesh/CGALTriangulation2Ddefs.H @@ -0,0 +1,76 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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/foamyHex2DMesh/CV2D.C b/applications/utilities/mesh/generation/foamyHex2DMesh/CV2D.C new file mode 100644 index 0000000000..deb0605aa0 --- /dev/null +++ b/applications/utilities/mesh/generation/foamyHex2DMesh/CV2D.C @@ -0,0 +1,1000 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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/foamyHex2DMesh/CV2D.H b/applications/utilities/mesh/generation/foamyHex2DMesh/CV2D.H new file mode 100644 index 0000000000..4948196fc2 --- /dev/null +++ b/applications/utilities/mesh/generation/foamyHex2DMesh/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