ENH: Adding HdotGradH interpolator to particleForces.

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -42,7 +42,7 @@ template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::preEvolve()
 {
     this->dispersion().cacheFields(true);
-    forces_.cacheFields(true);
+    forces_.cacheFields(true, interpolationSchemes_);
 }
 
 
@@ -152,7 +152,7 @@ void Foam::KinematicCloud<ParcelType>::postEvolve()
     }
 
     this->dispersion().cacheFields(false);
-    forces_.cacheFields(false);
+    forces_.cacheFields(false, interpolationSchemes_);
 
     this->postProcessing().post();
 }

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcel.C

@@ -162,6 +162,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
     // Momentum source due to particle forces
     const vector FCoupled = mass*td.cloud().forces().calcCoupled
     (
+        this->position(),
         cellI,
         dt,
         rhoc_,
@@ -173,6 +174,7 @@ const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 
     const vector FNonCoupled = mass*td.cloud().forces().calcNonCoupled
     (
+        this->position(),
         cellI,
         dt,
         rhoc_,

src/lagrangian/intermediate/particleForces/particleForces.C

@@ -30,6 +30,24 @@ License
 #include "mathematicalConstants.H"
 #include "electromagneticConstants.H"
 
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::particleForces::deleteFields()
+{
+    if (gradUPtr_)
+    {
+        delete gradUPtr_;
+        gradUPtr_ = NULL;
+    }
+
+    if (HdotGradHInterPtr_)
+    {
+        delete HdotGradHInterPtr_;
+        HdotGradHInterPtr_ = NULL;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::particleForces::particleForces
@@ -85,7 +103,7 @@ Foam::particleForces::particleForces(const particleForces& f)
 
 Foam::particleForces::~particleForces()
 {
-    cacheFields(false);
+    deleteFields();
 }
 
 
@@ -151,26 +169,46 @@ const Foam::word& Foam::particleForces::HdotGradHName() const
 }
 
 
-void Foam::particleForces::cacheFields(const bool store)
+void Foam::particleForces::cacheFields
+(
+    const bool store,
+    const dictionary& interpolationSchemes
+)
 {
-    if (store && pressureGradient_)
+    if (store)
     {
-        const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
+        if (pressureGradient_)
-        gradUPtr_ = fvc::grad(U).ptr();
+        {
+            const volVectorField& U =
+                mesh_.lookupObject<volVectorField>(UName_);
+
+            gradUPtr_ = fvc::grad(U).ptr();
+        }
+
+        if (paramagnetic_)
+        {
+            const volVectorField& HdotGradH = mesh_.lookupObject<volVectorField>
+            (
+                HdotGradHName_
+            );
+
+            HdotGradHInterPtr_ = interpolation<vector>::New
+            (
+                interpolationSchemes,
+                HdotGradH
+            ).ptr();
+        }
     }
     else
     {
-        if (gradUPtr_)
+        deleteFields();
-        {
-            delete gradUPtr_;
-            gradUPtr_ = NULL;
-        }
     }
 }
 
 
 Foam::vector Foam::particleForces::calcCoupled
 (
+    const vector& position,
     const label cellI,
     const scalar dt,
     const scalar rhoc,
@@ -205,6 +243,7 @@ Foam::vector Foam::particleForces::calcCoupled
 
 Foam::vector Foam::particleForces::calcNonCoupled
 (
+    const vector& position,
     const label cellI,
     const scalar dt,
     const scalar rhoc,
@@ -226,15 +265,12 @@ Foam::vector Foam::particleForces::calcNonCoupled
 
     if (paramagnetic_)
     {
-        const volVectorField& HdotGradH = mesh_.lookupObject<volVectorField>
+        const interpolation<vector>& HdotGradHInter = *HdotGradHInterPtr_;
-        (
-            HdotGradHName_
-        );
 
         accelTot +=
             3.0*constant::electromagnetic::mu0.value()/rho
            *magneticSusceptibility_/(magneticSusceptibility_ + 3)
-           *HdotGradH[cellI];
+           *HdotGradHInter.interpolate(position, cellI);
 
         // force is:
 

src/lagrangian/intermediate/particleForces/particleForces.H

@@ -40,6 +40,7 @@ SourceFiles
 #include "Switch.H"
 #include "vector.H"
 #include "volFieldsFwd.H"
+#include "interpolation.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,6 +70,9 @@ class particleForces
         //- Velocity gradient field
         const volTensorField* gradUPtr_;
 
+        //- HdotGradH interpolator
+        const interpolation<vector>* HdotGradHInterPtr_;
+
 
         // Forces to include in particle motion evaluation
 
@@ -101,6 +105,12 @@ class particleForces
             const word HdotGradHName_;
 
 
+    // Private Member Functions
+
+        //- Delete cached carrier fields
+        void deleteFields();
+
+
 public:
 
     // Constructors
@@ -159,11 +169,16 @@ public:
        // Evaluation
 
             //- Cache carrier fields
-            void cacheFields(const bool store);
+            void cacheFields
+            (
+                const bool store,
+                const dictionary& interpolationSchemes
+            );
 
             //- Calculate action/reaction forces between carrier and particles
             vector calcCoupled
             (
+                const vector& position,
                 const label cellI,
                 const scalar dt,
                 const scalar rhoc,
@@ -176,6 +191,7 @@ public:
             //- Calculate external forces applied to the particles
             vector calcNonCoupled
             (
+                const vector& position,
                 const label cellI,
                 const scalar dt,
                 const scalar rhoc,