ENH: Latest version of cvMesh. Squash of commits.

applications/utilities/mesh/generation/cvMesh/Allwmake

@@ -4,7 +4,8 @@ set -x
 
 wmake libso conformalVoronoiMesh
 wmake
-wmake cvMeshBackgroundMesh
+#wmake cvMeshBackgroundMesh
 (cd cvMeshSurfaceSimplify && ./Allwmake)
+wmake cellSizeAndAlignmentGrid
 
 # ----------------------------------------------------------------- end-of-file

applications/utilities/mesh/generation/cvMesh/Make/options

@@ -2,7 +2,7 @@ EXE_DEBUG = -DFULLDEBUG -g -O0
 EXE_FROUNDING_MATH = -frounding-math
 EXE_NDEBUG = -DNDEBUG
 
-CGAL_EXACT =
+CGAL_EXACT = /*-DCGAL_DONT_USE_LAZY_KERNEL*/
 CGAL_INEXACT = -DCGAL_INEXACT
 
 include $(GENERAL_RULES)/CGAL
@@ -19,7 +19,9 @@ EXE_INC = \
     -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -IvectorTools
 
 EXE_LIBS = \
     $(CGAL_LIBS) \
@@ -32,4 +34,5 @@ EXE_LIBS = \
     -ledgeMesh \
     -lfileFormats \
     -ltriSurface \
-    -ldynamicMesh
+    -ldynamicMesh \
+    -lsampling

applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/Make/files

+cellSizeAndAlignmentGrid.C
+EXE = $(FOAM_USER_APPBIN)/cellSizeAndAlignmentGrid

applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/Make/options

+EXE_DEBUG = -DFULLDEBUG -g -O0
+EXE_FROUNDING_MATH = -frounding-math
+EXE_NDEBUG = -DNDEBUG
+
+CGAL_EXACT = /*-DCGAL_DONT_USE_LAZY_KERNEL*/
+CGAL_INEXACT = -DCGAL_INEXACT
+
+include $(GENERAL_RULES)/CGAL
+
+
+EXE_INC =  \
+    ${EXE_FROUNDING_MATH} \
+    ${EXE_NDEBUG} \
+    ${CGAL_INEXACT} \
+    ${CGAL_INC} \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(HOME)/OpenFOAM/OpenFOAM-dev/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/lnInclude \
+    -I$(HOME)/OpenFOAM/OpenFOAM-dev/applications/utilities/mesh/generation/cvMesh/vectorTools
+
+EXE_LIBS =  \
+    $(CGAL_LIBS) \
+    -lmpfr \
+    -lboost_thread \
+    -lconformalVoronoiMesh \
+    -lfiniteVolume    \
+    -lmeshTools \
+    -ldecompositionMethods \
+    -L$(FOAM_LIBBIN)/dummy -lptscotchDecomp \
+    -ledgeMesh \
+    -ltriSurface \
+    -ldynamicMesh \
+    -lsampling \
+    -lfileFormats

applications/utilities/mesh/generation/cvMesh/cellSizeAndAlignmentGrid/cellSizeAndAlignmentGrid.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-distributedDelaunayMesh
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "CGALTriangulation3DKernel.H"
+
+#include "indexedVertex.H"
+#include "indexedCell.H"
+
+#include "argList.H"
+#include "Time.H"
+#include "DistributedDelaunayMesh.H"
+#include "backgroundMeshDecomposition.H"
+#include "searchableSurfaces.H"
+#include "conformationSurfaces.H"
+#include "PrintTable.H"
+#include "Random.H"
+#include "boundBox.H"
+#include "point.H"
+#include "cellShapeControlMesh.H"
+#include "triadField.H"
+#include "scalarIOField.H"
+#include "pointIOField.H"
+#include "triadIOField.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+template <class Triangulation, class Type>
+Foam::tmp<Foam::Field<Type> > filterFarPoints
+(
+    const Triangulation& mesh,
+    const Field<Type>& field
+)
+{
+    tmp<Field<Type> > tNewField(new Field<Type>(field.size()));
+    Field<Type>& newField = tNewField();
+
+    label added = 0;
+    label count = 0;
+
+    for
+    (
+        typename Triangulation::Finite_vertices_iterator vit =
+            mesh.finite_vertices_begin();
+        vit != mesh.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->real())
+        {
+            newField[added++] = field[count];
+        }
+
+        count++;
+    }
+
+    newField.resize(added);
+
+    return tNewField;
+}
+
+
+template <class T>
+autoPtr<mapDistribute> buildMap
+(
+    const T& mesh,
+    labelListList& pointPoints
+)
+{
+    pointPoints.setSize(mesh.vertexCount());
+
+    globalIndex globalIndexing(mesh.vertexCount());
+
+    for
+    (
+        typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
+        vit != mesh.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (!vit->real())
+        {
+            continue;
+        }
+
+        std::list<typename T::Vertex_handle> adjVerts;
+        mesh.finite_adjacent_vertices(vit, std::back_inserter(adjVerts));
+
+        DynamicList<label> indices(adjVerts.size());
+
+        for
+        (
+            typename std::list<typename T::Vertex_handle>::const_iterator
+                adjVertI = adjVerts.begin();
+            adjVertI != adjVerts.end();
+            ++adjVertI
+        )
+        {
+            typename T::Vertex_handle vh = *adjVertI;
+
+            if (!vh->farPoint())
+            {
+                indices.append
+                (
+                    globalIndexing.toGlobal(vh->procIndex(), vh->index())
+                );
+            }
+        }
+
+        pointPoints[vit->index()].transfer(indices);
+    }
+
+    List<Map<label> > compactMap;
+
+    return autoPtr<mapDistribute>
+    (
+        new mapDistribute
+        (
+            globalIndexing,
+            pointPoints,
+            compactMap
+        )
+    );
+}
+
+
+template <class T>
+Foam::tmp<Foam::triadField> buildAlignmentField(const T& mesh)
+{
+    tmp<triadField> tAlignments
+    (
+        new triadField(mesh.vertexCount(), triad::unset)
+    );
+    triadField& alignments = tAlignments();
+
+    for
+    (
+        typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
+        vit != mesh.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (!vit->real())
+        {
+            continue;
+        }
+
+        alignments[vit->index()] = triad
+        (
+            vit->alignment().x(),
+            vit->alignment().y(),
+            vit->alignment().z()
+        );
+    }
+
+    return tAlignments;
+}
+
+
+template <class T>
+Foam::tmp<Foam::pointField> buildPointField(const T& mesh)
+{
+    tmp<pointField> tPoints
+    (
+        new pointField(mesh.vertexCount(), point(GREAT, GREAT, GREAT))
+    );
+    pointField& points = tPoints();
+
+    for
+    (
+        typename T::Finite_vertices_iterator vit = mesh.finite_vertices_begin();
+        vit != mesh.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (!vit->real())
+        {
+            continue;
+        }
+
+        points[vit->index()] = topoint(vit->point());
+    }
+
+    return tPoints;
+}
+
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    label maxRefinementIterations = 0;
+    label maxSmoothingIterations = 200;
+    scalar minResidual = 0;
+    scalar defaultCellSize = 0.0004;
+    scalar nearFeatDistSqrCoeff = 1e-8;
+
+
+    // Need to decouple vertex and cell type from this class?
+    // Vertex must have:
+    // + index
+    // + procIndex
+    // - type should be optional
+    cellShapeControlMesh mesh(runTime);
+
+    IOdictionary cvMeshDict
+    (
+        IOobject
+        (
+            "cvMeshDict",
+            runTime.system(),
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    Random rndGen(64293*Pstream::myProcNo());
+
+    searchableSurfaces allGeometry
+    (
+        IOobject
+        (
+            "cvSearchableSurfaces",
+            runTime.constant(),
+            "triSurface",
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        cvMeshDict.subDict("geometry")
+    );
+
+    conformationSurfaces geometryToConformTo
+    (
+        runTime,
+        rndGen,
+        allGeometry,
+        cvMeshDict.subDict("surfaceConformation")
+    );
+
+    autoPtr<backgroundMeshDecomposition> bMesh;
+    if (Pstream::parRun())
+    {
+        bMesh.set
+        (
+            new backgroundMeshDecomposition
+            (
+                runTime,
+                rndGen,
+                geometryToConformTo,
+                cvMeshDict.subDict("backgroundMeshDecomposition")
+            )
+        );
+    }
+
+    // Nice to have IO for the delaunay mesh
+    // IO depend on vertex type.
+    //
+    // Define a delaunay mesh as:
+    // + list of points of the triangulation
+    // + optionally a list of cells
+
+    Info<< nl << "Loop over surfaces" << endl;
+
+    forAll(geometryToConformTo.surfaces(), sI)
+    {
+        const label surfI = geometryToConformTo.surfaces()[sI];
+
+        const searchableSurface& surface =
+            geometryToConformTo.geometry()[surfI];
+
+        Info<< nl << "Inserting points from surface " << surface.name()
+            << " (" << surface.type() << ")" << endl;
+
+        const tmp<pointField> tpoints = surface.points();
+        const pointField& points = tpoints();
+
+        Info<< "    Number of points = " << points.size() << endl;
+
+        forAll(points, pI)
+        {
+            // Is the point in the extendedFeatureEdgeMesh? If so get the
+            // point normal, otherwise get the surface normal from
+            // searchableSurface
+
+            pointIndexHit info;
+            label infoFeature;
+            geometryToConformTo.findFeaturePointNearest
+            (
+                points[pI],
+                nearFeatDistSqrCoeff,
+                info,
+                infoFeature
+            );
+
+
+            autoPtr<triad> pointAlignment;
+
+            if (info.hit())
+            {
+                const extendedFeatureEdgeMesh& features =
+                    geometryToConformTo.features()[infoFeature];
+
+                vectorField norms = features.featurePointNormals(info.index());
+
+                // Create a triad from these norms.
+                pointAlignment.set(new triad());
+                forAll(norms, nI)
+                {
+                    pointAlignment() += norms[nI];
+                }
+
+                pointAlignment().normalize();
+                pointAlignment().orthogonalize();
+            }
+            else
+            {
+                geometryToConformTo.findEdgeNearest
+                (
+                    points[pI],
+                    nearFeatDistSqrCoeff,
+                    info,
+                    infoFeature
+                );
+
+                if (info.hit())
+                {
+                    const extendedFeatureEdgeMesh& features =
+                        geometryToConformTo.features()[infoFeature];
+
+                    vectorField norms = features.edgeNormals(info.index());
+
+                    // Create a triad from these norms.
+                    pointAlignment.set(new triad());
+                    forAll(norms, nI)
+                    {
+                        pointAlignment() += norms[nI];
+                    }
+
+                    pointAlignment().normalize();
+                    pointAlignment().orthogonalize();
+                }
+                else
+                {
+                    pointField ptField(1, points[pI]);
+                    scalarField distField(1, nearFeatDistSqrCoeff);
+                    List<pointIndexHit> infoList(1, pointIndexHit());
+
+                    surface.findNearest(ptField, distField, infoList);
+
+                    vectorField normals(1);
+                    surface.getNormal(infoList, normals);
+
+                    pointAlignment.set(new triad(normals[0]));
+                }
+            }
+
+            if (Pstream::parRun())
+            {
+                if (bMesh().positionOnThisProcessor(points[pI]))
+                {
+                    CellSizeDelaunay::Vertex_handle vh = mesh.insert
+                    (
+                        points[pI],
+                        defaultCellSize,
+                        pointAlignment()
+                    );
+                }
+            }
+            else
+            {
+                CellSizeDelaunay::Vertex_handle vh = mesh.insert
+                (
+                    points[pI],
+                    defaultCellSize,
+                    pointAlignment()
+                );
+            }
+        }
+    }
+
+
+    for (label iter = 0; iter < maxRefinementIterations; ++iter)
+    {
+        DynamicList<point> 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<label>())
+            << endl;
+
+        forAll(ptsToInsert, ptI)
+        {
+            mesh.insert
+            (
+                ptsToInsert[ptI],
+                defaultCellSize,
+                triad::unset
+            );
+        }
+    }
+
+
+    if (Pstream::parRun())
+    {
+        mesh.distribute(bMesh);
+    }
+
+    labelListList pointPoints;
+    autoPtr<mapDistribute> meshDistributor = buildMap(mesh, pointPoints);
+
+    triadField alignments = buildAlignmentField(mesh);
+    pointField points = buildPointField(mesh);
+
+    mesh.printInfo(Info);
+
+
+    // Setup the sizes and alignments on each point
+    triadField fixedAlignments(mesh.vertexCount(), triad::unset);
+
+    for
+    (
+        CellSizeDelaunay::Finite_vertices_iterator vit =
+            mesh.finite_vertices_begin();
+        vit != mesh.finite_vertices_end();
+        ++vit
+    )
+    {
+        if (vit->real())
+        {
+            const tensor& alignment = vit->alignment();
+
+            fixedAlignments[vit->index()] = triad
+            (
+                alignment.x(),
+                alignment.y(),
+                alignment.z()
+            );
+        }
+    }
+
+