Merge pull request #124 from EXP-code/fixAccelFunc

Fix accel func and add center manipulation tools. See thread for full discussion.

Author: michael-petersen

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features
 
 project(
   EXP
-  VERSION "7.8.1"
+  VERSION "7.8.2"
   HOMEPAGE_URL https://github.com/EXP-code/EXP
   LANGUAGES C CXX Fortran)
 

expui/BasisFactory.H

@@ -263,6 +263,8 @@ namespace BasisClasses
     std::vector<std::string> getFieldLabels(void)
     { return getFieldLabels(coordinates); }
 
+    //! Get the basis expansion center
+    std::vector<double> getCenter() { return coefctr; }
   };
 
   using BasisPtr = std::shared_ptr<Basis>;

expui/BiorthBasis.cc

@@ -3685,11 +3685,17 @@ namespace BasisClasses
     //
     auto basis = std::get<0>(mod);
 
+    // Get expansion center
+    //
+    auto ctr = basis->getCenter();
+
     // Get fields
     //
     int rows = accel.rows();
     for (int n=0; n<rows; n++) {
-      auto v = basis->getFields(ps(n, 0), ps(n, 1), ps(n, 2));
+      auto v = basis->getFields(ps(n, 0) - ctr[0],
+				ps(n, 1) - ctr[1],
+				ps(n, 2) - ctr[2]);
       // First 6 fields are density and potential, follewed by acceleration
       for (int k=0; k<3; k++) accel(n, k) += v[6+k] - basis->pseudo(k);
     }

expui/CoefStruct.H

@@ -27,7 +27,7 @@ namespace CoefClasses
     Eigen::VectorXcd store;
 
     //! Center data
-    std::vector<double> ctr;
+    std::vector<double> ctr= {0.0, 0.0, 0.0};
 
     //! Destructor
     virtual ~CoefStruct() {}
@@ -62,6 +62,26 @@ namespace CoefClasses
     //! Read-only access to coefficient data (no copy)
     Eigen::Ref<const Eigen::VectorXcd> getCoefs() { return store; }
 
+    //! Set new center data (no copy)
+    void setCenter(std::vector<double>& STORE)
+    {
+      if (STORE.size() != ctr.size())
+        throw std::invalid_argument("CoefStruct::setCenter: center vector size does not match");
+      ctr = STORE;
+    }
+
+    //! Read-only access to center (no copy)
+    std::vector<double> getCenter() { return ctr; }
+
+    //! Set coefficient time (no copy)
+    void setTime(double& STORE)
+    {
+      time = STORE;
+    }
+
+    //! Read-only access to center (no copy)
+    double getTime() { return time; }
+
   };
 
   //! Specialization of CoefStruct for spheres

pyEXP/CoefWrappers.cc

@@ -712,11 +712,86 @@ void CoefficientClasses(py::module &m) {
                   str
                       geometry type
                   )")
-    .def_readwrite("time", &CoefStruct::time,
+    .def_readonly("time", &CoefStruct::time,
 		   R"(
                    float
                        data's time stamp
                    )")
+    .def_readonly("center", &CoefStruct::ctr,
+    R"(
+                float
+                    data's center value
+                )") 
+    .def("getCoefTime", &CoefStruct::getTime,
+        R"(
+        Read-only access to the coefficient time
+
+        Returns
+        -------
+        numpy.ndarray
+            vector of center data
+
+        See also
+        --------
+        setCoefTime : read-write access to the coefficient time
+        )")
+    .def("setCoefTime",
+        static_cast<void (CoefStruct::*)(double&)>(&CoefStruct::setTime),
+        py::arg("tval"),
+        R"(
+        Set the coefficient time
+
+        Parameters
+        ----------
+        tval  : float
+                time value
+
+        Returns
+        -------
+        None
+
+        Notes
+        -----
+
+        See also
+        --------
+        getCenter : read-only access to coefficient time
+        )")
+    .def("getCoefCenter", &CoefStruct::getCenter,
+          R"(
+          Read-only access to the center data
+  
+          Returns
+          -------
+          numpy.ndarray
+              vector of center data
+  
+          See also
+          --------
+          setCenter : read-write access to the center data
+          )")
+    .def("setCoefCenter",
+          static_cast<void (CoefStruct::*)(std::vector<double>&)>(&CoefStruct::setCenter),
+          py::arg("mat"),
+          R"(
+          Set the center vector
+  
+          Parameters
+          ----------
+          mat  : numpy.ndarray
+                center vector
+  
+          Returns
+          -------
+          None
+  
+          Notes
+          -----
+  
+          See also
+          --------
+          getCenter : read-only access to center data
+          )")
     .def("getCoefs", &CoefStruct::getCoefs,
         R"(
         Read-only access to the underlying data store