Skip to content
Snippets Groups Projects
autoSnapDriver.C 39.9 KiB
Newer Older
  • Learn to ignore specific revisions
  • /*---------------------------------------------------------------------------*\
      =========                 |
      \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       \\    /   O peration     |
    
        \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
    
         \\/     M anipulation  |
    -------------------------------------------------------------------------------
    License
        This file is part of OpenFOAM.
    
    
        OpenFOAM is free software: you can redistribute it and/or modify it
        under the terms of the GNU General Public License as published by
        the Free Software Foundation, either version 3 of the License, or
        (at your option) any later version.
    
    
        OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
        ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
        FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
        for more details.
    
        You should have received a copy of the GNU General Public License
    
        along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
    
    
    Description
        All to do with snapping to the surface
    
    \*----------------------------------------------------------------------------*/
    
    #include "autoSnapDriver.H"
    #include "motionSmoother.H"
    #include "polyTopoChange.H"
    #include "OFstream.H"
    #include "syncTools.H"
    #include "fvMesh.H"
    #include "Time.H"
    #include "OFstream.H"
    #include "mapPolyMesh.H"
    #include "pointEdgePoint.H"
    #include "PointEdgeWave.H"
    #include "mergePoints.H"
    #include "snapParameters.H"
    #include "refinementSurfaces.H"
    
    #include "unitConversion.H"
    
    
    // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
    
    
    defineTypeNameAndDebug(autoSnapDriver, 0);
    
    // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
    
    
    // Calculate geometrically collocated points, Requires PackedList to be
    
    Foam::label Foam::autoSnapDriver::getCollocatedPoints
    (
        const scalar tol,
        const pointField& points,
    
        PackedBoolList& isCollocatedPoint
    
    )
    {
        labelList pointMap;
        pointField newPoints;
        bool hasMerged = mergePoints
        (
            points,                         // points
            tol,                            // mergeTol
            false,                          // verbose
            pointMap,
            newPoints
        );
    
        if (!returnReduce(hasMerged, orOp<bool>()))
        {
            return 0;
        }
    
        // Determine which newPoints are referenced more than once
        label nCollocated = 0;
    
        // Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
        // twice)
        labelList firstOldPoint(newPoints.size(), -1);
        forAll(pointMap, oldPointI)
        {
            label newPointI = pointMap[oldPointI];
    
            if (firstOldPoint[newPointI] == -1)
            {
                // First use of oldPointI. Store.
                firstOldPoint[newPointI] = oldPointI;
            }
            else if (firstOldPoint[newPointI] == -2)
            {
                // Third or more reference of oldPointI -> non-manifold
                isCollocatedPoint.set(oldPointI, 1u);
                nCollocated++;
            }
            else
            {
                // Second reference of oldPointI -> non-manifold
                isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
                nCollocated++;
    
                isCollocatedPoint.set(oldPointI, 1u);
                nCollocated++;
    
                // Mark with special value to save checking next time round
                firstOldPoint[newPointI] = -2;
            }
        }
        return returnReduce(nCollocated, sumOp<label>());
    }
    
    
    // Calculate displacement as average of patch points.
    Foam::pointField Foam::autoSnapDriver::smoothPatchDisplacement
    (
    
        const motionSmoother& meshMover,
        const List<labelPair>& baffles
    
    ) const
    {
        const indirectPrimitivePatch& pp = meshMover.patch();
    
        // Calculate geometrically non-manifold points on the patch to be moved.
    
        PackedBoolList nonManifoldPoint(pp.nPoints());
    
        label nNonManifoldPoints = getCollocatedPoints
        (
            SMALL,
            pp.localPoints(),
            nonManifoldPoint
        );
        Info<< "Found " << nNonManifoldPoints << " non-mainfold point(s)."
            << endl;
    
    
        // Average points
        // ~~~~~~~~~~~~~~
    
        // We determine three points:
        // - average of (centres of) connected patch faces
        // - average of (centres of) connected internal mesh faces
        // - as fallback: centre of any connected cell
        // so we can do something moderately sensible for non/manifold points.
    
        // Note: the averages are calculated properly parallel. This is
        // necessary to get the points shared by processors correct.
    
    
        const labelListList& pointFaces = pp.pointFaces();
        const labelList& meshPoints = pp.meshPoints();
        const pointField& points = pp.points();
        const polyMesh& mesh = meshMover.mesh();
    
    
        // Get labels of faces to count (master of coupled faces and baffle pairs)
    
        PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
    
    
        {
            forAll(baffles, i)
            {
                label f0 = baffles[i].first();
                label f1 = baffles[i].second();
    
    
                if (isMasterFace.get(f0))
    
                {
                    // Make f1 a slave
    
                    isMasterFace.unset(f1);
    
                else if (isMasterFace.get(f1))
    
                    isMasterFace.unset(f0);
    
                }
                else
                {
                    FatalErrorIn("autoSnapDriver::smoothPatchDisplacement(..)")
                        << "Both sides of baffle consisting of faces " << f0
                        << " and " << f1 << " are already slave faces."
                        << abort(FatalError);
                }
            }
        }
    
    
    
        // Get average position of boundary face centres
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
        vectorField avgBoundary(pointFaces.size(), vector::zero);
        labelList nBoundary(pointFaces.size(), 0);
    
        forAll(pointFaces, patchPointI)
        {
            const labelList& pFaces = pointFaces[patchPointI];
    
            forAll(pFaces, pfI)
            {
    
                label faceI = pFaces[pfI];
    
    
                if (isMasterFace.get(pp.addressing()[faceI]))
    
                {
                    avgBoundary[patchPointI] += pp[faceI].centre(points);
                    nBoundary[patchPointI]++;
                }
    
            }
        }
    
        syncTools::syncPointList
        (
            mesh,
            pp.meshPoints(),
            avgBoundary,
            plusEqOp<point>(),  // combine op
    
    mattijs's avatar
    mattijs committed
            vector::zero        // null value
    
        );
        syncTools::syncPointList
        (
            mesh,
            pp.meshPoints(),
            nBoundary,
            plusEqOp<label>(),  // combine op
    
    mattijs's avatar
    mattijs committed
            0                   // null value
    
        );
    
        forAll(avgBoundary, i)
        {
            avgBoundary[i] /= nBoundary[i];
        }
    
    
        // Get average position of internal face centres
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
        vectorField avgInternal;
        labelList nInternal;
        {
            vectorField globalSum(mesh.nPoints(), vector::zero);
            labelList globalNum(mesh.nPoints(), 0);
    
            // Note: no use of pointFaces
            const faceList& faces = mesh.faces();
    
            for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
            {
                const face& f = faces[faceI];
                const point& fc = mesh.faceCentres()[faceI];
    
                forAll(f, fp)
                {
                    globalSum[f[fp]] += fc;
                    globalNum[f[fp]]++;
                }
            }
    
            // Count coupled faces as internal ones (but only once)
            const polyBoundaryMesh& patches = mesh.boundaryMesh();
    
            forAll(patches, patchI)
            {
    
                if
                (
                    patches[patchI].coupled()
                 && refCast<const coupledPolyPatch>(patches[patchI]).owner()
                )
    
                    const coupledPolyPatch& pp =
                        refCast<const coupledPolyPatch>(patches[patchI]);
    
    
                    const vectorField::subField faceCentres = pp.faceCentres();
    
    
                    forAll(pp, i)
    
                    {
                        const face& f = pp[i];
                        const point& fc = faceCentres[i];
    
                        forAll(f, fp)
                        {
                            globalSum[f[fp]] += fc;
                            globalNum[f[fp]]++;
                        }
                    }
                }
            }
    
            syncTools::syncPointList
            (
                mesh,
                globalSum,
                plusEqOp<vector>(), // combine op
    
    mattijs's avatar
    mattijs committed
                vector::zero        // null value
    
            );
            syncTools::syncPointList
            (
                mesh,
                globalNum,
                plusEqOp<label>(),  // combine op
    
    mattijs's avatar
    mattijs committed
                0                   // null value
    
            );
    
            avgInternal.setSize(meshPoints.size());
            nInternal.setSize(meshPoints.size());
    
            forAll(avgInternal, patchPointI)
            {
                label meshPointI = meshPoints[patchPointI];
    
                nInternal[patchPointI] = globalNum[meshPointI];
    
                if (nInternal[patchPointI] == 0)
                {
                    avgInternal[patchPointI] = globalSum[meshPointI];
                }
                else
                {
                    avgInternal[patchPointI] =
                        globalSum[meshPointI]
                      / nInternal[patchPointI];
                }
            }
        }
    
    
        // Precalculate any cell using mesh point (replacement of pointCells()[])
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
        labelList anyCell(mesh.nPoints(), -1);
        forAll(mesh.faceNeighbour(), faceI)
        {
            label own = mesh.faceOwner()[faceI];
            const face& f = mesh.faces()[faceI];
    
            forAll(f, fp)
            {
                anyCell[f[fp]] = own;
            }
        }
        for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
        {
            label own = mesh.faceOwner()[faceI];
    
            const face& f = mesh.faces()[faceI];
    
            forAll(f, fp)
            {
                anyCell[f[fp]] = own;
            }
        }
    
    
        // Displacement to calculate.
        pointField patchDisp(meshPoints.size(), vector::zero);
    
        forAll(pointFaces, i)
        {
            label meshPointI = meshPoints[i];
            const point& currentPos = pp.points()[meshPointI];
    
            // Now we have the two average points: avgBoundary and avgInternal
            // and how many boundary/internal faces connect to the point
            // (nBoundary, nInternal)
            // Do some blending between the two.
            // Note: the following section has some reasoning behind it but the
            // blending factors can be experimented with.
    
            point newPos;
    
    
            if (!nonManifoldPoint.get(i))
    
            {
                // Points that are manifold. Weight the internal and boundary
                // by their number of faces and blend with
                scalar internalBlend = 0.1;
                scalar blend = 0.1;
    
                point avgPos =
                    (
                       internalBlend*nInternal[i]*avgInternal[i]
                      +(1-internalBlend)*nBoundary[i]*avgBoundary[i]
                    )
                  / (internalBlend*nInternal[i]+(1-internalBlend)*nBoundary[i]);
    
                newPos = (1-blend)*avgPos + blend*currentPos;
            }
            else if (nInternal[i] == 0)
            {
                // Non-manifold without internal faces. Use any connected cell
                // as internal point instead. Use precalculated any cell to avoid
                // e.g. pointCells()[meshPointI][0]
    
                const point& cc = mesh.cellCentres()[anyCell[meshPointI]];
    
                scalar cellCBlend = 0.8;
                scalar blend = 0.1;
    
                point avgPos = (1-cellCBlend)*avgBoundary[i] + cellCBlend*cc;
    
                newPos = (1-blend)*avgPos + blend*currentPos;
            }
            else
            {
                // Non-manifold point with internal faces connected to them
                scalar internalBlend = 0.9;
                scalar blend = 0.1;
    
                point avgPos =
                    internalBlend*avgInternal[i]
                  + (1-internalBlend)*avgBoundary[i];
    
                newPos = (1-blend)*avgPos + blend*currentPos;
            }
    
            patchDisp[i] = newPos - currentPos;
        }
    
        return patchDisp;
    }
    
    
    Foam::tmp<Foam::scalarField> Foam::autoSnapDriver::edgePatchDist
    (
        const pointMesh& pMesh,
        const indirectPrimitivePatch& pp
    )
    {
        const polyMesh& mesh = pMesh();
    
        // Set initial changed points to all the patch points
        List<pointEdgePoint> wallInfo(pp.nPoints());
    
        forAll(pp.localPoints(), ppI)
        {
            wallInfo[ppI] = pointEdgePoint(pp.localPoints()[ppI], 0.0);
        }
    
        // Current info on points
        List<pointEdgePoint> allPointInfo(mesh.nPoints());
    
        // Current info on edges
        List<pointEdgePoint> allEdgeInfo(mesh.nEdges());
    
        PointEdgeWave<pointEdgePoint> wallCalc
        (
    
            pp.meshPoints(),
            wallInfo,
    
            allPointInfo,
            allEdgeInfo,
            mesh.globalData().nTotalPoints()  // max iterations
        );
    
        // Copy edge values into scalarField
        tmp<scalarField> tedgeDist(new scalarField(mesh.nEdges()));
        scalarField& edgeDist = tedgeDist();
    
        forAll(allEdgeInfo, edgeI)
        {
            edgeDist[edgeI] = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
        }
    
    
        //{
        //    // For debugging: dump to file
        //    pointScalarField pointDist
        //    (
        //        IOobject
        //        (
        //            "pointDist",
    
        //            meshRefiner_.timeName(),
    
        //            mesh.DB(),
        //            IOobject::NO_READ,
        //            IOobject::AUTO_WRITE
        //        ),
        //        pMesh,
        //        dimensionedScalar("pointDist", dimless, 0.0)
        //    );
        //
        //    forAll(allEdgeInfo, edgeI)
        //    {
        //        scalar d = Foam::sqrt(allEdgeInfo[edgeI].distSqr());
        //
        //        const edge& e = mesh.edges()[edgeI];
        //
        //        pointDist[e[0]] += d;
        //        pointDist[e[1]] += d;
        //    }
        //    forAll(pointDist, pointI)
        //    {
        //        pointDist[pointI] /= mesh.pointEdges()[pointI].size();
        //    }
        //    Info<< "Writing patch distance to " << pointDist.name()
    
        //        << " at time " << meshRefiner_.timeName() << endl;
    
        //
        //    pointDist.write();
        //}
    
        return tedgeDist;
    }
    
    
    void Foam::autoSnapDriver::dumpMove
    (
        const fileName& fName,
        const pointField& meshPts,
        const pointField& surfPts
    )
    {
        // Dump direction of growth into file
    
        Pout<< nl << "Dumping move direction to " << fName << endl;
    
    
        OFstream nearestStream(fName);
    
        label vertI = 0;
    
        forAll(meshPts, ptI)
        {
            meshTools::writeOBJ(nearestStream, meshPts[ptI]);
            vertI++;
    
            meshTools::writeOBJ(nearestStream, surfPts[ptI]);
            vertI++;
    
            nearestStream<< "l " << vertI-1 << ' ' << vertI << nl;
        }
    }
    
    
    // Check whether all displacement vectors point outwards of patch. Return true
    // if so.
    bool Foam::autoSnapDriver::outwardsDisplacement
    (
        const indirectPrimitivePatch& pp,
        const vectorField& patchDisp
    )
    {
        const vectorField& faceNormals = pp.faceNormals();
        const labelListList& pointFaces = pp.pointFaces();
    
        forAll(pointFaces, pointI)
        {
            const labelList& pFaces = pointFaces[pointI];
    
            vector disp(patchDisp[pointI]);
    
            scalar magDisp = mag(disp);
    
            if (magDisp > SMALL)
            {
                disp /= magDisp;
    
                bool outwards = meshTools::visNormal(disp, faceNormals, pFaces);
    
                if (!outwards)
                {
                    Warning<< "Displacement " << patchDisp[pointI]
                        << " at mesh point " << pp.meshPoints()[pointI]
                        << " coord " << pp.points()[pp.meshPoints()[pointI]]
                        << " points through the surrounding patch faces" << endl;
                    return false;
                }
            }
            else
            {
                //? Displacement small but in wrong direction. Would probably be ok.
            }
        }
        return true;
    }
    
    
    // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
    
    Foam::autoSnapDriver::autoSnapDriver
    (
        meshRefinement& meshRefiner,
        const labelList& globalToPatch
    )
    :
        meshRefiner_(meshRefiner),
        globalToPatch_(globalToPatch)
    {}
    
    
    // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
    
    Foam::autoPtr<Foam::mapPolyMesh> Foam::autoSnapDriver::mergeZoneBaffles
    (
        const List<labelPair>& baffles
    )
    {
    
    mattijs's avatar
    mattijs committed
        labelList zonedSurfaces = meshRefiner_.surfaces().getNamedSurfaces();
    
    
        autoPtr<mapPolyMesh> map;
    
        // No need to sync; all processors will have all same zonedSurfaces.
        label nBaffles = returnReduce(baffles.size(), sumOp<label>());
    
        if (zonedSurfaces.size() && nBaffles > 0)
    
        {
            // Merge any baffles
            Info<< "Converting " << nBaffles << " baffles back into zoned faces ..."
                << endl;
    
            map = meshRefiner_.mergeBaffles(baffles);
    
            Info<< "Converted baffles in = "
                << meshRefiner_.mesh().time().cpuTimeIncrement()
                << " s\n" << nl << endl;
        }
    
        return map;
    }
    
    
    Foam::scalarField Foam::autoSnapDriver::calcSnapDistance
    (
        const snapParameters& snapParams,
        const indirectPrimitivePatch& pp
    ) const
    {
        const edgeList& edges = pp.edges();
        const labelListList& pointEdges = pp.pointEdges();
        const pointField& localPoints = pp.localPoints();
        const fvMesh& mesh = meshRefiner_.mesh();
    
        scalarField maxEdgeLen(localPoints.size(), -GREAT);
    
        forAll(pointEdges, pointI)
        {
            const labelList& pEdges = pointEdges[pointI];
    
            forAll(pEdges, pEdgeI)
            {
                const edge& e = edges[pEdges[pEdgeI]];
    
                scalar len = e.mag(localPoints);
    
                maxEdgeLen[pointI] = max(maxEdgeLen[pointI], len);
            }
        }
    
        syncTools::syncPointList
        (
            mesh,
            pp.meshPoints(),
            maxEdgeLen,
            maxEqOp<scalar>(),  // combine op
    
    mattijs's avatar
    mattijs committed
            -GREAT              // null value
    
        return scalarField(snapParams.snapTol()*maxEdgeLen);
    
    }
    
    
    void Foam::autoSnapDriver::preSmoothPatch
    (
        const snapParameters& snapParams,
        const label nInitErrors,
        const List<labelPair>& baffles,
        motionSmoother& meshMover
    ) const
    {
        const fvMesh& mesh = meshRefiner_.mesh();
    
        labelList checkFaces;
    
        Info<< "Smoothing patch points ..." << endl;
        for
        (
            label smoothIter = 0;
            smoothIter < snapParams.nSmoothPatch();
            smoothIter++
        )
        {
            Info<< "Smoothing iteration " << smoothIter << endl;
            checkFaces.setSize(mesh.nFaces());
            forAll(checkFaces, faceI)
            {
                checkFaces[faceI] = faceI;
            }
    
    
            pointField patchDisp(smoothPatchDisplacement(meshMover, baffles));
    
    
            // The current mesh is the starting mesh to smooth from.
            meshMover.setDisplacement(patchDisp);
    
            meshMover.correct();
    
            scalar oldErrorReduction = -1;
    
            for (label snapIter = 0; snapIter < 2*snapParams.nSnap(); snapIter++)
            {
                Info<< nl << "Scaling iteration " << snapIter << endl;
    
                if (snapIter == snapParams.nSnap())
                {
                    Info<< "Displacement scaling for error reduction set to 0."
                        << endl;
                    oldErrorReduction = meshMover.setErrorReduction(0.0);
                }
    
                // Try to adapt mesh to obtain displacement by smoothly
                // decreasing displacement at error locations.
                if (meshMover.scaleMesh(checkFaces, baffles, true, nInitErrors))
                {
                    Info<< "Successfully moved mesh" << endl;
                    break;
                }
            }
    
            if (oldErrorReduction >= 0)
            {
                meshMover.setErrorReduction(oldErrorReduction);
            }
            Info<< endl;
        }
    
    
        // The current mesh is the starting mesh to smooth from.
        meshMover.correct();
    
        if (debug)
        {
            const_cast<Time&>(mesh.time())++;
    
            Pout<< "Writing patch smoothed mesh to time "
                << meshRefiner_.timeName() << '.' << endl;
            meshRefiner_.write
            (
                debug,
                mesh.time().path()/meshRefiner_.timeName()
            );
            Pout<< "Dumped mesh in = "
                << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
    
        }
    
        Info<< "Patch points smoothed in = "
            << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
    }
    
    
    // Get (pp-local) indices of points that are both on zone and on patched surface
    Foam::labelList Foam::autoSnapDriver::getZoneSurfacePoints
    (
        const indirectPrimitivePatch& pp,
        const word& zoneName
    ) const
    {
        const fvMesh& mesh = meshRefiner_.mesh();
    
        label zoneI = mesh.faceZones().findZoneID(zoneName);
    
        if (zoneI == -1)
        {
            FatalErrorIn
            (
                "autoSnapDriver::getZoneSurfacePoints"
                "(const indirectPrimitivePatch&, const word&)"
            )   << "Cannot find zone " << zoneName
                << exit(FatalError);
        }
    
        const faceZone& fZone = mesh.faceZones()[zoneI];
    
    
        // Could use PrimitivePatch & localFaces to extract points but might just
        // as well do it ourselves.
    
        boolList pointOnZone(pp.nPoints(), false);
    
        forAll(fZone, i)
        {
            const face& f = mesh.faces()[fZone[i]];
    
            forAll(f, fp)
            {
                label meshPointI = f[fp];
    
                Map<label>::const_iterator iter =
                    pp.meshPointMap().find(meshPointI);
    
                if (iter != pp.meshPointMap().end())
                {
                    label pointI = iter();
                    pointOnZone[pointI] = true;
                }
            }
        }
    
        return findIndices(pointOnZone, true);
    }
    
    
    Foam::vectorField Foam::autoSnapDriver::calcNearestSurface
    (
        const scalarField& snapDist,
        motionSmoother& meshMover
    ) const
    {
        Info<< "Calculating patchDisplacement as distance to nearest surface"
            << " point ..." << endl;
    
        const indirectPrimitivePatch& pp = meshMover.patch();
        const pointField& localPoints = pp.localPoints();
        const refinementSurfaces& surfaces = meshRefiner_.surfaces();
        const fvMesh& mesh = meshRefiner_.mesh();
    
        // Displacement per patch point
        vectorField patchDisp(localPoints.size(), vector::zero);
    
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
        if (returnReduce(localPoints.size(), sumOp<label>()) > 0)
    
            // Current surface snapped to
            labelList snapSurf(localPoints.size(), -1);
    
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
            // Divide surfaces into zoned and unzoned
    
    mattijs's avatar
    mattijs committed
            labelList zonedSurfaces =
                meshRefiner_.surfaces().getNamedSurfaces();
            labelList unzonedSurfaces =
                meshRefiner_.surfaces().getUnnamedSurfaces();
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
    
    
            // 1. All points to non-interface surfaces
            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                List<pointIndexHit> hitInfo;
                labelList hitSurface;
                surfaces.findNearest
                (
                    unzonedSurfaces,
                    localPoints,
    
                    sqr(snapDist),        // sqr of attract distance
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                    hitSurface,
                    hitInfo
                );
    
                forAll(hitInfo, pointI)
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                    if (hitInfo[pointI].hit())
                    {
                        patchDisp[pointI] =
                            hitInfo[pointI].hitPoint()
                          - localPoints[pointI];
    
    
                        snapSurf[pointI] = hitSurface[pointI];
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                    }
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
            // 2. All points on zones to their respective surface
            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
            // Surfaces with zone information
            const wordList& faceZoneNames = surfaces.faceZoneNames();
    
            // Current best snap distance
    
    mattijs's avatar
    mattijs committed
            scalarField minSnapDist(snapDist);
    
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
            forAll(zonedSurfaces, i)
            {
                label zoneSurfI = zonedSurfaces[i];
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                const labelList surfacesToTest(1, zoneSurfI);
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                // Get indices of points both on faceZone and on pp.
                labelList zonePointIndices
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                    getZoneSurfacePoints
                    (
                        pp,
                        faceZoneNames[zoneSurfI]
                    )
                );
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                // Find nearest for points both on faceZone and pp.
                List<pointIndexHit> hitInfo;
                labelList hitSurface;
                surfaces.findNearest
                (
                    labelList(1, zoneSurfI),
    
    mattijs's avatar
    mattijs committed
                    pointField(localPoints, zonePointIndices),
    
                    sqr(scalarField(minSnapDist, zonePointIndices)),
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                    hitSurface,
                    hitInfo
                );
    
    
    mattijs's avatar
    mattijs committed
                forAll(hitInfo, i)
    
                    label pointI = zonePointIndices[i];
    
    mattijs's avatar
    mattijs committed
    
    
    mattijs's avatar
    mattijs committed
                    if (hitInfo[i].hit())
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                    {
                        patchDisp[pointI] =
    
    mattijs's avatar
    mattijs committed
                            hitInfo[i].hitPoint()
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                          - localPoints[pointI];
    
    mattijs's avatar
    mattijs committed
    
                        minSnapDist[pointI] = min
                        (
                            minSnapDist[pointI],
                            mag(patchDisp[pointI])
                        );
    
    
                        snapSurf[pointI] = zoneSurfI;
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                    }
    
            // Check if all points are being snapped
            forAll(snapSurf, pointI)
            {
                if (snapSurf[pointI] == -1)
                {
                    WarningIn("autoSnapDriver::calcNearestSurface(..)")
                        << "For point:" << pointI
                        << " coordinate:" << localPoints[pointI]
                        << " did not find any surface within:"
                        << minSnapDist[pointI]
                        << " meter." << endl;
                }
            }
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
            {
                scalarField magDisp(mag(patchDisp));
    
    Mattijs Janssens's avatar
    Mattijs Janssens committed
                Info<< "Wanted displacement : average:"
                    << gSum(magDisp)/returnReduce(patchDisp.size(), sumOp<label>())
                    << " min:" << gMin(magDisp)
                    << " max:" << gMax(magDisp) << endl;
            }
    
        }
    
        Info<< "Calculated surface displacement in = "
            << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
    
    
        // Limit amount of movement.
        forAll(patchDisp, patchPointI)
        {
            scalar magDisp = mag(patchDisp[patchPointI]);
    
            if (magDisp > snapDist[patchPointI])
            {
                patchDisp[patchPointI] *= snapDist[patchPointI] / magDisp;
    
                Pout<< "Limiting displacement for " << patchPointI
                    << " from " << magDisp << " to " << snapDist[patchPointI]
                    << endl;
            }
        }
    
        // Points on zones in one domain but only present as point on other
        // will not do condition 2 on all. Sync explicitly.
        syncTools::syncPointList
        (
            mesh,
            pp.meshPoints(),
            patchDisp,
    
            minMagSqrEqOp<point>(),         // combine op
            vector(GREAT, GREAT, GREAT)     // null value (note: cannot use VGREAT)
    
        );
    
        return patchDisp;
    }
    
    
    void Foam::autoSnapDriver::smoothDisplacement
    (
        const snapParameters& snapParams,
        motionSmoother& meshMover
    ) const
    {
        const fvMesh& mesh = meshRefiner_.mesh();
        const indirectPrimitivePatch& pp = meshMover.patch();
    
        Info<< "Smoothing displacement ..." << endl;
    
        // Set edge diffusivity as inverse of distance to patch
    
        scalarField edgeGamma(1.0/(edgePatchDist(meshMover.pMesh(), pp) + SMALL));
    
        //scalarField edgeGamma(mesh.nEdges(), 1.0);
        //scalarField edgeGamma(wallGamma(mesh, pp, 10, 1));
    
        // Get displacement field
        pointVectorField& disp = meshMover.displacement();
    
        for (label iter = 0; iter < snapParams.nSmoothDispl(); iter++)
        {
            if ((iter % 10) == 0)
            {
                Info<< "Iteration " << iter << endl;
            }
            pointVectorField oldDisp(disp);
    
    mattijs's avatar
    mattijs committed
            meshMover.smooth(oldDisp, edgeGamma, disp);
    
        }
        Info<< "Displacement smoothed in = "
            << mesh.time().cpuTimeIncrement() << " s\n" << nl << endl;
    
        if (debug)
        {