BUG: probes - updated to allow sampling from fixed locations for moving mesh cases - mantis #1090

src/sampling/probes/probes.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,6 +28,7 @@ License
 #include "dictionary.H"
 #include "Time.H"
 #include "IOmanip.H"
+#include "mapPolyMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -41,6 +42,11 @@ defineTypeNameAndDebug(probes, 0);
 
 void Foam::probes::findElements(const fvMesh& mesh)
 {
+    if (debug)
+    {
+        Info<< "probes: resetting sample locations" << endl;
+    }
+
     elementList_.clear();
     elementList_.setSize(size());
 
@@ -92,7 +98,7 @@ void Foam::probes::findElements(const fvMesh& mesh)
     {
         const vector& location = operator[](probeI);
         label cellI = elementList_[probeI];
-        label faceI =  faceList_[probeI];
+        label faceI = faceList_[probeI];
 
         // Check at least one processor with cell.
         reduce(cellI, maxOp<label>());
@@ -174,8 +180,8 @@ Foam::label Foam::probes::prepare()
 
         if (debug)
         {
-            Info<< "Probing fields:" << currentFields << nl
-                << "Probing locations:" << *this << nl
+            Info<< "Probing fields: " << currentFields << nl
+                << "Probing locations: " << *this << nl
                 << endl;
         }
 
@@ -267,7 +273,10 @@ Foam::probes::probes
     pointField(0),
     name_(name),
     mesh_(refCast<const fvMesh>(obr)),
-    loadFromFiles_(loadFromFiles)
+    loadFromFiles_(loadFromFiles),
+    fieldSelection_(),
+    fixedLocations_(true),
+    interpolationScheme_("cell")
 {
     read(dict);
 }
@@ -323,10 +332,114 @@ void Foam::probes::read(const dictionary& dict)
     dict.lookup("probeLocations") >> *this;
     dict.lookup("fields") >> fieldSelection_;
 
-    // redetermined all cell locations
+    dict.readIfPresent("fixedLocations", fixedLocations_);
+    if (dict.readIfPresent("interpolationScheme", interpolationScheme_))
+    {
+        if (!fixedLocations_ && interpolationScheme_ != "cell")
+        {
+            WarningIn("void Foam::probes::read(const dictionary&)")
+                << "Only cell interpolation can be applied when "
+                << "not using fixedLocations.  InterpolationScheme "
+                << "entry will be ignored";
+        }
+    }
+
+    // Initialise cells to sample from supplied locations
     findElements(mesh_);
+
     prepare();
 }
 
 
+void Foam::probes::updateMesh(const mapPolyMesh& mpm)
+{
+    if (debug)
+    {
+        Info<< "probes: updateMesh" << endl;
+    }
+
+    if (fixedLocations_)
+    {
+        findElements(mesh_);
+    }
+    else
+    {
+        if (debug)
+        {
+            Info<< "probes: remapping sample locations" << endl;
+        }
+
+        // 1. Update cells
+        {
+            DynamicList<label> elems(elementList_.size());
+
+            const labelList& reverseMap = mpm.reverseCellMap();
+            forAll(elementList_, i)
+            {
+                label cellI = elementList_[i];
+                label newCellI = reverseMap[cellI];
+                if (newCellI == -1)
+                {
+                    // cell removed
+                }
+                else if (newCellI < -1)
+                {
+                    // cell merged
+                    elems.append(-newCellI - 2);
+                }
+                else
+                {
+                    // valid new cell
+                    elems.append(newCellI);
+                }
+            }
+
+            elementList_.transfer(elems);
+        }
+
+        // 2. Update faces
+        {
+            DynamicList<label> elems(faceList_.size());
+
+            const labelList& reverseMap = mpm.reverseFaceMap();
+            forAll(faceList_, i)
+            {
+                label faceI = faceList_[i];
+                label newFaceI = reverseMap[faceI];
+                if (newFaceI == -1)
+                {
+                    // face removed
+                }
+                else if (newFaceI < -1)
+                {
+                    // face merged
+                    elems.append(-newFaceI - 2);
+                }
+                else
+                {
+                    // valid new face
+                    elems.append(newFaceI);
+                }
+            }
+
+            faceList_.transfer(elems);
+        }
+    }
+}
+
+
+void Foam::probes::movePoints(const polyMesh&)
+{
+    if (debug)
+    {
+        Info<< "probes: movePoints" << endl;
+    }
+
+    if (fixedLocations_)
+    {
+        findElements(mesh_);
+    }
+}
+
+
 // ************************************************************************* //

src/sampling/probes/probes.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -84,7 +84,6 @@ protected:
             :
                 DynamicList<word>(0)
             {}
-
         };
 
 
@@ -106,6 +105,16 @@ protected:
             //- Names of fields to probe
             wordReList fieldSelection_;
 
+            //- Fixed locations, default = yes
+            //  Note: set to false for moving mesh calations where locations
+            //        should move with the mesh
+            bool fixedLocations_;
+
+            //- Interpolation scheme name
+            /// Note: only possible when fixedLocations_ is true
+            word interpolationScheme_;
+
+
         // Calculated
 
             //- Categorized scalar/vector/tensor vol fields
@@ -143,13 +152,14 @@ protected:
         //- Classify field types, returns the number of fields
         label classifyFields();
 
-        //- Find cells containing probes
+        //- Find cells and faces containing probes
         virtual void findElements(const fvMesh&);
 
         //- Classify field type and Open/close file streams,
-        //  returns number of fields
+        //  returns number of fields to sample
         label prepare();
 
+
 private:
 
         //- Sample and write a particular volume field
@@ -253,12 +263,10 @@ public:
         virtual void read(const dictionary&);
 
         //- Update for changes of mesh
-        virtual void updateMesh(const mapPolyMesh&)
-        {}
+        virtual void updateMesh(const mapPolyMesh&);
 
         //- Update for changes of mesh
-        virtual void movePoints(const polyMesh&)
-        {}
+        virtual void movePoints(const polyMesh&);
 
         //- Update for changes of mesh due to readUpdate
         virtual void readUpdate(const polyMesh::readUpdateState state)
@@ -285,7 +293,6 @@ public:
         (
             const GeometricField<Type, fvsPatchField, surfaceMesh>&
         ) const;
-
 };
 
 

src/sampling/probes/probesTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,6 +27,7 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 #include "IOmanip.H"
+#include "interpolation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -226,11 +227,36 @@ Foam::probes::sample
 
     Field<Type>& values = tValues();
 
-    forAll(*this, probeI)
+    if (fixedLocations_)
     {
-        if (elementList_[probeI] >= 0)
+        autoPtr<interpolation<Type> > interpolator
+        (
+            interpolation<Type>::New(interpolationScheme_, vField)
+        );
+
+        forAll(*this, probeI)
         {
-            values[probeI] = vField[elementList_[probeI]];
+            if (elementList_[probeI] >= 0)
+            {
+                const vector& position = operator[](probeI);
+
+                values[probeI] = interpolator().interpolate
+                (
+                    position,
+                    elementList_[probeI],
+                    -1
+                );
+            }
+        }
+    }
+    else
+    {
+        forAll(*this, probeI)
+        {
+            if (elementList_[probeI] >= 0)
+            {
+                values[probeI] = vField[elementList_[probeI]];
+            }
         }
     }