ENH: Adding surface fields to probes and patchProbes

src/sampling/probes/patchProbes.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -202,6 +202,7 @@ Foam::patchProbes::patchProbes
     // Not easy to workaround (apart from feeding through flag into constructor)
     // so clear out any cells found for now.
     elementList_.clear();
+    faceList_.clear();
 
     read(dict);
 }
@@ -222,6 +223,12 @@ void Foam::patchProbes::write()
         sampleAndWrite(sphericalTensorFields_);
         sampleAndWrite(symmTensorFields_);
         sampleAndWrite(tensorFields_);
+
+        sampleAndWriteSurfaceFields(surfaceScalarFields_);
+        sampleAndWriteSurfaceFields(surfaceVectorFields_);
+        sampleAndWriteSurfaceFields(surfaceSphericalTensorFields_);
+        sampleAndWriteSurfaceFields(surfaceSymmTensorFields_);
+        sampleAndWriteSurfaceFields(surfaceTensorFields_);
     }
 }
 

src/sampling/probes/patchProbes.H

@@ -74,11 +74,24 @@ class patchProbes
         );
 
 
+        //- Sample and write a particular surface field
+        template<class Type>
+        void sampleAndWrite
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        );
+
+
         //- Sample and write all the fields of the given type
         template <class Type>
         void sampleAndWrite(const fieldGroup<Type>&);
 
 
+         //- Sample and write all the surface fields of the given type
+        template<class Type>
+        void sampleAndWriteSurfaceFields(const fieldGroup<Type>&);
+
+
         //- Sample a volume field at all locations
         template<class Type>
         tmp<Field<Type> > sample
@@ -87,6 +100,14 @@ class patchProbes
         ) const;
 
 
+        //- Sample a surface field at all locations
+        template<class Type>
+        tmp<Field<Type> > sample
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        ) const;
+
+
         //- Sample a single field on all sample locations
         template <class Type>
         tmp<Field<Type> > sample(const word& fieldName) const;

src/sampling/probes/patchProbesTemplates.C

@@ -54,6 +54,30 @@ void Foam::patchProbes::sampleAndWrite
 }
 
 
+template<class Type>
+void Foam::patchProbes::sampleAndWrite
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
+)
+{
+    Field<Type> values(sample(sField));
+
+    if (Pstream::master())
+    {
+        unsigned int w = IOstream::defaultPrecision() + 7;
+        OFstream& probeStream = *probeFilePtrs_[sField.name()];
+
+        probeStream << setw(w) << sField.time().value();
+
+        forAll(values, probeI)
+        {
+            probeStream << ' ' << setw(w) << values[probeI];
+        }
+        probeStream << endl;
+    }
+}
+
+
 template <class Type>
 void Foam::patchProbes::sampleAndWrite
 (
@@ -106,6 +130,58 @@ void Foam::patchProbes::sampleAndWrite
 }
 
 
+template<class Type>
+void Foam::patchProbes::sampleAndWriteSurfaceFields
+(
+    const fieldGroup<Type>& fields
+)
+{
+    forAll(fields, fieldI)
+    {
+        if (loadFromFiles_)
+        {
+            sampleAndWrite
+            (
+                GeometricField<Type, fvsPatchField, surfaceMesh>
+                (
+                    IOobject
+                    (
+                        fields[fieldI],
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE,
+                        false
+                    ),
+                    mesh_
+                )
+            );
+        }
+        else
+        {
+            objectRegistry::const_iterator iter = mesh_.find(fields[fieldI]);
+
+            if
+            (
+                iter != objectRegistry::end()
+             && iter()->type()
+             == GeometricField<Type, fvsPatchField, surfaceMesh>::typeName
+            )
+            {
+                sampleAndWrite
+                (
+                    mesh_.lookupObject
+                    <GeometricField<Type, fvsPatchField, surfaceMesh> >
+                    (
+                        fields[fieldI]
+                    )
+                );
+            }
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -159,4 +235,39 @@ Foam::patchProbes::sample(const word& fieldName) const
 }
 
 
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::patchProbes::sample
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& sField
+) const
+{
+    const Type unsetVal(-VGREAT*pTraits<Type>::one);
+
+    tmp<Field<Type> > tValues
+    (
+        new Field<Type>(this->size(), unsetVal)
+    );
+
+    Field<Type>& values = tValues();
+
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    forAll(*this, probeI)
+    {
+        label faceI = elementList_[probeI];
+
+        if (faceI >= 0)
+        {
+            label patchI = patches.whichPatch(faceI);
+            label localFaceI = patches[patchI].whichFace(faceI);
+            values[probeI] = sField.boundaryField()[patchI][localFaceI];
+        }
+    }
+
+    Pstream::listCombineGather(values, isNotEqOp<Type>());
+    Pstream::listCombineScatter(values);
+
+    return tValues;
+}
 // ************************************************************************* //

src/sampling/probes/probes.C

@@ -41,16 +41,45 @@ void Foam::probes::findElements(const fvMesh& mesh)
     elementList_.clear();
     elementList_.setSize(size());
 
+    faceList_.clear();
+    faceList_.setSize(size());
+
     forAll(*this, probeI)
     {
         const vector& location = operator[](probeI);
 
-        elementList_[probeI] = mesh.findCell(location);
+        const label cellI = mesh.findCell(location);
 
-        if (debug && elementList_[probeI] != -1)
+        elementList_[probeI] = cellI;
+
+        if (cellI != -1)
+        {
+            const labelList& cellFaces = mesh.cells()[cellI];
+            const vector& cellCentre = mesh.cellCentres()[cellI];
+            scalar minDistance = GREAT;
+            label minFaceID = -1;
+            forAll (cellFaces, i)
+            {
+                label faceI = cellFaces[i];
+                vector dist = mesh.faceCentres()[faceI] - cellCentre;
+                if (mag(dist) < minDistance)
+                {
+                    minDistance = mag(dist);
+                    minFaceID = faceI;
+                }
+            }
+            faceList_[probeI] = minFaceID;
+        }
+        else
+        {
+            faceList_[probeI] = -1;
+        }
+
+        if (debug && (elementList_[probeI] != -1 || faceList_[probeI] != -1))
         {
             Pout<< "probes : found point " << location
-                << " in cell " << elementList_[probeI] << endl;
+                << " in cell " << elementList_[probeI]
+                << " and face " << faceList_[probeI] << endl;
         }
     }
 
@@ -60,25 +89,36 @@ void Foam::probes::findElements(const fvMesh& mesh)
     {
         const vector& location = operator[](probeI);
         label cellI = elementList_[probeI];
+        label faceI =  faceList_[probeI];
 
         // Check at least one processor with cell.
         reduce(cellI, maxOp<label>());
+        reduce(faceI, maxOp<label>());
 
         if (cellI == -1)
         {
             if (Pstream::master())
             {
-                WarningIn("probes::read()")
+                WarningIn("findElements::findElements(const fvMesh&)")
                     << "Did not find location " << location
                     << " in any cell. Skipping location." << endl;
             }
         }
+        else if (faceI == -1)
+        {
+            if (Pstream::master())
+            {
+                WarningIn("probes::findElements(const fvMesh&)")
+                    << "Did not find location " << location
+                    << " in any face. Skipping location." << endl;
+            }
+        }
         else
         {
             // Make sure location not on two domains.
             if (elementList_[probeI] != -1 && elementList_[probeI] != cellI)
             {
-                WarningIn("probes::read()")
+                WarningIn("probes::findElements(const fvMesh&)")
                     << "Location " << location
                     << " seems to be on multiple domains:"
                     << " cell " << elementList_[probeI]
@@ -89,6 +129,20 @@ void Foam::probes::findElements(const fvMesh& mesh)
                     << " a processor patch. Change the location slightly"
                     << " to prevent this." << endl;
             }
+
+            if (faceList_[probeI] != -1 && faceList_[probeI] != faceI)
+            {
+                WarningIn("probes::findElements(const fvMesh&)")
+                    << "Location " << location
+                    << " seems to be on multiple domains:"
+                    << " cell " << faceList_[probeI]
+                    << " on my domain " << Pstream::myProcNo()
+                        << " and face " << faceI << " on some other domain."
+                    << endl
+                    << "This might happen if the probe location is on"
+                    << " a processor patch. Change the location slightly"
+                    << " to prevent this." << endl;
+            }
         }
     }
 }
@@ -109,6 +163,12 @@ Foam::label Foam::probes::prepare()
         currentFields.insert(symmTensorFields_);
         currentFields.insert(tensorFields_);
 
+        currentFields.insert(surfaceScalarFields_);
+        currentFields.insert(surfaceVectorFields_);
+        currentFields.insert(surfaceSphericalTensorFields_);
+        currentFields.insert(surfaceSymmTensorFields_);
+        currentFields.insert(surfaceTensorFields_);
+
         if (debug)
         {
             Info<< "Probing fields:" << currentFields << nl
@@ -239,6 +299,12 @@ void Foam::probes::write()
         sampleAndWrite(sphericalTensorFields_);
         sampleAndWrite(symmTensorFields_);
         sampleAndWrite(tensorFields_);
+
+        sampleAndWriteSurfaceFields(surfaceScalarFields_);
+        sampleAndWriteSurfaceFields(surfaceVectorFields_);
+        sampleAndWriteSurfaceFields(surfaceSphericalTensorFields_);
+        sampleAndWriteSurfaceFields(surfaceSymmTensorFields_);
+        sampleAndWriteSurfaceFields(surfaceTensorFields_);
     }
 }
 

src/sampling/probes/probes.H

@@ -42,7 +42,8 @@ SourceFiles
 #include "polyMesh.H"
 #include "pointField.H"
 #include "volFieldsFwd.H"
-
+#include "surfaceFieldsFwd.H"
+#include "surfaceMesh.H"
 #include "wordReList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -104,16 +105,26 @@ protected:
 
         // Calculated
 
-            //- Categorized scalar/vector/tensor fields
+            //- Categorized scalar/vector/tensor vol fields
             fieldGroup<scalar> scalarFields_;
             fieldGroup<vector> vectorFields_;
             fieldGroup<sphericalTensor> sphericalTensorFields_;
             fieldGroup<symmTensor> symmTensorFields_;
             fieldGroup<tensor> tensorFields_;
 
+            //- Categorized scalar/vector/tensor surf fields
+            fieldGroup<scalar> surfaceScalarFields_;
+            fieldGroup<vector> surfaceVectorFields_;
+            fieldGroup<sphericalTensor> surfaceSphericalTensorFields_;
+            fieldGroup<symmTensor> surfaceSymmTensorFields_;
+            fieldGroup<tensor> surfaceTensorFields_;
+
             // Cells to be probed (obtained from the locations)
             labelList elementList_;
 
+            // Faces to be probed
+            labelList faceList_;
+
             //- Current open files
             HashPtrTable<OFstream> probeFilePtrs_;
 
@@ -145,10 +156,22 @@ private:
             const GeometricField<Type, fvPatchField, volMesh>&
         );
 
+
+        //- Sample and write a particular surface field
+        template<class Type>
+        void sampleAndWrite
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        );
+
         //- Sample and write all the fields of the given type
         template<class Type>
         void sampleAndWrite(const fieldGroup<Type>&);
 
+        //- Sample and write all the surface fields of the given type
+        template<class Type>
+        void sampleAndWriteSurfaceFields(const fieldGroup<Type>&);
+
         //- Disallow default bitwise copy construct
         probes(const probes&);
 
@@ -242,9 +265,21 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>&
         ) const;
 
-        //- Sample a single field on all sample locations
+        //- Sample a single vol field on all sample locations
         template <class Type>
         tmp<Field<Type> > sample(const word& fieldName) const;
+
+        //- Sample a single scalar field on all sample locations
+        template <class Type>
+        tmp<Field<Type> > sampleSurfaceFields(const word& fieldName) const;
+
+        //- Sample a surface field at all locations
+        template<class Type>
+        tmp<Field<Type> > sample
+        (
+            const GeometricField<Type, fvsPatchField, surfaceMesh>&
+        ) const;
+
 };
 
 

src/sampling/probes/probesGrouping.C

@@ -25,6 +25,7 @@ License
 
 #include "probes.H"
 #include "volFields.H"
+#include "surfaceFields.H"
 #include "IOobjectList.H"
 #include "stringListOps.H"
 
@@ -37,6 +38,12 @@ void Foam::probes::clearFieldGroups()
     sphericalTensorFields_.clear();
     symmTensorFields_.clear();
     tensorFields_.clear();
+
+    surfaceScalarFields_.clear();
+    surfaceVectorFields_.clear();
+    surfaceSphericalTensorFields_.clear();
+    surfaceSymmTensorFields_.clear();
+    surfaceTensorFields_.clear();
 }
 
 
@@ -71,6 +78,31 @@ Foam::label Foam::probes::appendFieldGroup
         tensorFields_.append(fieldName);
         return 1;
     }
+    else if (fieldType == surfaceScalarField::typeName)
+    {
+        surfaceScalarFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == surfaceVectorField::typeName)
+    {
+        surfaceVectorFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == surfaceSphericalTensorField::typeName)
+    {
+        surfaceSphericalTensorFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == surfaceSymmTensorField::typeName)
+    {
+        surfaceSymmTensorFields_.append(fieldName);
+        return 1;
+    }
+    else if (fieldType == surfaceTensorField::typeName)
+    {
+        surfaceTensorFields_.append(fieldName);
+        return 1;
+    }
 
     return 0;
 }

src/sampling/probes/probesTemplates.C

@@ -25,6 +25,7 @@ License
 
 #include "probes.H"
 #include "volFields.H"
+#include "surfaceFields.H"
 #include "IOmanip.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -86,11 +87,32 @@ void Foam::probes::sampleAndWrite
 }
 
 
-template <class Type>
+template<class Type>
 void Foam::probes::sampleAndWrite
 (
-    const fieldGroup<Type>& fields
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& vField
 )
+{
+    Field<Type> values(sample(vField));
+
+    if (Pstream::master())
+    {
+        unsigned int w = IOstream::defaultPrecision() + 7;
+        OFstream& os = *probeFilePtrs_[vField.name()];
+
+        os  << setw(w) << vField.time().value();
+
+        forAll(values, probeI)
+        {
+            os  << ' ' << setw(w) << values[probeI];
+        }
+        os  << endl;
+    }
+}
+
+
+template <class Type>
+void Foam::probes::sampleAndWrite(const fieldGroup<Type>& fields)
 {
     forAll(fields, fieldI)
     {
@@ -138,6 +160,54 @@ void Foam::probes::sampleAndWrite
 }
 
 
+template<class Type>
+void Foam::probes::sampleAndWriteSurfaceFields(const fieldGroup<Type>& fields)
+{
+    forAll(fields, fieldI)
+    {
+        if (loadFromFiles_)
+        {
+            sampleAndWrite
+            (
+                GeometricField<Type, fvsPatchField, surfaceMesh>
+                (
+                    IOobject
+                    (
+                        fields[fieldI],
+                        mesh_.time().timeName(),
+                        mesh_,
+                        IOobject::MUST_READ,
+                        IOobject::NO_WRITE,
+                        false
+                    ),
+                    mesh_
+                )
+            );
+        }
+        else
+        {
+            objectRegistry::const_iterator iter = mesh_.find(fields[fieldI]);
+
+            if
+            (
+                iter != objectRegistry::end()
+             && iter()->type()
+             == GeometricField<Type, fvsPatchField, surfaceMesh>::typeName
+            )
+            {
+                sampleAndWrite
+                (
+                    mesh_.lookupObject
+                    <GeometricField<Type, fvsPatchField, surfaceMesh> >
+                    (
+                        fields[fieldI]
+                    )
+                );
+            }
+        }
+    }
+}
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Type>
@@ -185,4 +255,48 @@ Foam::probes::sample(const word& fieldName) const
 }
 
 
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::probes::sample
+(
+    const GeometricField<Type, fvsPatchField, surfaceMesh>& vField
+) const
+{
+    const Type unsetVal(-VGREAT*pTraits<Type>::one);
+
+    tmp<Field<Type> > tValues
+    (
+        new Field<Type>(this->size(), unsetVal)
+    );
+
+    Field<Type>& values = tValues();
+
+    forAll(*this, probeI)
+    {
+        if (faceList_[probeI] >= 0)
+        {
+            values[probeI] = vField[faceList_[probeI]];
+        }
+    }
+
+    Pstream::listCombineGather(values, isNotEqOp<Type>());
+    Pstream::listCombineScatter(values);
+
+    return tValues;
+}
+
+
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::probes::sampleSurfaceFields(const word& fieldName) const
+{
+    return sample
+    (
+        mesh_.lookupObject<GeometricField<Type, fvsPatchField, surfaceMesh> >
+        (
+            fieldName
+        )
+    );
+}
+
 // ************************************************************************* //