Use syncTools::getMasterEdges to get robust parallel summation

src/autoMesh/Make/files

@@ -10,7 +10,6 @@ $(autoHexMesh)/meshRefinement/meshRefinement.C
 $(autoHexMesh)/meshRefinement/meshRefinementMerge.C
 $(autoHexMesh)/meshRefinement/meshRefinementRefine.C
 $(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
-$(autoHexMesh)/offsetTriSurfaceMesh/offsetTriSurfaceMesh.C
 $(autoHexMesh)/trackedParticle/trackedParticle.C
 $(autoHexMesh)/trackedParticle/trackedParticleCloud.C
 

src/autoMesh/autoHexMesh/offsetTriSurfaceMesh/offsetTriSurfaceMesh.C

@@ -1,199 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "offsetTriSurfaceMesh.H"
-#include "Random.H"
-#include "addToRunTimeSelectionTable.H"
-#include "triSurfaceTools.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-defineTypeNameAndDebug(offsetTriSurfaceMesh, 0);
-addToRunTimeSelectionTable(searchableSurface, offsetTriSurfaceMesh, dict);
-
-}
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::offsetTriSurfaceMesh::offsetTriSurfaceMesh
-(
-    const IOobject& io,
-    const triSurface& s,
-    const scalar offset)
-:
-    triSurfaceMesh(io, s),
-    offset_(offset)
-{}
-
-
-Foam::offsetTriSurfaceMesh::offsetTriSurfaceMesh
-(
-    const IOobject& io,
-    const scalar offset
-)
-:
-    triSurfaceMesh(io),
-    offset_(offset)
-{}
-
-
-Foam::offsetTriSurfaceMesh::offsetTriSurfaceMesh
-(
-    const word& name,
-    const objectRegistry& obj,
-    const dictionary& dict
-)
-:
-    triSurfaceMesh(name, obj, dict),
-    offset_(readScalar(dict.lookup("offset")))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::offsetTriSurfaceMesh::~offsetTriSurfaceMesh()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::pointIndexHit Foam::offsetTriSurfaceMesh::findNearest
-(
-    const point& sample,
-    const scalar nearestDistSqr
-) const
-{
-    // Find nearest (add offset to search span)
-    pointIndexHit surfNearest = triSurfaceMesh::findNearest
-    (
-        sample,
-        nearestDistSqr + Foam::sqr(offset_)
-    );
-
-    // Shift back onto surface
-    if (surfNearest.hit())
-    {
-        vector n(sample-surfNearest.hitPoint());
-        n /= mag(n)+VSMALL;
-        surfNearest.setPoint(surfNearest.hitPoint() + offset_*n);
-    }
-    return surfNearest;
-}
-
-
-Foam::pointIndexHit Foam::offsetTriSurfaceMesh::findNearestOnEdge
-(
-    const point& sample,
-    const scalar nearestDistSqr
-) const
-{
-    // Find nearest (add offset to search span)
-    pointIndexHit surfNearest = triSurfaceMesh::findNearestOnEdge
-    (
-        sample,
-        nearestDistSqr + Foam::sqr(offset_)
-    );
-
-    // Shift back onto surface
-    if (surfNearest.hit())
-    {
-        vector n = sample-surfNearest.hitPoint();
-        n /= mag(n)+VSMALL;
-        surfNearest.setPoint(surfNearest.hitPoint() + offset_*n);
-    }
-    return surfNearest;
-}
-
-
-Foam::searchableSurface::volumeType Foam::offsetTriSurfaceMesh::getVolumeType
-(
-    const point& sample
-) const
-{
-    // Find the nearest point on background surface
-    pointIndexHit surfNearest = triSurfaceMesh::findNearest
-    (
-        sample,
-        Foam::sqr(GREAT)
-    );
-
-    if (!surfNearest.hit())
-    {
-        FatalErrorIn("offsetTriSurfaceMesh::getVolumeType(const point&)")
-            << "treeBb:" << tree().bb()
-            << " sample:" << sample
-            << " surfNearest:" << surfNearest
-            << abort(FatalError);
-    }
-
-    // Offset sample to the point.
-    vector n(surfNearest.hitPoint()-sample);
-    n /= mag(n)+VSMALL;
-
-    triSurfaceTools::sideType t = triSurfaceTools::surfaceSide
-    (
-        *this,
-        sample+offset_*n,
-        surfNearest.index(),
-        surfNearest.hitPoint()
-    );
-
-    if (t == triSurfaceTools::UNKNOWN)
-    {
-        return searchableSurface::UNKNOWN;
-    }
-    else if (t == triSurfaceTools::INSIDE)
-    {
-        return searchableSurface::INSIDE;
-    }
-    else if (t == triSurfaceTools::OUTSIDE)
-    {
-        return searchableSurface::OUTSIDE;
-    }
-    else
-    {
-        FatalErrorIn("offsetTriSurfaceMesh::getVolumeType(const point&)")
-            << "problem" << abort(FatalError);
-        return searchableSurface::UNKNOWN;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
-// ************************************************************************* //

src/autoMesh/autoHexMesh/offsetTriSurfaceMesh/offsetTriSurfaceMesh.H

@@ -1,167 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2007 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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    offsetTriSurfaceMesh
-
-Description
-    triSurfaceMesh with offset. Used to define refinement boxes as a region
-    within certain distance to the refinement surface.
-    Note: reloads surface.
-
-SourceFiles
-    offsetTriSurfaceMesh.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef offsetTriSurfaceMesh_H
-#define offsetTriSurfaceMesh_H
-
-#include "triSurfaceMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class offsetTriSurfaceMesh Declaration
-\*---------------------------------------------------------------------------*/
-
-class offsetTriSurfaceMesh
-:
-    public triSurfaceMesh
-{
-    // Private data
-
-        const scalar offset_;
-
-public:
-
-    //- Runtime type information
-    TypeName("offsetTriSurfaceMesh");
-
-
-    // Constructors
-
-        //- Construct from triSurface and offset
-        offsetTriSurfaceMesh(const IOobject&, const triSurface&, const scalar);
-
-        //- Construct read and offset
-        offsetTriSurfaceMesh(const IOobject& io, const scalar);
-
-        //- Construct as searchableSurface
-        offsetTriSurfaceMesh
-        (
-            const word& name,
-            const objectRegistry& obj,
-            const dictionary& dict
-        );
-
-
-    // Destructor
-
-        virtual ~offsetTriSurfaceMesh();
-
-
-    // Member Functions
-
-        // searchableSurface implementation
-
-            //- Calculate nearest point on surface. Returns
-            //  - bool : any point found nearer than nearestDistSqr
-            //  - label: relevant index in surface
-            //  - point: actual nearest point found
-            virtual pointIndexHit findNearest
-            (
-                const point& sample,
-                const scalar nearestDistSqr
-            ) const;
-
-            //- Calculate nearest point on edge. Returns
-            //  - bool : any point found nearer than nearestDistSqr
-            //  - label: relevant index in surface
-            //  - point: actual nearest point found
-            virtual pointIndexHit findNearestOnEdge
-            (
-                const point& sample,
-                const scalar nearestDistSqr
-            ) const;
-
-            //- Find nearest to line. Returns
-            //  - bool : any point found?
-            //  - label: relevant index in shapes
-            //  - point: actual nearest point found
-            //  sets:
-            //  - tightest  : bounding box
-            //  - linePoint : corresponding nearest point on line
-            virtual pointIndexHit findNearest
-            (
-                const linePointRef& ln,
-                treeBoundBox& tightest,
-                point& linePoint
-            ) const
-            {
-                notImplemented("offsetTriSurfaceMesh::findNearest(..)");
-                return pointIndexHit();
-            }
-
-            //- Find nearest intersection of line between start and end.
-            virtual pointIndexHit findLine
-            (
-                const point& start,
-                const point& end
-            ) const
-            {
-                notImplemented("offsetTriSurfaceMesh::findLine(..)");
-                return pointIndexHit();
-            }
-
-            //- Find any intersection of line between start and end.
-            virtual pointIndexHit findLineAny
-            (
-                const point& start,
-                const point& end
-            ) const
-            {
-                notImplemented("offsetTriSurfaceMesh::findLine(..)");
-                return pointIndexHit();
-            }
-
-            //- Determine type (inside/outside/mixed) for point. unknown if
-            //  cannot be determined (e.g. non-manifold surface)
-            virtual volumeType getVolumeType(const point&) const;
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/dynamicMesh/motionSmoother/motionSmoother.C

@@ -1096,6 +1096,9 @@ void Foam::motionSmoother::updateMesh()
         isInternalPoint_.set(meshPoints[i], 0);
     }
 
+    // Calculate master edge addressing
+    isMasterEdge_ = syncTools::getMasterEdges(mesh_);
+
     makePatchPatchAddressing();
 }
 

src/dynamicMesh/motionSmoother/motionSmoother.H

@@ -160,6 +160,10 @@ class motionSmoother
         //- Is mesh point on boundary or not
         PackedList<1> isInternalPoint_;
 
+        //- Is edge master (always except if on coupled boundary and on
+        //  lower processor)
+        PackedList<1> isMasterEdge_;
+
         //- 2-D motion corrector
         twoDPointCorrector twoDCorrector_;
 
@@ -171,31 +175,6 @@ class motionSmoother
 
     // Private Member Functions
 
-        ////- (unweighted) average of all values on points connected via edges
-        ////  to pointI
-        //template <class Type>
-        //static Type avg
-        //(
-        //    const GeometricField<Type, pointPatchField, pointMesh>&,
-        //    const edgeList& edges,
-        //    const pointField& points,
-        //    const labelList& edgeLabels,
-        //    const label pointI
-        //);
-        //
-        ////- Distance weighted average. Is average with inverse distance
-        ////  weighting and varying edge-diffusivity.
-        //template <class Type>
-        //static Type avg
-        //(
-        //    const GeometricField<Type, pointPatchField, pointMesh>&,
-        //    const scalarField& edgeGamma,
-        //    const edgeList& edges,
-        //    const pointField& points,
-        //    const labelList& edgeLabels,
-        //    const label pointI
-        //);
-
         //- Average of connected points.
         template <class Type>
         tmp<GeometricField<Type, pointPatchField, pointMesh> > avg
@@ -481,13 +460,6 @@ public:
 
             // Helper functions to manipulate displacement vector.
 
-                ////- Point-jacobi smoothing of internal points
-                //template <class Type>
-                //void smooth
-                //(
-                //    GeometricField<Type, pointPatchField, pointMesh>&
-                //) const;
-
                 //- Fully explicit smoothing of internal points with varying
                 //  diffusivity.
                 template <class Type>

src/dynamicMesh/motionSmoother/motionSmootherCheck.C

@@ -227,7 +227,7 @@ bool Foam::motionSmoother::checkMesh
         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 
         Info<< "    faces with interpolation weights (0..1)  < "
-            << setw(4) << minWeight
+            << setw(5) << minWeight
             << "       : "
             << nNewWrongFaces-nWrongFaces << endl;
 
@@ -250,7 +250,7 @@ bool Foam::motionSmoother::checkMesh
         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 
         Info<< "    faces with volume ratio of neighbour cells < "
-            << setw(4) << minVolRatio
+            << setw(5) << minVolRatio
             << "     : "
             << nNewWrongFaces-nWrongFaces << endl;
 
@@ -519,7 +519,7 @@ bool Foam::motionSmoother::checkMesh
         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 
         Info<< "    faces with interpolation weights (0..1)  < "
-            << setw(4) << minWeight
+            << setw(5) << minWeight
             << "       : "
             << nNewWrongFaces-nWrongFaces << endl;
 
@@ -540,7 +540,7 @@ bool Foam::motionSmoother::checkMesh
         label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
 
         Info<< "    faces with volume ratio of neighbour cells < "
-            << setw(4) << minVolRatio
+            << setw(5) << minVolRatio
             << "     : "
             << nNewWrongFaces-nWrongFaces << endl;
 

src/dynamicMesh/motionSmoother/motionSmootherTemplates.C

@@ -182,18 +182,25 @@ Foam::tmp<Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh> >
     );
     GeometricField<Type, pointPatchField, pointMesh>& res = tres();
 
-
     const polyMesh& mesh = fld.mesh()();
 
-    scalarField sumWeight(mesh.nPoints(), 0.0);
 
+    // Sum local weighted values and weights
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Note: on coupled edges use only one edge (through isMasterEdge)
+    // This is done so coupled edges do not get counted double.
+
+    scalarField sumWeight(mesh.nPoints(), 0.0);
 
     const edgeList& edges = mesh.edges();
 
     forAll(edges, edgeI)
+    {
+        if (isMasterEdge_.get(edgeI) == 1)
         {
             const edge& e = edges[edgeI];
-        const scalar& w = edgeWeight[edgeI];
+            const scalar w = edgeWeight[edgeI];
 
             res[e[0]] += w*fld[e[1]];
             sumWeight[e[0]] += w;
@@ -201,88 +208,11 @@ Foam::tmp<Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh> >
             res[e[1]] += w*fld[e[0]];
             sumWeight[e[1]] += w;
         }
-
-
-    // Correct for coupled edges counted twice after summing below.
-
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    forAll(patches, patchI)
-    {
-        if (Pstream::parRun() && isA<processorPolyPatch>(patches[patchI]))
-        {
-            const processorPolyPatch& pp =
-                refCast<const processorPolyPatch>(patches[patchI]);
-
-            if (pp.myProcNo() < pp.neighbProcNo())
-            {
-                // Remove contribution of my edges to sum.
-
-                const labelList& meshEdges = pp.meshEdges();
-
-                forAll(meshEdges, i)
-                {
-                    label edgeI = meshEdges[i];
-                    const edge& e = edges[edgeI];
-                    const scalar& w = edgeWeight[edgeI];
-
-                    res[e[0]] -= w*fld[e[1]];
-                    sumWeight[e[0]] -= w;
-
-                    res[e[1]] -= w*fld[e[0]];
-                    sumWeight[e[1]] -= w;
-                }
-            }
-        }
-        else if (isA<cyclicPolyPatch>(patches[patchI]))
-        {
-            // Remove first half contribution.
-
-            const cyclicPolyPatch& pp =
-                refCast<const cyclicPolyPatch>(patches[patchI]);
-
-            SubList<face> half0Faces
-            (
-                mesh.faces(),
-                pp.size()/2,
-                pp.start()
-            );
-
-            labelList::subList half0Cells
-            (
-                mesh.faceOwner(),
-                pp.size()/2,
-                pp.start()
-            );
-
-            labelList meshEdges
-            (
-                primitivePatch(half0Faces, mesh.points()).meshEdges
-                (
-                    edges,
-                    mesh.cellEdges(),
-                    half0Cells
-                )
-            );
-
-            forAll(meshEdges, i)
-            {
-                label edgeI = meshEdges[i];
-                const edge& e = edges[edgeI];
-                const scalar& w = edgeWeight[edgeI];
-
-                res[e[0]] -= w*fld[e[1]];
-                sumWeight[e[0]] -= w;
-
-                res[e[1]] -= w*fld[e[0]];
-                sumWeight[e[1]] -= w;
-            }
-        }
     }
 
-    // Still to be done: multiply shared edges.
+
-    // (mesh.globalData().sharedEdgeLabels()) However needs for us to keep
+    // Add coupled contributions
-    // count how many have already been corrected above.
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
     syncTools::syncPointList
     (
@@ -303,6 +233,9 @@ Foam::tmp<Foam::GeometricField<Type, Foam::pointPatchField, Foam::pointMesh> >
     );
 
 
+    // Average
+    // ~~~~~~~
+
     forAll(res, pointI)
     {
         if (mag(sumWeight[pointI]) < VSMALL)
@@ -330,7 +263,6 @@ void Foam::motionSmoother::smooth
     const GeometricField<Type, pointPatchField, pointMesh>& fld,
     const scalarField& edgeWeight,
     const bool separation,
-
     GeometricField<Type, pointPatchField, pointMesh>& newFld
 ) const
 {
@@ -349,6 +281,7 @@ void Foam::motionSmoother::smooth
             newFld[pointI] = 0.5*fld[pointI] + 0.5*avgFld[pointI];
         }
     }
+
     newFld.correctBoundaryConditions();
     applyCornerConstraints(newFld);
 }