Updated sampling to reflect Mark's current revision and reorganized file locations

src/sampling/sampledSurface/patch/sampledPatch.C

@@ -37,7 +37,6 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(sampledPatch, 0);
-    // addToRunTimeSelectionTable(sampledSurface, sampledPatch, word);
     addNamedToRunTimeSelectionTable(sampledSurface, sampledPatch, word, patch);
 }
 
@@ -140,6 +139,8 @@ Foam::sampledPatch::sampledPatch
     faces_(0),
     patchFaceLabels_(0)
 {
+    // default: non-triangulated
+    triangulate() = getBool(dict, "triangulate", false);
     createGeometry();
 }
 

src/sampling/sampledSurface/patch/sampledPatch.H

@@ -26,6 +26,7 @@ Class
     Foam::sampledPatch
 
 Description
+    A sampledSurface on a patch. Non-triangulated by default.
 
 SourceFiles
     sampledPatch.C
@@ -61,7 +62,7 @@ class sampledPatch
         //- Zero size or copy of patch.localPoints()
         pointField points_;
 
-        //- Triangulated faces
+        //- Faces (triangulated or non-triangulated)
         faceList faces_;
 
         //- Local patch face labels
@@ -99,7 +100,7 @@ public:
             const word& name,
             const polyMesh& mesh,
             const word& patchName,
-            const bool triangulate = true
+            const bool triangulate = false
         );
 
         //- Construct from dictionary

src/sampling/sampledSurface/plane/sampledPlane.H

@@ -26,6 +26,7 @@ Class
     Foam::sampledPlane
 
 Description
+    A sampledSurface defined by a cuttingPlane. Triangulated by default.
 
 SourceFiles
     sampledPlane.C
@@ -57,7 +58,8 @@ class sampledPlane
         //- zone name (if restricted to zones)
         word zoneName_;
 
-        //- Triangulated faces in terms of intersection points
+        //- Triangulated faces in terms of intersection points.
+        //  Non-triangulated faces are obtained directly from cuttingPlane
         faceList faces_;
 
         //- For every triangulated face, the original cell in mesh

src/sampling/sampledSurface/sampledSurface/sampledSurfaceTemplates.C

@@ -22,25 +22,19 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "sampledSurface.H"
 
 template<class Type>
-void Foam::sampledSurface::checkFieldSize(const Field<Type>& field) const
+bool Foam::sampledSurface::checkFieldSize(const Field<Type>& field) const
 {
-    if (!faces().size())
+    if (faces().size() == 0 || field.size() == 0)
     {
-        FatalErrorIn
+        return false;
-        (
-            "sampledSurface::checkFieldSize(const Field<Type>&) const"
-        )
-            << "empty surface"
-            << exit(FatalError);
     }
-    else if (field.size() != faces().size())
+
+    if (field.size() != faces().size())
     {
         FatalErrorIn
         (
@@ -51,34 +45,65 @@ void Foam::sampledSurface::checkFieldSize(const Field<Type>& field) const
             << ") != surface (" << faces().size() << ")"
             << exit(FatalError);
     }
+
+    return true;
 }
 
+
 template<class Type>
 Type Foam::sampledSurface::integrate(const Field<Type>& field) const
 {
-    checkFieldSize(field);
+    Type value = pTraits<Type>::zero;
-    return sum(field * magSf());
+
+    if (checkFieldSize(field))
+    {
+        value = sum(field * magSf());
+    }
+
+    reduce(value, sumOp<Type>());
+    return value;
 }
 
+
 template<class Type>
 Type Foam::sampledSurface::integrate(const tmp<Field<Type> >& field) const
 {
-    checkFieldSize(field());
+    Type value = integrate(field());
-    return sum(field * magSf());
+    field.clear();
+    return value;
 }
 
+
 template<class Type>
 Type Foam::sampledSurface::average(const Field<Type>& field) const
 {
-    checkFieldSize(field);
+    Type value = pTraits<Type>::zero;
-    return sum(field * magSf()) / sum(magSf());
+
+    if (checkFieldSize(field))
+    {
+        value = sum(field * magSf());
+    }
+
+    reduce(value, sumOp<Type>());
+
+    // avoid divide-by-zero
+    if (area())
+    {
+        return value / area();
+    }
+    else
+    {
+        return pTraits<Type>::zero;
+    }
 }
 
+
 template<class Type>
 Type Foam::sampledSurface::average(const tmp<Field<Type> >& field) const
 {
-    checkFieldSize(field());
+    Type value = average(field());
-    return sum(field * magSf()) / sum(magSf());
+    field.clear();
+    return value;
 }
 
 
@@ -89,17 +114,22 @@ void Foam::sampledSurface::project
     const Field<Type>& field
 ) const
 {
-    checkFieldSize(field);
+    if (checkFieldSize(field))
-
-    const faceList& sampleFaces = faces();
-    const vectorField& norm = Sf();
-
-    forAll(sampleFaces, faceI)
     {
-        res[faceI] = field[faceI] & (norm[faceI] / mag(norm[faceI]));
+        const vectorField& norm = Sf();
+
+        forAll(norm, faceI)
+        {
+            res[faceI] = field[faceI] & (norm[faceI] / mag(norm[faceI]));
+        }
+    }
+    else
+    {
+        res.clear();
     }
 }
 
+
 template<class ReturnType, class Type>
 void Foam::sampledSurface::project
 (

src/sampling/sampledSurface/surface/sampledSurface.C

@@ -0,0 +1,335 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "sampledSurface.H"
+#include "polyMesh.H"
+#include "demandDrivenData.H"
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sampledSurface, 0);
+    defineRunTimeSelectionTable(sampledSurface, word);
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::sampledSurface::clearGeom() const
+{
+    deleteDemandDrivenData(SfPtr_);
+    deleteDemandDrivenData(magSfPtr_);
+    deleteDemandDrivenData(CfPtr_);
+    area_ = -1;
+}
+
+
+void Foam::sampledSurface::makeSf() const
+{
+    // It is an error to recalculate if the pointer is already set
+    if (SfPtr_)
+    {
+        FatalErrorIn("Foam::sampledSurface::makeSf()")
+            << "face area vectors already exist"
+            << abort(FatalError);
+    }
+
+    const faceList& theFaces = faces();
+    SfPtr_ = new vectorField(theFaces.size());
+
+    vectorField& values = *SfPtr_;
+    forAll(theFaces, faceI)
+    {
+        values[faceI] = theFaces[faceI].normal(points());
+    }
+}
+
+
+void Foam::sampledSurface::makeMagSf() const
+{
+    // It is an error to recalculate if the pointer is already set
+    if (magSfPtr_)
+    {
+        FatalErrorIn("Foam::sampledSurface::makeMagSf()")
+            << "mag face areas already exist"
+            << abort(FatalError);
+    }
+
+    const faceList& theFaces = faces();
+    magSfPtr_ = new scalarField(theFaces.size());
+
+    scalarField& values = *magSfPtr_;
+    forAll(theFaces, faceI)
+    {
+        values[faceI] = theFaces[faceI].mag(points());
+    }
+}
+
+
+void Foam::sampledSurface::makeCf() const
+{
+    // It is an error to recalculate if the pointer is already set
+    if (CfPtr_)
+    {
+        FatalErrorIn("Foam::sampledSurface::makeCf()")
+            << "face centres already exist"
+            << abort(FatalError);
+    }
+
+    const faceList& theFaces = faces();
+    CfPtr_ = new vectorField(theFaces.size());
+
+    vectorField& values = *CfPtr_;
+    forAll(theFaces, faceI)
+    {
+        values[faceI] = theFaces[faceI].centre(points());
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+
+Foam::autoPtr<Foam::sampledSurface>
+Foam::sampledSurface::New
+(
+    const word& sampleType,
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+{
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_
+            ->find(sampleType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "sampledSurface::New(const word&, "
+            "const polyMesh&, const dictionary&)"
+        )   << "Unknown sample type " << sampleType
+            << endl << endl
+            << "Valid sample types : " << endl
+            << wordConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<sampledSurface>
+    (
+        cstrIter()
+        (
+            mesh,
+            dict
+        )
+    );
+}
+
+
+bool Foam::sampledSurface::getBool
+(
+    const dictionary& dict,
+    const word& key,
+    const bool defaultVal
+)
+{
+    if (dict.found(key))
+    {
+        return readBool(dict.lookup(key));
+    }
+    else
+    {
+        return defaultVal;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from mesh, name
+Foam::sampledSurface::sampledSurface
+(
+    const polyMesh& mesh,
+    const word& name,
+    const bool triangulate
+)
+:
+    mesh_(mesh),
+    name_(name),
+    triangulate_(triangulate),
+    interpolate_(false),
+    SfPtr_(NULL),
+    magSfPtr_(NULL),
+    CfPtr_(NULL),
+    area_(-1)
+{}
+
+
+// Construct from dictionary
+Foam::sampledSurface::sampledSurface
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    mesh_(mesh),
+    name_(type()),
+    triangulate_(getBool(dict, "triangulate", true)),
+    interpolate_(getBool(dict, "interpolate", false)),
+    SfPtr_(NULL),
+    magSfPtr_(NULL),
+    CfPtr_(NULL),
+    area_(-1)
+{
+    if (dict.found("name"))
+    {
+        dict.lookup("name") >> name_;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledSurface::~sampledSurface()
+{
+    clearGeom();
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::vectorField& Foam::sampledSurface::Sf() const
+{
+    if (!SfPtr_)
+    {
+        makeSf();
+    }
+
+    return *SfPtr_;
+}
+
+
+const Foam::scalarField& Foam::sampledSurface::magSf() const
+{
+    if (!magSfPtr_)
+    {
+        makeMagSf();
+    }
+
+    return *magSfPtr_;
+}
+
+
+const Foam::vectorField& Foam::sampledSurface::Cf() const
+{
+    if (!CfPtr_)
+    {
+        makeCf();
+    }
+
+    return *CfPtr_;
+}
+
+
+Foam::scalar Foam::sampledSurface::area() const
+{
+    if (area_ < 0)
+    {
+        area_ = sum(magSf());
+        reduce(area_, sumOp<scalar>());
+    }
+
+    return area_;
+}
+
+
+// do not project scalar - just copy values
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::sampledSurface::project(const Field<scalar>& field) const
+{
+    tmp<Field<scalar> > tRes(new Field<scalar>(faces().size()));
+    Field<scalar>& res = tRes();
+
+    forAll(faces(), faceI)
+    {
+        res[faceI] = field[faceI];
+    }
+
+    return tRes;
+}
+
+
+Foam::tmp<Foam::Field<Foam::scalar> >
+Foam::sampledSurface::project(const Field<vector>& field) const
+{
+    tmp<Field<scalar> > tRes(new Field<scalar>(faces().size()));
+    project(tRes(), field);
+    return tRes;
+}
+
+
+Foam::tmp<Foam::Field<Foam::vector> >
+Foam::sampledSurface::project(const Field<sphericalTensor>& field) const
+{
+    tmp<Field<vector> > tRes(new Field<vector>(faces().size()));
+    project(tRes(), field);
+    return tRes;
+}
+
+
+Foam::tmp<Foam::Field<Foam::vector> >
+Foam::sampledSurface::project(const Field<symmTensor>& field) const
+{
+    tmp<Field<vector> > tRes(new Field<vector>(faces().size()));
+    project(tRes(), field);
+    return tRes;
+}
+
+
+Foam::tmp<Foam::Field<Foam::vector> >
+Foam::sampledSurface::project(const Field<tensor>& field) const
+{
+    tmp<Field<vector> > tRes(new Field<vector>(faces().size()));
+    project(tRes(), field);
+    return tRes;
+}
+
+
+void Foam::sampledSurface::print(Ostream& os) const
+{
+    os << type();
+}
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream &os, const sampledSurface& s)
+{
+    s.print(os);
+    os.check("Ostream& operator<<(Ostream&, const sampledSurface&");
+    return os;
+}
+
+// ************************************************************************* //

src/sampling/sampledSurface/surface/sampledSurface.H

@@ -0,0 +1,414 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    Foam::sampledSurface
+
+Description
+    An abstract class for surfaces with sampling.
+
+SourceFiles
+    sampledSurface.C
+    sampledSurfaceTemplates.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledSurface_H
+#define sampledSurface_H
+
+#include "pointField.H"
+#include "word.H"
+#include "labelList.H"
+#include "faceList.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "volFieldsFwd.H"
+#include "polyMesh.H"
+#include "coordinateSystems.H"
+#include "interpolation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class sampledSurface Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledSurface
+{
+    // Private data
+
+        //- Reference to mesh
+        const polyMesh& mesh_;
+
+        //- Name of sample surface
+        word name_;
+
+        //- Triangulate faces or keep faces as it
+        bool triangulate_;
+
+        //- Do we intend to interpolate the information?
+        bool interpolate_;
+
+    // Demand-driven data
+
+        //- Face area vectors
+        mutable vectorField* SfPtr_;
+
+        //- Mag face area vectors
+        mutable scalarField* magSfPtr_;
+
+        //- Face centres
+        mutable vectorField* CfPtr_;
+
+        //- Total surface area
+        mutable scalar area_;
+
+    // Make geometric data
+
+        //- make Sf
+        void makeSf() const;
+
+        //- make magSf
+        void makeMagSf() const;
+
+        //- make Cf
+        void makeCf() const;
+
+    // Service methods
+
+        //- Check field size matches surface size
+        template<class Type>
+        bool checkFieldSize(const Field<Type>&) const;
+
+        //- project field onto surface
+        template<class ReturnType, class Type>
+        void project
+        (
+            Field<ReturnType>&,
+            const Field<Type>&
+        ) const;
+
+        //- project field onto surface
+        template<class ReturnType, class Type>
+        void project
+        (
+            Field<ReturnType>&,
+            const tmp<Field<Type> >&
+        ) const;
+
+        //- project field onto surface
+        template<class ReturnType, class Type>
+        tmp<Field<ReturnType> > project(const tmp<Field<Type> >&) const;
+
+protected:
+
+    // Protected static functions
+
+        //- Read bool from dictionary. Return provided value if not found
+        static bool getBool(const dictionary&, const word&, const bool);
+
+    // Protected Member functions
+
+        virtual void clearGeom() const;
+
+        //- Non-const access to triangulation
+        bool& triangulate()
+        {
+            return triangulate_;
+        }
+
+public:
+
+    //- Runtime type information
+      TypeName("sampledSurface");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            sampledSurface,
+            word,
+            (
+                const polyMesh& mesh,
+                const dictionary& dict
+            ),
+            (mesh, dict)
+        );
+
+
+        //- Class used for the PtrLists read-construction
+        //  Has the ability to rewrite coordinate systems as required
+        class iNew
+        :
+            public coordinateSystems
+        {
+            //- Reference to the volume mesh
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                coordinateSystems(mesh),
+                mesh_(mesh)
+            {}
+
+            iNew
+            (
+                const polyMesh& mesh,
+                const coordinateSystems& cs
+            )
+            :
+                coordinateSystems(cs),
+                mesh_(mesh)
+            {}
+
+            autoPtr<sampledSurface> operator()(Istream& is) const
+            {
+                word sampleType(is);
+                dictionary dict(is);
+                rewriteDict(dict, true);
+
+                return sampledSurface::New(sampleType, mesh_, dict);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from mesh, name
+        sampledSurface
+        (
+            const polyMesh& mesh,
+            const word& name,
+            const bool triangulate = true
+        );
+
+        //- Construct from dictionary
+        sampledSurface
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Clone
+        autoPtr<sampledSurface> clone() const
+        {
+            notImplemented("autoPtr<sampledSurface> clone() const");
+            return autoPtr<sampledSurface>(NULL);
+        }
+
+
+    // Selectors
+
+        //- Return a reference to the selected surface
+        static autoPtr<sampledSurface> New
+        (
+            const word& sampleType,
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~sampledSurface();
+
+
+    // Member Functions
+
+      // Access
+
+        //- Access to the underlying mesh
+        const polyMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        //- Name of surface
+        const word& name() const
+        {
+            return name_;
+        }
+
+        //- interpolation requested for surface
+        bool interpolate() const
+        {
+            return interpolate_;
+        }
+
+        //- triangulation requested for surface
+        bool triangulate() const
+        {
+            return triangulate_;
+        }
+
+        //- Points of surface
+        virtual const pointField& points() const = 0;
+
+        //- Faces of surface
+        virtual const faceList& faces() const = 0;
+
+        //- Correct for mesh movement and/or field changes
+        virtual void correct(const bool meshChanged) = 0;
+
+        //- Return face area vectors
+        virtual const vectorField& Sf() const;
+
+        //- Return face area magnitudes
+        virtual const scalarField& magSf() const;
+
+        //- Return face centres as vectorField
+        virtual const vectorField& Cf() const;
+
+        //- The total surface area
+        scalar area() const;
+
+        //- integration of a field across the surface
+        template<class Type>
+        Type integrate(const Field<Type>&) const;
+
+        //- integration of a field across the surface
+        template<class Type>
+        Type integrate(const tmp<Field<Type> >&) const;
+
+        //- area-averaged value of a field across the surface
+        template<class Type>
+        Type average(const Field<Type>&) const;
+
+        //- area-averaged value of a field across the surface
+        template<class Type>
+        Type average(const tmp<Field<Type> >&) const;
+
+        //- project field onto surface
+        tmp<Field<scalar> > project(const Field<scalar>&) const;
+        tmp<Field<scalar> > project(const Field<vector>&) const;
+        tmp<Field<vector> > project(const Field<sphericalTensor>&) const;
+        tmp<Field<vector> > project(const Field<symmTensor>&) const;
+        tmp<Field<vector> > project(const Field<tensor>&) const;
+
+        //- sample field on surface
+        virtual tmp<scalarField> sample
+        (
+            const volScalarField&
+        ) const = 0;
+
+        //- sample field on surface
+        virtual tmp<vectorField> sample
+        (
+            const volVectorField&
+        ) const = 0;
+
+        //- sample field on surface
+        virtual tmp<sphericalTensorField> sample
+        (
+            const volSphericalTensorField&
+        ) const = 0;
+
+        //- sample field on surface
+        virtual tmp<symmTensorField> sample
+        (
+            const volSymmTensorField&
+        ) const = 0;
+
+        //- sample field on surface
+        virtual tmp<tensorField> sample
+        (
+            const volTensorField&
+        ) const = 0;
+
+
+        //- interpolate field on surface
+        virtual tmp<scalarField> interpolate
+        (
+            const interpolation<scalar>&
+        ) const = 0;
+
+
+        //- interpolate field on surface
+        virtual tmp<vectorField> interpolate
+        (
+            const interpolation<vector>&
+        ) const = 0;
+
+        //- interpolate field on surface
+        virtual tmp<sphericalTensorField> interpolate
+        (
+            const interpolation<sphericalTensor>&
+        ) const = 0;
+
+        //- interpolate field on surface
+        virtual tmp<symmTensorField> interpolate
+        (
+            const interpolation<symmTensor>&
+        ) const = 0;
+
+        //- interpolate field on surface
+        virtual tmp<tensorField> interpolate
+        (
+            const interpolation<tensor>&
+        ) const = 0;
+
+
+      // Edit
+
+        //- Rename
+        virtual void rename(const word& newName)
+        {
+            name_ = newName;
+        }
+
+      // Write
+
+        //- Write
+        virtual void print(Ostream&) const;
+
+       // IOstream operators
+
+        friend Ostream& operator<<(Ostream&, const sampledSurface&);
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "sampledSurfaceTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSurface/surfaces/sampledSurfaces.C

@@ -0,0 +1,390 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "sampledSurfaces.H"
+#include "volFields.H"
+#include "dictionary.H"
+#include "Time.H"
+#include "IOmanip.H"
+#include "ListListOps.H"
+#include "mergePoints.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    //- Used to offset faces in Pstream::combineOffset
+    template <>
+    class offsetOp<face>
+    {
+
+    public:
+
+        face operator()
+        (
+            const face& x,
+            const label offset
+        ) const
+        {
+            face result(x.size());
+
+            forAll(x, xI)
+            {
+                result[xI] = x[xI] + offset;
+            }
+            return result;
+        }
+    };
+
+    defineTypeNameAndDebug(sampledSurfaces, 0);
+}
+
+bool Foam::sampledSurfaces::verbose_ = false;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+bool Foam::sampledSurfaces::checkFieldTypes()
+{
+    wordList fieldTypes(fieldNames_.size());
+
+    // check files for a particular time
+    if (loadFromFiles_)
+    {
+        forAll(fieldNames_, fieldI)
+        {
+            IOobject io
+            (
+                fieldNames_[fieldI],
+                obr_.time().timeName(),
+                refCast<const polyMesh>(obr_),
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            );
+
+            if (io.headerOk())
+            {
+                fieldTypes[fieldI] = io.headerClassName();
+            }
+            else
+            {
+                fieldTypes[fieldI] = "(notFound)";
+            }
+        }
+    }
+    else
+    {
+        // check objectRegistry
+        forAll(fieldNames_, fieldI)
+        {
+            objectRegistry::const_iterator iter =
+                obr_.find(fieldNames_[fieldI]);
+
+            if (iter != obr_.end())
+            {
+                fieldTypes[fieldI] = iter()->type();
+            }
+            else
+            {
+                fieldTypes[fieldI] = "(notFound)";
+            }
+        }
+    }
+
+
+    label nFields = 0;
+
+    // classify fieldTypes
+    nFields += grep(scalarFields_, fieldTypes);
+    nFields += grep(vectorFields_, fieldTypes);
+    nFields += grep(sphericalTensorFields_, fieldTypes);
+    nFields += grep(symmTensorFields_, fieldTypes);
+    nFields += grep(tensorFields_, fieldTypes);
+
+    if (Pstream::master)
+    {
+        if (debug)
+        {
+            Pout<< "timeName = " << obr_.time().timeName() << nl
+                << "scalarFields    " << scalarFields_ << nl
+                << "vectorFields    " << vectorFields_ << nl
+                << "sphTensorFields " << sphericalTensorFields_ << nl
+                << "symTensorFields " << symmTensorFields_ <<nl
+                << "tensorFields    " << tensorFields_ <<nl;
+        }
+
+        if (nFields > 0)
+        {
+            Pout<< "Creating directory " << outputPath_/obr_.time().timeName()
+                << nl << endl;
+
+            mkDir(outputPath_/obr_.time().timeName());
+        }
+    }
+
+    return nFields > 0;
+}
+
+
+void Foam::sampledSurfaces::mergeSurfaces()
+{
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
+    // Merge close points (1E-10 of mesh bounding box)
+    const scalar mergeTol = 1e-10;
+
+    const polyMesh& mesh = refCast<const polyMesh>(obr_);
+    const boundBox& bb = mesh.globalData().bb();
+
+    scalar mergeDim = mergeTol * mag(bb.max() - bb.min());
+
+    if (Pstream::master() && debug)
+    {
+        Pout<< nl << "Merging all points within "
+            << mergeDim << " meter" << endl;
+    }
+
+    mergeList_.setSize(size());
+    forAll(*this, surfI)
+    {
+        sampledSurface& s = operator[](surfI);
+
+        // Collect points from all processors
+        List<pointField> gatheredPoints(Pstream::nProcs());
+        gatheredPoints[Pstream::myProcNo()] = s.points();
+        Pstream::gatherList(gatheredPoints);
+
+        if (Pstream::master())
+        {
+            mergeList_[surfI].points = ListListOps::combine<pointField>
+            (
+                gatheredPoints,
+                accessOp<pointField>()
+            );
+        }
+
+        // Collect faces from all processors and renumber using sizes of
+        // gathered points
+        List<faceList> gatheredFaces(Pstream::nProcs());
+        gatheredFaces[Pstream::myProcNo()] = s.faces();
+        Pstream::gatherList(gatheredFaces);
+
+        if (Pstream::master())
+        {
+            mergeList_[surfI].faces = static_cast<const faceList&>
+            (
+                ListListOps::combineOffset<faceList>
+                (
+                    gatheredFaces,
+                    ListListOps::subSizes
+                    (
+                        gatheredPoints,
+                        accessOp<pointField>()
+                    ),
+                    accessOp<faceList>(),
+                    offsetOp<face>()
+                )
+            );
+        }
+
+        pointField newPoints;
+        labelList oldToNew;
+
+        bool hasMerged = mergePoints
+        (
+            mergeList_[surfI].points,
+            mergeDim,
+            false,                  // verbosity
+            oldToNew,
+            newPoints
+        );
+
+        if (hasMerged)
+        {
+            // Store point mapping
+            mergeList_[surfI].pointsMap.transfer(oldToNew);
+
+            // Copy points
+            mergeList_[surfI].points.transfer(newPoints);
+
+            // Relabel faces
+            faceList& faces = mergeList_[surfI].faces;
+
+            forAll(faces, faceI)
+            {
+                inplaceRenumber(mergeList_[surfI].pointsMap, faces[faceI]);
+            }
+
+            if (Pstream::master() && debug)
+            {
+                Pout<< "For surface " << surfI << " merged from "
+                    << mergeList_[surfI].pointsMap.size() << " points down to "
+                    << mergeList_[surfI].points.size()    << " points" << endl;
+            }
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sampledSurfaces::sampledSurfaces
+(
+    const word& name,
+    const objectRegistry& obr,
+    const dictionary& dict,
+    const bool loadFromFiles
+)
+:
+    PtrList<sampledSurface>(),
+    name_(name),
+    obr_(obr),
+    loadFromFiles_(loadFromFiles_),
+    outputPath_(fileName::null),
+    pMeshPtr_(NULL),
+    pInterpPtr_(NULL),
+    fieldNames_(),
+    interpolationScheme_(word::null),
+    writeFormat_(word::null),
+    mergeList_(),
+    scalarFields_(),
+    vectorFields_(),
+    sphericalTensorFields_(),
+    symmTensorFields_(),
+    tensorFields_()
+{
+    if (Pstream::parRun())
+    {
+        outputPath_ = obr_.time().path()/".."/name_;
+    }
+    else
+    {
+        outputPath_ = obr_.time().path()/name_;
+    }
+
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sampledSurfaces::~sampledSurfaces()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::sampledSurfaces::verbose(const bool verbosity)
+{
+    verbose_ = verbosity;
+}
+
+
+void Foam::sampledSurfaces::write()
+{
+    if (size() && checkFieldTypes())
+    {
+        sampleAndWrite(scalarFields_);
+        sampleAndWrite(vectorFields_);
+        sampleAndWrite(sphericalTensorFields_);
+        sampleAndWrite(symmTensorFields_);
+        sampleAndWrite(tensorFields_);
+    }
+}
+
+
+void Foam::sampledSurfaces::read(const dictionary& dict)
+{
+    fieldNames_ = wordList(dict.lookup("fields"));
+
+    interpolationScheme_ = "cell";
+    if (dict.found("interpolationScheme"))
+    {
+        dict.lookup("interpolationScheme") >> interpolationScheme_;
+    }
+
+    writeFormat_ = "null";
+    if (dict.found("surfaceFormat"))
+    {
+        dict.lookup("surfaceFormat") >> writeFormat_;
+    }
+
+
+    PtrList<sampledSurface> newList
+    (
+        dict.lookup("surfaces"),
+        sampledSurface::iNew(refCast<const polyMesh>(obr_))
+    );
+
+    transfer(newList);
+    mergeSurfaces();
+
+    if (Pstream::master() && debug)
+    {
+        Pout<< "sample fields:" << fieldNames_ << nl
+            << "sample surfaces:" << nl << "(" << nl;
+
+        forAll(*this, surfI)
+        {
+            Pout << "  " << operator[](surfI) << endl;
+        }
+        Pout << ")" << endl;
+    }
+}
+
+
+void Foam::sampledSurfaces::correct()
+{
+    forAll(*this, surfI)
+    {
+        operator[](surfI).correct(true);
+    }
+
+    // reset interpolation for later
+    pMeshPtr_.clear();
+    pInterpPtr_.clear();
+
+    mergeSurfaces();
+}
+
+
+void Foam::sampledSurfaces::updateMesh(const mapPolyMesh&)
+{
+    correct();
+}
+
+
+void Foam::sampledSurfaces::movePoints(const pointField&)
+{
+    correct();
+}
+
+
+// ************************************************************************* //

src/sampling/sampledSurface/surfaces/sampledSurfaces.H

@@ -0,0 +1,277 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    Foam::sampledSurfaces
+
+Description
+    Set of surfaces to sample.
+
+    The write() method is used to sample and write files.
+
+SourceFiles
+    sampledSurfaces.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sampledSurfaces_H
+#define sampledSurfaces_H
+
+#include "sampledSurface.H"
+#include "surfaceWriter.H"
+#include "volFields.H"
+#include "volPointInterpolation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class objectRegistry;
+class dictionary;
+
+/*---------------------------------------------------------------------------*\
+                      Class sampledSurfaces Declaration
+\*---------------------------------------------------------------------------*/
+
+class sampledSurfaces
+:
+    public PtrList<sampledSurface>
+{
+    // Private classes
+
+        //- Class used for grouping field types
+        template<class Type>
+        class fieldGroup
+        :
+            public wordList
+        {
+        public:
+            //- surface formatter
+            autoPtr<surfaceWriter<Type> > formatter;
+
+            //- Construct null
+            fieldGroup()
+            :
+                wordList(0),
+                formatter(NULL)
+            {}
+
+            //- Construct for a particular surface format
+            fieldGroup(const word& writeFormat)
+            :
+                wordList(0),
+                formatter(surfaceWriter<Type>::New(writeFormat))
+            {}
+
+            //- Construct for a particular surface format and a list of field names
+            fieldGroup
+            (
+                const word& writeFormat,
+                const wordList& fieldNames
+            )
+            :
+                wordList(fieldNames),
+                formatter(surfaceWriter<Type>::New(writeFormat))
+            {}
+
+            void operator=(const word& writeFormat)
+            {
+                formatter = surfaceWriter<Type>::New(writeFormat);
+            }
+
+            void operator=(const wordList& fieldNames)
+            {
+                wordList::operator=(fieldNames);
+            }
+
+        };
+
+        //- Class used for surface merging information
+        class mergeInfo
+        {
+        public:
+            pointField points;
+            faceList   faces;
+            labelList  pointsMap;
+        };
+
+
+    // Static data members
+
+            //- output verbosity
+            static bool verbose_;
+
+
+    // Private data
+
+        //- Name of this set of surfaces,
+        //  Also used as the name of the sampledSurfaces directory.
+        word name_;
+
+        //- Const reference to objectRegistry
+        const objectRegistry& obr_;
+
+        //- Load fields from files (not from objectRegistry)
+        bool loadFromFiles_;
+
+        //- output path
+        fileName outputPath_;
+
+        //- pointMesh for interpolation
+        autoPtr<pointMesh> pMeshPtr_;
+
+        //- volPointInterpolation for interpolation
+        autoPtr<volPointInterpolation> pInterpPtr_;
+
+
+        // Read from dictonary
+
+            //- Names of fields to sample
+            wordList fieldNames_;
+
+            //- Interpolation scheme to use
+            word interpolationScheme_;
+
+            //- Output format to use
+            word writeFormat_;
+
+
+        // surfaces
+
+            //- information for merging surfaces
+            List<mergeInfo> mergeList_;
+
+
+        // Calculated
+
+            //- Categorized scalar/vector/tensor fields
+            fieldGroup<scalar> scalarFields_;
+            fieldGroup<vector> vectorFields_;
+            fieldGroup<sphericalTensor> sphericalTensorFields_;
+            fieldGroup<symmTensor> symmTensorFields_;
+            fieldGroup<tensor> tensorFields_;
+
+
+    // Private Member Functions
+
+        //- classify field types, return true if nFields > 0
+        bool checkFieldTypes();
+
+        //- merge points on surfaces
+        void mergeSurfaces();
+
+        //- Correct for mesh changes
+        void correct();
+
+        //- Find the fields in the list of the given type, return count
+        template<class Type>
+        label grep
+        (
+            fieldGroup<Type>& fieldList,
+            const wordList& fieldTypes
+        ) const;
+
+        //- set interpolator for the field
+        template<class Type>
+        autoPtr<interpolation<Type> > setInterpolator
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        );
+
+        //- Sample and write a particular volume field
+        template<class Type>
+        void sampleAndWrite
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&,
+            const surfaceWriter<Type>& formatter
+        );
+
+        //- Sample and write all the fields of the given type
+        template <class Type>
+        void sampleAndWrite(fieldGroup<Type>&);
+
+        //- Disallow default bitwise copy construct and assignment
+        sampledSurfaces(const sampledSurfaces&);
+        void operator=(const sampledSurfaces&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("surfaces");
+
+
+    // Constructors
+
+        //- Construct for given objectRegistry and dictionary
+        //  allow the possibility to load fields from files
+        sampledSurfaces
+        (
+            const word& name,
+            const objectRegistry&,
+            const dictionary&,
+            const bool loadFromFiles = false
+        );
+
+
+    // Destructor
+
+        virtual ~sampledSurfaces();
+
+
+    // Member Functions
+
+       //- set verbosity level
+        void verbose(const bool verbosity = true);
+
+        //- Sample and write
+        virtual void write();
+
+        //- Read the sampledSurfaces
+        virtual void read(const dictionary&);
+
+        //- Update for changes of mesh
+        virtual void updateMesh(const mapPolyMesh&);
+
+        //- Update for mesh point-motion
+        virtual void movePoints(const pointField&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "sampledSurfacesTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/sampledSurface/surfaces/surfacesFunctionObject.C

@@ -0,0 +1,43 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "surfacesFunctionObject.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineNamedTemplateTypeNameAndDebug(surfacesFunctionObject, 0);
+
+    addToRunTimeSelectionTable
+    (
+        functionObject,
+        surfacesFunctionObject,
+        dictionary
+    );
+}
+
+// ************************************************************************* //

src/sampling/sampledSurface/surfaces/surfacesFunctionObject.H

@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Typedef
+    surfacesFunctionObject
+
+Description
+    FunctionObject wrapper around surfaces to allow them to be created via the
+    functions list within controlDict.
+
+SourceFiles
+    surfacesFunctionObject.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfacesFunctionObject_H
+#define surfacesFunctionObject_H
+
+#include "sampledSurfaces.H"
+#include "OutputFilterFunctionObject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef OutputFilterFunctionObject<sampledSurfaces> surfacesFunctionObject;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //