View Issue Details

IDProjectCategoryView StatusLast Update
0003176OpenFOAMBugpublic2019-02-22 16:34
Reporterphilipc Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityalways
Status closedResolutionno change required 
PlatformLinuxOSUbuntuOS Version16.04
Summary0003176: PrimitivePatch::pointNormals() returns inconsistent point normals in parallel
DescriptionExamining the source code for PrimitivePatch::calcPointNormals() (PrimitivePatchMeshData.C:189): for points on the edge/boundary of a patch, the point normal does not take the neighbour-processor-face-normal data into account (it is only a function of the face-normals on the current processor). This means (for a non-flat patch) that a point may/will have different pointNormals on different processors.

I am unaware of classes/solvers/utilities directly affected by this (though there may be many) but I found this issue when creating a mesh surface smoothing procedure.
Additional InformationI have attached a workaround function called pointNormalsConsistent.

One possible solution going forward is to update the pointNormals function or alternatively to create a new function in the PrimitivePatch class called consistentPointNormals.
TagsNo tags attached.

Activities

philipc

2019-02-22 15:41

reporter  

pointNormalsConsistent.H (2,031 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     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
    pointNormalsConsistent

Description
    The default PrimitvePatch pointNormals function does not consistently
    calculate pointNormals at parallel boundaries i.e. points have different
    normals on each patch.

    This function is calculates the patch pointNormals consistently in parallel.

Author
    Philip Cardiff, UCD. All rights reserved.

\*---------------------------------------------------------------------------*/

#ifndef pointNormalsConsistent_H
#define pointNormalsConsistent_H

#include "polyPatch.H"

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

namespace Foam
{

    tmp<pointField> pointNormalsConsistent(const polyPatch& ppatch);

} // End namespace Foam

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

#endif


// ************************************************************************* //
pointNormalsConsistent.H (2,031 bytes)   
pointNormalsConsistent.C (4,957 bytes)   
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright held by original author
     \\/     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 "pointNormalsConsistent.H"
#include "fvMesh.H"
#include "pointMesh.H"
#include "pointFields.H"

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


Foam::tmp<Foam::pointField> Foam::pointNormalsConsistent
(
    const polyPatch& ppatch
)
{
    tmp<pointField> tresult(new pointField(ppatch.nPoints(), vector::zero));

    // Calculate weight for each point on the patch
    // Initialise point weights

    const labelListList& pf = ppatch.pointFaces();
    const labelList& meshPoints = ppatch.meshPoints();

    scalarListList weights(ppatch.nPoints());
    forAll(weights, pointI)
    {
        weights[pointI].setSize(pf[pointI].size());
    }

    pointScalarField sumWeights
    (
        IOobject
        (
            "sumWeights_" + ppatch.name(),
            ppatch.boundaryMesh().mesh().polyMesh::instance(),
            ppatch.boundaryMesh().mesh()
        ),
        pointMesh::New(ppatch.boundaryMesh().mesh()),
        dimensionedScalar("zero", dimless, 0.0)
    );

    // Set un-normalised weights to 1.0 i.e. arithmetic average
    forAll(weights, pI)
    {
        scalarList& pw = weights[pI];
        const labelList& curPointFaces = pf[pI];

        forAll(curPointFaces, pfI)
        {
            // Set the un-normalised point weight
            pw[pfI] = 1.0;

            // Increment the weights sum
            sumWeights[meshPoints[pI]] += 1.0;
        }
    }

    // Sync weights in parallel
    forAll(sumWeights.boundaryField(), patchI)
    {
        if (sumWeights.boundaryField()[patchI].coupled())
        {
            sumWeights.boundaryField()[patchI].initAddField();
        }
    }

    forAll(sumWeights.boundaryField(), patchI)
    {
        if (sumWeights.boundaryField()[patchI].coupled())
        {
            sumWeights.boundaryField()[patchI].addField
            (
                sumWeights.internalField()
            );
        }
    }

    // Normalise the weights
    forAll(weights, pI)
    {
        scalarList& pw = weights[pI];

        forAll(pw, pfI)
        {
            pw[pfI] /= sumWeights[meshPoints[pI]];
        }
    }

    // Calculate consistent point normals

    const pointField& fn = ppatch.faceNormals();

    pointVectorField pointNormalField
    (
        IOobject
        (
            "pointNormalField_" + ppatch.name(),
            ppatch.boundaryMesh().mesh().polyMesh::instance(),
            ppatch.boundaryMesh().mesh()
        ),
        pointMesh::New(ppatch.boundaryMesh().mesh()),
        dimensionedVector("zero", dimless, vector::zero)
    );

    forAll(weights, pointI)
    {
        point& curNormal = pointNormalField[meshPoints[pointI]];
        const labelList& curFaces = pf[pointI];

        forAll(curFaces, faceI)
        {
            curNormal += weights[pointI][faceI]*fn[curFaces[faceI]];
        }
    }

    // Add contributions from other processors
    forAll(pointNormalField.boundaryField(), patchI)
    {
        if (pointNormalField.boundaryField()[patchI].coupled())
        {
            pointNormalField.boundaryField()[patchI].initAddField();
        }
    }

    forAll(pointNormalField.boundaryField(), patchI)
    {
        if (pointNormalField.boundaryField()[patchI].coupled())
        {
            pointNormalField.boundaryField()[patchI].addField
            (
                pointNormalField.internalField()
            );
        }
    }

    // Update coupled and constrained boundaries
    pointNormalField.correctBoundaryConditions();

    // Copy result into tmp field
    tresult() =
        pointNormalField.boundaryField()[ppatch.index()].patchInternalField();

    return tresult;
}


// ************************************************************************* //
pointNormalsConsistent.C (4,957 bytes)   

philipc

2019-02-22 15:43

reporter   ~0010322

I forgot to mention: the attached pointNormalsConsistent.H/C files compile with foam-extend-4.0 (as opposed to OpenFOAM-6); however, the issue is present across all main versions of OpenFOAM.

MattijsJ

2019-02-22 16:12

reporter   ~0010323

You can use PatchToools::pointNormals instead. It 'knows' about polyMesh and hence processor connectivity. See e.g. $FOAM_APP/Test/PatchTools

Issue History

Date Modified Username Field Change
2019-02-22 15:41 philipc New Issue
2019-02-22 15:41 philipc File Added: pointNormalsConsistent.H
2019-02-22 15:41 philipc File Added: pointNormalsConsistent.C
2019-02-22 15:43 philipc Note Added: 0010322
2019-02-22 16:12 MattijsJ Note Added: 0010323
2019-02-22 16:34 henry Assigned To => henry
2019-02-22 16:34 henry Status new => closed
2019-02-22 16:34 henry Resolution open => no change required