Merge branch 'WrapperFix' into devel

Author: Martin D. Weinberg

expui/BasisFactory.H

@@ -259,6 +259,10 @@ namespace BasisClasses
       if (Naccel > 0) pseudo = currentAccel(time);
     }
 
+    //! Get the field label vector
+    std::vector<std::string> getFieldLabels(void)
+    { return getFieldLabels(coordinates); }
+
   };
 
   using BasisPtr = std::shared_ptr<Basis>;

expui/Coefficients.cc

@@ -2792,7 +2792,7 @@ namespace CoefClasses
 
     for (int t=0; t<ntim; t++) {
       auto & cof = *(coefs[roundTime(times[t])]->coefs);
-      for (int i=0; i<4; i++) {
+      for (int i=0; i<Nfld; i++) {
 	for (int l=0; l<(Lmax+2)*(Lmax+1)/2; l++) {
 	  for (int n=0; n<Nmax; n++) {
 	    ret(i, l, n, t) = cof(i, l, n);

expui/FieldBasis.H

@@ -98,9 +98,13 @@ namespace BasisClasses
     virtual std::vector<double>
     crt_eval(double x, double y, double z);
 
+    //! Get the field labels
+    std::vector<std::string> getFieldLabels(const Coord ctype)
+    { return fieldLabels; }
+  
   public:
 
-    //! Constructor from YAML node
+   //! Constructor from YAML node
     FieldBasis(const YAML::Node& conf,
 	       const std::string& name="FieldBasis") : Basis(conf, name)
     { configure(); }
@@ -158,10 +162,6 @@ namespace BasisClasses
     {
     }
 
-    //! Get the field labels
-    std::vector<std::string> getFieldLabels(const Coord ctype)
-    { return fieldLabels; }
-  
     //! Return current maximum harmonic order in expansion
     int getLmax() { return lmax; }
 

expui/FieldBasis.cc

@@ -731,8 +731,8 @@ namespace BasisClasses
     double r   = sqrt(R*R + z*z);
 
     double vr  = (u*x + v*y + w*z)/r;
-    double vt  = (u*z*x + v*z*y - w*R)/R/r;
-    double vp  = (u*y - v*x)/R;
+    double vt  = (u*z*x + v*z*y - w*R*R)/R/r;
+    double vp  = (v*x - u*y)/R;
 
     return {vr, vt, vp, vr*vr, vt*vt, vp*vp};
   }

pyEXP/BasisWrappers.cc

@@ -928,8 +928,25 @@ void BasisFactoryClasses(py::module &m)
          Returns
          -------
          None
-         )");
+         )")
+    .def("getFieldLabels",
+	 [](BasisClasses::Basis& A)
+	 {
+	   return A.getFieldLabels();
+	 },
+	 R"(
+         Provide the field labels for the basis functions
 
+         Parameters
+         ----------
+         None
+
+         Returns
+         -------
+         list: str
+           list of basis function labels
+      )"
+    );
 
   py::class_<BasisClasses::BiorthBasis, std::shared_ptr<BasisClasses::BiorthBasis>, PyBiorthBasis, BasisClasses::Basis>
     (m, "BiorthBasis")
@@ -1279,9 +1296,9 @@ void BasisFactoryClasses(py::module &m)
     // orthoCheck is not in the base class and needs to have different
     // parameters depending on the basis type.  Here, the quadrature
     // is determined by the scale of the meridional grid.
-    .def("orthoCheck", [](BasisClasses::Cylindrical& A)
+    .def("orthoCheck", [](BasisClasses::Cylindrical& A, int knots)
 	 {
-	   return A.orthoCheck();
+	   return A.orthoCheck(knots);
 	 },
 	R"(
         Check orthgonality of basis functions by quadrature
@@ -1298,7 +1315,7 @@ void BasisFactoryClasses(py::module &m)
         -------
         list(numpy.ndarray)
 	    list of numpy.ndarrays from [0, ... , Mmax]
-        )")
+        )", py::arg("knots")=400)
     .def_static("cacheInfo", [](std::string cachefile)
     {
       return BasisClasses::Cylindrical::cacheInfo(cachefile);
@@ -2044,11 +2061,11 @@ void BasisFactoryClasses(py::module &m)
       // orthoCheck is not in the base class and needs to have
       // different parameters depending on the basis type.  Here the
       // user can and will often need to specify a quadrature value.
-      .def("orthoCheck", [](BasisClasses::FieldBasis& A)
+    .def("orthoCheck", [](BasisClasses::FieldBasis& A)
       {
 	return A.orthoCheck();
       },
-	R"(
+      R"(
         Check orthgonality of basis functions by quadrature
 
         Inner-product matrix of orthogonal functions
@@ -2062,10 +2079,10 @@ void BasisFactoryClasses(py::module &m)
         numpy.ndarray
 	    orthogonality matrix
         )"
-	);
+      );
 
   py::class_<BasisClasses::VelocityBasis, std::shared_ptr<BasisClasses::VelocityBasis>, BasisClasses::FieldBasis>(m, "VelocityBasis")
-      .def(py::init<const std::string&>(),
+    .def(py::init<const std::string&>(),
 	 R"(
          Create a orthogonal velocity-field basis
 

pyEXP/TensorToArray.H

@@ -11,7 +11,7 @@ py::array_t<T> make_ndarray3(Eigen::Tensor<T, 3>& mat)
   // Check rank
   if (dims.size() != 3) {
     std::ostringstream sout;
-    sout << "make_ndarray: tensor rank must be 3, found " << dims.size();
+    sout << "make_ndarray3: tensor rank must be 3, found " << dims.size();
     throw std::runtime_error(sout.str());
   }
 
@@ -37,10 +37,17 @@ py::array_t<T> make_ndarray4(Eigen::Tensor<T, 4>& mat)
   // Check rank
   if (dims.size() != 4) {
     std::ostringstream sout;
-    sout << "make_ndarray: tensor rank must be 4, found " << dims.size();
+    sout << "make_ndarray4: tensor rank must be 4, found " << dims.size();
     throw std::runtime_error(sout.str());
   }
 
+  // Sanity check
+  for (int i=0; i<mat.size(); i++) {
+    if (isnan(std::abs(mat.data()[i]))) {
+      throw std::runtime_error("make_ndarray4: NaN encountered");
+    }
+  }
+
   // Make the memory mapping
   return py::array_t<T>
     (
@@ -107,11 +114,11 @@ Eigen::Tensor<T, 4> make_tensor4(py::array_t<T> array)
 
   // Build result tensor with col-major ordering
   Eigen::Tensor<T, 4> tensor(shape[0], shape[1], shape[2], shape[3]);
-  for (int i=0, l=0; i < shape[0]; i++) {
+  for (int i=0, c=0; i < shape[0]; i++) {
     for (int j=0; j < shape[1]; j++) {
       for (int k=0; k < shape[2]; k++) {
-	for (int l=0; l < shape[3]; k++) {
-	  tensor(i, j, k, l) = data[l++];
+	for (int l=0; l < shape[3]; l++, c++) {
+	  tensor(i, j, k, l) = data[c];
 	}
       }
     }