Changeset 2122 for XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate_coordinate.cpp

Timestamp:

04/01/21 11:13:01 (3 years ago)

Author:

ymipsl

Message:

Improve axis interpolate feature :
introduce : coordinate_src and coordinate_dst feature :

• level -> pressure
• pressure -> pressure

YM

File:

• XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate_coordinate.cpp (modified) (6 diffs)

XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate_coordinate.cpp

@@ -55 +55 @@
   vector<string> CAxisAlgorithmInterpolateCoordinate::getAuxFieldId(void)
   {
-    if (hasCoordinate_) return vector<string>(1,coordinate_) ;
+    if (hasCoordinateSrc_ && hasCoordinateDest_) return {coordinateSrc_,coordinateDest_} ;
+    else if (hasCoordinateSrc_) return {coordinateSrc_} ;
+    else if (hasCoordinateDest_) return {coordinateDest_} ;
     else return vector<string>() ;
+  }
+  
+  bool CAxisAlgorithmInterpolateCoordinate::transformAuxField(int pos)
+  {
+    if (pos==0)
+    {
+      if (hasCoordinateSrc_) return true ;
+      else if(hasCoordinateDest_) return false ;
+    }
+    if (pos==1) return false ;
   }
 
@@ -69 +81 @@
     if (!interpAxis->coordinate.isEmpty())
     {
-      coordinate_ = interpAxis->coordinate.getValue();
+      coordinateSrc_ = interpAxis->coordinate.getValue();
       hasCoordinate_=true ;
-    }
-  
-    ngloSrc_=axisSource->n_glo ;
+      hasCoordinateSrc_=true ;
+    }
+
+    if (!interpAxis->coordinate_src.isEmpty())
+    {
+      coordinateSrc_ = interpAxis->coordinate_src.getValue();
+      hasCoordinate_=true ;
+      hasCoordinateSrc_=true ;
+    }
+
+    if (!interpAxis->coordinate_dst.isEmpty())
+    {
+      coordinateDest_ = interpAxis->coordinate_dst.getValue();
+      hasCoordinate_=true ;
+      hasCoordinateDest_=true ;
+    }
+
+    ngloSrc_=axisSource->n_glo ; 
     nDest_ =  axisDest_-> getLocalView(CElementView::WORKFLOW)->getSize() ;
-    CArray<double,1> coord ;
-    CLocalConnector destConnector(axisDest_->getLocalView(CElementView::FULL), axisDest_->getLocalView(CElementView::WORKFLOW)) ;
-    destConnector.computeConnector() ;
-    destConnector.transfer(axisDest_->value, coord) ;
-    destCoordinate_ = vector<double>(coord.dataFirst(), coord.dataFirst()+nDest_) ;
     
+    if (!hasCoordinateDest_)
+    {
+      CArray<double,1> coord ;
+      CLocalConnector destConnector(axisDest_->getLocalView(CElementView::FULL), axisDest_->getLocalView(CElementView::WORKFLOW)) ;
+      destConnector.computeConnector() ;
+      destConnector.transfer(axisDest_->value, coord) ;
+      destCoordinate_ = vector<double>(coord.dataFirst(), coord.dataFirst()+nDest_) ;
+    }
+
     CArray<size_t,1> globalIndex(ngloSrc_) ;
     for(int i=0;i<ngloSrc_;i++) 
@@ -93 +124 @@
 
     this->computeAlgorithm(axisSource->getLocalView(CElementView::WORKFLOW), axisSourceGlo.getView(CElementView::FULL)) ;
+
+    if (!hasCoordinateSrc_)
+    {
+      CArray<double,1> coord ;
+      CArray<double,1> coordGlo ;
+      CLocalConnector srcConnector(axisSrc_->getLocalView(CElementView::FULL), axisSrc_->getLocalView(CElementView::WORKFLOW)) ;
+      srcConnector.computeConnector() ;
+      srcConnector.transfer(axisSrc_->value, coord) ; // full view value -> workflow value
+      transferTransformConnector_ -> transfer(coord, coordGlo) ; // workflow view -> full global view
+      srcCoordinate_ = vector<double>(coordGlo.dataFirst(), coordGlo.dataFirst()+ngloSrc_) ;
+    }
   }
   CATCH
@@ -98 +140 @@
   CTransformFilter* CAxisAlgorithmInterpolateCoordinate::createTransformFilter(CGarbageCollector& gc, CGridAlgorithm* algo, bool detectMissingValues, double defaultValue)
   {
-    return new CTransformFilter(gc, 2, algo, detectMissingValues, defaultValue) ;  
+    if (hasCoordinateSrc_ && hasCoordinateDest_) return new CTransformFilter(gc, 3, algo, detectMissingValues, defaultValue) ;
+    else return new CTransformFilter(gc, 2, algo, detectMissingValues, defaultValue) ;  
   }
 
@@ -105 +148 @@
   {
     CArray<double,1> dataInTmp;
-    CArray<double,1> auxDataInTmp ;
+    CArray<double,1> auxDataInSrc ;
     transferTransformConnector_ -> transfer(dimBefore, dimAfter, dataIn, dataInTmp) ;
-    transferTransformConnector_ -> transfer(dimBefore, dimAfter, auxDataIn[0], auxDataInTmp) ;
-  
+    if (hasCoordinateSrc_)  transferTransformConnector_ -> transfer(dimBefore, dimAfter, auxDataIn[0], auxDataInSrc) ;
+
+   
     dataOut.resize(dimBefore*dimAfter*nDest_) ;
-    const double* pressure = auxDataInTmp.dataFirst() ;
+    const double* pressureSrc ;
+    const double* pressureDest ;
+    if (hasCoordinateSrc_)  pressureSrc = auxDataInSrc.dataFirst() ;
+    if (hasCoordinateDest_ && hasCoordinateSrc_) pressureDest = auxDataIn[1].dataFirst() ;
+    else if (hasCoordinateDest_ && !hasCoordinateSrc_ ) pressureDest = auxDataIn[0].dataFirst() ;
+    
     const double* in = dataInTmp.dataFirst() ;
     double* out = dataOut.dataFirst() ;
 
-    size_t sliceSrc = dimAfter*ngloSrc_ ;
+    size_t sliceSrc  = dimAfter*ngloSrc_ ;
     size_t sliceDest = dimAfter*nDest_ ;
     vector<double> srcCoordinate(ngloSrc_) ;
@@ -123 +172 @@
     std::vector<int> destIndex(nDest_);
 
-    size_t nsrc ;
+    size_t nsrc, nDest ;
+
     for(size_t j=0, posJsrc=0, posJdest=0 ;  j<dimBefore ; j++, posJsrc+=sliceSrc, posJdest+=sliceDest )
       for(size_t k=0, posKsrc=posJsrc, posKdest=posJdest ; k<dimAfter ; k++, posKsrc++,posKdest++)
       {
-        nsrc=0 ;
-        for(size_t i=0, posIsrc=posKsrc, posIdest=posKdest ; i<ngloSrc_ ; i++, posIsrc+=dimAfter, posIdest+=dimAfter)
-        {
-          if ( !( std::isnan(pressure[posIsrc]) || std::isnan(in[posIsrc]) ) )
-          {
-            srcCoordinate[nsrc]=pressure[posIsrc] ;
-            srcValue[nsrc] = in[posIsrc] ;
-            nsrc++ ;
-          }
-        }
-        destCoordinate = destCoordinate_ ;    
-        computeInterp(nsrc, srcCoordinate, srcValue, srcIndex, nDest_, destCoordinate, destValue, destIndex) ;
-
-        for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)  out[posIdest] = destValue[i] ;
+        if (hasCoordinateSrc_)
+        {
+          nsrc=0 ;
+          for(size_t i=0, posIsrc=posKsrc, posIdest=posKdest ; i<ngloSrc_ ; i++, posIsrc+=dimAfter, posIdest+=dimAfter)
+          {
+            if ( !( std::isnan(pressureSrc[posIsrc]) || std::isnan(in[posIsrc]) ) )
+            {
+              srcCoordinate[nsrc]=pressureSrc[posIsrc] ;
+              srcValue[nsrc] = in[posIsrc] ;
+              nsrc++ ;
+            }
+          }
+        }
+        else
+        {
+          nsrc=0 ;
+          for(size_t i=0, posIsrc=posKsrc ; i<ngloSrc_ ; i++, posIsrc+=dimAfter)
+          {
+            if ( !( std::isnan(srcCoordinate_[i]) || std::isnan(in[posIsrc]) ) )
+            {
+              srcCoordinate[nsrc]=srcCoordinate_[i] ;
+              srcValue[nsrc] = in[posIsrc] ;
+              nsrc++ ;
+            }
+          }
+        }  
+
+        if (hasCoordinateDest_)
+        {
+          nDest=0 ;
+          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)
+          {
+            if ( !( std::isnan(pressureDest[posIdest]) ) )
+            {
+              destCoordinate[nDest]=pressureDest[posIdest] ;
+              nDest++ ;
+            }
+          }
+        }
+        else
+        {
+          nDest=0 ;
+          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)
+          {
+            if ( !( std::isnan(destCoordinate[i]) ) )
+            {
+              destCoordinate[nDest]=destCoordinate_[i] ;
+              nDest++ ;
+            }
+          }
+        }
+  
+        computeInterp(nsrc, srcCoordinate, srcValue, srcIndex, nDest, destCoordinate, destValue, destIndex) ;
+        
+        if (hasCoordinateDest_)
+        {
+          nDest=0 ;
+          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)  
+          {
+            if ( !( std::isnan(pressureDest[posIdest]) ) )  
+            {
+              out[posIdest] = destValue[nDest] ;
+              nDest++ ;
+            }
+            else out[posIdest] = std::numeric_limits<double>::quiet_NaN() ;
+          }
+        }
+        else
+        {
+          nDest=0 ;
+          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)  
+          {
+            if ( !( std::isnan(destCoordinate[i]) ) )  
+            {
+              out[posIdest] = destValue[nDest] ;
+              nDest++ ;
+            }
+            else out[posIdest] = std::numeric_limits<double>::quiet_NaN() ;
+          }
+        }
       }
-
   }