conformalVoronoiMesh.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2012-2016 OpenFOAM Foundation
    Copyright (C) 2020-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
\*---------------------------------------------------------------------------*/
 
#include "conformalVoronoiMesh.H"
#include "initialPointsMethod.H"
#include "relaxationModel.H"
#include "faceAreaWeightModel.H"
#include "meshSearch.H"
#include "vectorTools.H"
#include "IOmanip.H"
#include "indexedCellChecks.H"
#include "controlMeshRefinement.H"
#include "smoothAlignmentSolver.H"
#include "OBJstream.H"
#include "indexedVertexOps.H"
#include "DelaunayMeshTools.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    defineTypeNameAndDebug(conformalVoronoiMesh, 0);
}
 
const Foam::Enum
<
    Foam::conformalVoronoiMesh::dualMeshPointType
>
Foam::conformalVoronoiMesh::dualMeshPointTypeNames_
({
    { dualMeshPointType::internal, "internal" },
    { dualMeshPointType::surface, "surface" },
    { dualMeshPointType::featureEdge, "featureEdge" },
    { dualMeshPointType::featurePoint, "featurePoint" },
    { dualMeshPointType::constrained, "constrained" },
});
 
 
// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
void Foam::conformalVoronoiMesh::cellSizeMeshOverlapsBackground() const
{
    const cellShapeControlMesh& cellSizeMesh =
        cellShapeControl_.shapeControlMesh();
 
    DynamicList<Foam::point> pts(number_of_vertices());
 
    for
    (
        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
        vit != finite_vertices_end();
        ++vit
    )
    {
        if (vit->internalOrBoundaryPoint() && !vit->referred())
        {
            pts.append(topoint(vit->point()));
        }
    }
 
    boundBox bb(pts);
 
    boundBox cellSizeMeshBb = cellSizeMesh.bounds();
 
    bool fullyContained = cellSizeMeshBb.contains(bb);
 
    if (!fullyContained)
    {
        Pout<< "Triangulation not fully contained in cell size mesh." << endl
            << "Cell Size Mesh Bounds = " << cellSizeMeshBb << endl
            << "foamyHexMesh Bounds         = " << bb << endl;
    }
 
    Info<< "Triangulation is "
        << (returnReduceAnd(fullyContained) ? "fully" : "not fully")
        << " contained in the cell size mesh"
        << endl;
}
 
 
void Foam::conformalVoronoiMesh::insertInternalPoints
(
    List<Point>& points,
    bool distribute
)
{
    const label nPoints = returnReduce(points.size(), sumOp<label>());
 
    Info<< "    " << nPoints << " points to insert..." << endl;
 
    if (Pstream::parRun() && distribute)
    {
        List<Foam::point> transferPoints(points.size());
 
        forAll(points, pI)
        {
            transferPoints[pI] = topoint(points[pI]);
        }
 
        // Send the points that are not on this processor to the appropriate
        // place
        Foam::autoPtr<Foam::mapDistribute> map
        (
            decomposition_().distributePoints(transferPoints)
        );
 
        transferPoints.clear();
 
        map().distribute(points);
    }
 
    label preReinsertionSize(number_of_vertices());
 
    insert(points.begin(), points.end());
 
    const label nInserted = returnReduce
    (
        label(number_of_vertices()) - preReinsertionSize,
        sumOp<label>()
    );
 
    Info<< "    " << nInserted << " points inserted"
        << ", failed to insert " << nPoints - nInserted
        << " ("
        << 100.0*(nPoints - nInserted)/(nInserted + SMALL)
        << " %)"<< endl;
 
    for
    (
        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
        vit != finite_vertices_end();
        ++vit
    )
    {
        if (CGAL::indexedVertexOps::uninitialised(vit))
        {
            vit->index() = getNewVertexIndex();
            vit->type() = Vb::vtInternal;
        }
    }
}
 
 
Foam::Map<Foam::label> Foam::conformalVoronoiMesh::insertPointPairs
(
    List<Vb>& vertices,
    bool distribute,
    bool reIndex
)
{
    if (Pstream::parRun() && distribute)
    {
        autoPtr<mapDistribute> mapDist =
            decomposition_().distributePoints(vertices);
 
        // Re-index the point pairs if one or both have been distributed.
        // If both, remove
 
        // If added a point, then need to know its point pair
        // If one moved, then update procIndex locally
 
        forAll(vertices, vI)
        {
            vertices[vI].procIndex() = Pstream::myProcNo();
        }
    }
 
    label preReinsertionSize(number_of_vertices());
 
    Map<label> oldToNewIndices =
        this->DelaunayMesh<Delaunay>::insertPoints(vertices, reIndex);
 
    const label nReinserted = returnReduce
    (
        label(number_of_vertices()) - preReinsertionSize,
        sumOp<label>()
    );
 
    Info<< "    Reinserted " << nReinserted << " vertices out of "
        << returnReduce(vertices.size(), sumOp<label>())
        << endl;
 
    return oldToNewIndices;
}
 
 
void Foam::conformalVoronoiMesh::insertSurfacePointPairs
(
    const pointIndexHitAndFeatureList& surfaceHits,
    const fileName fName,
    DynamicList<Vb>& pts
)
{
    forAll(surfaceHits, i)
    {
        vectorField norm(1);
 
        const pointIndexHit surfaceHit = surfaceHits[i].first();
        const label featureIndex = surfaceHits[i].second();
 
        allGeometry_[featureIndex].getNormal
        (
            List<pointIndexHit>(1, surfaceHit),
            norm
        );
 
        const vector& normal = norm[0];
 
        const Foam::point& surfacePt = surfaceHit.hitPoint();
 
        extendedFeatureEdgeMesh::sideVolumeType meshableSide =
            geometryToConformTo_.meshableSide(featureIndex, surfaceHit);
 
        if (meshableSide == extendedFeatureEdgeMesh::BOTH)
        {
            createBafflePointPair
            (
                pointPairDistance(surfacePt),
                surfacePt,
                normal,
                true,
                pts
            );
        }
        else if (meshableSide == extendedFeatureEdgeMesh::INSIDE)
        {
            createPointPair
            (
                pointPairDistance(surfacePt),
                surfacePt,
                normal,
                true,
                pts
            );
        }
        else if (meshableSide == extendedFeatureEdgeMesh::OUTSIDE)
        {
            createPointPair
            (
                pointPairDistance(surfacePt),
                surfacePt,
                -normal,
                true,
                pts
            );
        }
        else
        {
            WarningInFunction
                << meshableSide << ", bad"
                << endl;
        }
    }
 
    if (foamyHexMeshControls().objOutput() && !fName.empty())
    {
        DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
    }
}
 
 
void Foam::conformalVoronoiMesh::insertEdgePointGroups
(
    const pointIndexHitAndFeatureList& edgeHits,
    const fileName fName,
    DynamicList<Vb>& pts
)
{
    forAll(edgeHits, i)
    {
        if (edgeHits[i].first().hit())
        {
            const extendedFeatureEdgeMesh& feMesh
            (
                geometryToConformTo_.features()[edgeHits[i].second()]
            );
 
//            const bool isBaffle =
//                geometryToConformTo_.isBaffleFeature(edgeHits[i].second());
 
            createEdgePointGroup
            (
                feMesh,
                edgeHits[i].first(),
                pts
            );
        }
    }
 
    if (foamyHexMeshControls().objOutput() && !fName.empty())
    {
        DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
    }
}
 
 
bool Foam::conformalVoronoiMesh::nearFeaturePt(const Foam::point& pt) const
{
    scalar exclusionRangeSqr = featurePointExclusionDistanceSqr(pt);
 
    pointIndexHit info;
    label featureHit;
 
    geometryToConformTo_.findFeaturePointNearest
    (
        pt,
        exclusionRangeSqr,
        info,
        featureHit
    );
 
    return info.hit();
}
 
 
bool Foam::conformalVoronoiMesh::surfacePtNearFeatureEdge
(
    const Foam::point& pt
) const
{
    scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
 
    pointIndexHit info;
    label featureHit;
 
    geometryToConformTo_.findEdgeNearest
    (
        pt,
        exclusionRangeSqr,
        info,
        featureHit
    );
 
    return info.hit();
}
 
 
void Foam::conformalVoronoiMesh::insertInitialPoints()
{
    Info<< nl << "Inserting initial points" << endl;
 
    timeCheck("Before initial points call");
 
    List<Point> initPts = initialPointsMethod_->initialPoints();
 
    timeCheck("After initial points call");
 
    // Assume that the initial points method made the correct decision for
    // which processor each point should be on, so give distribute = false
    insertInternalPoints(initPts);
 
    if (initialPointsMethod_->fixInitialPoints())
    {
        for
        (
            Finite_vertices_iterator vit = finite_vertices_begin();
            vit != finite_vertices_end();
            ++vit
        )
        {
            vit->fixed() = true;
        }
    }
 
    if (foamyHexMeshControls().objOutput())
    {
        DelaunayMeshTools::writeOBJ
        (
            time().path()/"initialPoints.obj",
            *this,
            Foam::indexedVertexEnum::vtInternal
        );
    }
}
 
 
void Foam::conformalVoronoiMesh::distribute()
{
    if (!Pstream::parRun())
    {
        return;
    }
 
    DynamicList<Foam::point> points(number_of_vertices());
    DynamicList<Foam::indexedVertexEnum::vertexType> types
    (
        number_of_vertices()
    );
    DynamicList<scalar> sizes(number_of_vertices());
    DynamicList<tensor> alignments(number_of_vertices());
 
    for
    (
        Finite_vertices_iterator vit = finite_vertices_begin();
        vit != finite_vertices_end();
        ++vit
    )
    {
        if (vit->real())
        {
            points.append(topoint(vit->point()));
            types.append(vit->type());
            sizes.append(vit->targetCellSize());
            alignments.append(vit->alignment());
        }
    }
 
    autoPtr<mapDistribute> mapDist =
        DistributedDelaunayMesh<Delaunay>::distribute(decomposition_(), points);
 
    mapDist().distribute(types);
    mapDist().distribute(sizes);
    mapDist().distribute(alignments);
 
    // Reset the entire tessellation
    DelaunayMesh<Delaunay>::reset();
 
    Info<< nl << "    Inserting distributed tessellation" << endl;
 
    // Internal points have to be inserted first
 
    DynamicList<Vb> verticesToInsert(points.size());
 
    forAll(points, pI)
    {
        verticesToInsert.append
        (
            Vb
            (
                toPoint(points[pI]),
                -1,
                types[pI],
                Pstream::myProcNo()
            )
        );
 
        verticesToInsert.last().targetCellSize() = sizes[pI];
        verticesToInsert.last().alignment() = alignments[pI];
    }
 
    this->rangeInsertWithInfo
    (
        verticesToInsert.begin(),
        verticesToInsert.end(),
        true
    );
 
    Info<< "    Total number of vertices after redistribution "
        << returnReduce
           (
               label(number_of_vertices()), sumOp<label>()
           )
        << endl;
}
 
 
void Foam::conformalVoronoiMesh::buildCellSizeAndAlignmentMesh()
{
    controlMeshRefinement meshRefinement
    (
        cellShapeControl_
    );
 
    smoothAlignmentSolver meshAlignmentSmoother
    (
        cellShapeControl_.shapeControlMesh()
    );
 
    meshRefinement.initialMeshPopulation(decomposition_);
 
    cellShapeControlMesh& cellSizeMesh = cellShapeControl_.shapeControlMesh();
 
    if (Pstream::parRun())
    {
        if (!distributeBackground(cellSizeMesh))
        {
            // Synchronise the cell size mesh if it has not been distributed
            cellSizeMesh.distribute(decomposition_());
        }
    }
 
    const dictionary& motionControlDict
        = foamyHexMeshControls().foamyHexMeshDict().subDict("motionControl");
 
    const label nMaxIter =
        motionControlDict.get<label>("maxRefinementIterations");
 
    Info<< "Maximum number of refinement iterations : " << nMaxIter << endl;
 
    for (label i = 0; i < nMaxIter; ++i)
    {
        label nAdded = meshRefinement.refineMesh(decomposition_);
        //cellShapeControl_.refineMesh(decomposition_);
        reduce(nAdded, sumOp<label>());
 
        if (Pstream::parRun())
        {
            cellSizeMesh.distribute(decomposition_());
        }
 
        Info<< "    Iteration " << i
            << " Added = " << nAdded << " points"
            << endl;
 
        if (nAdded == 0)
        {
            break;
        }
    }
 
    if (Pstream::parRun())
    {
        // Need to distribute the cell size mesh to cover the background mesh
        if (!distributeBackground(cellSizeMesh))
        {
            cellSizeMesh.distribute(decomposition_());
        }
    }
 
    meshAlignmentSmoother.smoothAlignments
    (
        motionControlDict.get<label>("maxSmoothingIterations")
    );
 
    Info<< "Background cell size and alignment mesh:" << endl;
    cellSizeMesh.printInfo(Info);
 
    Info<< "Triangulation is "
        << (cellSizeMesh.is_valid() ? "valid" : "not valid!" )
        << endl;
 
    if (foamyHexMeshControls().writeCellShapeControlMesh())
    {
        //cellSizeMesh.writeTriangulation();
        cellSizeMesh.write();
    }
 
    if (foamyHexMeshControls().printVertexInfo())
    {
        cellSizeMesh.printVertexInfo(Info);
    }
 
//    Info<< "Estimated number of cells in final mesh = "
//        << returnReduce
//           (
//               cellSizeMesh.estimateCellCount(decomposition_),
//               sumOp<label>()
//           )
//        << endl;
}
 
 
void Foam::conformalVoronoiMesh::setVertexSizeAndAlignment()
{
    Info<< nl << "Calculating target cell alignment and size" << endl;
 
    for
    (
        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
        vit != finite_vertices_end();
        vit++
    )
    {
        if (vit->internalOrBoundaryPoint())
        {
            pointFromPoint pt = topoint(vit->point());
 
            cellShapeControls().cellSizeAndAlignment
            (
                pt,
                vit->targetCellSize(),
                vit->alignment()
            );
        }
    }
}
 
 
Foam::face Foam::conformalVoronoiMesh::buildDualFace
(
    const Delaunay::Finite_edges_iterator& eit
) const
{
    Cell_circulator ccStart = incident_cells(*eit);
    Cell_circulator cc1 = ccStart;
    Cell_circulator cc2 = cc1;
 
    // Advance the second circulator so that it always stays on the next
    // cell around the edge;
    cc2++;
 
    DynamicList<label> verticesOnFace;
 
    label nUniqueVertices = 0;
 
    do
    {
        if
        (
            cc1->hasFarPoint() || cc2->hasFarPoint()
         || is_infinite(cc1) || is_infinite(cc2)
        )
        {
            Cell_handle c = eit->first;
            Vertex_handle vA = c->vertex(eit->second);
            Vertex_handle vB = c->vertex(eit->third);
 
//            DelaunayMeshTools::drawDelaunayCell(Pout, cc1);
//            DelaunayMeshTools::drawDelaunayCell(Pout, cc2);
 
            WarningInFunction
                << "Dual face uses circumcenter defined by a "
                << "Delaunay tetrahedron with no internal "
                << "or boundary points.  Defining Delaunay edge ends: "
                << vA->info() << " "
                << vB->info() << nl
                <<endl;//<< exit(FatalError);
        }
        else
        {
            label cc1I = cc1->cellIndex();
            label cc2I = cc2->cellIndex();
 
            if (cc1I != cc2I)
            {
                if (!verticesOnFace.found(cc1I))
                {
                    nUniqueVertices++;
                }
 
                verticesOnFace.append(cc1I);
            }
        }
 
        cc1++;
        cc2++;
 
    } while (cc1 != ccStart);
 
    verticesOnFace.shrink();
 
    if (verticesOnFace.size() >= 3 && nUniqueVertices < 3)
    {
        // There are not enough unique vertices on this face to
        // justify its size, it may have a form like:
 
        // Vertices:
        // A                                  B
        // A                                  B
 
        // Face:
        // ABAB
 
        // Setting the size to be below 3, so that it will not be
        // created
 
        verticesOnFace.setSize(nUniqueVertices);
    }
 
    return face(verticesOnFace);
}
 
 
Foam::label Foam::conformalVoronoiMesh::maxFilterCount
(
    const Delaunay::Finite_edges_iterator& eit
) const
{
    Cell_circulator ccStart = incident_cells(*eit);
    Cell_circulator cc = ccStart;
 
    label maxFC = 0;
 
    do
    {
        if (cc->hasFarPoint())
        {
            Cell_handle c = eit->first;
            Vertex_handle vA = c->vertex(eit->second);
            Vertex_handle vB = c->vertex(eit->third);
 
            FatalErrorInFunction
                << "Dual face uses circumcenter defined by a "
                << "Delaunay tetrahedron with no internal "
                << "or boundary points.  Defining Delaunay edge ends: "
                << topoint(vA->point()) << " "
                << topoint(vB->point()) << nl
                << exit(FatalError);
        }
 
        if (cc->filterCount() > maxFC)
        {
            maxFC = cc->filterCount();
        }
 
        cc++;
 
    } while (cc != ccStart);
 
    return maxFC;
}
 
 
bool Foam::conformalVoronoiMesh::ownerAndNeighbour
(
    Vertex_handle vA,
    Vertex_handle vB,
    label& owner,
    label& neighbour
) const
{
    bool reverse = false;
 
    owner = -1;
    neighbour = -1;
 
    label dualCellIndexA = vA->index();
 
    if (!vA->internalOrBoundaryPoint() || vA->referred())
    {
        if (!vA->constrained())
        {
            dualCellIndexA = -1;
        }
    }
 
    label dualCellIndexB = vB->index();
 
    if (!vB->internalOrBoundaryPoint() || vB->referred())
    {
        if (!vB->constrained())
        {
            dualCellIndexB = -1;
        }
    }
 
    if (dualCellIndexA == -1 && dualCellIndexB == -1)
    {
        FatalErrorInFunction
            << "Attempting to create a face joining "
            << "two unindexed dual cells "
            << exit(FatalError);
    }
    else if (dualCellIndexA == -1 || dualCellIndexB == -1)
    {
        // boundary face, find which is the owner
 
        if (dualCellIndexA == -1)
        {
            owner = dualCellIndexB;
 
            reverse = true;
        }
        else
        {
            owner = dualCellIndexA;
        }
    }
    else
    {
        // internal face, find the lower cell to be the owner
 
        if (dualCellIndexB > dualCellIndexA)
        {
            owner = dualCellIndexA;
            neighbour = dualCellIndexB;
        }
        else
        {
            owner = dualCellIndexB;
            neighbour = dualCellIndexA;
 
            // reverse face order to correctly orientate normal
            reverse = true;
        }
    }
 
    return reverse;
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::conformalVoronoiMesh::conformalVoronoiMesh
(
    const Time& runTime,
    const dictionary& foamyHexMeshDict,
    const fileName& decompDictFile
)
:
    DistributedDelaunayMesh<Delaunay>(runTime),
    runTime_(runTime),
    rndGen_(64293*Pstream::myProcNo()),
    foamyHexMeshControls_(foamyHexMeshDict),
    allGeometry_
    (
        IOobject
        (
            "cvSearchableSurfaces",
            runTime_.constant(),
            "triSurface",
            runTime_,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        ),
        foamyHexMeshDict.subDict("geometry"),
        foamyHexMeshDict.getOrDefault("singleRegionName", true)
    ),
    geometryToConformTo_
    (
        runTime_,
        rndGen_,
        allGeometry_,
        foamyHexMeshDict.subDict("surfaceConformation")
    ),
    decomposition_
    (
        Pstream::parRun()
      ? new backgroundMeshDecomposition
        (
            runTime_,
            rndGen_,
            geometryToConformTo_,
            foamyHexMeshControls().foamyHexMeshDict().subDict
            (
                "backgroundMeshDecomposition"
            ),
            decompDictFile
        )
      : nullptr
    ),
    cellShapeControl_
    (
        runTime_,
        foamyHexMeshControls_,
        allGeometry_,
        geometryToConformTo_
    ),
    limitBounds_(),
    ptPairs_(*this),
    ftPtConformer_(*this),
    edgeLocationTreePtr_(),
    surfacePtLocationTreePtr_(),
    surfaceConformationVertices_(),
    initialPointsMethod_
    (
        initialPointsMethod::New
        (
            foamyHexMeshDict.subDict("initialPoints"),
            runTime_,
            rndGen_,
            geometryToConformTo_,
            cellShapeControl_,
            decomposition_
        )
    ),
    relaxationModel_
    (
        relaxationModel::New
        (
            foamyHexMeshDict.subDict("motionControl"),
            runTime_
        )
    ),
    faceAreaWeightModel_
    (
        faceAreaWeightModel::New
        (
            foamyHexMeshDict.subDict("motionControl")
        )
    )
{}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
Foam::conformalVoronoiMesh::~conformalVoronoiMesh()
{}
 
 
// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
void Foam::conformalVoronoiMesh::initialiseForMotion()
{
    if (foamyHexMeshControls().objOutput())
    {
        geometryToConformTo_.writeFeatureObj("foamyHexMesh");
    }
 
    buildCellSizeAndAlignmentMesh();
 
    insertInitialPoints();
 
    insertFeaturePoints(true);
 
    setVertexSizeAndAlignment();
 
    cellSizeMeshOverlapsBackground();
 
    // Improve the guess that the backgroundMeshDecomposition makes with the
    // initial positions.  Use before building the surface conformation to
    // better balance the surface conformation load.
    distributeBackground(*this);
 
    buildSurfaceConformation();
 
    // The introduction of the surface conformation may have distorted the
    // balance of vertices, distribute if necessary.
    distributeBackground(*this);
 
    if (Pstream::parRun())
    {
        sync(decomposition_().procBounds());
    }
 
    // Do not store the surface conformation until after it has been
    // (potentially) redistributed.
    storeSurfaceConformation();
 
    // Report any Delaunay vertices that do not think that they are in the
    // domain the processor they are on.
    // reportProcessorOccupancy();
 
    cellSizeMeshOverlapsBackground();
 
    if (foamyHexMeshControls().printVertexInfo())
    {
        printVertexInfo(Info);
    }
 
    if (foamyHexMeshControls().objOutput())
    {
        DelaunayMeshTools::writeOBJ
        (
            time().path()/"internalPoints_" + time().timeName() + ".obj",
            *this,
            Foam::indexedVertexEnum::vtUnassigned,
            Foam::indexedVertexEnum::vtExternalFeaturePoint
        );
    }
}
 
 
void Foam::conformalVoronoiMesh::initialiseForConformation()
{
    if (Pstream::parRun())
    {
        decomposition_.reset
        (
            new backgroundMeshDecomposition
            (
                runTime_,
                rndGen_,
                geometryToConformTo_,
                foamyHexMeshControls().foamyHexMeshDict().subDict
                (
                    "backgroundMeshDecomposition"
                )
            )
        );
    }
 
    insertInitialPoints();
 
    insertFeaturePoints();
 
    // Improve the guess that the backgroundMeshDecomposition makes with the
    // initial positions.  Use before building the surface conformation to
    // better balance the surface conformation load.
    distributeBackground(*this);
 
    buildSurfaceConformation();
 
    // The introduction of the surface conformation may have distorted the
    // balance of vertices, distribute if necessary.
    distributeBackground(*this);
 
    if (Pstream::parRun())
    {
        sync(decomposition_().procBounds());
    }
 
    cellSizeMeshOverlapsBackground();
 
    if (foamyHexMeshControls().printVertexInfo())
    {
        printVertexInfo(Info);
    }
}
 
 
void Foam::conformalVoronoiMesh::move()
{
    timeCheck("Start of move");
 
    scalar relaxation = relaxationModel_->relaxation();
 
    Info<< nl << "Relaxation = " << relaxation << endl;
 
    pointField dualVertices(number_of_finite_cells());
 
    this->resetCellCount();
 
    // Find the dual point of each tetrahedron and assign it an index.
    for
    (
        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
        cit != finite_cells_end();
        ++cit
    )
    {
        cit->cellIndex() = Cb::ctUnassigned;
 
        if (cit->anyInternalOrBoundaryDualVertex())
        {
            cit->cellIndex() = getNewCellIndex();
 
            dualVertices[cit->cellIndex()] = cit->dual();
        }
 
        if (cit->hasFarPoint())
        {
            cit->cellIndex() = Cb::ctFar;
        }
    }
 
    dualVertices.setSize(cellCount());
 
    setVertexSizeAndAlignment();
 
    timeCheck("Determined sizes and alignments");
 
    Info<< nl << "Determining vertex displacements" << endl;
 
    vectorField cartesianDirections(3);
 
    cartesianDirections[0] = vector(1, 0, 0);
    cartesianDirections[1] = vector(0, 1, 0);
    cartesianDirections[2] = vector(0, 0, 1);
 
    vectorField displacementAccumulator
    (
        number_of_vertices(),
        Zero
    );
 
    bitSet pointToBeRetained(number_of_vertices(), true);
 
    DynamicList<Point> pointsToInsert(number_of_vertices());
 
    for
    (
        Delaunay::Finite_edges_iterator eit = finite_edges_begin();
        eit != finite_edges_end();
        ++eit
    )
    {
        Cell_handle c = eit->first;
        Vertex_handle vA = c->vertex(eit->second);
        Vertex_handle vB = c->vertex(eit->third);
 
        if
        (
            (
                vA->internalPoint() && !vA->referred()
             && vB->internalOrBoundaryPoint()
            )
         || (
                vB->internalPoint() && !vB->referred()
             && vA->internalOrBoundaryPoint()
            )
        )
        {
            pointFromPoint dVA = topoint(vA->point());
            pointFromPoint dVB = topoint(vB->point());
 
            Field<vector> alignmentDirsA
            (
                vA->alignment().T() & cartesianDirections
            );
            Field<vector> alignmentDirsB
            (
                vB->alignment().T() & cartesianDirections
            );
 
            Field<vector> alignmentDirs(alignmentDirsA);
 
            forAll(alignmentDirsA, aA)
            {
                const vector& a = alignmentDirsA[aA];
 
                scalar maxDotProduct = 0.0;
 
                forAll(alignmentDirsB, aB)
                {
                    const vector& b = alignmentDirsB[aB];
 
                    const scalar dotProduct = a & b;
 
                    if (mag(dotProduct) > maxDotProduct)
                    {
                        maxDotProduct = mag(dotProduct);
 
                        alignmentDirs[aA] = a + sign(dotProduct)*b;
 
                        alignmentDirs[aA].normalise();
                    }
                }
            }
 
            vector rAB = dVA - dVB;
 
            scalar rABMag = mag(rAB);
 
            if (rABMag < SMALL)
            {
                // Removal of close points
 
                if
                (
                    vA->internalPoint() && !vA->referred() && !vA->fixed()
                 && vB->internalPoint() && !vB->referred() && !vB->fixed()
                )
                {
                    // Only insert a point at the midpoint of
                    // the short edge if neither attached
                    // point has already been identified to be
                    // removed.
 
                    if
                    (
                        pointToBeRetained.test(vA->index())
                     && pointToBeRetained.test(vB->index())
                    )
                    {
                        const Foam::point pt(0.5*(dVA + dVB));
 
                        if (internalPointIsInside(pt))
                        {
                            pointsToInsert.append(toPoint(pt));
                        }
                    }
                }
 
                if (vA->internalPoint() && !vA->referred() && !vA->fixed())
                {
                    pointToBeRetained.unset(vA->index());
                }
 
                if (vB->internalPoint() && !vB->referred() && !vB->fixed())
                {
                    pointToBeRetained.unset(vB->index());
                }
 
                // Do not consider this Delaunay edge any further
 
                continue;
            }
 
            forAll(alignmentDirs, aD)
            {
                vector& alignmentDir = alignmentDirs[aD];
 
                scalar dotProd = rAB & alignmentDir;
 
                if (dotProd < 0)
                {
                    // swap the direction of the alignment so that has the
                    // same sense as rAB
                    alignmentDir *= -1;
                    dotProd *= -1;
                }
 
                const scalar alignmentDotProd = dotProd/rABMag;
 
                if
                (
                    alignmentDotProd
                  > foamyHexMeshControls().cosAlignmentAcceptanceAngle()
                )
                {
                    scalar targetCellSize =
                        CGAL::indexedVertexOps::averageCellSize(vA, vB);
 
                    scalar targetFaceArea = sqr(targetCellSize);
 
                    const vector originalAlignmentDir = alignmentDir;
 
                    // Update cell size and face area
                    cellShapeControls().aspectRatio().updateCellSizeAndFaceArea
                    (
                        alignmentDir,
                        targetFaceArea,
                        targetCellSize
                    );
 
                    // Vector to move end points around middle of vector
                    // to align edge (i.e. dual face normal) with alignment
                    // directions.
                    vector delta = alignmentDir - 0.5*rAB;
 
                    face dualFace = buildDualFace(eit);
 
//                    Pout<< dualFace << endl;
//                    Pout<< "    " << vA->info() << endl;
//                    Pout<< "    " << vB->info() << endl;
 
                    const scalar faceArea = dualFace.mag(dualVertices);
 
                    // Update delta vector
                    cellShapeControls().aspectRatio().updateDeltaVector
                    (
                        originalAlignmentDir,
                        targetCellSize,
                        rABMag,
                        delta
                    );
 
//                    if (targetFaceArea == 0)
//                    {
//                        Pout<< vA->info() << vB->info();
//
//                        Cell_handle ch = locate(vA->point());
//                        if (is_infinite(ch))
//                        {
//                            Pout<< "vA " << vA->targetCellSize() << endl;
//                        }
//
//                        ch = locate(vB->point());
//                        if (is_infinite(ch))
//                        {
//                            Pout<< "vB " << vB->targetCellSize() << endl;
//                        }
//                    }
 
                    delta *= faceAreaWeightModel_->faceAreaWeight
                    (
                        faceArea/targetFaceArea
                    );
 
                    if
                    (
                        (
                            (vA->internalPoint() && vB->internalPoint())
                         && (!vA->referred() || !vB->referred())
//                         ||
//                            (
//                                vA->referredInternalPoint()
//                             && vB->referredInternalPoint()
//                            )
                        )
                     && rABMag
                      > foamyHexMeshControls().insertionDistCoeff()
                       *targetCellSize
                     && faceArea
                      > foamyHexMeshControls().faceAreaRatioCoeff()
                       *targetFaceArea
                     && alignmentDotProd
                      > foamyHexMeshControls().cosInsertionAcceptanceAngle()
                    )
                    {
                        // Point insertion
                        if
                        (
                            !geometryToConformTo_.findSurfaceAnyIntersection
                            (
                                dVA,
                                dVB
                            )
                        )
                        {
                            const Foam::point newPt(0.5*(dVA + dVB));
 
                            // Prevent insertions spanning surfaces
                            if (internalPointIsInside(newPt))
                            {
                                if (Pstream::parRun())
                                {
                                    if
                                    (
                                        decomposition().positionOnThisProcessor
                                        (
                                            newPt
                                        )
                                    )
                                    {
                                        pointsToInsert.append(toPoint(newPt));
                                    }
                                }
                                else
                                {
                                    pointsToInsert.append(toPoint(newPt));
                                }
                            }
                        }
                    }
                    else if
                    (
                        (
                            (vA->internalPoint() && !vA->referred())
                         || (vB->internalPoint() && !vB->referred())
                        )
                     && rABMag
                      < foamyHexMeshControls().removalDistCoeff()
                       *targetCellSize
                    )
                    {
                        // Point removal
                        if
                        (
                            (
                                vA->internalPoint()
                             && !vA->referred()
                             && !vA->fixed()
                            )
                         &&
                            (
                                vB->internalPoint()
                             && !vB->referred()
                             && !vB->fixed()
                            )
                        )
                        {
                            // Only insert a point at the midpoint of
                            // the short edge if neither attached
                            // point has already been identified to be
                            // removed.
                            if
                            (
                                pointToBeRetained.test(vA->index())
                             && pointToBeRetained.test(vB->index())
                            )
                            {
                                const Foam::point pt(0.5*(dVA + dVB));
 
                                if (internalPointIsInside(pt))
                                {
                                    pointsToInsert.append(toPoint(pt));
                                }
                            }
                        }
 
                        if
                        (
                            vA->internalPoint()
                         && !vA->referred()
                         && !vA->fixed()
                        )
                        {
                            pointToBeRetained.unset(vA->index());
                        }
 
                        if
                        (
                            vB->internalPoint()
                         && !vB->referred()
                         && !vB->fixed()
                        )
                        {
                            pointToBeRetained.unset(vB->index());
                        }
                    }
                    else
                    {
                        if
                        (
                            vA->internalPoint()
                         && !vA->referred()
                         && !vA->fixed()
                        )
                        {
                            if (vB->fixed())
                            {
                                displacementAccumulator[vA->index()] += 2*delta;
                            }
                            else
                            {
                                displacementAccumulator[vA->index()] += delta;
                            }
                        }
 
                        if
                        (
                            vB->internalPoint()
                         && !vB->referred()
                         && !vB->fixed()
                        )
                        {
                            if (vA->fixed())
                            {
                                displacementAccumulator[vB->index()] -= 2*delta;
                            }
                            else
                            {
                                displacementAccumulator[vB->index()] -= delta;
                            }
                        }
                    }
                }
            }
        }
    }
 
    Info<< "Limit displacements" << endl;
 
    // Limit displacements that pierce, or get too close to the surface
    for
    (
        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
        vit != finite_vertices_end();
        ++vit
    )
    {
        if (vit->internalPoint() && !vit->referred() && !vit->fixed())
        {
            if (pointToBeRetained.test(vit->index()))
            {
                limitDisplacement
                (
                    vit,
                    displacementAccumulator[vit->index()]
                );
            }
        }
    }
 
    vector totalDisp = gSum(displacementAccumulator);
    scalar totalDist = gSum(mag(displacementAccumulator));
 
    displacementAccumulator *= relaxation;
 
    Info<< "Sum displacements" << endl;
 
    label nPointsToRetain = 0;
    label nPointsToRemove = 0;
 
    for
    (
        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
        vit != finite_vertices_end();
        ++vit
    )
    {
        if (vit->internalPoint() && !vit->referred() && !vit->fixed())
        {
            if (pointToBeRetained.test(vit->index()))
            {
                // Convert vit->point() to FOAM vector (double) to do addition,
                // avoids memory increase because a record of the constructions
                // would be kept otherwise.
                // See cgal-discuss@lists-sop.inria.fr:
                // "Memory issue with openSUSE 11.3, exact kernel, adding
                //  points/vectors"
                // 14/1/2011.
                // Only necessary if using an exact constructions kernel
                // (extended precision)
                Foam::point pt
                (
                    topoint(vit->point())
                  + displacementAccumulator[vit->index()]
                );
 
                if (internalPointIsInside(pt))
                {
                    pointsToInsert.append(toPoint(pt));
                    nPointsToRemove++;
                }
 
                nPointsToRetain++;
            }
        }
    }
 
    pointsToInsert.shrink();
 
    Info<< returnReduce
           (
               nPointsToRetain - nPointsToRemove,
               sumOp<label>()
           )
        << " internal points are outside the domain. "
        << "They will not be inserted." << endl;
 
    // Save displacements to file.
    if (foamyHexMeshControls().objOutput() && time().writeTime())
    {
        Info<< "Writing point displacement vectors to file." << endl;
        OFstream str
        (
            time().path()/"displacements_" + runTime_.timeName() + ".obj"
        );
 
        for
        (
            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
            vit != finite_vertices_end();
            ++vit
        )
        {
            if (vit->internalPoint() && !vit->referred())
            {
                if (pointToBeRetained.test(vit->index()))
                {
                    meshTools::writeOBJ(str, topoint(vit->point()));
 
                    str << "vn "
                        << displacementAccumulator[vit->index()][0] << " "
                        << displacementAccumulator[vit->index()][1] << " "
                        << displacementAccumulator[vit->index()][2] << " "
                        << endl;
                }
            }
        }
    }
 
    // Remove the entire tessellation
    DelaunayMesh<Delaunay>::reset();
 
    timeCheck("Displacement calculated");
 
    Info<< nl<< "Inserting displaced tessellation" << endl;
 
    insertFeaturePoints(true);
 
    insertInternalPoints(pointsToInsert, true);
 
    // Fix points that have not been significantly displaced
//    for
//    (
//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
//        vit != finite_vertices_end();
//        ++vit
//    )
//    {
//        if (vit->internalPoint())
//        {
//            if
//            (
//                mag(displacementAccumulator[vit->index()])
//              < 0.1*targetCellSize(topoint(vit->point()))
//            )
//            {
//                vit->setVertexFixed();
//            }
//        }
//    }
 
    timeCheck("Internal points inserted");
 
    {
        // Check that no index is shared between any of the local points
        labelHashSet usedIndices;
        for
        (
            Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
            vit != finite_vertices_end();
            ++vit
        )
        {
            if (!vit->referred() && !usedIndices.insert(vit->index()))
            {
                FatalErrorInFunction
                    << "Index already used! Could not insert: " << nl
                    << vit->info()
                    << abort(FatalError);
            }
        }
    }
 
    conformToSurface();
 
    if (foamyHexMeshControls().objOutput())
    {
        DelaunayMeshTools::writeOBJ
        (
            time().path()/"internalPoints_" + time().timeName() + ".obj",
            *this,
            Foam::indexedVertexEnum::vtInternal
        );
 
        if (reconformToSurface())
        {
            DelaunayMeshTools::writeBoundaryPoints
            (
                time().path()/"boundaryPoints_" + time().timeName() + ".obj",
                *this
            );
 
            DelaunayMeshTools::writeOBJ
            (
                time().path()/"internalBoundaryPoints_" + time().timeName()
              + ".obj",
                *this,
                Foam::indexedVertexEnum::vtInternalSurface,
                Foam::indexedVertexEnum::vtInternalFeaturePoint
            );
 
            DelaunayMeshTools::writeOBJ
            (
                time().path()/"externalBoundaryPoints_" + time().timeName()
              + ".obj",
                *this,
                Foam::indexedVertexEnum::vtExternalSurface,
                Foam::indexedVertexEnum::vtExternalFeaturePoint
            );
 
            OBJstream multipleIntersections
            (
                "multipleIntersections_"
              + time().timeName()
              + ".obj"
            );
 
            for
            (
                Delaunay::Finite_edges_iterator eit = finite_edges_begin();
                eit != finite_edges_end();
                ++eit
            )
            {
                Cell_handle c = eit->first;
                Vertex_handle vA = c->vertex(eit->second);
                Vertex_handle vB = c->vertex(eit->third);
 
                Foam::point ptA(topoint(vA->point()));
                Foam::point ptB(topoint(vB->point()));
 
                List<pointIndexHit> surfHits;
                labelList hitSurfaces;
 
                geometryToConformTo().findSurfaceAllIntersections
                (
                    ptA,
                    ptB,
                    surfHits,
                    hitSurfaces
                );
 
                if
                (
                    surfHits.size() >= 2
                 || (
                     surfHits.size() == 0
                  && (
                          (vA->externalBoundaryPoint()
                       && vB->internalBoundaryPoint())
                       || (vB->externalBoundaryPoint()
                       && vA->internalBoundaryPoint())
                     )
                    )
                )
                {
                    multipleIntersections.writeLine(ptA, ptB);
                }
            }
        }
    }
 
    timeCheck("After conformToSurface");
 
    if (foamyHexMeshControls().printVertexInfo())
    {
        printVertexInfo(Info);
    }
 
    if (time().writeTime())
    {
        writeMesh(time().timeName());
    }
 
    Info<< nl
        << "Total displacement = " << totalDisp << nl
        << "Total distance = " << totalDist << nl
        << endl;
}
 
 
void Foam::conformalVoronoiMesh::checkCoPlanarCells() const
{
    typedef CGAL::Exact_predicates_exact_constructions_kernel   Kexact;
    typedef CGAL::Point_3<Kexact>                               PointExact;
 
    if (!is_valid())
    {
        Pout<< "Triangulation is invalid!" << endl;
    }
 
    OFstream str("badCells.obj");
 
    label badCells = 0;
 
    for
    (
        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
        cit != finite_cells_end();
        ++cit
    )
    {
        const scalar quality = foamyHexMeshChecks::coplanarTet(cit, 1e-16);
 
        if (quality == 0)
        {
            Pout<< "COPLANAR: " << cit->info() << nl
                << "    quality = " << quality << nl
                << "    dual    = " << topoint(cit->dual()) << endl;
 
            DelaunayMeshTools::drawDelaunayCell(str, cit, badCells++);
 
            FixedList<PointExact, 4> cellVerticesExact(PointExact(0,0,0));
            forAll(cellVerticesExact, vI)
            {
                cellVerticesExact[vI] = PointExact
                (
                    cit->vertex(vI)->point().x(),
                    cit->vertex(vI)->point().y(),
                    cit->vertex(vI)->point().z()
                );
            }
 
            PointExact synchronisedDual = CGAL::circumcenter<Kexact>
            (
                cellVerticesExact[0],
                cellVerticesExact[1],
                cellVerticesExact[2],
                cellVerticesExact[3]
            );
 
            Foam::point exactPt
            (
                CGAL::to_double(synchronisedDual.x()),
                CGAL::to_double(synchronisedDual.y()),
                CGAL::to_double(synchronisedDual.z())
            );
 
            Info<< "inexact = " << cit->dual() << nl
                << "exact   = " << exactPt << endl;
        }
    }
 
    Pout<< "There are " << badCells << " bad cells out of "
        << number_of_finite_cells() << endl;
 
 
    label nNonGabriel = 0;
    for
    (
        Delaunay::Finite_facets_iterator fit = finite_facets_begin();
        fit != finite_facets_end();
        ++fit
    )
    {
        if (!is_Gabriel(*fit))
        {
            nNonGabriel++;//Pout<< "Non-gabriel face" << endl;
        }
    }
 
    Pout<< "There are " << nNonGabriel << " non-Gabriel faces out of "
        << number_of_finite_facets() << endl;
}
 
 
// ************************************************************************* //