Changeset 937

Timestamp:

09/22/16 10:59:04 (7 years ago)

Author:

mhnguyen

Message:

Smalll improvements for interpolation axis

+) Strictly check order of interpolation
+) Change default order to 1
+) Add precision check

Test
+) On Curie
+) Correct

Location:

XIOS/trunk/src

Files:

• node/interpolate_axis.cpp (modified) (2 diffs)
• transformation/axis_algorithm_interpolate.cpp (modified) (4 diffs)

XIOS/trunk/src/node/interpolate_axis.cpp

@@ -40 +40 @@
   void CInterpolateAxis::checkValid(CAxis* axisSrc)
   {
-    if (this->order.isEmpty()) this->order.setValue(2);
+    if (this->order.isEmpty()) this->order.setValue(1);
     int order = this->order.getValue();
     if (order >= axisSrc->n_glo.getValue())
@@ -46 +46 @@
       ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)",
              << "Order of interpolation is greater than global size of axis source"
+             << "Size of axis source " <<axisSrc->getId() << " is " << axisSrc->n_glo.getValue()  << std::endl
+             << "Order of interpolation is " << order );
+    }
+
+    if (order < 1)
+    {
+      ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)",
+             << "Order of interpolation is smaller than 1"
              << "Size of axis source " <<axisSrc->getId() << " is " << axisSrc->n_glo.getValue()  << std::endl
              << "Order of interpolation is " << order );

XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp

@@ -97 +97 @@
 {
   std::vector<double>::const_iterator itb = axisValue.begin(), ite = axisValue.end();
-  std::vector<double>::const_iterator itLowerBound, itUpperBound, it, iteRange;
+  std::vector<double>::const_iterator itLowerBound, itUpperBound, it, iteRange, itfirst, itsecond;
   const double sfmax = NumTraits<double>::sfmax();
+  const double precision = NumTraits<double>::dummy_precision();
 
   int ibegin = axisDest_->begin.getValue();
@@ -116 +117 @@
 
     if ((ite == itLowerBound) || (ite == itUpperBound)) outOfRange = true;
-    // If the value is not in the range, that means we'll do extra-polation
-//    if (ite == itLowerBound) // extra-polation
-//    {
-//      itLowerBound = itb;
-//      itUpperBound = itb + order_+1;
-//    }
-//    else if (ite == itUpperBound) // extra-polation
-//    {
-//      itLowerBound = itUpperBound - order_-1;
-//    }
-//else
+
     // We don't do extrapolation FOR NOW, maybe in the future
     if (!outOfRange)
     {
       if ((itLowerBound == itUpperBound) && (itb != itLowerBound)) --itLowerBound;
+      double distanceToLower = destValue - *itLowerBound;
+      double distanceToUpper = *itUpperBound - destValue;
       int order = (order_ + 1) - 2;
-      bool down = true;
+      bool down = (distanceToLower < distanceToUpper) ? true : false;
       for (int k = 0; k < order; ++k)
       {
@@ -138 +131 @@
         {
           --itLowerBound;
-          down = false;
+          distanceToLower = destValue - *itLowerBound;
+          down = (distanceToLower < distanceToUpper) ? true : false;
           continue;
         }
@@ -144 +138 @@
         {
           ++itUpperBound;
-          down = true;
-        }
-      }
-
+          distanceToUpper = *itUpperBound - destValue;
+          down = (distanceToLower < distanceToUpper) ? true : false;
+
+        }
+      }
 
       iteRange = (ite == itUpperBound) ? itUpperBound : itUpperBound + 1;
-      for (it = itLowerBound; it != iteRange; ++it)
-      {
+      itsecond = it = itLowerBound; ++itsecond;
+      while (it < iteRange)
+      {
+        while (((*itsecond -*it) < precision) && itsecond < ite)
+        { ++itsecond; ++it; }
         int index = std::distance(itb, it);
         interpolatingIndexValues[idx+ibegin].push_back(make_pair(indexVec[index],*it));
-      }
+        ++it; ++itsecond;
+      }
+
     }
   }