Changeset 918 for XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp

Timestamp:

08/12/16 19:57:50 (8 years ago)

Author:

mhnguyen

Message:

Modifying vertical interpolation

+) Only process interpolation, the extrapolation will be ignored (value will be unidentified (masked))
+) Make sure even unidenfitified from one transformation can be known in the next transformation

Test
+) On Curie
+) Vertical interpolation: Correct
+) Vertical interpolation + horizontal interpolation: Correct

File:

• XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp (modified) (4 diffs)

XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp

@@ -69 +69 @@
 {
   std::vector<double>::const_iterator itb = axisValue.begin(), ite = axisValue.end();
-  std::vector<double>::const_iterator itLowerBound, itUpperBound, it;
+  std::vector<double>::const_iterator itLowerBound, itUpperBound, it, iteRange;
   const double sfmax = NumTraits<double>::sfmax();
 
@@ -79 +79 @@
   for (int idx = 0; idx < numValue; ++idx)
   {
+    bool outOfRange = false;
     double destValue = axisDestValue(idx);
+    if (destValue < *itb) outOfRange = true;
+
     itLowerBound = std::lower_bound(itb, ite, destValue);
     itUpperBound = std::upper_bound(itb, ite, destValue);
     if ((ite != itUpperBound) && (sfmax == *itUpperBound)) itUpperBound = ite;
 
+    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
-    {
-      if (itb != itLowerBound) --itLowerBound;
-      if (ite != itUpperBound) ++itUpperBound;
+//    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;
       int order = (order_ + 1) - 2;
       bool down = true;
@@ -114 +119 @@
         }
       }
-    }
-
-    for (it = itLowerBound; it != itUpperBound; ++it)
-    {
-      int index = std::distance(itb, it);
-      interpolatingIndexValues[idx+ibegin].push_back(make_pair(indexVec[index],*it));
+
+
+      iteRange = (ite == itUpperBound) ? itUpperBound : itUpperBound + 1;
+      for (it = itLowerBound; it != iteRange; ++it)
+      {
+        int index = std::distance(itb, it);
+        interpolatingIndexValues[idx+ibegin].push_back(make_pair(indexVec[index],*it));
+      }
     }
   }
@@ -163 +170 @@
     }
   }
+  if (!transPosition_.empty() && this->transformationPosition_[transPos].empty())
+    (this->transformationPosition_[transPos])[0] = transPosition_[transPos];
+
 }