ENH: Updated Curle function object

The function object now computes the acoustic pressure at a list of user
specified locations, or from the face centres of a user-supplied surface.
When operating on an input surface, the output can be written back to the
surface or as a list of point values.

Example of function object specification:

    Curle1
    {
        type            Curle;
        libs            ("libfieldFunctionObjects.so");
        ...
        patches         (surface1 surface2);
        c0              330;

        // Input - either points or surface

        input           points;
        observerPositions ((0 0 0)(1 0 0));

        //input           surface;
        //surface         "inputSurface.obj"

        // Output - either points or surface
        output          points;

        //output          surface;
        //surfaceType     ensight;
    }

    Where the entries comprise:
        Property     | Description                  | Required | Default value
        type         | Type name: Curle             | yes      |
        p            | Pressure field name          | no       | p
        patches      | Sound generation patch names | yes      |
        c0           | Reference speed of sound     | yes      |
        input        | Input type                   | yes      |
        observerPositions | List of observer positions (x y z) | no      |
        surface      | Input surface file name      | no       |
        output       | Output type                  | yes      |
        surfaceType  | Output surface type          | no       |

src/functionObjects/field/Curle/Curle.C

@@ -43,57 +43,12 @@ namespace functionObjects
 }
 }
 
-
+const Foam::Enum<Foam::functionObjects::Curle::modeType>
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+Foam::functionObjects::Curle::modeTypeNames_
-
+({
-bool Foam::functionObjects::Curle::calc()
+    { modeType::point, "point" },
-{
+    { modeType::surface, "surface" },
-    if (foundObject<volScalarField>(fieldName_))
+});
-    {
-        // Evaluate pressure force time derivative
-
-        const volScalarField& p = lookupObject<volScalarField>(fieldName_);
-        const volScalarField dpdt(scopedName("dpdt"), fvc::ddt(p));
-        const volScalarField::Boundary& dpdtBf = dpdt.boundaryField();
-        const surfaceVectorField::Boundary& SfBf = mesh_.Sf().boundaryField();
-
-        dimensionedVector dfdt("dfdt", p.dimensions()*dimArea/dimTime, Zero);
-
-        for (const label patchi : patchSet_)
-        {
-            dfdt.value() += sum(dpdtBf[patchi]*SfBf[patchi]);
-        }
-
-        reduce(dfdt.value(), sumOp<vector>());
-
-
-        // Construct and store Curle acoustic pressure
-
-        const volVectorField& C = mesh_.C();
-
-        auto tpDash = tmp<volScalarField>::New
-        (
-            IOobject
-            (
-                resultName_,
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar(p.dimensions(), Zero)
-        );
-        auto& pDash = tpDash.ref();
-
-        const volVectorField d(scopedName("d"), C - x0_);
-        pDash = (d/magSqr(d) & dfdt)/(4.0*mathematical::pi*c0_);
-
-        return store(resultName_, tpDash);
-    }
-
-    return false;
-}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -105,14 +60,17 @@ Foam::functionObjects::Curle::Curle
     const dictionary& dict
 )
 :
-    fieldExpression(name, runTime, dict, "p"),
+    fvMeshFunctionObject(name, runTime, dict),
+    writeFile(mesh_, name),
+    pName_("p"),
     patchSet_(),
-    x0_("x0", dimLength,  Zero),
+    observerPositions_(),
-    c0_("c0", dimVelocity, Zero)
+    c0_(0),
+    rawFilePtrs_(),
+    inputSurface_(),
+    surfaceWriterPtr_(nullptr)
 {
     read(dict);
-
-    setResultName(typeName, fieldName_);
 }
 
 
@@ -120,57 +78,192 @@ Foam::functionObjects::Curle::Curle
 
 bool Foam::functionObjects::Curle::read(const dictionary& dict)
 {
-    if (fieldExpression::read(dict))
+    if (!(fvMeshFunctionObject::read(dict) && writeFile::read(dict)))
     {
-        patchSet_ =
+        return false;
-            mesh_.boundaryMesh().patchSet
-            (
-                dict.get<wordRes>("patches")
-            );
-
-        if (patchSet_.empty())
-        {
-            WarningInFunction
-                << "No patches defined"
-                << endl;
-
-            return false;
-        }
-
-        // Read the reference speed of sound
-        dict.readEntry("c0", c0_);
-
-        if (c0_.value() < VSMALL)
-        {
-            FatalErrorInFunction
-                << "Reference speed of sound = " << c0_
-                << " cannot be negative or zero."
-                << abort(FatalError);
-        }
-
-        // Set the location of the effective point source to the area-average
-        // of the patch face centres
-        const volVectorField::Boundary& Cbf = mesh_.C().boundaryField();
-        const surfaceScalarField::Boundary& magSfBf =
-            mesh_.magSf().boundaryField();
-
-        x0_.value() = vector::zero;
-        scalar sumMagSf = 0;
-        for (auto patchi : patchSet_)
-        {
-            x0_.value() += sum(Cbf[patchi]*magSfBf[patchi]);
-            sumMagSf += sum(magSfBf[patchi]);
-        }
-
-        reduce(x0_.value(), sumOp<vector>());
-        reduce(sumMagSf, sumOp<scalar>());
-
-        x0_.value() /= sumMagSf + ROOTVSMALL;
-
-        return true;
     }
 
-    return false;
+    dict.readIfPresent("p", pName_);
+
+    patchSet_ = mesh_.boundaryMesh().patchSet(dict.get<wordRes>("patches"));
+
+    if (patchSet_.empty())
+    {
+        WarningInFunction
+            << "No patches defined"
+            << endl;
+
+        return false;
+    }
+
+    const modeType inputMode = modeTypeNames_.get("input", dict);
+
+    switch (inputMode)
+    {
+        case modeType::point:
+        {
+            observerPositions_ = dict.get<List<point>>("observerPositions");
+            break;
+        }
+        case modeType::surface:
+        {
+            const fileName fName = (dict.get<fileName>("surface")).expand();
+            inputSurface_ = MeshedSurface<face>(fName);
+
+            observerPositions_ = inputSurface_.Cf();
+        }
+    }
+
+    if (observerPositions_.empty())
+    {
+        WarningInFunction
+            << "No observer positions defined"
+            << endl;
+
+        return false;
+    }
+
+    const auto outputMode = modeTypeNames_.get("output", dict);
+
+    switch (outputMode)
+    {
+        case modeType::point:
+        {
+            rawFilePtrs_.setSize(observerPositions_.size());
+            forAll(rawFilePtrs_, filei)
+            {
+                rawFilePtrs_.set
+                (
+                    filei,
+                    createFile("observer" + Foam::name(filei))
+                );
+
+                if (rawFilePtrs_.set(filei))
+                {
+                    OFstream& os = rawFilePtrs_[filei];
+
+                    writeHeaderValue(os, "Observer", filei);
+                    writeHeaderValue(os, "Position", observerPositions_[filei]);
+                    writeCommented(os, "Time");
+                    writeTabbed(os, "p(Curle)");
+                    os  << endl;
+                }
+            }
+            break;
+        }
+        case modeType::surface:
+        {
+            if (inputMode != modeType::surface)
+            {
+                FatalIOErrorInFunction(dict)
+                    << "Surface output is only available when input mode is "
+                    << "type '" << modeTypeNames_[modeType::surface] << "'"
+                    << abort(FatalIOError);
+            }
+
+            const word surfaceType(dict.get<word>("surfaceType"));
+            surfaceWriterPtr_ = surfaceWriter::New
+            (
+                surfaceType,
+                dict.subOrEmptyDict("formatOptions").subOrEmptyDict(surfaceType)
+            );
+
+            // Use outputDir/TIME/surface-name
+            surfaceWriterPtr_->useTimeDir() = true;
+        }
+    }
+
+    // Read the reference speed of sound
+    dict.readEntry("c0", c0_);
+
+    if (c0_ < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Reference speed of sound = " << c0_
+            << " cannot be negative or zero."
+            << abort(FatalError);
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::Curle::execute()
+{
+    if (!foundObject<volScalarField>(pName_))
+    {
+        return false;
+    }
+
+    const volScalarField& p = lookupObject<volScalarField>(pName_);
+    const volScalarField::Boundary& pBf = p.boundaryField();
+    const volScalarField dpdt(scopedName("dpdt"), fvc::ddt(p));
+    const volScalarField::Boundary& dpdtBf = dpdt.boundaryField();
+    const surfaceVectorField::Boundary& SfBf = mesh_.Sf().boundaryField();
+    const surfaceVectorField::Boundary& CfBf = mesh_.Cf().boundaryField();
+
+    scalarField pDash(observerPositions_.size(), 0);
+
+    for (auto patchi : patchSet_)
+    {
+        const scalarField& pp = pBf[patchi];
+        const scalarField& dpdtp = dpdtBf[patchi];
+        const vectorField& Sfp = SfBf[patchi];
+        const vectorField& Cfp = CfBf[patchi];
+
+        forAll(observerPositions_, pointi)
+        {
+            const vectorField r(Cfp - observerPositions_[pointi]);
+            const scalarField invMagR(1/(mag(r) + ROOTVSMALL));
+
+            pDash[pointi] +=
+                sum((pp*sqr(invMagR) + invMagR/c0_*dpdtp)*(Sfp & r));
+        }
+    }
+    Pstream::listCombineGather(pDash, plusEqOp<scalar>());
+    Pstream::listCombineScatter(pDash);
+
+    if (surfaceWriterPtr_)
+    {
+        if (Pstream::master())
+        {
+            // Time-aware, with time spliced into the output path
+            surfaceWriterPtr_->beginTime(time_);
+
+            surfaceWriterPtr_->open
+            (
+                inputSurface_.points(),
+                inputSurface_.surfFaces(),
+                (baseFileDir()/name()/"surface"),
+                false  // serial - already merged
+            );
+
+            surfaceWriterPtr_->write("p(Curle)", pDash);
+
+            surfaceWriterPtr_->endTime();
+            surfaceWriterPtr_->clear();
+        }
+    }
+    else
+    {
+        forAll(observerPositions_, pointi)
+        {
+            if (rawFilePtrs_.set(pointi))
+            {
+                OFstream& os = rawFilePtrs_[pointi];
+                writeCurrentTime(os);
+                os  << pDash[pointi] << endl;
+            }
+        }
+    }
+
+    return true;
+}
+
+
+bool Foam::functionObjects::Curle::write()
+{
+    return true;
 }
 
 

src/functionObjects/field/Curle/Curle.H

@@ -35,15 +35,21 @@ Description
     Curle's analogy is implemented as:
 
     \f[
-        p' = \frac{1}{4 \pi c_0}\frac{\vec d}{|\vec d|^2}\cdot\frac{d\vec F}{dt}
+        p' = \frac{1}{4 \pi}
+             \frac{\vec{r}}{| \vec{r} | ^2} \cdot
+             \left(
+                \frac{\vec{F}}{| \vec{r} |}
+              + \frac{1}{c_0}\frac{d(\vec{F})}{d(t)}
+             \right)
     \f]
 
     where
     \vartable
-      p'        | Curle's acoustic pressure     [Pa] or [Pa \f$(m^3/\rho)\f$]
+        p'          | Curle's acoustic pressure [Pa] or [Pa (m3/kg)]
-      c_0       | Reference speed of sound              [m/s]
+        c_0         | Reference speed of sound [m/s]
-      \vec d    | Distance vector to observer locations [m]
+        \vec{r}     | Distance vector to observer locations [m]
-      \vec F    | Force                         [N] or [N (\f$m^3/\rho\f$)]
+        \vec{F}     | Force [N] or [N (m3/kg)]
+        t           | time [s]
     \endvartable
 
     Operands:
@@ -62,16 +68,27 @@ Usage
     \verbatim
     Curle1
     {
-        // Mandatory entries (unmodifiable)
         type            Curle;
-        libs            (fieldFunctionObjects);
+        libs            ("libfieldFunctionObjects.so");
-
-        // Mandatory entries (runtime modifiable)
-        patches         (<patch1> <patch2> ... <patchN>)
-        c0              343;
-
-        // Optional (inherited) entries
         ...
+        patches         (surface1 surface2);
+        c0              330;
+
+
+        // Input - either points or surface
+
+        input           points;
+        observerPositions ((0 0 0)(1 0 0));
+
+        //input           surface;
+        //surface         "inputSurface.obj"
+
+
+        // Output - either points or surface
+        output          points;
+
+        //output          surface;
+        //surfaceType     ensight;
     }
     \endverbatim
 
@@ -80,8 +97,14 @@ Usage
       Property     | Description                        | Type | Req'd | Dflt
       type         | Type name: Curle                   | word |  yes  | -
       libs         | Library name: fieldFunctionObjects | word |  yes  | -
-      patches      | Names of the operand patches       | wordList | yes | -
+      p            | Pressure field name                | word |  no   | p
-      c0           | Reference speed of sound [m/s]     | scalar | yes | -
+      patches      | Sound generation patch names       | wordList | yes | -
+      c0           | Reference speed of sound           | scalar | yes | -
+      input        | Input type                         | word | yes   | -
+      observerPositions | List of observer positions (x y z) | pointList | no |-
+      surface      | Input surface file name            | word | no    | -
+      output       | Output type                        | word | yes   | -
+      surfaceType  | Output surface type                | word | no    | -
     \endtable
 
     The inherited entries are elaborated in:
@@ -93,7 +116,7 @@ Usage
 See also
     - Foam::functionObject
     - Foam::functionObjects::fvMeshFunctionObject
-    - Foam::functionObjects::fieldExpression
+    - Foam::functionObjects::writeFile
     - ExtendedCodeGuide::functionObjects::field::Curle
 
 SourceFiles
@@ -104,9 +127,11 @@ SourceFiles
 #ifndef functionObjects_Curle_H
 #define functionObjects_Curle_H
 
-#include "fieldExpression.H"
+#include "fvMeshFunctionObject.H"
-#include "dimensionedScalar.H"
+#include "writeFile.H"
-#include "dimensionedVector.H"
+#include "Enum.H"
+#include "MeshedSurface.H"
+#include "surfaceWriter.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -121,26 +146,40 @@ namespace functionObjects
 
 class Curle
 :
-    public fieldExpression
+    public fvMeshFunctionObject,
+    public writeFile
 {
+    enum class modeType
+    {
+        point,
+        surface
+    };
+
+    static const Enum<modeType> modeTypeNames_;
+
+
     // Private Data
 
-        // Read from dictionary
+        //- Name of pressure field; default = p
+        word pName_;
 
-            //- Patches to integrate forces over
+        //- Patches to integrate forces over
-            labelHashSet patchSet_;
+        labelHashSet patchSet_;
 
-            //- Area-averaged centre of patch faces
+        //- Observer positions
-            dimensionedVector x0_;
+        List<point> observerPositions_;
 
-            //- Reference speed of sound
+        //- Reference speed of sound
-            dimensionedScalar c0_;
+        scalar c0_;
 
+        //- Output files when sampling points
+        PtrList<OFstream> rawFilePtrs_;
 
-protected:
+        //- Input surface when sampling a surface
+        MeshedSurface<face> inputSurface_;
 
-        //- Calculate the acoustic pressure field and return true if successful
+        //- Ouput surface when sampling a surface
-        virtual bool calc();
+        autoPtr<surfaceWriter> surfaceWriterPtr_;
 
 
 public:
@@ -174,6 +213,10 @@ public:
 
         //- Read the Curle data
         virtual bool read(const dictionary&);
+
+        virtual bool execute();
+
+        virtual bool write();
 };