Merge branch 'surface-area' into 'develop'

Provide common face area/normal support in PrimitivePatch

For polyPatch, both the faceAreas() and faceCentres() are masked by their subField equivalents.
Since there is no polyPatch method for magFaceAreas(), the PrimitivePatch method will be seen.

See merge request !74

applications/utilities/surface/surfaceMeshConvertTesting/surfaceMeshConvertTesting.C

@@ -161,7 +161,8 @@ int main(int argc, char *argv[])
 
         Info<< "Read surface:" << endl;
         surf.writeStats(Info);
-        Info<< endl;
+        Info<< "Area         : " << sum(surf.magSf()) << nl
+            << endl;
 
         // check: output to ostream, construct from istream
         {
@@ -205,7 +206,8 @@ int main(int argc, char *argv[])
             Info<< " with scaling " << scaleFactor << endl;
             surf.scalePoints(scaleFactor);
             surf.writeStats(Info);
-            Info<< endl;
+            Info<< "Area         : " << sum(surf.magSf()) << nl
+                << endl;
         }
 
         if (optStdout)
@@ -224,7 +226,8 @@ int main(int argc, char *argv[])
 
         Info<< "Read surface:" << endl;
         surf.writeStats(Info);
-        Info<< endl;
+        Info<< "Area         : " << sum(surf.magSf()) << nl
+            << endl;
 
         // check: output to ostream, construct from istream
         {
@@ -268,7 +271,8 @@ int main(int argc, char *argv[])
             Info<< " with scaling " << scaleFactor << endl;
             surf.scalePoints(scaleFactor);
             surf.writeStats(Info);
-            Info<< endl;
+            Info<< "Area         : " << sum(surf.magSf()) << nl
+                << endl;
         }
 
         if (optStdout)
@@ -286,7 +290,8 @@ int main(int argc, char *argv[])
 
         Info<< "Read surface:" << endl;
         surf.writeStats(Info);
-        Info<< endl;
+        Info<< "Area         : " << sum(surf.magSf()) << nl
+            << endl;
 
         // check: output to ostream, construct from istream
         {
@@ -330,7 +335,8 @@ int main(int argc, char *argv[])
             Info<< " with scaling " << scaleFactor << endl;
             surf.scalePoints(scaleFactor);
             surf.writeStats(Info);
-            Info<< endl;
+            Info<< "Area         : " << sum(surf.magSf()) << nl
+                << endl;
         }
 
         if (optStdout)
@@ -348,7 +354,8 @@ int main(int argc, char *argv[])
 
         Info<< "Read surface:" << endl;
         surf.writeStats(Info);
-        Info<< endl;
+        Info<< "Area         : " << sum(surf.magSf()) << nl
+            << endl;
 
         // check: output to ostream, construct from istream
         {
@@ -392,7 +399,8 @@ int main(int argc, char *argv[])
             Info<< " with scaling " << scaleFactor << endl;
             surf.scalePoints(scaleFactor);
             surf.writeStats(Info);
-            Info<< endl;
+            Info<< "Area         : " << sum(surf.magSf()) << nl
+                << endl;
         }
 
         if (optStdout)

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C

@@ -58,6 +58,8 @@ PrimitivePatch
     localPointsPtr_(nullptr),
     localPointOrderPtr_(nullptr),
     faceCentresPtr_(nullptr),
+    faceAreasPtr_(nullptr),
+    magFaceAreasPtr_(nullptr),
     faceNormalsPtr_(nullptr),
     pointNormalsPtr_(nullptr)
 {}
@@ -94,6 +96,8 @@ PrimitivePatch
     localPointsPtr_(nullptr),
     localPointOrderPtr_(nullptr),
     faceCentresPtr_(nullptr),
+    faceAreasPtr_(nullptr),
+    magFaceAreasPtr_(nullptr),
     faceNormalsPtr_(nullptr),
     pointNormalsPtr_(nullptr)
 {}
@@ -131,6 +135,8 @@ PrimitivePatch
     localPointsPtr_(nullptr),
     localPointOrderPtr_(nullptr),
     faceCentresPtr_(nullptr),
+    faceAreasPtr_(nullptr),
+    magFaceAreasPtr_(nullptr),
     faceNormalsPtr_(nullptr),
     pointNormalsPtr_(nullptr)
 {}
@@ -167,6 +173,8 @@ PrimitivePatch
     localPointsPtr_(nullptr),
     localPointOrderPtr_(nullptr),
     faceCentresPtr_(nullptr),
+    faceAreasPtr_(nullptr),
+    magFaceAreasPtr_(nullptr),
     faceNormalsPtr_(nullptr),
     pointNormalsPtr_(nullptr)
 {}
@@ -524,6 +532,46 @@ faceCentres() const
 }
 
 
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+const Foam::Field<PointType>&
+Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
+faceAreas() const
+{
+    if (!faceAreasPtr_)
+    {
+        calcFaceAreas();
+    }
+
+    return *faceAreasPtr_;
+}
+
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+const Foam::Field<Foam::scalar>&
+Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
+magFaceAreas() const
+{
+    if (!magFaceAreasPtr_)
+    {
+        calcMagFaceAreas();
+    }
+
+    return *magFaceAreasPtr_;
+}
+
+
 template
 <
     class Face,

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H

@@ -168,6 +168,12 @@ private:
         //- Face centres
         mutable Field<PointType>* faceCentresPtr_;
 
+        //- Face area vectors
+        mutable Field<PointType>* faceAreasPtr_;
+
+        //- Mag face area
+        mutable Field<scalar>* magFaceAreasPtr_;
+
         //- Face unit normals
         mutable Field<PointType>* faceNormalsPtr_;
 
@@ -210,6 +216,12 @@ private:
         //- Calculate face centres
         void calcFaceCentres() const;
 
+        //- Calculate face area vectors
+        void calcFaceAreas() const;
+
+        //- Calculate face area magnitudes
+        void calcMagFaceAreas() const;
+
         //- Calculate unit face normals
         void calcFaceNormals() const;
 
@@ -381,7 +393,13 @@ public:
             //- Return face centres for patch
             const Field<PointType>& faceCentres() const;
 
-            //- Return face normals for patch
+            //- Return face area vectors for patch
+            const Field<PointType>& faceAreas() const;
+
+            //- Return face area magnitudes for patch
+            const Field<scalar>& magFaceAreas() const;
+
+            //- Return face unit normals for patch
             const Field<PointType>& faceNormals() const;
 
             //- Return point normals for patch

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C

@@ -47,6 +47,8 @@ clearGeom()
 
     deleteDemandDrivenData(localPointsPtr_);
     deleteDemandDrivenData(faceCentresPtr_);
+    deleteDemandDrivenData(faceAreasPtr_);
+    deleteDemandDrivenData(magFaceAreasPtr_);
     deleteDemandDrivenData(faceNormalsPtr_);
     deleteDemandDrivenData(pointNormalsPtr_);
 }

src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C

@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -52,7 +52,7 @@ calcMeshData() const
     if (meshPointsPtr_ || localFacesPtr_)
     {
         FatalErrorInFunction
-            << "meshPointsPtr_ or localFacesPtr_already allocated"
+            << "meshPointsPtr_ or localFacesPtr_ already allocated"
             << abort(FatalError);
     }
 
@@ -210,7 +210,7 @@ calcLocalPoints() const
     if (localPointsPtr_)
     {
         FatalErrorInFunction
-            << "localPointsPtr_already allocated"
+            << "localPointsPtr_ already allocated"
             << abort(FatalError);
     }
 
@@ -259,7 +259,7 @@ calcPointNormals() const
     if (pointNormalsPtr_)
     {
         FatalErrorInFunction
-            << "pointNormalsPtr_already allocated"
+            << "pointNormalsPtr_ already allocated"
             << abort(FatalError);
     }
 
@@ -323,7 +323,7 @@ calcFaceCentres() const
     if (faceCentresPtr_)
     {
         FatalErrorInFunction
-            << "faceCentresPtr_already allocated"
+            << "faceCentresPtr_ already allocated"
             << abort(FatalError);
     }
 
@@ -346,6 +346,98 @@ calcFaceCentres() const
 }
 
 
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+void
+Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
+calcMagFaceAreas() const
+{
+    if (debug)
+    {
+        Pout<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
+               "calcMagFaceAreas() : "
+               "calculating magFaceAreas in PrimitivePatch"
+            << endl;
+    }
+
+    // It is an error to calculate these more than once.
+    if (magFaceAreasPtr_)
+    {
+        FatalErrorInFunction
+            << "magFaceAreasPtr_ already allocated"
+            << abort(FatalError);
+    }
+
+    magFaceAreasPtr_ = new Field<scalar>(this->size());
+    Field<scalar>& a = *magFaceAreasPtr_;
+
+    forAll(a, facei)
+    {
+        a[facei] = this->operator[](facei).mag(points_);
+    }
+
+    if (debug)
+    {
+        Pout<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
+               "calcMagFaceAreas() : "
+               "finished calculating magFaceAreas in PrimitivePatch"
+            << endl;
+    }
+}
+
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+void
+Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
+calcFaceAreas() const
+{
+    if (debug)
+    {
+        Pout<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
+               "calcFaceAreas() : "
+               "calculating faceAreas in PrimitivePatch"
+            << endl;
+    }
+
+    // It is considered an error to attempt to recalculate faceNormals
+    // if they have already been calculated.
+    if (faceAreasPtr_)
+    {
+        FatalErrorInFunction
+            << "faceAreasPtr_ already allocated"
+            << abort(FatalError);
+    }
+
+    faceAreasPtr_ = new Field<PointType>(this->size());
+
+    Field<PointType>& n = *faceAreasPtr_;
+
+    forAll(n, facei)
+    {
+        n[facei] = this->operator[](facei).normal(points_);
+    }
+
+    if (debug)
+    {
+        Pout<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
+               "calcFaceAreas() : "
+               "finished calculating faceAreas in PrimitivePatch"
+            << endl;
+    }
+}
+
+
 template
 <
     class Face,
@@ -370,7 +462,7 @@ calcFaceNormals() const
     if (faceNormalsPtr_)
     {
         FatalErrorInFunction
-            << "faceNormalsPtr_already allocated"
+            << "faceNormalsPtr_ already allocated"
             << abort(FatalError);
     }
 

src/randomProcesses/noise/noiseModels/surfaceNoise/surfaceNoise.C

@@ -301,14 +301,8 @@ Foam::scalar surfaceNoise::writeSurfaceData
                 false
             );
 
-            // TODO: Move faceAreas to demand-driven function in MeshedSurface
-            // scalarField faceAreas(surf.faces().size());
-            // forAll(faceAreas, i)
-            // {
-            //     faceAreas[i] = surf.faces()[i].mag(surf.points());
-            // }
-            //
-            // areaAverage = sum(allData*faceAreas)/sum(faceAreas);
+            // TO BE VERIFIED: area-averaged values
+            // areaAverage = sum(allData*surf.magSf())/sum(surf.magSf());
             areaAverage = sum(allData)/allData.size();
         }
         Pstream::scatter(areaAverage);
@@ -330,14 +324,8 @@ Foam::scalar surfaceNoise::writeSurfaceData
             false
         );
 
-        // TODO: Move faceAreas to demand-driven function in MeshedSurface
-        // scalarField faceAreas(surf.faces().size());
-        // forAll(faceAreas, i)
-        // {
-        //     faceAreas[i] = surf.faces()[i].mag(surf.points());
-        // }
-        //
-        // return sum(data*faceAreas)/sum(faceAreas);
+        // TO BE VERIFIED: area-averaged values
+        // return sum(data*surf.magSf())/sum(surf.magSf());
         return sum(data)/data.size();
     }
 }
@@ -387,14 +375,8 @@ Foam::scalar surfaceNoise::surfaceAverage
                 }
             }
 
-            // TODO: Move faceAreas to demand-driven function in MeshedSurface
-            scalarField faceAreas(surf.faces().size());
-            forAll(faceAreas, i)
-            {
-                faceAreas[i] = surf.faces()[i].mag(surf.points());
-            }
-
-//            areaAverage = sum(allData*faceAreas)/sum(faceAreas);
+            // TO BE VERIFIED: area-averaged values
+            // areaAverage = sum(allData*surf.magSf())/sum(surf.magSf());
             areaAverage = sum(allData)/allData.size();
         }
         Pstream::scatter(areaAverage);
@@ -405,14 +387,8 @@ Foam::scalar surfaceNoise::surfaceAverage
     {
         const meshedSurface& surf = readerPtr_->geometry();
 
-        // TODO: Move faceAreas to demand-driven function in MeshedSurface
-        scalarField faceAreas(surf.faces().size());
-        forAll(faceAreas, i)
-        {
-            faceAreas[i] = surf.faces()[i].mag(surf.points());
-        }
-
-//        return sum(data*faceAreas)/sum(faceAreas);
+        // TO BE VERIFIED: area-averaged values
+        // return sum(data*surf.magSf())/sum(surf.magSf());
         return sum(data)/data.size();
     }
 }
@@ -550,9 +526,6 @@ void surfaceNoise::calculate()
             surfPSD13f[freqI][faceI] = PSD13f.y()[freqI];
             surfPrms13f2[freqI][faceI] = Prms13f2.y()[freqI];
         }
-
-        // Free the storage for p
-//        p.clear();
     }
 
     // Output directory for graphs

src/sampling/sampledSurface/distanceSurface/distanceSurface.H

@@ -193,6 +193,25 @@ public:
             return facesPtr_;
         }
 
+        //- Face area vectors
+        virtual const vectorField& Sf() const
+        {
+            return surface().Sf();
+        }
+
+        //- Face area magnitudes
+        virtual const scalarField& magSf() const
+        {
+            return surface().magSf();
+        }
+
+        //- Face centres
+        virtual const vectorField& Cf() const
+        {
+            return surface().Cf();
+        }
+
+
         const triSurface& surface() const
         {
             if (cell_)

src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H

@@ -158,6 +158,11 @@ public:
 
     // Member Functions
 
+        const isoSurface& surface() const
+        {
+            return surfPtr_();
+        }
+
         //- Does the surface need an update?
         virtual bool needsUpdate() const;
 
@@ -194,16 +199,32 @@ public:
             return facesPtr_;
         }
 
+        //- Face area magnitudes
+        virtual const vectorField& Sf() const
+        {
+            return surface().Sf();
+        }
 
-        const isoSurface& surface() const
+        //- Face area magnitudes
+        virtual const scalarField& magSf() const
         {
-            return surfPtr_();
+            return surface().magSf();
         }
 
-        //- Lookup or read isoField. Sets volFieldPtr_ and pointFieldPtr_.
+        //- Face centres
+        virtual const vectorField& Cf() const
+        {
+            return surface().Cf();
+        }
+
+
+        //- Lookup or read isoField.
+        //  Sets volFieldPtr_ and pointFieldPtr_.
         void getIsoField();
 
 
+    // Sample
+
         //- Sample field on surface
         virtual tmp<scalarField> sample
         (
@@ -235,6 +256,8 @@ public:
         ) const;
 
 
+    // Interpolate
+
         //- Interpolate field on surface
         virtual tmp<scalarField> interpolate
         (
@@ -265,6 +288,9 @@ public:
             const interpolation<tensor>&
         ) const;
 
+
+    // Output
+
         //- Write
         virtual void print(Ostream&) const;
 };

src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H

@@ -54,6 +54,9 @@ class sampledIsoSurfaceCell
     public sampledSurface,
     public triSurface
 {
+    //- Private typedef for convenience
+    typedef triSurface MeshStorage;
+
     // Private data
 
         //- Field to get isoSurface of
@@ -145,7 +148,7 @@ public:
         //- Points of surface
         virtual const pointField& points() const
         {
-            return triSurface::points();
+            return MeshStorage::points();
         }
 
         //- Faces of surface
@@ -165,6 +168,24 @@ public:
             return facesPtr_;
         }
 
+        //- Face area magnitudes
+        virtual const vectorField& Sf() const
+        {
+            return MeshStorage::Sf();
+        }
+
+        //- Face area magnitudes
+        virtual const scalarField& magSf() const
+        {
+            return MeshStorage::magSf();
+        }
+
+        //- Face centres
+        virtual const vectorField& Cf() const
+        {
+            return MeshStorage::Cf();
+        }
+
 
         //- Sample field on surface
         virtual tmp<scalarField> sample

src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H

@@ -141,6 +141,12 @@ public:
 
     // Member Functions
 
+        //const isoSurfaceCell& surface() const
+        const isoSurface& surface() const
+        {
+            return isoSurfPtr_();
+        }
+
         //- Does the surface need an update?
         virtual bool needsUpdate() const;
 
@@ -176,13 +182,27 @@ public:
             return facesPtr_;
         }
 
+        //- Face area magnitudes
+        virtual const vectorField& Sf() const
+        {
+            return surface().Sf();
+        }
 
-        //const isoSurfaceCell& surface() const
-        const isoSurface& surface() const
+        //- Face area magnitudes
+        virtual const scalarField& magSf() const