Adding more flexibility to the IdealH user override functor in the AS…

…PH object

src/PYB11/Utilities/SpheralFunctor.py

@@ -12,7 +12,7 @@ def pyinit(self):
 
     @PYB11pure_virtual
     @PYB11const
-    def __call__(self, x="%(argT)s"):
+    def __call__(self, x1="const %(argT)s&"):
         "Required operator() to map %(argT)s --> %(retT)s"
         return "%(retT)s"
 
@@ -26,8 +26,8 @@ def pyinit(self):
     @PYB11pure_virtual
     @PYB11const
     def __call__(self,
-                 x = "%(argT1)s",
-                 y = "%(argT2)s"):
+                 x1 = "const %(argT1)s&",
+                 x2 = "const %(argT2)s&"):
         "Required operator() to map %(argT1)s %(argT2)s --> %(retT)s"
         return "%(retT)s"
 
@@ -41,8 +41,25 @@ def pyinit(self):
     @PYB11pure_virtual
     @PYB11const
     def __call__(self,
-                 x = "%(argT1)s",
-                 y = "%(argT2)s",
-                 z = "%(argT3)s"):
+                 x1 = "const %(argT1)s&",
+                 x2 = "const %(argT2)s&",
+                 x3 = "const %(argT3)s&"):
         "Required operator() to map %(argT1)s %(argT2)s %(argT3)s --> %(retT)s"
         return "%(retT)s"
+
+@PYB11namespace("Spheral::PythonBoundFunctors")
+@PYB11holder("std::shared_ptr")
+@PYB11template("argT1", "argT2", "argT3", "argT4", "retT")
+class Spheral4ArgFunctor:
+    def pyinit(self):
+        return
+
+    @PYB11pure_virtual
+    @PYB11const
+    def __call__(self,
+                 x1 = "const %(argT1)s&",
+                 x2 = "const %(argT2)s&",
+                 x3 = "const %(argT3)s&",
+                 x4 = "const %(argT4)s&"):
+        "Required operator() to map %(argT1)s %(argT2)s %(argT3)s %(argT4)s --> %(retT)s"
+        return "%(retT)s"

src/PYB11/Utilities/Utilities_PYB11.py

@@ -272,7 +272,7 @@ def computeShepardsInterpolation(fieldList = "const FieldList<%(Dimension)s, %(D
 VectorScalarFunctor%(ndim)id = PYB11TemplateClass(SpheralFunctor, template_parameters=("%(Vector)s", "double"))
 VectorVectorFunctor%(ndim)id = PYB11TemplateClass(SpheralFunctor, template_parameters=("%(Vector)s", "%(Vector)s"))
 VectorPairScalarFunctor%(ndim)id = PYB11TemplateClass(SpheralFunctor, template_parameters=("%(Vector)s", "std::pair<double,double>"))
-SizetSizetSymTensorSymTensorFunctor%(ndim)id = PYB11TemplateClass(Spheral3ArgFunctor, template_parameters=("size_t", "size_t", "%(SymTensor)s", "%(SymTensor)s"))
+SizetSizetSymTensorSymTensorSymTensorFunctor%(ndim)id = PYB11TemplateClass(Spheral4ArgFunctor, template_parameters=("size_t", "size_t", "%(SymTensor)s", "%(SymTensor)s", "%(SymTensor)s"))
 
 # boundingVolumes
 boundingBoxVec%(ndim)id = PYB11TemplateFunction(boundingBoxVec, template_parameters="%(Vector)s", pyname="boundingBox")

src/SmoothingScale/ASPHSmoothingScale.cc

@@ -150,15 +150,17 @@ polySecondMoment(const Dim<3>::FacetedVolume& poly,
   return result;
 }
 
+//------------------------------------------------------------------------------
 // A default no-op functor for the Hideal filter
+//------------------------------------------------------------------------------
 template<typename Dimension>
 class HidealPassthrough:
-    public PythonBoundFunctors::Spheral3ArgFunctor<size_t, size_t, typename Dimension::SymTensor, typename Dimension::SymTensor> {
+    public PythonBoundFunctors::Spheral4ArgFunctor<size_t, size_t, typename Dimension::SymTensor, typename Dimension::SymTensor, typename Dimension::SymTensor> {
 public:
   using SymTensor = typename Dimension::SymTensor;
-  HidealPassthrough(): PythonBoundFunctors::Spheral3ArgFunctor<size_t, size_t, typename Dimension::SymTensor, typename Dimension::SymTensor>() {}
+  HidealPassthrough(): PythonBoundFunctors::Spheral4ArgFunctor<size_t, size_t, SymTensor, SymTensor, SymTensor>() {}
   virtual ~HidealPassthrough() {}
-  virtual SymTensor __call__(const size_t nodeListi, const size_t i, const SymTensor Hideal) const override { return Hideal; }
+  virtual SymTensor __call__(const size_t& nodeListi, const size_t& i, const SymTensor& H0, const SymTensor& Hideal) const override { return Hideal; }
 };
 
 }
@@ -670,8 +672,8 @@ finalize(const Scalar time,
         // Build the new H tensor
         // Hi = constructSymTensorWithBoundedDiagonal(fscale*eigenT.eigenValues, hmaxInv, hminInv);
         // Hi.rotationalTransform(eigenT.eigenVectors);
-        Hi = (*mHidealFilterPtr)(k, i, T.Inverse());
-        Hideali = Hi;
+        Hideali = (*mHidealFilterPtr)(k, i, Hi, T.Inverse());
+        Hi = Hideali;     // Since this is the after all our regular state update gotta update the actual H
 
         // // If requested, move toward the cell centroid
         // if (mfHourGlass > 0.0 and surfacePoint(k,i) == 0) {

src/SmoothingScale/ASPHSmoothingScale.hh

@@ -25,7 +25,7 @@ public:
   using Tensor = typename Dimension::Tensor;
   using SymTensor = typename Dimension::SymTensor;
   using FacetedVolume = typename Dimension::FacetedVolume;
-  using HidealFilterType = PythonBoundFunctors::Spheral3ArgFunctor<size_t, size_t, SymTensor, SymTensor>;
+  using HidealFilterType = PythonBoundFunctors::Spheral4ArgFunctor<size_t, size_t, SymTensor, SymTensor, SymTensor>;
 
   // Constructors, destructor.
   ASPHSmoothingScale(const HEvolutionType HUpdate,

src/Utilities/Functors.hh

@@ -19,8 +19,8 @@ class SpheralFunctor {
 public:
   SpheralFunctor() {};
   virtual ~SpheralFunctor() {};
-  virtual retT operator()(const argT x) const { return __call__(x); }
-  virtual retT __call__(const argT x) const = 0;
+  virtual retT operator()(const argT& x1) const { return __call__(x1); }
+  virtual retT __call__(const argT& x1) const = 0;
 };
 
 // retT F(argT1, argT2) 
@@ -29,8 +29,8 @@ class Spheral2ArgFunctor {
 public:
   Spheral2ArgFunctor() {};
   virtual ~Spheral2ArgFunctor() {};
-  virtual retT operator()(const argT1 x, const argT2 y) const { return __call__(x, y); }
-  virtual retT __call__(const argT1 x, const argT2 y) const = 0;
+  virtual retT operator()(const argT1& x1, const argT2& x2) const { return __call__(x1, x2); }
+  virtual retT __call__(const argT1& x1, const argT2& x2) const = 0;
 };
 
 // retT F(argT1, argT2, argT3) 
@@ -39,8 +39,18 @@ class Spheral3ArgFunctor {
 public:
   Spheral3ArgFunctor() {};
   virtual ~Spheral3ArgFunctor() {};
-  virtual retT operator()(const argT1 x, const argT2 y, const argT3 z) const { return __call__(x, y, z); }
-  virtual retT __call__(const argT1 x, const argT2 y, const argT3 z) const = 0;
+  virtual retT operator()(const argT1& x1, const argT2& x2, const argT3& x3) const { return __call__(x1, x2, x3); }
+  virtual retT __call__(const argT1& x1, const argT2& x2, const argT3& x3) const = 0;
+};
+
+// retT F(argT1, argT2, argT3, argT4) 
+template<typename argT1, typename argT2, typename argT3, typename argT4, typename retT>
+class Spheral4ArgFunctor {
+public:
+  Spheral4ArgFunctor() {};
+  virtual ~Spheral4ArgFunctor() {};
+  virtual retT operator()(const argT1& x1, const argT2& x2, const argT3& x3, const argT3& x4) const { return __call__(x1, x2, x3, x4); }
+  virtual retT __call__(const argT1& x1, const argT2& x2, const argT3& x3, const argT4& x4) const = 0;
 };
 
 }