Changeset 2033

Timestamp:

02/11/21 13:20:03 (3 years ago)

Author:

ymipsl

Message:

Add extrapolation to vertical interpolation.
YM

Location:

XIOS/trunk/src

Files:

• config/interpolate_axis_attribute.conf (modified) (1 diff)
• transformation/axis_algorithm_interpolate.cpp (modified) (5 diffs)
• transformation/axis_algorithm_interpolate.hpp (modified) (3 diffs)

XIOS/trunk/src/config/interpolate_axis_attribute.conf

@@ -5 +5 @@
 DECLARE_ATTRIBUTE(StdString, coordinate_src)
 DECLARE_ATTRIBUTE(StdString, coordinate_dst)
+DECLARE_ATTRIBUTE(bool, extrapolate)
 

XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp

@@ -55 +55 @@
   interpAxis->checkValid(axisSource);
   order_ = interpAxis->order.getValue();
+  if (interpAxis->extrapolate.isEmpty()) extrapolate_=false ;
+  else extrapolate_=interpAxis->extrapolate ;
+  
   this->idAuxInputs_.clear();
   if (!interpAxis->coordinate.isEmpty())
@@ -132 +135 @@
   if(!coordinateDST_.empty())
   {
-    int nDomPoint = (*dataAuxInputs[0]).numElements()/numValue ;
     int dst_position_in_data = dataAuxInputs.size()-1;
+    int nDomPoint = (*dataAuxInputs[dst_position_in_data]).numElements()/numValue ;
     for(int ii=0; ii<numValue; ii++)
     {
@@ -192 +195 @@
 
     }
-  }
-  computeWeightedValueAndMapping(interpolatingIndexValues, transPos);
+    else
+    {
+      it=itb ;
+      if (destValue <= *it) 
+      {
+        int numVal=0 ;
+        while(numVal <= order_ && it!=ite)
+        {
+          if (*it != sfmax)
+          { 
+            interpolatingIndexValues[idx+ibegin].push_back(make_pair(indexVec[std::distance(itb, it)],*it));
+            ++numVal ;
+          }
+          ++it ;
+        }
+      }
+      
+      it=ite ;
+      --it ;
+      if (destValue >= *it) 
+      {
+        int numVal=0 ;
+        do
+        {
+          if (*it != sfmax)
+          {
+            interpolatingIndexValues[idx+ibegin].push_back(make_pair(indexVec[std::distance(itb, it)],*it));
+            ++numVal ;
+          }
+          --it ;
+        } while(it!=itb && numVal<=order_) ;
+      }
+    }
+  }
+  
+  computeWeightedValueAndMapping(axisDestValue, interpolatingIndexValues, transPos);
 }
 CATCH
@@ -203 +240 @@
   \param [in] interpolatingIndexValues the necessary axis value to calculate the coeffs
 */
-void CAxisAlgorithmInterpolate::computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos)
+void CAxisAlgorithmInterpolate::computeWeightedValueAndMapping(CArray<double,1>& axisDestValue, const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos)
 TRY
 {
@@ -214 +251 @@
   {
     int globalIndexDest = it->first;
-    double localValue = axisDest_->value(globalIndexDest - ibegin);
+//    double localValue = axisDest_->value(globalIndexDest - ibegin);
+    double localValue = axisDestValue(globalIndexDest - ibegin);
     const std::vector<std::pair<int,double> >& interpVal = it->second;
     int interpSize = interpVal.size();

XIOS/trunk/src/transformation/axis_algorithm_interpolate.hpp

@@ -18 +18 @@
 class CGrid;
 class CInterpolateAxis;
-
+ 
 /*!
   \class CAxisAlgorithmInterpolate
@@ -45 +45 @@
   void computeInterpolantPoint(const std::vector<double>& recvBuff, const std::vector<int>&, 
                                const std::vector<CArray<double,1>* >& dataAuxInputs, int transPos = 0);
-  void computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos = 0);
+  void computeWeightedValueAndMapping(CArray<double,1>& axisDestValue, const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos = 0);
   void fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue,
                        const std::vector<CArray<double,1>* >& dataAuxInputs);
@@ -52 +52 @@
   // Interpolation order
   int order_;
+  bool extrapolate_ ;
   StdString coordinate_;
   StdString coordinateDST_;