meshRefinementGapRefine.C

    List<volumeType>& shellGapMode
) const
{
    const labelList& cellLevel = meshCutter_.cellLevel();
    const pointField& cellCentres = mesh_.cellCentres();

    // Collect cells to test
    cellMap.setSize(cellLevel.size()-nRefine);
    label compactI = 0;

    forAll(cellLevel, cellI)
    {
        if (refineCell[cellI] == -1)
        {
            cellMap[compactI++] = cellI;
        }
    }
    Info<< "Selected " << returnReduce(compactI, sumOp<label>())
        << " unmarked cells out of "
        << mesh_.globalData().nTotalCells() << endl;
    cellMap.setSize(compactI);

    // Do test to see whether cells are inside/outside shell with
    // applicable specification (minLevel <= celllevel < maxLevel)
    shells_.findHigherGapLevel
    (
        pointField(cellCentres, cellMap),
        UIndirectList<label>(cellLevel, cellMap)(),
        shellGapInfo,
        shellGapMode
    );

    // Compact out hits

    labelList map(shellGapInfo.size());
    compactI = 0;
    forAll(shellGapInfo, i)
    {
        if (shellGapInfo[i][2] > 0)
        {
            map[compactI++] = i;
        }
    }

    Info<< "Selected " << returnReduce(compactI, sumOp<label>())
        << " cells inside gap shells out of "
        << mesh_.globalData().nTotalCells() << endl;

    map.setSize(compactI);
    cellMap = UIndirectList<label>(cellMap, map)();
    shellGapInfo = UIndirectList<FixedList<label, 3>>(shellGapInfo, map)();
    shellGapMode = UIndirectList<volumeType>(shellGapMode, map)();
}


void Foam::meshRefinement::mergeGapInfo
(
    const FixedList<label, 3>& shellGapInfo,
    const volumeType shellGapMode,
    const FixedList<label, 3>& surfGapInfo,
    const volumeType surfGapMode,
    FixedList<label, 3>& gapInfo,
    volumeType& gapMode
) const
{
    if (surfGapInfo[0] == 0)
    {
        gapInfo = shellGapInfo;
        gapMode = shellGapMode;
    }
    else if (shellGapInfo[0] == 0)
    {
        gapInfo = surfGapInfo;
        gapMode = surfGapMode;
    }
    else
    {
        // Both specify a level. Does surface level win? Or does information
        // need to be merged?

        //gapInfo[0] = max(surfGapInfo[0], shellGapInfo[0]);
        //gapInfo[1] = min(surfGapInfo[1], shellGapInfo[1]);
        //gapInfo[2] = max(surfGapInfo[2], shellGapInfo[2]);
        gapInfo = surfGapInfo;
        gapMode = surfGapMode;
    }
}


Foam::label Foam::meshRefinement::markInternalGapRefinement
(
    const scalar planarCos,
    const bool spreadGapSize,
    const label nAllowRefine,

    labelList& refineCell,
    label& nRefine,
    labelList& numGapCells,
    scalarField& detectedGapSize
) const
{
    detectedGapSize.setSize(mesh_.nCells());
    detectedGapSize = GREAT;
    numGapCells.setSize(mesh_.nCells());
    numGapCells = -1;

    const labelList& cellLevel = meshCutter_.cellLevel();
    const pointField& cellCentres = mesh_.cellCentres();
    const scalar edge0Len = meshCutter_.level0EdgeLength();

    const List<FixedList<label, 3>>& extendedGapLevel =
        surfaces_.extendedGapLevel();
    const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();

    // Get the gap level for the shells
    const labelList maxLevel(shells_.maxGapLevel());

    label oldNRefine = nRefine;

    if (max(maxLevel) > 0)
    {
        // Collect cells to test
        labelList cellMap;
        List<FixedList<label, 3>> shellGapInfo;
        List<volumeType> shellGapMode;
        selectGapCandidates
        (
            refineCell,
            nRefine,

            cellMap,
            shellGapInfo,
            shellGapMode
        );

        // Find nearest point and normal on the surfaces
        List<pointIndexHit> nearInfo;
        vectorField nearNormal;
        labelList nearSurface;
        labelList nearRegion;
        {
            // Now we have both the cell-level and the gap size information. Use
            // this to calculate the gap size
            scalarField gapSize(cellMap.size());
            forAll(cellMap, i)
            {
                label cellI = cellMap[i];
                scalar cellSize = edge0Len/pow(2.0, cellLevel[cellI]);
                gapSize[i] = shellGapInfo[i][0]*cellSize;
            }

            surfaces_.findNearestRegion
            (
                identity(surfaces_.surfaces().size()),
                pointField(cellCentres, cellMap),
                sqr(gapSize),
                nearSurface,
                nearInfo,
                nearRegion,
                nearNormal
            );
        }



        DynamicList<label> map(nearInfo.size());
        DynamicField<point> rayStart(nearInfo.size());
        DynamicField<point> rayEnd(nearInfo.size());
        DynamicField<scalar> gapSize(nearInfo.size());

        DynamicField<point> rayStart2(nearInfo.size());
        DynamicField<point> rayEnd2(nearInfo.size());
        DynamicField<scalar> gapSize2(nearInfo.size());

        label nTestCells = 0;

        forAll(nearInfo, i)
        {
            if (nearInfo[i].hit())
            {
                label globalRegionI = surfaces_.globalRegion
                (
                    nearSurface[i],
                    nearRegion[i]
                );

                // Combine info from shell and surface
                FixedList<label, 3> gapInfo;
                volumeType gapMode;
                mergeGapInfo
                (
                    shellGapInfo[i],
                    shellGapMode[i],
                    extendedGapLevel[globalRegionI],
                    extendedGapMode[globalRegionI],
                    gapInfo,
                    gapMode
                );

                // Store wanted number of cells in gap
                label cellI = cellMap[i];
                label cLevel = cellLevel[cellI];
                if (cLevel >= gapInfo[1] && cLevel < gapInfo[2])
                {
                    numGapCells[cellI] = max(numGapCells[cellI], gapInfo[0]);
                }

                // Construct one or more rays to test for oppositeness
                label nRays = generateRays
                (
                    false,
                    nearInfo[i].hitPoint(),
                    nearNormal[i],
                    gapInfo,
                    gapMode,

                    cellCentres[cellI],
                    cLevel,

                    rayStart,
                    rayEnd,
                    gapSize,

                    rayStart2,
                    rayEnd2,
                    gapSize2
                );
                if (nRays > 0)
                {
                    nTestCells++;
                    for (label j = 0; j < nRays; j++)
                    {
                        map.append(i);
                    }
                }
            }
        }

        Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
            << " cells for testing out of "
            << mesh_.globalData().nTotalCells() << endl;
        map.shrink();
        rayStart.shrink();
        rayEnd.shrink();
        gapSize.shrink();

        rayStart2.shrink();
        rayEnd2.shrink();
        gapSize2.shrink();

        cellMap = UIndirectList<label>(cellMap, map)();
        nearNormal = UIndirectList<vector>(nearNormal, map)();
        shellGapInfo.clear();
        shellGapMode.clear();
        nearInfo.clear();
        nearSurface.clear();
        nearRegion.clear();


        // Do intersections in pairs
        labelList surf1;
        List<pointIndexHit> hit1;
        vectorField normal1;
        surfaces_.findNearestIntersection
        (
            rayStart,
            rayEnd,
            surf1,
            hit1,
            normal1
        );

        labelList surf2;
        List<pointIndexHit> hit2;
        vectorField normal2;
        surfaces_.findNearestIntersection
        (
            rayStart2,
            rayEnd2,
            surf2,
            hit2,
            normal2
        );

        // Extract cell based gap size
        forAll(surf1, i)
        {
            if (surf1[i] != -1 && surf2[i] != -1)
            {
                // Found intersections with surface. Check for
                // - small gap
                // - coplanar normals

                label cellI = cellMap[i];

                scalar d2 = magSqr(hit1[i].hitPoint()-hit2[i].hitPoint());

                if
                (
                    cellI != -1
                 && (mag(normal1[i]&normal2[i]) > planarCos)
                 && (d2 < Foam::sqr(gapSize[i]))
                )
                {
                    detectedGapSize[cellI] = min
                    (
                        detectedGapSize[cellI],
                        Foam::sqrt(d2)
                    );
                }
            }
        }

        // Spread it
        if (spreadGapSize)
        {
            // Field on cells and faces
            List<transportData> cellData(mesh_.nCells());
            List<transportData> faceData(mesh_.nFaces());

            // Start of walk
            const pointField& faceCentres = mesh_.faceCentres();

            DynamicList<label> frontFaces(mesh_.nFaces());
            DynamicList<transportData> frontData(mesh_.nFaces());
            for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
            {
                label own = mesh_.faceOwner()[faceI];
                label nei = mesh_.faceNeighbour()[faceI];

                scalar minSize = min
                (
                    detectedGapSize[own],
                    detectedGapSize[nei]
                );

                if (minSize < GREAT)
                {
                    frontFaces.append(faceI);
                    frontData.append
                    (
                        transportData
                        (
                            faceCentres[faceI],
                            minSize,
                            0.0
                        )
                    );
                }
            }
            for
            (
                label faceI = mesh_.nInternalFaces();
                faceI < mesh_.nFaces();
                faceI++
            )
            {
                label own = mesh_.faceOwner()[faceI];

                if (detectedGapSize[own] < GREAT)
                {
                    frontFaces.append(faceI);
                    frontData.append
                    (
                        transportData
                        (
                            faceCentres[faceI],
                            detectedGapSize[own],
                            0.0
                        )
                    );
                }
            }

            Info<< "Selected "
                << returnReduce(frontFaces.size(), sumOp<label>())
                << " faces for spreading gap size out of "
                << mesh_.globalData().nTotalFaces() << endl;


            transportData::trackData td(surfaceIndex_);

            FaceCellWave<transportData, transportData::trackData> deltaCalc
            (
                mesh_,
                frontFaces,
                frontData,
                faceData,
                cellData,
                mesh_.globalData().nTotalCells()+1,
                td
            );


            forAll(cellMap, i)
            {
                label cellI = cellMap[i];
                if
                (
                    cellI != -1
                 && cellData[cellI].valid(deltaCalc.data())
                 && numGapCells[cellI] != -1
                )
                {
                    // Update transported gap size
                    detectedGapSize[cellI] = min
                    (
                        detectedGapSize[cellI],
                        cellData[cellI].data()
                    );
                }
            }
        }


        // Use it
        forAll(cellMap, i)
        {
            label cellI = cellMap[i];

            if (cellI != -1 && numGapCells[cellI] != -1)
            {
                // Needed gap size
                label cLevel = cellLevel[cellI];
                scalar cellSize =
                    meshCutter_.level0EdgeLength()/pow(2.0, cLevel);
                scalar neededGapSize = numGapCells[cellI]*cellSize;

                if (neededGapSize > detectedGapSize[cellI])
                {
                    if
                    (
                       !markForRefine
                        (
                            123,
                            nAllowRefine,
                            refineCell[cellI],
                            nRefine
                        )
                    )
                    {
                        break;
                    }
                }
            }
        }


        if
        (
            returnReduce(nRefine, sumOp<label>())
          > returnReduce(nAllowRefine, sumOp<label>())
        )
        {
            Info<< "Reached refinement limit." << endl;
        }
    }

    return returnReduce(nRefine-oldNRefine, sumOp<label>());
}


Foam::label Foam::meshRefinement::markSmallFeatureRefinement
(
    const scalar planarCos,
    const label nAllowRefine,
    const labelList& neiLevel,
    const pointField& neiCc,

    labelList& refineCell,
    label& nRefine
) const
{
    const labelList& cellLevel = meshCutter_.cellLevel();
    const labelList& surfaceIndices = surfaces_.surfaces();
    const List<FixedList<label, 3>>& extendedGapLevel =
        surfaces_.extendedGapLevel();
    const List<volumeType>& extendedGapMode = surfaces_.extendedGapMode();

    label oldNRefine = nRefine;

    // Check that we're using any gap refinement
    labelList shellMaxLevel(shells_.maxGapLevel());

    if (max(shellMaxLevel) == 0)
    {
        return 0;
    }

    //- Force calculation of tetBasePt
    (void)mesh_.tetBasePtIs();
    (void)mesh_.cellTree();


    forAll(surfaceIndices, surfI)
    {
        label geomI = surfaceIndices[surfI];
        const searchableSurface& geom = surfaces_.geometry()[geomI];


        // Get the element index in a roundabout way. Problem is e.g.
        // distributed surface where local indices differ from global
        // ones (needed for getRegion call)

        pointField ctrs;
        labelList region;
        vectorField normal;
        {
            // Representative local coordinates and bounding sphere
            scalarField radiusSqr;
            geom.boundingSpheres(ctrs, radiusSqr);

            List<pointIndexHit> info;
            geom.findNearest(ctrs, radiusSqr, info);

            forAll(info, i)
            {
                if (!info[i].hit())
                {
                    FatalErrorInFunction
                        << "fc:" << ctrs[i]
                        << " radius:" << radiusSqr[i]
                        << exit(FatalError);
                }
            }

            geom.getRegion(info, region);
            geom.getNormal(info, normal);
        }

        // Do test to see whether triangles are inside/outside shell with
        // applicable specification (minLevel <= celllevel < maxLevel)
        List<FixedList<label, 3>> shellGapInfo;
        List<volumeType> shellGapMode;
        shells_.findHigherGapLevel
        (
            ctrs,
            labelList(ctrs.size(), 0),
            shellGapInfo,
            shellGapMode
        );


        DynamicList<label> map(ctrs.size());
        DynamicList<label> cellMap(ctrs.size());

        DynamicField<point> rayStart(ctrs.size());
        DynamicField<point> rayEnd(ctrs.size());
        DynamicField<scalar> gapSize(ctrs.size());

        label nTestCells = 0;

        forAll(ctrs, i)
        {
            if (shellGapInfo[i][2] > 0)
            {
                label globalRegionI = surfaces_.globalRegion(surfI, region[i]);

                // Combine info from shell and surface
                FixedList<label, 3> gapInfo;
                volumeType gapMode;
                mergeGapInfo
                (
                    shellGapInfo[i],
                    shellGapMode[i],
                    extendedGapLevel[globalRegionI],
                    extendedGapMode[globalRegionI],
                    gapInfo,
                    gapMode
                );

                //- Option 1: use octree nearest searching inside polyMesh
                //label cellI = mesh_.findCell(pt);

                //- Option 2: use octree 'inside' searching inside polyMesh. Is
                //            much faster.
                label cellI = -1;
                const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
                if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
                {
                    cellI = tree.findInside(ctrs[i]);
                }

                if (cellI != -1 && refineCell[cellI] == -1)
                {
                    // Construct one or two rays to test for oppositeness
                    // Note that we always want to use the surface normal
                    // and not the vector from cell centre to surface point

                    label nRays = generateRays
                    (
                        ctrs[i],
                        normal[i],
                        gapInfo,
                        gapMode,

                        cellLevel[cellI],

                        rayStart,
                        rayEnd
                    );

                    if (nRays > 0)
                    {
                        nTestCells++;
                        for (label j = 0; j < nRays; j++)
                        {
                            cellMap.append(cellI);
                            map.append(i);
                        }
                    }
                }
            }
        }

        Info<< "Selected " << returnReduce(nTestCells, sumOp<label>())
            << " cells containing triangle centres out of "
            << mesh_.globalData().nTotalCells() << endl;
        map.shrink();
        cellMap.shrink();
        rayStart.shrink();
        rayEnd.shrink();

        ctrs.clear();
        region.clear();
        shellGapInfo.clear();
        shellGapMode.clear();
        normal = UIndirectList<vector>(normal, map)();

        // Do intersections.
        labelList surfaceHit;
        vectorField surfaceNormal;
        surfaces_.findNearestIntersection
        (
            rayStart,
            rayEnd,
            surfaceHit,
            surfaceNormal
        );


        label nOldRefine = 0;


        forAll(surfaceHit, i)
        {
            if (surfaceHit[i] != -1) // && surf2[i] != -1)
            {
                // Found intersection with surface. Check coplanar normals.
                label cellI = cellMap[i];

                if (mag(normal[i]&surfaceNormal[i]) > planarCos)
                {
                    if
                    (
                       !markForRefine
                        (
                            surfaceHit[i],
                            nAllowRefine,
                            refineCell[cellI],
                            nRefine
                        )
                    )
                    {
                        break;
                    }
                }
            }
        }

        Info<< "For surface " << geom.name() << " found "
            << returnReduce(nRefine-nOldRefine, sumOp<label>())
            << " cells in small gaps" << endl;

        if
        (
            returnReduce(nRefine, sumOp<label>())
          > returnReduce(nAllowRefine, sumOp<label>())
        )
        {
            Info<< "Reached refinement limit." << endl;
        }
    }

    return returnReduce(nRefine-oldNRefine, sumOp<label>());
}


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