controlMeshRefinement.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) 2013-2015 OpenFOAM Foundation
    Copyright (C) 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 "controlMeshRefinement.H"
#include "cellSizeAndAlignmentControl.H"
#include "OFstream.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(controlMeshRefinement, 0);
}

// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

Foam::scalar Foam::controlMeshRefinement::calcFirstDerivative
(
    const Foam::point& a,
    const scalar& cellSizeA,
    const Foam::point& b,
    const scalar& cellSizeB
) const
{
    return (cellSizeA - cellSizeB)/mag(a - b);
}


//Foam::scalar Foam::controlMeshRefinement::calcSecondDerivative
//(
//    const Foam::point& a,
//    const scalar& cellSizeA,
//    const Foam::point& midPoint,
//    const scalar& cellSizeMid,
//    const Foam::point& b,
//    const scalar& cellSizeB
//) const
//{
//    return (cellSizeA - 2*cellSizeMid + cellSizeB)/magSqr((a - b)/2);
//}


bool Foam::controlMeshRefinement::detectEdge
(
    const Foam::point& startPt,
    const Foam::point& endPt,
    pointHit& pointFound,
    const scalar tolSqr,
    const scalar secondDerivTolSqr
) const
{
    Foam::point a(startPt);
    Foam::point b(endPt);

    Foam::point midPoint = (a + b)/2.0;

    //label nIterations = 0;

    while (true)
    {
        //++nIterations;

        if
        (
            magSqr(a - b) < tolSqr
        )
        {
            pointFound.hitPoint(midPoint);
            return true;
        }

        // Split into two regions

        scalar cellSizeA = sizeControls_.cellSize(a);
        scalar cellSizeB = sizeControls_.cellSize(b);

//        if (magSqr(cellSizeA - cellSizeB) < 1e-6)
//        {
//            return false;
//        }

        scalar cellSizeMid = sizeControls_.cellSize(midPoint);

        // Region 1
        Foam::point midPoint1 = (a + midPoint)/2.0;
        const scalar cellSizeMid1 = sizeControls_.cellSize(midPoint1);

//        scalar firstDerivative1 =
//            calcFirstDerivative(cellSizeA, cellSizeMid);

        scalar secondDerivative1 =
            calcSecondDerivative
            (
                a,
                cellSizeA,
                midPoint1,
                cellSizeMid1,
                midPoint,
                cellSizeMid
            );

        // Region 2
        Foam::point midPoint2 = (midPoint + b)/2.0;
        const scalar cellSizeMid2 = sizeControls_.cellSize(midPoint2);

//        scalar firstDerivative2 =
//            calcFirstDerivative(f, cellSizeMid, cellSizeB);

        scalar secondDerivative2 =
            calcSecondDerivative
            (
                midPoint,
                cellSizeMid,
                midPoint2,
                cellSizeMid2,
                b,
                cellSizeB
            );

        // Neither region appears to have an inflection
        // To be sure should use higher order derivatives
        if
        (
            magSqr(secondDerivative1) < secondDerivTolSqr
         && magSqr(secondDerivative2) < secondDerivTolSqr
        )
        {
            return false;
        }

        // Pick region with greatest second derivative
        if (magSqr(secondDerivative1) > magSqr(secondDerivative2))
        {
            b = midPoint;
            midPoint = midPoint1;
        }
        else
        {
            a = midPoint;
            midPoint = midPoint2;
        }
    }
}


Foam::pointHit Foam::controlMeshRefinement::findDiscontinuities
(
    const linePointRef& l
) const
{
    pointHit p(point::max);

    const scalar tolSqr = sqr(1e-3);
    const scalar secondDerivTolSqr = sqr(1e-3);

    detectEdge
    (
        l.start(),
        l.end(),
        p,
        tolSqr,
        secondDerivTolSqr
    );

    return p;
}


// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::controlMeshRefinement::controlMeshRefinement
(
    cellShapeControl& shapeController
)
:
    shapeController_(shapeController),
    mesh_(shapeController.shapeControlMesh()),
    sizeControls_(shapeController.sizeAndAlignment()),
    geometryToConformTo_(sizeControls_.geometryToConformTo())
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::controlMeshRefinement::~controlMeshRefinement()
{}


// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

void Foam::controlMeshRefinement::initialMeshPopulation
(
    const autoPtr<backgroundMeshDecomposition>& decomposition
)
{
    if (shapeController_.shapeControlMesh().vertexCount() > 0)
    {
        // Mesh already populated.
        Info<< "Cell size and alignment mesh already populated." << endl;
        return;
    }

    autoPtr<boundBox> overallBoundBox;

    // Need to pass in the background mesh decomposition so that can test if
    // a point to insert is on the processor.
    if (Pstream::parRun())
    {
//        overallBoundBox.set(new boundBox(decomposition().procBounds()));
    }
    else
    {
//        overallBoundBox.set
//        (
//            new boundBox(geometryToConformTo_.geometry().bounds())
//        );
//
//        mesh_.insertBoundingPoints
//        (
//            overallBoundBox(),
//            sizeControls_
//        );
    }

    Map<label> priorityMap;

    const PtrList<cellSizeAndAlignmentControl>& controlFunctions =
        sizeControls_.controlFunctions();

    forAll(controlFunctions, fI)
    {
        const cellSizeAndAlignmentControl& controlFunction =
            controlFunctions[fI];

        const Switch forceInsertion =
            controlFunction.forceInitialPointInsertion();

        Info<< "Inserting points from " << controlFunction.name()
            << " (" << controlFunction.type() << ")" << endl;
        Info<< "    Force insertion is " << forceInsertion.c_str() << endl;

        pointField pts;
        scalarField sizes;
        triadField alignments;

        controlFunction.initialVertices(pts, sizes, alignments);

        Info<< "    Got initial vertices list of size " << pts.size() << endl;

        List<Vb> vertices(pts.size());

        // Clip the minimum size
        for (label vI = 0; vI < pts.size(); ++vI)
        {
            vertices[vI] = Vb(pts[vI], Vb::vtInternalNearBoundary);

            label maxPriority = -1;
            scalar size = sizeControls_.cellSize(pts[vI], maxPriority);

            if (maxPriority > controlFunction.maxPriority())
            {
                vertices[vI].targetCellSize() = max
                (
                    size,
                    shapeController_.minimumCellSize()
                );
            }
//            else if (maxPriority == controlFunction.maxPriority())
//            {
//                vertices[vI].targetCellSize() = max
//                (
//                    min(sizes[vI], size),
//                    shapeController_.minimumCellSize()
//                );
//            }
            else
            {
                vertices[vI].targetCellSize() = max
                (
                    sizes[vI],
                    shapeController_.minimumCellSize()
                );
            }

            vertices[vI].alignment() = alignments[vI];
        }

        Info<< "    Clipped minimum size" << endl;

        pts.clear();
        sizes.clear();
        alignments.clear();

        bitSet keepVertex(vertices.size(), true);

        forAll(vertices, vI)
        {
            bool keep = true;

            pointFromPoint pt = topoint(vertices[vI].point());

            if (Pstream::parRun())
            {
                keep = decomposition().positionOnThisProcessor(pt);
            }

            if (keep && geometryToConformTo_.wellOutside(pt, SMALL))
            {
                keep = false;
            }

            if (!keep)
            {
                keepVertex.unset(vI);
            }
        }

        inplaceSubset(keepVertex, vertices);

        const label preInsertedSize = mesh_.number_of_vertices();

        Info<< "    Check sizes" << endl;

        forAll(vertices, vI)
        {
            bool insertPoint = false;

            pointFromPoint pt(topoint(vertices[vI].point()));

            if
            (
                mesh_.dimension() < 3
             || mesh_.is_infinite
                (
                    mesh_.locate(vertices[vI].point())
                )
            )
            {
                insertPoint = true;
            }

            const scalar interpolatedCellSize = shapeController_.cellSize(pt);
            const triad interpolatedAlignment =
                shapeController_.cellAlignment(pt);
            const scalar calculatedCellSize = vertices[vI].targetCellSize();
            const triad calculatedAlignment = vertices[vI].alignment();

            if (debug)
            {
                Info<< "Point = " << pt << nl
                    << "  Size(interp) = " << interpolatedCellSize << nl
                    << "    Size(calc) = " << calculatedCellSize << nl
                    << " Align(interp) = " << interpolatedAlignment << nl
                    << "   Align(calc) = " << calculatedAlignment << nl
                    << endl;
            }

            const scalar sizeDiff =
                mag(interpolatedCellSize - calculatedCellSize);
            const scalar alignmentDiff =
                diff(interpolatedAlignment, calculatedAlignment);

            if (debug)
            {
                Info<< "    size difference = " << sizeDiff << nl
                    << ", alignment difference = " << alignmentDiff << endl;
            }

            // TODO: Also need to base it on the alignments
            if
            (
                sizeDiff/interpolatedCellSize > 0.1
             || alignmentDiff > 0.15
            )
            {
                insertPoint = true;
            }

            if (forceInsertion || insertPoint)
            {
                const label oldSize = mesh_.vertexCount();

                cellShapeControlMesh::Vertex_handle insertedVert = mesh_.insert
                (
                    pt,
                    calculatedCellSize,
                    vertices[vI].alignment(),
                    Vb::vtInternalNearBoundary
                );

                if (oldSize == mesh_.vertexCount() - 1)
                {
                    priorityMap.insert
                    (
                        insertedVert->index(),
                        controlFunction.maxPriority()
                    );
                }
            }
        }

        //mesh_.rangeInsertWithInfo(vertices.begin(), vertices.end());

        Info<< "    Inserted "
            << returnReduce
               (
                   label(mesh_.number_of_vertices()) - preInsertedSize,
                   sumOp<label>()
               )
            << "/" << returnReduce(vertices.size(), sumOp<label>())
            << endl;
    }



    forAll(controlFunctions, fI)
    {
        const cellSizeAndAlignmentControl& controlFunction =
            controlFunctions[fI];

        const Switch forceInsertion =
            controlFunction.forceInitialPointInsertion();

        Info<< "Inserting points from " << controlFunction.name()
            << " (" << controlFunction.type() << ")" << endl;
        Info<< "    Force insertion is " << forceInsertion.c_str() << endl;

        DynamicList<Foam::point> extraPts;
        DynamicList<scalar> extraSizes;

        controlFunction.cellSizeFunctionVertices(extraPts, extraSizes);

        List<Vb> vertices(extraPts.size());

        // Clip the minimum size
        for (label vI = 0; vI < extraPts.size(); ++vI)
        {
            vertices[vI] = Vb(extraPts[vI], Vb::vtUnassigned);

            label maxPriority = -1;
            scalar size = sizeControls_.cellSize(extraPts[vI], maxPriority);

            if (maxPriority > controlFunction.maxPriority())
            {
                vertices[vI].targetCellSize() = max
                (
                    size,
                    shapeController_.minimumCellSize()
                );
            }
            else if (maxPriority == controlFunction.maxPriority())
            {
                vertices[vI].targetCellSize() = max
                (
                    min(extraSizes[vI], size),
                    shapeController_.minimumCellSize()
                );
            }
            else
            {
                vertices[vI].targetCellSize() = max
                (
                    extraSizes[vI],
                    shapeController_.minimumCellSize()
                );
            }
        }

        bitSet keepVertex(vertices.size(), true);

        forAll(vertices, vI)
        {
            bool keep = true;

            pointFromPoint pt = topoint(vertices[vI].point());

            if (Pstream::parRun())
            {
                keep = decomposition().positionOnThisProcessor(pt);
            }

            if (keep && geometryToConformTo_.wellOutside(pt, SMALL))
            {
                keep = false;
            }

            if (!keep)
            {
                keepVertex.unset(vI);
            }
        }

        inplaceSubset(keepVertex, vertices);

        const label preInsertedSize = mesh_.number_of_vertices();

        forAll(vertices, vI)
        {
            bool insertPoint = false;

            pointFromPoint pt(topoint(vertices[vI].point()));

            if
            (
                mesh_.dimension() < 3
             || mesh_.is_infinite
                (
                    mesh_.locate(vertices[vI].point())
                )
            )
            {
                insertPoint = true;
            }

            const scalar interpolatedCellSize = shapeController_.cellSize(pt);
            const scalar calculatedCellSize = vertices[vI].targetCellSize();

            if (debug)
            {
                Info<< "Point = " << pt << nl
                    << "  Size(interp) = " << interpolatedCellSize << nl
                    << "    Size(calc) = " << calculatedCellSize << nl
                    << endl;
            }

            const scalar sizeDiff =
                mag(interpolatedCellSize - calculatedCellSize);

            if (debug)
            {
                Info<< "    size difference = " << sizeDiff << endl;
            }

            // TODO: Also need to base it on the alignments
            if (sizeDiff/interpolatedCellSize > 0.1)
            {
                insertPoint = true;
            }

            if (forceInsertion || insertPoint)
            {
                // Check the priority

//                cellShapeControlMesh::Cell_handle ch =
//                    mesh_.locate(toPoint<cellShapeControlMesh::Point>(pt));

//                if (mesh_.is_infinite(ch))
//                {
//                    continue;
//                }

//                const label newPtPriority = controlFunction.maxPriority();

//                label highestPriority = -1;
//                for (label cI = 0; cI < 4; ++cI)
//                {
//                    if (mesh_.is_infinite(ch->vertex(cI)))
//                    {
//                        continue;
//                    }

//                    const label vertPriority =
//                        priorityMap[ch->vertex(cI)->index()];

//                    if (vertPriority > highestPriority)
//                    {
//                        highestPriority = vertPriority;
//                    }
//                }

//                if (newPtPriority >= highestPriority)
//                {
//                    const label oldSize = mesh_.vertexCount();
//
//                    cellShapeControlMesh::Vertex_handle insertedVert =
                        mesh_.insert
                        (
                            pt,
                            calculatedCellSize,
                            vertices[vI].alignment(),
                            Vb::vtInternal
                        );

//                    if (oldSize == mesh_.vertexCount() - 1)
//                    {
//                        priorityMap.insert
//                        (
//                            insertedVert->index(),
//                            newPtPriority
//                        );
//                    }
//                }
            }
        }

       //mesh_.rangeInsertWithInfo(vertices.begin(), vertices.end());

        Info<< "    Inserted extra points "
            << returnReduce
               (
                   label(mesh_.number_of_vertices()) - preInsertedSize,
                   sumOp<label>()
               )
            << "/" << returnReduce(vertices.size(), sumOp<label>())
            << endl;
    }

    // Change cell size function of bounding points to be consistent
    // with their nearest neighbours
//    for
//    (
//        CellSizeDelaunay::Finite_vertices_iterator vit =
//            mesh_.finite_vertices_begin();
//        vit != mesh_.finite_vertices_end();
//        ++vit
//    )
//    {
//        if (vit->uninitialised())
//        {
//            // Get its adjacent vertices
//            std::list<CellSizeDelaunay::Vertex_handle> adjacentVertices;
//
//            mesh_.adjacent_vertices
//            (
//                vit,
//                std::back_inserter(adjacentVertices)
//            );
//
//            scalar totalCellSize = 0;
//            label nVerts = 0;
//
//            for
//            (
//                std::list<CellSizeDelaunay::Vertex_handle>::iterator avit =
//                    adjacentVertices.begin();
//                avit != adjacentVertices.end();
//                ++avit
//            )
//            {
//                if (!(*avit)->uninitialised())
//                {
//                    totalCellSize += (*avit)->targetCellSize();
//                    nVerts++;
//                }
//            }
//
//            Pout<< "Changing " << vit->info();
//
//            vit->targetCellSize() = totalCellSize/nVerts;
//            vit->type() = Vb::vtInternalNearBoundary;
//
//            Pout<< "to " << vit->info() << endl;
//        }
//    }
}


Foam::label Foam::controlMeshRefinement::refineMesh
(
    const autoPtr<backgroundMeshDecomposition>& decomposition
)
{
    Info<< "Iterate over "
        << returnReduce(label(mesh_.number_of_finite_edges()), sumOp<label>())
        << " cell size mesh edges" << endl;

    DynamicList<Vb> verts(mesh_.number_of_vertices());

    label count = 0;

    for
    (
        CellSizeDelaunay::Finite_edges_iterator eit =
            mesh_.finite_edges_begin();
        eit != mesh_.finite_edges_end();
        ++eit
    )
    {
        if (count % 10000 == 0)
        {
            Info<< count << " edges, inserted " << verts.size()
                << " Time = " << mesh_.time().elapsedCpuTime()
                << endl;
        }
        count++;

        CellSizeDelaunay::Cell_handle c = eit->first;
        CellSizeDelaunay::Vertex_handle vA = c->vertex(eit->second);
        CellSizeDelaunay::Vertex_handle vB = c->vertex(eit->third);

        if
        (
            mesh_.is_infinite(vA)
         || mesh_.is_infinite(vB)
         || (vA->referred() && vB->referred())
         || (vA->referred() && (vA->procIndex() > vB->procIndex()))
         || (vB->referred() && (vB->procIndex() > vA->procIndex()))
        )
        {
            continue;
        }

        pointFromPoint ptA(topoint(vA->point()));
        pointFromPoint ptB(topoint(vB->point()));

        linePointRef l(ptA, ptB);

        const pointHit hitPt = findDiscontinuities(l);

        if (hitPt.hit())
        {
            const Foam::point& pt = hitPt.point();

            if (!geometryToConformTo_.inside(pt))
            {
                continue;
            }

            if (Pstream::parRun())
            {
                if (!decomposition().positionOnThisProcessor(pt))
                {
                    continue;
                }
            }

            verts.append
            (
                Vb
                (
                    toPoint(pt),
                    Vb::vtInternal
                )
            );

            verts.last().targetCellSize() = sizeControls_.cellSize(pt);
            verts.last().alignment() = triad::unset;
        }
    }

    mesh_.insertPoints(verts, false);

    return verts.size();
}


// * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * * //


// * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * * //


// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //


// ************************************************************************* //