NDGrid dual and some cosmetic changes to polynomial patch

src/Classic/Fields/Interpolation/NDGrid.cpp

@@ -224,5 +224,22 @@ Mesh::Iterator NDGrid::getNearest(const double* position) const {
     return Mesh::Iterator(index, this);
 }
 
+
+Mesh* NDGrid::dual() const {
+    std::vector<std::vector<double> > coord(coord_m.size());
+    for (size_t i = 0; i < coord.size(); ++i) {
+        if (coord_m[i].size() <= 1) {
+            throw(OpalException("NDGrid::dual(...)",
+                                "ND Grid must be at least 2x2x...x2 grid"));
+        }
+        coord[i] = std::vector<double>(coord_m[i].size()-1);
+        for (size_t j = 0; j < coord[i].size(); ++j) {
+            coord[i][j] = (coord_m[i][j] + coord_m[i][j+1])/2.;
+        }
+    }
+    return new NDGrid(coord);
+}
+
+
 ////// NDGrid END ////////
 }

src/Classic/Fields/Interpolation/NDGrid.h

@@ -62,7 +62,7 @@ class NDGrid : public Mesh {
     inline NDGrid* clone();
 
     /** Not implemented */
-    inline Mesh* dual() const;
+    Mesh* dual() const;
 
     /** Build a default, empty grid with zero dimension */
     NDGrid();
@@ -348,10 +348,6 @@ NDGrid* NDGrid::clone() {
     return new NDGrid(*this);
 }
 
-Mesh* NDGrid::dual() const {
-    return NULL;
-}
-
 void NDGrid::setCoord(int dimension, int nCoords, double * x) {
     coord_m[dimension] = std::vector<double>(x, x+nCoords);
 }

src/Classic/Fields/Interpolation/PPSolveFactory.cpp

@@ -233,7 +233,7 @@ void PPSolveFactory::getDerivPoints() {
             }
         }
     }
-    
+
     int nPolyCoeffs = SquarePolynomialVector::NumberOfPolynomialCoefficients
                                                     (posDim, smoothingOrder_m);
     for (size_t i = 0; i < derivPoints_m.size(); ++i) {

src/Classic/Fields/Interpolation/PPSolveFactory.h

@@ -68,7 +68,8 @@ class PPSolveFactory {
      *
      *  \param points Set of points on which values are stored. Must be a
      *                rectangular grid. PPSolveFactory takes ownership of
-     *                points (will delete on exit).
+     *                points (will delete on exit). The solve will return a
+     *                wrong answer if the grid does not have regular spacing.
      *  \param values Set of values to which we fit. Must be one value per mesh
      *                point and each value must have the same size.
      *  \param polyPatchOrder The order of the fitted part of the polynomial.

src/Classic/Fields/Interpolation/SquarePolynomialVector.cpp

@@ -25,9 +25,11 @@
  *  POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include <math.h>
+
 #include <iomanip>
 #include <sstream>
-#include <math.h>
+#include <algorithm>
 
 #include "gsl/gsl_sf_gamma.h"
 
@@ -245,4 +247,11 @@ void SquarePolynomialVector::PrintContainer(std::ostream& out, const Container&
   if(pad_at_start) out << strstr2.str();
 }
 
+unsigned int SquarePolynomialVector::PolynomialOrder() const {
+    std::vector<int> index = IndexByPower(_polyCoeffs.num_col(), _pointDim);
+    size_t maxPower = *std::max_element(index.begin(), index.end());
+    return maxPower;
+}
+
+
 }
\ No newline at end of file

src/Classic/Fields/Interpolation/SquarePolynomialVector.h

@@ -63,6 +63,9 @@ namespace interpolation {
  *  n-dimensional square. The distinction should be made with a PolynomialVector
  *  where coefficients include all polynomial coefficients with 
  *  \f$ \Sigma(i_j) <= n. \f$
+ * 
+ *  Nb2: It is quite tricky to get right - mostly an excercise in indexing, but
+ *  that is very easy to get wrong.
  */
 
 class SquarePolynomialVector {

tests/classic_src/Fields/Interpolation/NDGridTest.cpp

@@ -62,22 +62,19 @@ public:
 
     interpolation::NDGrid* grid_m;
     interpolation::NDGrid* grid2_m;
+    OpalTestUtilities::SilenceTest silencer;
 
 private:
 
 };
 
 TEST_F(NDGridTest, DefaultConstructorTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     interpolation::NDGrid grid;
     EXPECT_EQ(grid.begin(), grid.end());
 }
 
 // also tests size
 TEST_F(NDGridTest, Constructor1Test) {
-    OpalTestUtilities::SilenceTest silencer;
-
     int size[] = {5, 6, 7, 8};
     double spacing[] = {1., 2., 3., 4.};
     double min[] = {-1., -2., -3., -4.};
@@ -92,8 +89,6 @@ TEST_F(NDGridTest, Constructor1Test) {
 }
 
 TEST_F(NDGridTest, Constructor2Test) {
-    OpalTestUtilities::SilenceTest silencer;
-
     std::vector<int> size(2);
     size[0] = 2;
     size[1] = 3;
@@ -112,8 +107,6 @@ TEST_F(NDGridTest, Constructor2Test) {
 }
 
 TEST_F(NDGridTest, Constructor3Test) {
-    OpalTestUtilities::SilenceTest silencer;
-
     std::vector< std::vector<double> > gridCoordinates(2);
     gridCoordinates[0] = std::vector<double>(2, 0.);
     gridCoordinates[1] = std::vector<double>(3, 1.);
@@ -129,8 +122,6 @@ TEST_F(NDGridTest, Constructor3Test) {
 }
 
 TEST_F(NDGridTest, CoordTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     std::vector< std::vector<double> > gridCoordinates(2);
     gridCoordinates[0] = std::vector<double>(2, 0.);
     gridCoordinates[1] = std::vector<double>(3, 1.);
@@ -144,8 +135,6 @@ TEST_F(NDGridTest, CoordTest) {
 }
 
 TEST_F(NDGridTest, CoordVectorTest) {  // and newCoordArray
-    OpalTestUtilities::SilenceTest silencer;
-
     std::vector< std::vector<double> > gridCoordinates(2);
     gridCoordinates[0] = std::vector<double>(2, 0.);
     gridCoordinates[1] = std::vector<double>(3, 1.);
@@ -164,8 +153,6 @@ TEST_F(NDGridTest, CoordVectorTest) {  // and newCoordArray
 }
 
 TEST_F(NDGridTest, CoordLowerBoundTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     // first dimension ... 0., 3.;
     int index = -1;
     grid_m->coordLowerBound(-1., 0, index);
@@ -202,8 +189,6 @@ TEST_F(NDGridTest, CoordLowerBoundTest) {
 }
 
 TEST_F(NDGridTest, LowerBoundTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     // first dimension ... 0., 3.;
     // second dimension ...  1., 5., 9.
     std::vector<int> index1(2, -2);
@@ -226,8 +211,6 @@ TEST_F(NDGridTest, LowerBoundTest) {
 }
 
 TEST_F(NDGridTest, MinMaxTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     EXPECT_EQ(grid_m->min(0), 0.);
     EXPECT_EQ(grid_m->min(1), 1.);
     EXPECT_EQ(grid_m->max(0), 3.);
@@ -235,8 +218,6 @@ TEST_F(NDGridTest, MinMaxTest) {
 }
 
 TEST_F(NDGridTest, SetCoordTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     double xNew[] = {5., 10., 12., 15.};
     grid_m->setCoord(0, 4, xNew);
     std::vector<double> xTest = grid_m->coordVector(0);
@@ -247,8 +228,6 @@ TEST_F(NDGridTest, SetCoordTest) {
 }
 
 TEST_F(NDGridTest, BeginEndTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     ASSERT_EQ(grid_m->begin().getState().size(), (unsigned int)2);
     EXPECT_EQ(grid_m->begin().getState()[0], 1);
     EXPECT_EQ(grid_m->begin().getState()[1], 1);
@@ -258,8 +237,6 @@ TEST_F(NDGridTest, BeginEndTest) {
 }
 
 TEST_F(NDGridTest, GetPositionTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     std::vector<double> position(3, -1);
     interpolation::Mesh::Iterator it = grid_m->begin();
     grid_m->getPosition(it, &position[0]);
@@ -274,8 +251,6 @@ TEST_F(NDGridTest, GetPositionTest) {
 }
 
 TEST_F(NDGridTest, GetSetConstantSpacingTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     EXPECT_TRUE(grid_m->getConstantSpacing());
     grid_m->setConstantSpacing(false);
     EXPECT_FALSE(grid_m->getConstantSpacing());
@@ -298,8 +273,6 @@ TEST_F(NDGridTest, GetSetConstantSpacingTest) {
 }
 
 TEST_F(NDGridTest, ToIntegerTest) {
-    OpalTestUtilities::SilenceTest silencer;
-
     interpolation::Mesh::Iterator it = grid_m->begin();
     EXPECT_EQ(grid_m->toInteger(it), 0);
     it[0] = 2;
@@ -308,7 +281,6 @@ TEST_F(NDGridTest, ToIntegerTest) {
 }
 
 TEST_F(NDGridTest, GetNearestTest) {
-    OpalTestUtilities::SilenceTest silencer;
 
     // first dimension ... 0., 3.;
     // second dimension ...  1., 5., 9.
@@ -335,4 +307,20 @@ TEST_F(NDGridTest, GetNearestTest) {
     EXPECT_EQ(it[1], 3);
 }
 
+TEST_F(NDGridTest, DualTest) {
+    using interpolation::NDGrid;
+    NDGrid* new_grid = dynamic_cast<NDGrid*>(grid2_m->dual());
+    std::vector<std::vector<double> > ref = {
+      {1.5},
+      {3.0, 7.5}
+    };
+    ASSERT_EQ(new_grid->getPositionDimension(), int(ref.size()));
+    for (size_t i = 0; i < ref.size(); ++i) {
+        ASSERT_EQ(new_grid->size(i), int(ref[i].size()));
+        for (size_t j = 0; j < ref[i].size(); ++j) {
+            EXPECT_NEAR(new_grid->coordVector(i)[j], ref[i][j], 1e-9) << i << " " << j;
+        }
+    }
+}
+
 } // namespace ndgridtest
\ No newline at end of file