meshOctreeAutomaticRefinementRef.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | cfMesh: A library for mesh generation
   \\    /   O peration     |
    \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
     \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
-------------------------------------------------------------------------------
License
    This file is part of cfMesh.
 
    cfMesh 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.
 
    cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
 
Description
 
\*---------------------------------------------------------------------------*/
 
#include "meshOctreeAutomaticRefinement.H"
#include "triSurface.H"
#include "demandDrivenData.H"
#include "triSurfacePartitioner.H"
#include "triSurfaceCurvatureEstimator.H"
#include "meshOctreeAddressing.H"
#include "triSurf.H"
#include "helperFunctions.H"
#include "meshOctreeModifier.H"
 
#include "Map.H"
 
# ifdef USE_OMP
#include <omp.h>
# endif
 
//#define DEBUGAutoRef
 
# ifdef DEBUGAutoRef
#include "pointSet.H"
#include "IOdictionary.H"
#include "objectRegistry.H"
#include "Time.H"
# endif
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
namespace Foam
{
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
void meshOctreeAutomaticRefinement::activateHexRefinement()
{
    hexRefinement_ = true;
}
 
void meshOctreeAutomaticRefinement::automaticRefinement()
{
    Info << "Performing automatic refinement" << endl;
 
    if( !maxRefLevel_ )
        return;
 
    curvatureRefinement();
 
    proximityRefinement();
 
    Info << "Finished with automatic refinement" << endl;
}
 
bool meshOctreeAutomaticRefinement::curvatureRefinement()
{
    labelList refineBox(octree_.numberOfLeaves(), direction(0));
    labelLongList refinementCandidates;
    forAll(refineBox, i)
        refinementCandidates.append(i);
    while( refineBasedOnCurvature(refineBox, refinementCandidates) )
    {
        refineSelectedBoxes(refineBox, refinementCandidates);
    }
 
    return false;
}
 
bool meshOctreeAutomaticRefinement::proximityRefinement()
{
    bool refine(false);
    labelList refineBox(octree_.numberOfLeaves(), direction(0));
    labelLongList refinementCandidates;
    forAll(refineBox, i)
        refinementCandidates.append(i);
    while( refineBasedOnContainedCorners(refineBox, refinementCandidates) )
    {
        refineSelectedBoxes(refineBox, refinementCandidates);
        refine = true;
    }
 
    refinementCandidates.clear();
    forAll(refineBox, i)
        refinementCandidates.append(i);
    while( refineBasedOnContainedPartitions(refineBox, refinementCandidates) )
    {
        refineSelectedBoxes(refineBox, refinementCandidates);
        refine = true;
    }
 
    refinementCandidates.clear();
    forAll(refineBox, i)
        refinementCandidates.append(i);
    while( refineBasedOnProximityTests(refineBox, refinementCandidates) )
    {
        refineSelectedBoxes(refineBox, refinementCandidates);
        refine = true;
    }
 
    return refine;
}
 
bool meshOctreeAutomaticRefinement::refineBasedOnContainedCorners
(
    labelList& refineBox,
    const labelLongList& refCandidates
)
{
    meshOctreeModifier octreeModifier(octree_);
    const LongList<meshOctreeCube*>& leaves = octreeModifier.leavesAccess();
    const boundBox& rootBox = octree_.rootBox();
    const triSurf& surface = octree_.surface();
    const pointField& points = surface.points();
    const triSurfacePartitioner& sPart = this->partitioner();
 
    //- find leaves which contains corner nodes
    labelList cornerInLeaf(refineBox.size(), -1);
    const labelList& corners = sPart.corners();
 
    label nMarked(0);
 
    forAll(corners, cornerI)
    {
        const label cLabel =
        octree_.findLeafContainingVertex(points[corners[cornerI]]);
 
        if( cLabel < 0 )
            continue;
 
        if( cornerInLeaf[cLabel] != -1 )
        {
            // refine this box because it already contains some other corner
            ++nMarked;
            refineBox[cLabel] = 1;
        }
        else
        {
            cornerInLeaf[cLabel] = corners[cornerI];
        }
    }
 
    DynList<label> leavesInBox;
    # ifdef USE_OMP
    # pragma omp parallel for if( refCandidates.size() > 1000 ) \
    private(leavesInBox) shared(cornerInLeaf) \
    reduction(+ : nMarked) schedule(dynamic, 20)
    # endif
    forAll(refCandidates, refI)
    {
        const label leafI = refCandidates[refI];
        if( leaves[leafI]->level() >= maxRefLevel_ )
            continue;
        if( cornerInLeaf[leafI] == -1 )
            continue;
 
        // check if there exist some corners in the neighbour boxes
        // refine the box if some corners are found in the neighbouring boxes
        const point c = leaves[leafI]->centre(rootBox);
        const scalar r = 1.732 * leaves[leafI]->size(rootBox);
 
        boundBox bb(c - point(r, r, r), c + point(r, r, r));
 
        leavesInBox.clear();
        octree_.findLeavesContainedInBox(bb, leavesInBox);
 
        forAll(leavesInBox, i)
        {
            const label nei = leavesInBox[i];
 
            if( nei < 0 )
                continue;
 
            if( nei == leafI )
                continue;
            if( cornerInLeaf[nei] == -1 )
                continue;
 
            if( mag(points[cornerInLeaf[nei]] - c) < r )
            {
                ++nMarked;
                refineBox[nei] = 1;
                refineBox[leafI] = 1;
                break;
            }
        }
    }
 
    reduce(nMarked, sumOp<label>());
    Info << nMarked << " boxes marked by the corner criteria" << endl;
 
    if( nMarked )
        return true;
 
    return false;
}
 
bool meshOctreeAutomaticRefinement::refineBasedOnContainedPartitions
(
    labelList& refineBox,
    const labelLongList& refCandidates
)
{
    const boundBox& rootBox = octree_.rootBox();
    const triSurfacePartitioner& sPart = this->partitioner();
 
    //- find leaves which contains corner nodes
    const List<labelHashSet>& pPatches = sPart.patchPatches();
    const labelList& edgeGroups = sPart.edgeGroups();
    const List<labelHashSet>& eNeiGroups = sPart.edgeGroupEdgeGroups();
 
    # ifdef DEBUGAutoRef
    Info << "pPart " << pPart << endl;
    # endif
 
    label nMarked(0);
 
    meshOctreeModifier octreeModifier(octree_);
    const triSurf& surf = octree_.surface();
    const LongList<meshOctreeCube*>& leaves = octreeModifier.leavesAccess();
 
    DynList<label> patches, eGroups, helper;
    # ifdef USE_OMP
    # pragma omp parallel for if( refCandidates.size() > 1000 ) \
    private(patches, eGroups, helper) \
    reduction(+ : nMarked) schedule(dynamic, 20)
    # endif
    forAll(refCandidates, refI)
    {
        const label leafI = refCandidates[refI];
        if( !leaves[leafI]->hasContainedElements() )
            continue;
        if( leaves[leafI]->level() >= maxRefLevel_ )
            continue;
 
        const meshOctreeCubeBasic& oc = *leaves[leafI];
        const point c = oc.centre(rootBox);
        const scalar s = 1.733 * oc.size(rootBox);
        const boundBox bb(c - point(s, s, s), c + point(s, s, s));
 
        //- find triangle patches contained in this box
        octree_.findTrianglesInBox(bb, helper);
        patches.clear();
        forAll(helper, i)
            patches.appendIfNotIn(surf[helper[i]].region());
 
        //- find edge partitions contained in this box
        helper.clear();
        octree_.findEdgesInBox(bb, helper);
        eGroups.clear();
        forAll(helper, i)
            eGroups.appendIfNotIn(edgeGroups[helper[i]]);
 
        # ifdef DEBUGAutoRef
        Info << "patches for leaf " << leafI << " are " << patches << endl;
        # endif
 
        bool refine(false);
        forAll(patches, patchI)
        {
            for(label patchJ=(patchI+1);patchJ<patches.size();++patchJ)
                if( !pPatches[patches[patchI]].found(patches[patchJ]) )
                {
                    # ifdef DEBUGAutoRef
                    Info << "2.Here" << endl;
                    # endif
 
                    refine = true;
                    break;
                }
        }
 
        forAll(eGroups, egI)
        {
            for(label egJ=egI+1;egJ<eGroups.size();++egJ)
                if( !eNeiGroups[eGroups[egI]].found(eGroups[egJ]) )
                {
                    refine = true;
                    break;
                }
        }
 
        if( refine )
        {
            # ifdef DEBUGAutoRef
            Info << "Selecting leaf " << leafI
                << " for auto-refinement" << endl;
            # endif
 
            ++nMarked;
            refineBox[leafI] = 1;
        }
    }
 
    reduce(nMarked, sumOp<label>());
    Info << nMarked << " boxed marked by partitioning criteria" << endl;
 
    if( nMarked )
        return true;
 
    return false;
}
 
bool meshOctreeAutomaticRefinement::refineBasedOnCurvature
(
    labelList& refineBox,
    const labelLongList& refCandidates
)
{
    const triSurfaceCurvatureEstimator& curv = curvature();
 
    const boundBox& rootBox = octree_.rootBox();
 
    label nMarked(0);
    DynList<label> containedTrias;
    # ifdef USE_OMP
    # pragma omp parallel for if( refCandidates.size() > 10000 ) \
    private(containedTrias) \
    reduction(+ : nMarked) schedule(dynamic, 100)
    # endif
    forAll(refCandidates, refI)
    {
        const label leafI = refCandidates[refI];
 
        if( !octree_.hasContainedTriangles(leafI) )
            continue;
 
        const meshOctreeCubeBasic& oc = octree_.returnLeaf(leafI);
        if( oc.level() >= maxRefLevel_ )
            continue;
 
        //- search for the minimum curvature radius at surface triangles
        octree_.containedTriangles(leafI, containedTrias);
 
        scalar maxCurv(0.0);
        forAll(containedTrias, i)
            maxCurv =
                Foam::max
                (
                    maxCurv,
                    mag(curv.meanCurvatureAtTriangle(containedTrias[i]))
                );
 
        //- check the edge curvature
        if( octree_.hasContainedEdges(leafI) )
        {
            octree_.containedEdges(leafI, containedTrias);
 
            forAll(containedTrias, i)
            {
                maxCurv =
                    Foam::max
                    (
                        maxCurv,
                        mag(curv.curvatureAtEdge(containedTrias[i]))
                    );
            }
        }
 
        if( oc.size(rootBox) > 0.2835 / (maxCurv + SMALL) )
        {
            refineBox[leafI] = 1;
            ++nMarked;
        }
    }
 
    reduce(nMarked, sumOp<label>());
    Info << nMarked << " boxes marked by curvature criteria!" << endl;
 
    if( nMarked )
        return true;
 
    return false;
}
 
bool meshOctreeAutomaticRefinement::refineBasedOnProximityTests
(
    labelList& refineBox,
    const labelLongList& refCandidates
)
{
    const boundBox& rootBox = octree_.rootBox();
    const triSurf& surf = octree_.surface();
 
    label nMarked(0);
    DynList<label> helper;
    # ifdef USE_OMP
    # pragma omp parallel for if( refCandidates.size() > 1000 ) \
    private(helper) reduction(+ : nMarked) schedule(dynamic, 20)
    # endif
    forAll(refCandidates, refI)
    {
        const label leafI = refCandidates[refI];
 
        if( !octree_.hasContainedTriangles(leafI) )
            continue;
 
        const meshOctreeCubeBasic& oc = octree_.returnLeaf(leafI);
        if( oc.level() >= maxRefLevel_ )
            continue;
 
        const point c = oc.centre(rootBox);
        const scalar s = 1.732 * oc.size(rootBox);
        boundBox bb(c - point(s, s, s), c + point(s, s, s));
 
        labelHashSet triaInRange(100), edgesInRange(100);
        //- find triangles in range
        helper.clear();
        octree_.findTrianglesInBox(bb, helper);
        forAll(helper, i)
            triaInRange.insert(helper[i]);
 
        //- find edges contained in the neighbourhood
        helper.clear();
        octree_.findEdgesInBox(bb, helper);
        forAll(helper, i)
            edgesInRange.insert(helper[i]);
 
        //- refine boxes with more than two face groups
        if
        (
            (help::numberOfFaceGroups(triaInRange, c, s, surf) > 1) ||
            (help::numberOfEdgeGroups(edgesInRange, c, s, surf) > 1)
        )
        {
            ++nMarked;
            refineBox[leafI] = 1;
        }
    }
 
    reduce(nMarked, sumOp<label>());
    Info << nMarked << " boxed marked by proximity criteria" << endl;
 
    if( nMarked != 0 )
        return true;
 
    return false;
}
 
void meshOctreeAutomaticRefinement::refineSelectedBoxes
(
    labelList& refineBox,
    labelLongList& refCandidates
)
{
    deleteDemandDrivenData(octreeAddressingPtr_);
 
    meshOctreeModifier octreeModifier(octree_);
    LongList<meshOctreeCube*> leaves = octreeModifier.leavesAccess();
 
    octreeModifier.markAdditionalLayers(refineBox, 1);
    octreeModifier.refineSelectedBoxes(refineBox, hexRefinement_);
 
    //- find the cubes which have been marked for refinement
    LongList<meshOctreeCubeCoordinates> refinedCubes;
    forAll(refineBox, i)
    {
        if( refineBox[i] )
            refinedCubes.append(leaves[i]->coordinates());
    }
    leaves.setSize(0);
 
    //- perform load distribution in case od parallel runs
    octreeModifier.loadDistribution();
 
    //- communicate the cubes selected for refinement with other processors
    LongList<meshOctreeCubeCoordinates> receivedCoordinates;
    if( Pstream::parRun() )
    {
        octree_.exchangeRequestsWithNeighbourProcessors
        (
            refinedCubes,
            receivedCoordinates
        );
    }
 
    forAll(refinedCubes, i)
        receivedCoordinates.append(refinedCubes[i]);
    refinedCubes.setSize(0);
 
    //- find the cubes which shall checked in the next iteration
    refCandidates.clear();
    forAll(receivedCoordinates, i)
    {
        const meshOctreeCubeCoordinates& cc = receivedCoordinates[i];
 
        for(label scI=0;scI<8;++scI)
        {
            const meshOctreeCubeCoordinates child = cc.refineForPosition(scI);
 
            meshOctreeCube* oc = octreeModifier.findCubeForPosition(child);
 
            if( !oc || !oc->isLeaf() )
                continue;
 
            refCandidates.append(oc->cubeLabel());
        }
    }
 
    refineBox.setSize(octree_.numberOfLeaves());
    refineBox = direction(0);
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
} // End namespace Foam
 
// ************************************************************************* //