Adding radixgeometry subpackage.

cmake/Dependencies.cmake

@@ -3,6 +3,7 @@ TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
      bug  radixbug      SS OPTIONAL
      io   radixio       SS OPTIONAL
      math radixmath     SS OPTIONAL
+     geometry radixgeometry SS OPTIONAL
      plot radixplot     SS OPTIONAL
      widgets radixwidgets SS OPTIONAL
 )

radixgeometry/CMakeLists.txt

+TRIBITS_SUBPACKAGE(geometry)
+
+
+SET(SOURCE
+grid.cc
+)
+
+SET(HEADERS
+grid.hh
+)
+
+#
+# Add library
+TRIBITS_ADD_LIBRARY(radixgeometry
+  SOURCES ${SOURCE}
+  NOINSTALLHEADERS ${HEADERS}
+)
+
+#
+# Add test directory
+TRIBITS_ADD_TEST_DIRECTORIES(tests)
+
+INSTALL(FILES ${HEADERS} DESTINATION "include/radixgeometry/")
+TRIBITS_SUBPACKAGE_POSTPROCESS()

radixgeometry/cmake/Dependencies.cmake

+TRIBITS_PACKAGE_DEFINE_DEPENDENCIES(
+LIB_REQUIRED_PACKAGES radixmath radixbug
+LIB_OPTIONAL_PACKAGES
+TEST_REQUIRED_PACKAGES testframework
+TEST_OPTIONAL_PACKAGES
+LIB_REQUIRED_TPLS
+LIB_OPTIONAL_TPLS
+TEST_REQUIRED_TPLS
+TEST_OPTIONAL_TPLS
+)

radixgeometry/grid.cc

+#include "radixgeometry/grid.hh"
+
+#include <algorithm>
+#include <sstream>
+
+namespace radix
+{
+const std::vector<Real>& Grid::xplanes() const
+{
+    return mXPlanes;
+}
+
+
+const std::vector<Real>& Grid::yplanes() const
+{
+    return mYPlanes;
+}
+
+const std::vector<Real>& Grid::zplanes() const
+{
+    return mZPlanes;
+}
+
+std::vector<Real>& Grid::xplanes()
+{
+    return mXPlanes;
+}
+
+
+std::vector<Real>& Grid::yplanes()
+{
+    return mYPlanes;
+}
+
+std::vector<Real>& Grid::zplanes()
+{
+    return mZPlanes;
+}
+
+
+Grid::Grid() {
+}
+Grid::Grid(const Grid& orig)
+    : mXPlanes(orig.mXPlanes)
+    , mYPlanes(orig.mYPlanes)
+    , mZPlanes(orig.mZPlanes) {
+}
+
+Grid::Grid(Point3D lowerLeftCorner
+         , Real deltaX, Real deltaY, Real deltaZ
+         , int numX, int numY, int numZ)
+        : mXPlanes(numX+1)
+        , mYPlanes(numY+1)
+        , mZPlanes(numZ+1)
+    {
+        radix_check(mXPlanes.size() >= 2);
+        radix_check(mYPlanes.size() >= 2);
+        radix_check(mZPlanes.size() >= 2);
+        //
+        // Initialize first planes with lower left corner
+        mXPlanes[0] = lowerLeftCorner.x;
+        mYPlanes[0] = lowerLeftCorner.y;
+        mZPlanes[0] = lowerLeftCorner.z;
+        //
+        // auto populate grid with delta#
+        for(size_t i = 1; i < mXPlanes.size(); ++i)
+        {
+            mXPlanes[i] = mXPlanes[i-1] + deltaX;
+        }
+        for(size_t i = 1; i < mYPlanes.size(); ++i)
+        {
+            mYPlanes[i] = mYPlanes[i-1] + deltaY;
+        }
+        for(size_t i = 1; i < mZPlanes.size(); ++i)
+        {
+            mZPlanes[i] = mZPlanes[i-1] + deltaZ;
+        }
+    }
+
+    Grid::Grid(const std::vector<Real>&xplanes
+            , const std::vector<Real>&yplanes
+            , const std::vector<Real>&zplanes)
+    : mXPlanes(xplanes)
+    , mYPlanes(yplanes)
+    , mZPlanes(zplanes) {
+        radix_check(xplanes.size() >= 2);
+        radix_check(yplanes.size() >= 2);
+        radix_check(zplanes.size() >= 2);
+    }
+
+    Grid::~Grid() {
+    }
+    void Grid::origin(int xi, int yi, int zi, Point3D &o) const
+    {
+        radix_require(!mXPlanes.empty());
+        radix_require(!mYPlanes.empty());
+        radix_require(!mZPlanes.empty());
+        radix_check(xi < int(mXPlanes.size())-1);
+        radix_check(yi < int(mYPlanes.size())-1);
+        radix_check(zi < int(mZPlanes.size())-1);
+
+        Real x1 = mXPlanes[xi+1];
+        Real x0 = mXPlanes[xi];
+        Real y1 = mYPlanes[yi+1];
+        Real y0 = mYPlanes[yi];
+        Real z1 = mZPlanes[zi+1];
+        Real z0 = mZPlanes[zi];
+        o.x = x0 + (0.5 * (x1 - x0));
+        o.y = y0 + (0.5 * (y1 - y0));
+        o.z = z0 + (0.5 * (z1 - z0));
+    }
+    std::string Grid::toString() const {
+        std::stringstream stream;
+        stream << "Grid x0="
+               << mXPlanes.front()
+               << ",x1=" << mXPlanes.back()
+               << " (" << mXPlanes.size() << ")"
+               << " y0=" << mYPlanes.front()
+               << ",y1=" << mYPlanes.back()
+               << " (" << mYPlanes.size() << ")"
+               << " z0=" << mZPlanes.front()
+               << ",z1=" << mZPlanes.back()
+               << " (" << mZPlanes.size() << ")";
+        return stream.str();
+    }
+
+    bool Grid::inside(Real x, Real y, Real z)const {
+        radix_require(mXPlanes.size() >= 2);
+        radix_require(mYPlanes.size() >= 2);
+        radix_require(mZPlanes.size() >= 2);
+        arrayCoordinates(x,y,z);
+
+        bool in = x + halfKEpsilon >= mXPlanes.front()
+                  && x - halfKEpsilon <= mXPlanes.back()
+                  && y + halfKEpsilon >= mYPlanes.front()
+                  && y - halfKEpsilon <= mYPlanes.back()
+                  && z + halfKEpsilon >= mZPlanes.front()
+                  && z - halfKEpsilon <= mZPlanes.back();
+        //        radix_tagged_line(toString()
+        //                <<", tested against x="<<x<<",y="<<y<<",z="<<z
+        //                <<", inside ? "<<(in?"yes":"no"));
+        return in;
+    }
+
+    /**
+     * Determine the starting array index for the given point
+     * @return true, if and only if the point starts inside the array
+     */
+    bool Grid::start(Real xo, Real yo, Real zo,
+            Real xd, Real yd, Real zd,
+            int& xi, int& yi, int&zi
+            ) const {
+        radix_tagged_line(Point3D(xo,yo,zo)<<" translated into by "<<translation<<"...");
+        arrayCoordinates(xo,yo,zo); // put world coordinates into array
+        radix_tagged_line(Point3D(xo,yo,zo));
+        Real tx = xd < 0 ? xo - halfKEpsilon : (xd > 0 ? xo + halfKEpsilon : xo);
+        Real ty = yd < 0 ? yo - halfKEpsilon : (yd > 0 ? yo + halfKEpsilon : yo);
+        Real tz = zd < 0 ? zo - halfKEpsilon : (zd > 0 ? zo + halfKEpsilon : zo);
+        std::vector<Real>::const_iterator xitr = std::lower_bound(mXPlanes.begin(), mXPlanes.end(), tx);
+        std::vector<Real>::const_iterator yitr = std::lower_bound(mYPlanes.begin(), mYPlanes.end(), ty);
+        std::vector<Real>::const_iterator zitr = std::lower_bound(mZPlanes.begin(), mZPlanes.end(), tz);
+        xi = xitr - mXPlanes.begin();
+        yi = yitr - mYPlanes.begin();
+        zi = zitr - mZPlanes.begin();
+        radix_tagged_line("Start lower_bound indices are "<<xi<<","<<yi<<","<<zi
+                <<", xo="<<xo<<",yo="<<yo<<",zo="<<zo
+                <<", xd="<<xd<<",yd="<<yd<<",zd="<<zd
+                <<", "<<toString());
+        bool foundx = xi != int(mXPlanes.size());
+        bool foundy = yi != int(mYPlanes.size());
+        bool foundz = zi != int(mZPlanes.size());
+
+        bool isxplane = foundx && std::abs(mXPlanes[xi] - xo) < kEpsilon;
+        bool isyplane = foundy && std::abs(mYPlanes[yi] - yo) < kEpsilon;
+        bool iszplane = foundz && std::abs(mZPlanes[zi] - zo) < kEpsilon;
+        //        radix_tagged_line("isx ? "<<(isxplane?"yes":"no")
+        //                        <<", isy ? "<<(isyplane?"yes":"no")
+        //                        <<", isz ? "<<(iszplane?"yes":"no"));
+        //
+        // determine if any of the axis are prior to their respective first plane
+        //
+        bool prior = (!isxplane && xi == 0)
+                     || (!isyplane && yi == 0)
+                     || (!iszplane && zi == 0);
+
+        bool isfirstxplane = isxplane && xi == 0;
+        bool isfirstyplane = isyplane && yi == 0;
+        bool isfirstzplane = iszplane && zi == 0;
+        bool islastxplane = isxplane && xi == int(mXPlanes.size()) - 1;
+        bool islastyplane = isyplane && yi == int(mYPlanes.size()) - 1;
+        bool islastzplane = iszplane && zi == int(mZPlanes.size()) - 1;
+
+
+        bool exiting = (isfirstxplane && xd < 0)// on first xplane heading away
+                       || (isfirstyplane && yd < 0)// on first yplane heading away
+                       || (isfirstzplane && zd < 0)// on first zplane heading away
+                       || (islastxplane && xd > 0)// on last xplane heading away
+                       || (islastyplane && yd > 0)// on last yplane heading away
+                       || (islastzplane && zd > 0)// on last zplane heading away
+                ;
+
+
+        if (!isxplane || xd) {
+            xi -= 1;
+        }
+        if (!isyplane || yd < 0) {
+            yi -= 1;
+        }
+        if (!iszplane || zd < 0) {
+            zi -= 1;
+        }
+//        radix_tagged_line("start - xi="<<xi<<",yi="<<yi<<",zi="<<zi
+//                <<" out of xcount="<<xplanes.size()<<",ycount="<<yplanes.size()
+//                <<",zcount="<<zplanes.size());
+
+        return foundx && foundy && foundz && !prior && !exiting;
+    } // end of start
+
+    /**
+     * Determine the next array index that the given ray will move into
+     * @param xo - the x origin
+     * @param yo - the y origin
+     * @param zo - the z origin
+     * @param xd - the x direction
+     * @param yd - the y direction
+     * @param zd - the z direction
+     * @param xi - the xi index that contains xo, and once complete the new x index
+     * @param yi - the yi index that contains yo, and once complete the new y index
+     * @param zi - the zi index that contains zo, and once complete the new z index
+     * @param t  - the track length to populate with the distance to the next index (if any)
+     * NOTE: It is assumed that the boundary conditions are managed by the caller
+     */
+    void Grid::next(Real xo, Real yo, Real zo,
+            Real xd, Real yd, Real zd,
+            int& xi, int& yi, int& zi,
+            Real& t
+            ) const {
+        // assume boundary conditions managed outside
+        arrayCoordinates(xo,yo,zo); // put world coordinates into array
+        Real tx = maxReal;
+        Real ty = maxReal;
+        Real tz = maxReal;
+        int cx = xi;
+        if (xd > 0) {
+            cx = xi + 1;
+            tx = (mXPlanes[cx] - xo) / xd;
+        } else if (xd < 0) {
+            cx = xi - 1;
+            tx = (mXPlanes[(mXPlanes[xi]<xo?xi:cx)] - xo) / xd;
+        }
+        int cy = yi;
+        if (yd > 0) {
+            cy = yi + 1;
+            ty = (mYPlanes[cy] - yo) / yd;
+        } else if (yd < 0) {
+            cy = yi - 1;
+            ty = (mYPlanes[(mYPlanes[yi]<yo?yi:cy)] - yo) / yd;
+        }
+        int cz = zi;
+        if (zd > 0) {
+            cz = zi + 1;
+            tz = (mZPlanes[cz] - zo) / zd;
+        } else if (zd < 0) {
+            cz = zi - 1;
+            tz = (mZPlanes[(mZPlanes[zi]<zo?zi:cz)] - zo) / zd;
+        }
+//        radix_tagged_line("xo=" << xo << ",yo=" << yo << ",zo=" << zo
+//                << ",xd=" << xd << ",yd=" << yd << ",zd=" << zd
+//                << ",tx=" << tx << ",ty=" << ty << ",tz=" << tz);
+        bool txety = std::abs(tx - ty) < halfKEpsilon;
+        bool txetz = std::abs(tx - tz) < halfKEpsilon;
+        bool tyetz = std::abs(ty - tz) < halfKEpsilon;
+
+        bool txlty = tx < ty;
+        bool txltz = tx < tz;
+        bool tyltz = ty < tz;
+        t = maxReal; // make sure t is prepared for our algorithm below
+
+        if ((txlty || txety) && (txltz || txetz)) {
+            xi = cx;
+            t = tx;
+        }
+        if ((txety || !txlty) && (tyltz || tyetz)) {
+            yi = cy;
+            t = std::min(ty, t);
+        }
+        if (!txltz && !tyltz) {
+            zi = cz;
+            t = std::min(tz, t);
+        }
+    }// end of next
+} // end of namespace

radixgeometry/grid.hh

+#ifndef RADIX_RADIXGEOMETRY_GRID_HH_
+#define RADIX_RADIXGEOMETRY_GRID_HH_
+#include <iostream>
+#include <sstream>
+#include <vector>
+#include <cmath>
+#include <limits>
+
+#include "radixbug/bug.hh"
+#include "radixmath/constants.hh"
+#include "radixmath/point3d.hh"
+#include "radixmath/vector3d.hh"
+
+namespace radix
+{
+class Grid {
+protected:
+    Vector3D mTranslation;
+    std::vector<Real> mXPlanes;
+    std::vector<Real> mYPlanes;
+    std::vector<Real> mZPlanes;
+public:
+    Grid();
+    /**
+         * @brief Grid - initialize the cuboidal array with the given array data
+         * @param arrayData
+         */
+    Grid(const Grid& orig);
+    Grid(Point3D lowerLeftCorner
+         , Real deltaX, Real deltaY, Real deltaZ
+         , int numX, int numY, int numZ);
+    Grid(const std::vector<Real>&xplanes
+         , const std::vector<Real>&yplanes
+         , const std::vector<Real>&zplanes);
+    ~Grid();
+    /**
+         * @brief origin - acquire the origin of the desired array index in array space coordinates
+         * @param xi - x index [0,xsize-1]
+         * @param yi - y index [0,ysize-1]
+         * @param zi - z index [0,zsize-1]
+         * @param o - the point to populate with the origin
+         */
+    void origin(int xi, int yi, int zi, Point3D &o) const;
+
+    inline size_t xcount() const {
+        return mXPlanes.size();
+    }
+
+    inline size_t ycount() const {
+        return mYPlanes.size();
+    }
+
+    inline size_t zcount() const {
+        return mZPlanes.size();
+    }
+    std::string toString() const;
+    /**
+         * Determine if the given x,y,z is inside (or on) the array
+         */
+    bool inside(Real x, Real y, Real z)const;
+    /**
+         * Determine the starting array index for the given point
+         * @param xo - the x origin
+         * @param yo - the y origin
+         * @param zo - the z origin
+         * @param xd - the x direction
+         * @param yd - the y direction
+         * @param zd - the z direction
+         * @param xi - the x index to be populated if xo starts inside the array
+         * @param yi - the y index to be populated if yo starts inside the array
+         * @param zi - the z index to be populated if zo starts inside the array
+         * @return true, if and only if the point starts inside the array
+         */
+    bool start(Real xo, Real yo, Real zo,
+               Real xd, Real yd, Real zd,
+               int& xi, int& yi, int&zi
+               ) const;
+    /**
+         * Determine the next array index that the given ray will move into
+         * @param xo - the x origin
+         * @param yo - the y origin
+         * @param zo - the z origin
+         * @param xd - the x direction
+         * @param yd - the y direction
+         * @param zd - the z direction
+         * @param xi - the xi index that contains xo, and once complete the new x index
+         * @param yi - the yi index that contains yo, and once complete the new y index
+         * @param zi - the zi index that contains zo, and once complete the new z index
+         * @param t  - the track length to populate with the distance to the next index (if any)
+         * NOTE: It is assumed that the boundary conditions are managed by the caller
+         */
+    void next(Real xo, Real yo, Real zo,
+              Real xd, Real yd, Real zd,
+              int& xi, int& yi, int& zi,
+              Real& t) const;
+    void setTranslation(const Vector3D & t){mTranslation = t;}
+
+    void arrayCoordinates(Real& xo, Real& yo, Real& zo)const{
+        xo+=mTranslation.x;
+        yo+=mTranslation.y;
+        zo+=mTranslation.z;
+    }
+
+
+    const std::vector<Real>& xplanes() const;
+    const std::vector<Real>& yplanes() const;
+    const std::vector<Real>& zplanes() const;
+    std::vector<Real>& xplanes();
+    std::vector<Real>& yplanes();
+    std::vector<Real>& zplanes();
+};
+
+} // end of namespace
+#endif /** RADIX_RADIXGEOMETRY_GRID_HH_ **/

radixgeometry/tests/CMakeLists.txt

+INCLUDE(GoogleTest)
+
+ADD_GOOGLE_TEST(tstGrid.cc NP 1)

radixgeometry/tests/tstGrid.cc

+/*
+ * File:   tstGrid.cc
+ * Author: Jordan P. Lefebvre, lefebvrejp@ornl.gov
+ *
+ * Created on November 23, 2013, 3:09 PM
+ */
+
+#include "radixgeometry/grid.hh"
+#include "gtest/gtest.h"
+#include <cstdlib>
+#include <algorithm>
+#include <cmath>
+
+using namespace std;
+using namespace radix;
+TEST(Grid, Constructor)
+{
+    Point3D point(-116.15, 37.16667, 0);
+    Grid grid(point, 0.05, 0.05, 1000, 40, 20, 5);
+    EXPECT_EQ(grid.xcount(), 41);
+    EXPECT_EQ(grid.ycount(), 21);
+    EXPECT_EQ(grid.zcount(), 6);
+
+    EXPECT_FLOAT_EQ(grid.xplanes()[0], point.x);
+    EXPECT_FLOAT_EQ(grid.yplanes()[0], point.y);
+    EXPECT_FLOAT_EQ(grid.zplanes()[0], point.z);
+
+    EXPECT_FLOAT_EQ(grid.xplanes()[40], point.x+(40*.05));
+    EXPECT_FLOAT_EQ(grid.yplanes()[20], point.y+(20*.05));
+    EXPECT_FLOAT_EQ(grid.zplanes()[5], 5000);
+
+    // override the zplanes
+
+    std::vector<Real>& zplanes = grid.zplanes();
+    zplanes[1] = 500;
+    zplanes[2] = 1500;
+    zplanes[3] = 2500;
+    zplanes[4] = 3500;
+    zplanes[5] = 4500;
+
+    EXPECT_EQ(zplanes, grid.zplanes());
+}
+
+TEST(Grid, Uniform)
+{ // uniform grid and copy constructor test
+    Real xo = 1;
+    Real yo = 1;
+    Real zo = 1;
+    Real xd = 1;
+    Real yd = 0;
+    Real zd = 0;
+    int xi = 0;
+    int yi = 0;
+    int zi = 0;
+    Real p [] = {1, 2, 3};
+    std::vector<Real> planes(p, p + 3);
+    Grid array(planes, planes, planes);
+    { // copy constructor
+        Grid cpy(array);
+        EXPECT_EQ(array.xcount(), 3);
+        EXPECT_EQ(array.ycount(), 3);
+        EXPECT_EQ(array.zcount(), 3);
+    }
+    //
+    // test inside for each scenario
+    //
+    xo = 1;
+    yo = 1;
+    zo = 1;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    EXPECT_TRUE(array.inside(xo, yo, zo));
+    zo = 0;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    radix_insist(!array.inside(xo, yo, zo), "Point result not as expected!");
+    yo = 0;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    radix_insist(!array.inside(xo, yo, zo), "Point result not as expected!");
+    xo = 0;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    radix_insist(!array.inside(xo, yo, zo), "Point result not as expected!");
+    xo = 1;
+    yo = 1;
+    zo = 3;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    radix_insist(array.inside(xo, yo, zo), "Point result not as expected!");
+    yo = 3;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    radix_insist(array.inside(xo, yo, zo), "Point result not as expected!");
+    xo = 3;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    radix_insist(array.inside(xo, yo, zo), "Point result not as expected!");
+    zo = 4;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    radix_insist(!array.inside(xo, yo, zo), "Point result not as expected!");
+    yo = 4;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    radix_insist(!array.inside(xo, yo, zo), "Point result not as expected!");
+    xo = 4;
+    radix_tagged_line(xo << "," << yo << "," << zo << " inside ? " << array.inside(xo, yo, zo));
+    radix_insist(!array.inside(xo, yo, zo), "Point result not as expected!");
+    //
+    // start
+    //
+    xo = 1;
+    yo = 1;
+    zo = 1;
+    xd = 1;
+    yd = 0;
+    zd = 0;
+    radix_tagged_line(xo << "," << yo << "," << zo << " start ? " << array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi));
+    radix_insist(array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi), "Start result not as expected!");
+    radix_tagged_line("xi=" << xi << ",yi=" << yi << ",zi=" << zi);
+    radix_insist(xi==0 && yi == 0 && zi == 0, "Indices not as expected!");
+    xo=1.2;
+    radix_tagged_line(xo << "," << yo << "," << zo << " start ? " << array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi));
+    radix_insist(array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi), "Start result not as expected!");
+    radix_tagged_line("xi=" << xi << ",yi=" << yi << ",zi=" << zi);
+    radix_insist(xi==0 && yi == 0 && zi == 0, "Indices not as expected!");
+    xo = 1;
+    zo = 0;
+    radix_tagged_line(xo << "," << yo << "," << zo << " start ? " << array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi));
+    radix_insist(!array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi), "Start result not as expected!");
+    radix_tagged_line("xi=" << xi << ",yi=" << yi << ",zi=" << zi);
+    yo = 0;
+    radix_tagged_line(xo << "," << yo << "," << zo << " start ? " << array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi));
+    radix_insist(!array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi), "Start result not as expected!");
+    radix_tagged_line("xi=" << xi << ",yi=" << yi << ",zi=" << zi);
+    xo = 0;
+    radix_tagged_line(xo << "," << yo << "," << zo << " start ? " << array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi));
+    radix_insist(!array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi), "Start result not as expected!");
+    radix_tagged_line("xi=" << xi << ",yi=" << yi << ",zi=" << zi);
+    xo = 1;
+    yo = 1;
+    zo = 3;
+    radix_tagged_line(xo << "," << yo << "," << zo << " start ? " << array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi));
+    radix_insist(array.start(xo, yo, zo, xd, yd, zd, xi, yi, zi), "Start result not as expected!");
+    radix_tagged_line("xi=" << xi << ",yi=" << yi << ",zi=" << zi);
+    radix_insist(xi == 0 && yi == 0 && zi == 2, "indices not as expected!");