[cig-commits] r21147 - in cs/spatialdata/trunk: libsrc/spatialdb tests/libtests/spatialdb

brad at geodynamics.org brad at geodynamics.org
Thu Dec 13 15:37:27 PST 2012


Author: brad
Date: 2012-12-13 15:37:26 -0800 (Thu, 13 Dec 2012)
New Revision: 21147

Modified:
   cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridAscii.cc
   cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.cc
   cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.hh
   cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.icc
   cs/spatialdata/trunk/tests/libtests/spatialdb/TestSimpleGridDB.cc
Log:
Fixed bugs in using gridded spatial data with lower dimension topology.

Modified: cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridAscii.cc
===================================================================
--- cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridAscii.cc	2012-12-13 22:38:03 UTC (rev 21146)
+++ cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridAscii.cc	2012-12-13 23:37:26 UTC (rev 21147)
@@ -282,6 +282,7 @@
   db->_cs->initialize();
 } // _readHeader
 
+#include <iostream>
 // ----------------------------------------------------------------------
 // Read data values.
 void
@@ -365,23 +366,8 @@
     for (int iDim=0; iDim < spaceDim; ++iDim) {
       buffer >> coords[iDim];
     } // for
-    
-    int indexX = 0;
-    int indexY = 0;
-    int indexZ = 0;
-    if (spaceDim > 2) {
-      indexX = int(std::floor(db->_search(coords[0], db->_x, numX)+0.5));
-      indexY = int(std::floor(db->_search(coords[1], db->_y, numY)+0.5));
-      indexZ = int(std::floor(db->_search(coords[2], db->_z, numZ)+0.5));
-    } else if (spaceDim > 1) {
-      indexX = int(std::floor(db->_search(coords[0], db->_x, numX)+0.5));
-      indexY = int(std::floor(db->_search(coords[1], db->_y, numY)+0.5));
-    } else {
-      assert(1 == spaceDim);
-      indexX = int(std::floor(db->_search(coords[0], db->_x, numX)+0.5));
-    } // if
-    
-    const int indexData = db->_dataIndex(indexX, indexY, indexZ);
+
+    const int indexData = db->_dataIndex(coords, spaceDim);
     for (int iVal=0; iVal < db->_numValues; ++iVal) {
       buffer >> db->_data[indexData+iVal];
     } // for
@@ -489,7 +475,7 @@
   for (int iZ=0; iZ < numZ; ++iZ) {
     for (int iY=0; iY < numY; ++iY) {
       for (int iX=0; iX < numX; ++iX) {
-	const int iD = db._dataIndex(iX, iY, iZ);
+	const int iD = db._dataIndex(iX, numX, iY, numY, iZ, numZ);
 	fileout 
 	  << std::setw(14) << db._x[iX]
 	  << std::setw(14) << db._y[iY]

Modified: cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.cc
===================================================================
--- cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.cc	2012-12-13 22:38:03 UTC (rev 21146)
+++ cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.cc	2012-12-13 23:37:26 UTC (rev 21147)
@@ -202,78 +202,47 @@
   int queryFlag = 0;
   const int spaceDim = _spaceDim;
 
-  double indexX = 0.0;
-  double indexY = 0.0;
-  double indexZ = 0.0;
-  int numX = 0;
-  int numY = 0;
-  int numZ = 0;
+  double index0 = 0.0;
+  double index1 = 0.0;
+  double index2 = 0.0;
+  int size0 = 0;
+  int size1 = 0;
+  int size2 = 0;
   if (3 == spaceDim) {
-    indexX = _search(_xyz[0], _x, _numX);
-    indexY = _search(_xyz[1], _y, _numY);
-    indexZ = _search(_xyz[2], _z, _numZ);
-    numX = _numX;
-    numY = _numY;
-    numZ = _numZ;
-    if (2 == _dataDim) {
-      if (1 == _numX) {
-	indexX = indexY;
-	numX = _numY;
-	indexY = indexZ;
-	numY = _numZ;
-	indexZ = 0;
-      } else if (1 == numY) {
-	indexY = indexZ;
-	numY = _numZ;
-	indexZ = 0;
-      } // if/else
-    } else if (1 == _dataDim) {
-      if (_numY > 1) {
-	indexX = indexY;
-	numX = _numY;
-	indexY = 0;
-      } else if (_numZ > 1) {
-	indexX = indexZ;
-	numX = _numZ;
-	indexZ = 0;
-      } // if
-    } // if
+    index0 = _search(_xyz[0], _x, _numX);
+    index1 = _search(_xyz[1], _y, _numY);
+    index2 = _search(_xyz[2], _z, _numZ);
+    _reindex3d(&index0, &size0, &index1, &size1, &index2, &size2);
   } else if (2 == spaceDim) {
-    indexX = _search(_xyz[0], _x, _numX);
-    indexY = _search(_xyz[1], _y, _numY);
-    numX = _numX;
-    numY = _numY;
-    if (1 == _dataDim && 1 == _numX) {
-      indexX = indexY;
-      numX = _numY;
-      indexY = 0;
-    } // if
+    index0 = _search(_xyz[0], _x, _numX);
+    index1 = _search(_xyz[1], _y, _numY);
+    _reindex2d(&index0, &size0, &index1, &size1);
   } else { // else
     assert(1 == spaceDim);
-    indexX = _search(_xyz[0], _x, _numX);
-    numX = _numX;
+    index0 = _search(_xyz[0], _x, _numX);
+    size0 = _numX;
   } // if/else
   
   switch (_queryType) {
   case LINEAR : 
-    if (indexX < 0.0 || (indexX > 0 && indexX > numX-1.0) ||
-	indexY < 0.0 || (indexY > 0 && indexY > numY-1.0) ||
-        indexZ < 0.0 || (indexZ > 0 && indexZ > numZ-1.0)) {
-  queryFlag = 1;
-  return queryFlag;
+    if (index0 < 0.0 || (index0 > 0 && index0 > size0-1.0) ||
+	index1 < 0.0 || (index1 > 0 && index1 > size1-1.0) ||
+        index2 < 0.0 || (index2 > 0 && index2 > size2-1.0)) {
+      queryFlag = 1;
+      return queryFlag;
     } // if
 
     switch (_dataDim) {
     case 1: {
-      _interpolate1D(vals, numVals, indexX, numX);
+      _interpolate1D(vals, numVals, index0, size0);
       break;
     } // case 1
     case 2: {
-      _interpolate2D(vals, numVals, indexX, numX, indexY, numY);
+      _interpolate2D(vals, numVals, index0, size0, index1, size1);
       break;
     } // case 2
     case 3 : {
-      _interpolate3D(vals, numVals, indexX, indexY, indexZ);
+      _interpolate3D(vals, numVals, index0, index1, index2);
       break;
     } // case 3
     default :
@@ -282,24 +251,25 @@
     } // switch
     break;
   case NEAREST : {
-    indexX = std::min(indexX, numX-1.0);
-    indexX = std::max(indexX, 0.0);
-    indexY = std::min(indexY, numY-1.0);
-    indexY = std::max(indexY, 0.0);
-    indexZ = std::min(indexZ, numZ-1.0);
-    indexZ = std::max(indexZ, 0.0);
-    const int indexNearestX = int(std::floor(indexX+0.5));
-    const int indexNearestY = int(std::floor(indexY+0.5));
-    const int indexNearestZ = int(std::floor(indexZ+0.5));
-    const int indexData = _dataIndex(indexNearestX, indexNearestY, indexNearestZ);
+    index0 = std::min(index0, size0-1.0);
+    index0 = std::max(index0, 0.0);
+    index1 = std::min(index1, size1-1.0);
+    index1 = std::max(index1, 0.0);
+    index2 = std::min(index2, size2-1.0);
+    index2 = std::max(index2, 0.0);
+    const int indexNearest0 = int(std::floor(index0+0.5));
+    const int indexNearest1 = int(std::floor(index1+0.5));
+    const int indexNearest2 = int(std::floor(index2+0.5));
+    const int indexData = _dataIndex(indexNearest0, size0, indexNearest1, size1, indexNearest2, size2);
+
     for (int iVal=0; iVal < querySize; ++iVal) {
       vals[iVal] = _data[indexData+_queryVals[iVal]];
 #if 0 // DEBUGGING
     std::cout << "val["<<iVal<<"]: " << vals[iVal]
 	      << ", indexData: " << indexData
-	      << ", indexX: " << indexX
-	      << ", indexY: " << indexY
-	      << ", indexZ: " << indexZ
+	      << ", index0: " << index0
+	      << ", index1: " << index1
+	      << ", index2: " << index2
 	      << std::endl;
 #endif
     } // for
@@ -447,26 +417,10 @@
 
   assert(_data);
   for (int iLoc=0; iLoc < numLocs; ++iLoc) {
-    int indexX = 0;
-    int indexY = 0;
-    int indexZ = 0;
-    const int ii = iLoc*spaceDim;
-    if (spaceDim > 2) {
-      indexX = int(std::floor(_search(coords[ii+0], _x, _numX)+0.5));
-      indexY = int(std::floor(_search(coords[ii+1], _y, _numY)+0.5));
-      indexZ = int(std::floor(_search(coords[ii+2], _z, _numZ)+0.5));
-    } else if (spaceDim > 1) {
-      indexX = int(std::floor(_search(coords[ii+0], _x, _numX)+0.5));
-      indexY = int(std::floor(_search(coords[ii+1], _y, _numY)+0.5));
-    } else {
-      assert(1 == spaceDim);
-      indexX = int(std::floor(_search(coords[ii+0], _x, _numX)+0.5));
-    } // if
-
-    const int iD = _dataIndex(indexX, indexY, indexZ);
+    const int indexData = _dataIndex(&coords[iLoc*spaceDim], spaceDim);
     const int jj = iLoc*numValues;
     for (int iV=0; iV < numValues; ++iV) {
-      _data[iD+iV] = values[jj+iV];
+      _data[indexData+iV] = values[jj+iV];
     } // for
   } // for
 } // data
@@ -591,7 +545,7 @@
 double
 spatialdata::spatialdb::SimpleGridDB::_search(const double target,
 					      const double* vals,
-					      const int nvals)
+					      const int nvals) const
 { // _search
   if (1 == nvals) {
     return 0.0;
@@ -643,8 +597,8 @@
   assert(0 <= indexX0 && indexX0 < numX);
   assert(0 <= indexX1 && indexX1 < numX);
 
-  const int index000 = _dataIndex(indexX0, 0, 0);
-  const int index100 = _dataIndex(indexX1, 0, 0);
+  const int index000 = _dataIndex(indexX0, numX, 0, 0, 0, 0);
+  const int index100 = _dataIndex(indexX1, numX, 0, 0, 0, 0);
 
   const double wt000 = wtX0;
   const double wt100 = wtX1;
@@ -692,10 +646,10 @@
   assert(0 <= indexY0 && indexY0 < numY);
   assert(0 <= indexY1 && indexY1 < numY);
   
-  const int index000 = _dataIndex(indexX0, indexY0, 0);
-  const int index010 = _dataIndex(indexX0, indexY1, 0);
-  const int index100 = _dataIndex(indexX1, indexY0, 0);
-  const int index110 = _dataIndex(indexX1, indexY1, 0);
+  const int index000 = _dataIndex(indexX0, numX, indexY0, numY, 0, 0);
+  const int index010 = _dataIndex(indexX0, numX, indexY1, numY, 0, 0);
+  const int index100 = _dataIndex(indexX1, numX, indexY0, numY, 0, 0);
+  const int index110 = _dataIndex(indexX1, numX, indexY1, numY, 0, 0);
 
   const double wt000 = wtX0 * wtY0;
   const double wt010 = wtX0 * wtY1;
@@ -760,14 +714,14 @@
   assert(0 <= indexZ0 && indexZ0 < numZ);
   assert(0 <= indexZ1 && indexZ1 < numZ);
 
-  const int index000 = _dataIndex(indexX0, indexY0, indexZ0);
-  const int index001 = _dataIndex(indexX0, indexY0, indexZ1);
-  const int index010 = _dataIndex(indexX0, indexY1, indexZ0);
-  const int index011 = _dataIndex(indexX0, indexY1, indexZ1);
-  const int index100 = _dataIndex(indexX1, indexY0, indexZ0);
-  const int index101 = _dataIndex(indexX1, indexY0, indexZ1);
-  const int index110 = _dataIndex(indexX1, indexY1, indexZ0);
-  const int index111 = _dataIndex(indexX1, indexY1, indexZ1);
+  const int index000 = _dataIndex(indexX0, numX, indexY0, numY, indexZ0, numZ);
+  const int index001 = _dataIndex(indexX0, numX, indexY0, numY, indexZ1, numZ);
+  const int index010 = _dataIndex(indexX0, numX, indexY1, numY, indexZ0, numZ);
+  const int index011 = _dataIndex(indexX0, numX, indexY1, numY, indexZ1, numZ);
+  const int index100 = _dataIndex(indexX1, numX, indexY0, numY, indexZ0, numZ);
+  const int index101 = _dataIndex(indexX1, numX, indexY0, numY, indexZ1, numZ);
+  const int index110 = _dataIndex(indexX1, numX, indexY1, numY, indexZ0, numZ);
+  const int index111 = _dataIndex(indexX1, numX, indexY1, numY, indexZ1, numZ);
 
   const double wt000 = wtX0 * wtY0 * wtZ0;
   const double wt001 = wtX0 * wtY0 * wtZ1;
@@ -807,4 +761,117 @@
 } // _interpolate3D
 
 
+// ----------------------------------------------------------------------
+// Adjust indices to account for optimizations for lower dimension
+// distribution.
+void
+spatialdata::spatialdb::SimpleGridDB::_reindex2d(double* const index0,
+						 int* const size0,
+						 double* const index1,
+						 int* const size1) const
+{ // _reindex2d
+  assert(index0);
+  assert(index1);
+  assert(size0);
+  assert(size1);
+
+  *size0 = _numX;
+  *size1 = _numY;
+  if (1 == _dataDim && 1 == _numX) {
+    *index0 = *index1;
+    *size0 = *size1;
+    *index1 = 0;
+    *size1 = 1;
+  } // if
+} // _reindex2d
+
+
+// ----------------------------------------------------------------------
+// Adjust indices to account for optimizations for lower dimension
+// distribution.
+void
+spatialdata::spatialdb::SimpleGridDB::_reindex3d(double* const index0,
+						 int* const size0,
+						 double* const index1,
+						 int* const size1,
+						 double* const index2,
+						 int* const size2) const
+{ // _reindex3d
+  assert(index0);
+  assert(index1);
+  assert(index2);
+  assert(size0);
+  assert(size1);
+  assert(size2);
+
+  *size0 = _numX;
+  *size1 = _numY;
+  *size2 = _numZ;
+  if (2 == _dataDim) {
+    if (1 == _numX) {
+      *index0 = *index1;
+      *size0 = *size1;
+      *index1 = *index2;
+      *size1 = *size2;
+      *index2 = 0;
+      *size2 = 1;
+    } else if (1 == _numY) {
+      *index1 = *index2;
+      *size1 = *size2;
+      *index2 = 0;
+      *size2 = 1;
+    } // if/else
+  } else if (1 == _dataDim) {
+    if (_numY > 1) {
+      *index0 = *index1;
+      *size0 = *size1;
+      *index1 = 0;
+      *size1 = 1;
+      *index2 = 0;
+      *size2 = 1;
+    } else if (_numZ > 1) {
+      *index0 = *index2;
+      *size0 = *size2;
+      *index1 = 0;
+      *size1 = 1;
+      *index2 = 0;
+      *size2 = 1;
+    } // if
+  } // if/else
+} // _reindex3d
+
+
+// ----------------------------------------------------------------------
+// Get index into data array.
+int
+spatialdata::spatialdb::SimpleGridDB::_dataIndex(const double* const coords,
+						 const int spaceDim) const
+{ // _dataIndex
+  assert(coords);
+
+  double index0 = 0;
+  double index1 = 0;
+  double index2 = 0;
+  int size0 = 0;
+  int size1 = 0;
+  int size2 = 0;
+  if (spaceDim > 2) {
+    index0 = std::floor(_search(coords[0], _x, _numX)+0.5);
+    index1 = std::floor(_search(coords[1], _y, _numY)+0.5);
+    index2 = std::floor(_search(coords[2], _z, _numZ)+0.5);
+    _reindex3d(&index0, &size0, &index1, &size1, &index2, &size2);
+  } else if (spaceDim > 1) {
+    index0 = std::floor(_search(coords[0], _x, _numX)+0.5);
+    index1 = std::floor(_search(coords[1], _y, _numY)+0.5);
+    _reindex2d(&index0, &size0, &index1, &size1);
+  } else {
+    assert(1 == spaceDim);
+    index0 = std::floor(_search(coords[0], _x, _numX)+0.5);
+  } // if
+
+  const int indexData = _dataIndex(int(index0), size0, int(index1), size1, int(index2), size2);
+  return indexData;
+} // _dataIndex
+
+
 // End of file

Modified: cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.hh
===================================================================
--- cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.hh	2012-12-13 22:38:03 UTC (rev 21146)
+++ cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.hh	2012-12-13 23:37:26 UTC (rev 21147)
@@ -196,7 +196,7 @@
    */
   double _search(const double target,
 		 const double* vals,
-		 const int nvals);
+		 const int nvals) const;
 
   /** Interpolate in 1-D to get values at target location defined by
    * indices.
@@ -246,18 +246,64 @@
 		      const double indexY,
 		      const double indexZ) const;
 
+  /** Adjust indices to account for optimizations for lower dimension
+   *  distribution.
+   *
+   * @param index0 Adjusted index for coordinate 0.
+   * @param size0 Adjusted size for coordinate 0.
+   * @param index1 Adjusted index for dimension 1.
+   * @param size1 Adjusted size for dimension 1.
+   */
+  void _reindex2d(double* const index0,
+		  int* const size0,
+		  double* const index1,
+		  int* const size1) const;
+  
+  /** Adjust indices to account for optimizations for lower dimension
+   *  distribution.
+   *
+   * @param index0 Adjusted index for coordinate 0.
+   * @param size0 Adjusted size for coordinate 0.
+   * @param index1 Adjusted index for dimension 1.
+   * @param size1 Adjusted size for dimension 1.
+   * @param index2 Adjusted index for dimension 2.
+   * @param size2 Adjusted size for dimension 2.
+   */
+  void _reindex3d(double* const index0,
+		  int* const size0,
+		  double* const index1,
+		  int* const size1,
+		  double* const index2,
+		  int* const size2) const;
+
   /** Get index into data array.
    *
-   * @param indexX Index along x dimension.
-   * @param indexY Index along y dimension.
-   * @param indexZ Index along z dimension.
+   * @param index0 Adjusted index for coordinate 0.
+   * @param size0 Adjusted size for coordinate 0.
+   * @param index1 Adjusted index for dimension 1.
+   * @param size1 Adjusted size for dimension 1.
+   * @param index2 Adjusted index for dimension 2.
+   * @param size2 Adjusted size for dimension 2.
    *
    * @returns Index into data array.
    */
-  int _dataIndex(const int indexX,
-		 const int indexY,
-		 const int indexZ) const;
+  int _dataIndex(const int index0,
+		 const int size0,
+		 const int index1,
+		 const int size1,
+		 const int index2,
+		 const int size2) const;
 
+  /** Get index into data array.
+   *
+   * @param Coordinates of point.
+   * @param spaceDim Number of coordinate dimensions.
+   *
+   * @returns Index into data array.
+   */
+  int _dataIndex(const double* const coords,
+		 const int spaceDim) const;
+
 // PRIVATE MEMBERS //////////////////////////////////////////////////////
 private :
   

Modified: cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.icc
===================================================================
--- cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.icc	2012-12-13 22:38:03 UTC (rev 21146)
+++ cs/spatialdata/trunk/libsrc/spatialdb/SimpleGridDB.icc	2012-12-13 23:37:26 UTC (rev 21147)
@@ -24,12 +24,15 @@
 // Get index into data array.
 inline
 int
-spatialdata::spatialdb::SimpleGridDB::_dataIndex(const int indexX,
-						 const int indexY,
-						 const int indexZ) const
+spatialdata::spatialdb::SimpleGridDB::_dataIndex(const int index0,
+						 const int size0,
+						 const int index1,
+						 const int size1,
+						 const int index2,
+						 const int size2) const
 { // _dataIndex
   // Order points so indexing works in any dimension.
-  const int locIndex = indexZ*_numY*_numX + indexY*_numX + indexX;
+  const int locIndex = index2*size1*size0 + index1*size0 + index0;
   return locIndex*_numValues;
 } // _dataIndex
 

Modified: cs/spatialdata/trunk/tests/libtests/spatialdb/TestSimpleGridDB.cc
===================================================================
--- cs/spatialdata/trunk/tests/libtests/spatialdb/TestSimpleGridDB.cc	2012-12-13 22:38:03 UTC (rev 21146)
+++ cs/spatialdata/trunk/tests/libtests/spatialdb/TestSimpleGridDB.cc	2012-12-13 23:37:26 UTC (rev 21147)
@@ -173,14 +173,14 @@
   db._numZ = 5;
   db._numValues = 10;
 
-  CPPUNIT_ASSERT_EQUAL(0, db._dataIndex(0, 0, 0));
-  CPPUNIT_ASSERT_EQUAL(1*3*4*10, db._dataIndex(0, 0, 1));
+  CPPUNIT_ASSERT_EQUAL(0, db._dataIndex(0, db._numX, 0, db._numY, 0, db._numZ));
+  CPPUNIT_ASSERT_EQUAL(1*3*4*10, db._dataIndex(0, db._numX, 0, db._numY, 1, db._numZ));
 
-  CPPUNIT_ASSERT_EQUAL(1*4*10, db._dataIndex(0, 1, 0));
-  CPPUNIT_ASSERT_EQUAL(4*3*4*10 + 1*4*10, db._dataIndex(0, 1, 4));
+  CPPUNIT_ASSERT_EQUAL(1*4*10, db._dataIndex(0, db._numX, 1, db._numY, 0, db._numZ));
+  CPPUNIT_ASSERT_EQUAL(4*3*4*10 + 1*4*10, db._dataIndex(0, db._numX, 1, db._numY, 4, db._numZ));
 
-  CPPUNIT_ASSERT_EQUAL(1*10, db._dataIndex(1, 0, 0));
-  CPPUNIT_ASSERT_EQUAL(3*4*3*10 + 1*4*10 + 2*10, db._dataIndex(2, 1, 3));
+  CPPUNIT_ASSERT_EQUAL(1*10, db._dataIndex(1, db._numX, 0, db._numY, 0, db._numZ));
+  CPPUNIT_ASSERT_EQUAL(3*4*3*10 + 1*4*10 + 2*10, db._dataIndex(2, db._numX, 1, db._numY, 3, db._numZ));
 } // testDataIndex
 
 // ----------------------------------------------------------------------
@@ -282,14 +282,14 @@
 void
 spatialdata::spatialdb::TestSimpleGridDB::testIO(void)
 { // testIO
-  const int numX = 3;
+  const int numX = 1;
   const int numY = 2;
   const int numZ = 3;
   const int spaceDim = 3;
   const int numVals = 2;
-  const int dataDim = 3;
+  const int dataDim = 2;
 
-  const double x[numX] = { -2.0, 0.0, 3.0 };
+  const double x[numX] = { -2.0 };
   const double y[numY] = { 0.0, 1.0 };
   const double z[numZ] = { -2.0, -1.0, 2.0 };
   
@@ -300,41 +300,17 @@
     -2.0,  1.0, -1.0,
     -2.0,  0.0,  2.0,
     -2.0,  1.0,  2.0,
-     0.0,  0.0, -2.0,
-     0.0,  1.0, -2.0,
-     0.0,  0.0, -1.0,
-     0.0,  1.0, -1.0,
-     0.0,  0.0,  2.0,
-     0.0,  1.0,  2.0,
-     3.0,  0.0, -2.0,
-     3.0,  1.0, -2.0,
-     3.0,  0.0, -1.0,
-     3.0,  1.0, -1.0,
-     3.0,  0.0,  2.0,
-     3.0,  1.0,  2.0,
   };
   const double data[numX*numY*numZ*numVals] = {
     6.6,  3.4,
     5.5,  6.7,
     2.3,  4.1,
     5.7,  2.0,
-    6.3,  6.7,
+    6.3,  6.9,
     3.4,  6.4,
-    7.2,  6.8,
-    5.7,  8.2,
-    3.4,  9.8,
-    5.7,  2.3,
-    9.4,  8.5,
-    7.2,  9.3,
-    4.8,  7.5,
-    9.2,  8.3,
-    5.8,  8.5,
-    4.7,  8.9,
-    7.8,  6.2,
-    2.9,  8.3,
   };
-  const char* names[] = { "One", "Two" };
-  const char* units[] = { "m", "m" };
+  const char* names[numVals] = { "One", "Two" };
+  const char* units[numVals] = { "m", "m" };
 
   geocoords::CSCart csOut;
   csOut.initialize();
@@ -371,18 +347,44 @@
     CPPUNIT_ASSERT_EQUAL(std::string(units[iVal]), dbIn._units[iVal]);
   } // for
 
+  // Check to make sure values were read in correctly
   CPPUNIT_ASSERT(dbIn._data);
   const double tolerance = 1.0e-06;
   for (int iX=0, i=0; iX < numX; ++iX) {
     for (int iZ=0; iZ < numZ; ++iZ) {
       for (int iY=0; iY < numY; ++iY) {
-	const int iD = dbIn._dataIndex(iX, iY, iZ);
+	const int iD = dbIn._dataIndex(iX, numX, iY, numY, iZ, numZ);
 	for (int iVal=0; iVal < numVals; ++iVal, ++i) {
 	  CPPUNIT_ASSERT_DOUBLES_EQUAL(dbIn._data[iD+iVal]/data[i], 1.0, tolerance);
 	} // for
       } // for
     } // for
   } // for
+
+  // Perform simple nearest query to ensure consistency of read/query
+  dbIn.queryVals(names, numVals);
+  const int numLocs = 3;
+  const double points[numLocs*spaceDim] = {
+    -2.0, 1.0, -2.0,
+    -5.0, 0.0,  2.0,
+    +6.0, 1.0, -1.0,
+  };
+  const double dataE[numLocs*numVals] = {
+    5.5, 6.7,
+    6.3, 6.9,
+    5.7, 2.0,
+  };
+  const int errE[numLocs] = { 0, 0 };
+  
+  for (int iLoc=0; iLoc < numLocs; ++iLoc) {
+    double data[numVals];
+    int err = dbIn.query(data, numVals, &points[iLoc*spaceDim], spaceDim, &csOut);
+    CPPUNIT_ASSERT_EQUAL(errE[iLoc], err);
+    for (int iVal=0; iVal < numVals; ++iVal) {
+      CPPUNIT_ASSERT_DOUBLES_EQUAL(data[iVal]/dataE[iLoc*numVals+iVal], 1.0, tolerance);
+    } // for
+  } // for
+
 } // testIO
 
 



More information about the CIG-COMMITS mailing list