source: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/axis_algorithm/axis_algorithm_interpolate_coordinate.cpp @ 2328

Last change on this file since 2328 was 2328, checked in by jderouillat, 2 years ago

Set extrapolation as an option (set by default to true) in axis interpolation using the coordinate attribute
  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 13.6 KB
Line 
1/*!
2   \file axis_algorithm_interpolate.cpp
3   \author Ha NGUYEN
4   \since 23 June 2015
5   \date 02 Jul 2015
6
7   \brief Algorithm for interpolation on an axis.
8 */
9#include "axis_algorithm_interpolate_coordinate.hpp"
10#include "axis.hpp"
11#include "interpolate_axis.hpp"
12#include <algorithm>
13#include <numeric>
14#include "context.hpp"
15#include "context_client.hpp"
16#include "utils.hpp"
17#include "grid.hpp"
18#include "grid_transformation_factory_impl.hpp"
19#include "distribution_client.hpp"
20#include "transform_filter.hpp"
21#include "timer.hpp"
22
23namespace xios
24{
25  shared_ptr<CGenericAlgorithmTransformation> CAxisAlgorithmInterpolateCoordinate::create(bool isSource, CGrid* gridDst, CGrid* gridSrc,
26                                                                     CTransformation<CAxis>* transformation,
27                                                                     int elementPositionInGrid,
28                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
29                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
30                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
31                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
32                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
33                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
34  TRY
35  {
36    std::vector<CAxis*> axisListDestP = gridDst->getAxis();
37    std::vector<CAxis*> axisListSrcP  = gridSrc->getAxis();
38
39    CInterpolateAxis* interpolateAxis = dynamic_cast<CInterpolateAxis*> (transformation);
40    int axisDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
41    int axisSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
42
43    return make_shared<CAxisAlgorithmInterpolateCoordinate>(isSource, axisListDestP[axisDstIndex], axisListSrcP[axisSrcIndex], interpolateAxis);
44  }
45  CATCH
46
47  bool CAxisAlgorithmInterpolateCoordinate::dummyRegistered_ = CAxisAlgorithmInterpolateCoordinate::registerTrans();
48  bool CAxisAlgorithmInterpolateCoordinate::registerTrans()
49  TRY
50  {
51    return CGridTransformationFactory<CAxis>::registerTransformation(TRANS_INTERPOLATE_AXIS, create);
52  }
53  CATCH
54
55  vector<string> CAxisAlgorithmInterpolateCoordinate::getAuxFieldId(void)
56  {
57    if (hasCoordinateSrc_ && hasCoordinateDest_) return {coordinateSrc_,coordinateDest_} ;
58    else if (hasCoordinateSrc_) return {coordinateSrc_} ;
59    else if (hasCoordinateDest_) return {coordinateDest_} ;
60    else return vector<string>() ;
61  }
62 
63  bool CAxisAlgorithmInterpolateCoordinate::transformAuxField(int pos)
64  {
65    if (pos==0)
66    {
67      if (hasCoordinateSrc_) return true ;
68      else if(hasCoordinateDest_) return false ;
69    }
70    if (pos==1) return false ;
71
72    MISSING_RETURN( "bool CAxisAlgorithmInterpolateCoordinate::transformAuxField(int pos)" );
73    return false;
74  }
75
76  CAxisAlgorithmInterpolateCoordinate::CAxisAlgorithmInterpolateCoordinate(bool isSource, CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis)
77  : CAlgorithmTransformationTransfer(isSource), axisSrc_(axisSource), axisDest_(axisDestination)
78  TRY
79  {
80    interpAxis->checkValid(axisSource);
81    axisDestination->checkAttributes() ;
82
83    order_ = interpAxis->order.getValue();
84   
85    if (interpAxis->extrapolate.isEmpty()) extrapolate_=true ;
86    else extrapolate_=interpAxis->extrapolate ;
87
88    if (!interpAxis->coordinate.isEmpty())
89    {
90      coordinateSrc_ = interpAxis->coordinate.getValue();
91      hasCoordinate_=true ;
92      hasCoordinateSrc_=true ;
93    }
94
95    if (!interpAxis->coordinate_src.isEmpty())
96    {
97      coordinateSrc_ = interpAxis->coordinate_src.getValue();
98      hasCoordinate_=true ;
99      hasCoordinateSrc_=true ;
100    }
101
102    if (!interpAxis->coordinate_dst.isEmpty())
103    {
104      coordinateDest_ = interpAxis->coordinate_dst.getValue();
105      hasCoordinate_=true ;
106      hasCoordinateDest_=true ;
107    }
108
109    ngloSrc_=axisSource->n_glo ; 
110    nDest_ =  axisDest_-> getLocalView(CElementView::WORKFLOW)->getSize() ;
111   
112    if (!hasCoordinateDest_)
113    {
114      CArray<double,1> coord ;
115      auto destConnector = make_shared<CLocalConnector>(axisDest_->getLocalView(CElementView::FULL), axisDest_->getLocalView(CElementView::WORKFLOW)) ;
116      destConnector->computeConnector() ;
117      destConnector->transfer(axisDest_->value, coord) ;
118      destCoordinate_ = vector<double>(coord.dataFirst(), coord.dataFirst()+nDest_) ;
119    }
120
121    CArray<size_t,1> globalIndex(ngloSrc_) ;
122    for(int i=0;i<ngloSrc_;i++) 
123    {
124      transformationMapping_[i] = i ;
125      globalIndex(i) = i ;
126    }
127   
128
129    shared_ptr<CLocalElement> axisSourceGlo = make_shared<CLocalElement>(CContext::getCurrent()->getIntraCommRank(), ngloSrc_, globalIndex) ;
130    axisSourceGlo->addFullView() ; 
131
132    this->computeAlgorithm(axisSource->getLocalView(CElementView::WORKFLOW), axisSourceGlo->getView(CElementView::FULL)) ;
133
134    if (!hasCoordinateSrc_)
135    {
136      CArray<double,1> coord ;
137      CArray<double,1> coordGlo ;
138      auto srcConnector = make_shared<CLocalConnector>(axisSrc_->getLocalView(CElementView::FULL), axisSrc_->getLocalView(CElementView::WORKFLOW)) ;
139      srcConnector->computeConnector() ;
140      srcConnector->transfer(axisSrc_->value, coord) ; // full view value -> workflow value
141      transferTransformConnector_ -> transfer(coord, coordGlo) ; // workflow view -> full global view
142      srcCoordinate_ = vector<double>(coordGlo.dataFirst(), coordGlo.dataFirst()+ngloSrc_) ;
143    }
144  }
145  CATCH
146
147  CTransformFilter* CAxisAlgorithmInterpolateCoordinate::createTransformFilter(CGarbageCollector& gc, shared_ptr<CGridAlgorithm> algo, bool detectMissingValues, double defaultValue)
148  {
149    if (hasCoordinateSrc_ && hasCoordinateDest_) return new CTransformFilter(gc, 3, algo, detectMissingValues, defaultValue) ;
150    else return new CTransformFilter(gc, 2, algo, detectMissingValues, defaultValue) ; 
151  }
152
153  void CAxisAlgorithmInterpolateCoordinate::apply(int dimBefore, int dimAfter, const CArray<double,1>& dataIn, 
154                                                   const vector<CArray<double,1>>& auxDataIn, CArray<double,1>& dataOut)
155  {
156    CArray<double,1> dataInTmp;
157    CArray<double,1> auxDataInSrc ;
158    transferTransformConnector_ -> transfer(dimBefore, dimAfter, dataIn, dataInTmp) ;
159    if (hasCoordinateSrc_)  transferTransformConnector_ -> transfer(dimBefore, dimAfter, auxDataIn[0], auxDataInSrc) ;
160
161   
162    dataOut.resize(dimBefore*dimAfter*nDest_) ;
163    const double* pressureSrc ;
164    const double* pressureDest ;
165    if (hasCoordinateSrc_)  pressureSrc = auxDataInSrc.dataFirst() ;
166    if (hasCoordinateDest_ && hasCoordinateSrc_) pressureDest = auxDataIn[1].dataFirst() ;
167    else if (hasCoordinateDest_ && !hasCoordinateSrc_ ) pressureDest = auxDataIn[0].dataFirst() ;
168   
169    const double* in = dataInTmp.dataFirst() ;
170    double* out = dataOut.dataFirst() ;
171
172    size_t sliceSrc  = dimAfter*ngloSrc_ ;
173    size_t sliceDest = dimAfter*nDest_ ;
174    vector<double> srcCoordinate(ngloSrc_) ;
175    vector<double> destCoordinate(nDest_) ;
176    std::vector<int> srcIndex(ngloSrc_);
177    vector<double> srcValue(ngloSrc_) ;
178    vector<double> destValue(nDest_) ;
179    std::vector<int> destIndex(nDest_);
180
181    size_t nsrc, nDest ;
182
183    for(size_t j=0, posJsrc=0, posJdest=0 ;  j<dimBefore ; j++, posJsrc+=sliceSrc, posJdest+=sliceDest )
184      for(size_t k=0, posKsrc=posJsrc, posKdest=posJdest ; k<dimAfter ; k++, posKsrc++,posKdest++)
185      {
186        if (hasCoordinateSrc_)
187        {
188          nsrc=0 ;
189          for(size_t i=0, posIsrc=posKsrc, posIdest=posKdest ; i<ngloSrc_ ; i++, posIsrc+=dimAfter, posIdest+=dimAfter)
190          {
191            if ( !( std::isnan(pressureSrc[posIsrc]) || std::isnan(in[posIsrc]) ) )
192            {
193              srcCoordinate[nsrc]=pressureSrc[posIsrc] ;
194              srcValue[nsrc] = in[posIsrc] ;
195              nsrc++ ;
196            }
197          }
198        }
199        else
200        {
201          nsrc=0 ;
202          for(size_t i=0, posIsrc=posKsrc ; i<ngloSrc_ ; i++, posIsrc+=dimAfter)
203          {
204            if ( !( std::isnan(srcCoordinate_[i]) || std::isnan(in[posIsrc]) ) )
205            {
206              srcCoordinate[nsrc]=srcCoordinate_[i] ;
207              srcValue[nsrc] = in[posIsrc] ;
208              nsrc++ ;
209            }
210          }
211        } 
212
213        if (hasCoordinateDest_)
214        {
215          nDest=0 ;
216          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)
217          {
218            if ( !( std::isnan(pressureDest[posIdest]) ) )
219            {
220              destCoordinate[nDest]=pressureDest[posIdest] ;
221              nDest++ ;
222            }
223          }
224        }
225        else
226        {
227          nDest=0 ;
228          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)
229          {
230            if ( !( std::isnan(destCoordinate[i]) ) )
231            {
232              destCoordinate[nDest]=destCoordinate_[i] ;
233              nDest++ ;
234            }
235          }
236        }
237 
238        computeInterp(nsrc, srcCoordinate, srcValue, srcIndex, nDest, destCoordinate, destValue, destIndex) ;
239       
240        if (hasCoordinateDest_)
241        {
242          nDest=0 ;
243          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter) 
244          {
245            if ( !( std::isnan(pressureDest[posIdest]) ) ) 
246            {
247              out[posIdest] = destValue[nDest] ;
248              nDest++ ;
249            }
250            else out[posIdest] = std::numeric_limits<double>::quiet_NaN() ;
251          }
252        }
253        else
254        {
255          nDest=0 ;
256          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter) 
257          {
258            if ( !( std::isnan(destCoordinate[i]) ) ) 
259            {
260              out[posIdest] = destValue[nDest] ;
261              nDest++ ;
262            }
263            else out[posIdest] = std::numeric_limits<double>::quiet_NaN() ;
264          }
265        }
266      }
267  }
268
269  void CAxisAlgorithmInterpolateCoordinate::computeInterp(int nsrc, vector<double>& srcCoordinate, vector<double>& srcValue, vector<int>& srcIndex,
270                                                          int ndst, vector<double>& dstCoordinate, vector<double>& dstValue, vector<int>& dstIndex)
271  {
272    double x,y ;
273    double d ;
274 
275    iota(srcIndex.data(), srcIndex.data()+nsrc, 0); // sort array and retrive sorted index
276    stable_sort(srcIndex.data(), srcIndex.data()+nsrc, [&srcCoordinate](size_t i1, size_t i2) {return srcCoordinate[i1] < srcCoordinate[i2];});
277
278    iota(dstIndex.data(), dstIndex.data()+ndst, 0);
279    stable_sort(dstIndex.data(), dstIndex.data()+ndst, [&dstCoordinate](size_t i1, size_t i2) {return dstCoordinate[i1] < dstCoordinate[i2];});
280
281    if (order_==1 || nsrc<=2)
282    {
283      if (nsrc<=1) dstValue.assign(ndst,std::numeric_limits<double>::quiet_NaN());
284      else
285      {
286
287        double x0,x1 ;
288        double y0,y1 ;
289        int lastj=0;
290     
291        for(int i=0; i < ndst;i++)
292        {
293          bool intrapolating( true );
294          x=dstCoordinate[dstIndex[i]] ;
295          if ( x<=srcCoordinate[srcIndex[0]])
296          {
297            lastj=0 ;
298            intrapolating = false;
299          }
300          else if (x>=srcCoordinate[srcIndex[nsrc-1]])
301          {
302            lastj=nsrc-2 ;
303            intrapolating = false;
304          }
305          else
306          {
307            for(int j=lastj; j<nsrc; j++)
308            { 
309              lastj=j ;
310              if (x >= srcCoordinate[srcIndex[j]] && x<srcCoordinate[srcIndex[j+1]]) break ;
311            } 
312          }
313          if ((!extrapolate_) && (!intrapolating))
314          {
315            dstValue[dstIndex[i]]= std::numeric_limits<double>::quiet_NaN() ;
316          }
317          else
318          {
319            x0=srcCoordinate[srcIndex[lastj]] ;
320            x1=srcCoordinate[srcIndex[lastj+1]] ;
321            y0=srcValue[srcIndex[lastj]] ;
322            y1=srcValue[srcIndex[lastj+1]] ;
323            y=((x-x1)*y0-(x-x0)*y1)/(x0-x1) ;
324            dstValue[dstIndex[i]]=y ;
325          }
326        }
327      }
328    }
329    else if (order_==2)
330    {
331      double x0,x1,x2 ;
332      double y0,y1,y2 ;
333      int lastj=0, cj ;
334     
335      for(int i=0; i < ndst;i++)
336      {
337        bool intrapolating( true );
338        x=dstCoordinate[dstIndex[i]] ;
339        if ( x<=srcCoordinate[srcIndex[0]])
340        {
341          lastj=0 ;
342          intrapolating = false;
343        }
344        else if (x>=srcCoordinate[srcIndex[nsrc-1]])
345        {
346          lastj=nsrc-2 ;
347          intrapolating = false;
348        }
349        else
350        {
351          for(int j=lastj; j<nsrc; j++)
352          { 
353            lastj=j ;
354            if (x >= srcCoordinate[srcIndex[j]] && x<srcCoordinate[srcIndex[j+1]])  break ;
355          } 
356        }
357       
358        if (lastj==0) cj=1 ;
359        else if (lastj>=nsrc-2) cj=nsrc-2 ;
360        else
361        {
362          if ( (x-srcCoordinate[srcIndex[lastj-1]]) > (srcCoordinate[srcIndex[lastj+2]]-x) ) cj=lastj ;
363          else cj=lastj+1 ;
364        }
365        if ((!extrapolate_) && (!intrapolating))
366        {
367          dstValue[dstIndex[i]]= std::numeric_limits<double>::quiet_NaN() ;
368        }
369        else
370        {
371          x0=srcCoordinate[srcIndex[cj-1]] ;
372          x1=srcCoordinate[srcIndex[cj]] ;
373          x2=srcCoordinate[srcIndex[cj+1]] ;
374          y0=srcValue[srcIndex[cj-1]] ;
375          y1=srcValue[srcIndex[cj]] ;
376          y2=srcValue[srcIndex[cj+1]] ;
377         
378          y=y0*(x-x1)*(x-x2)/((x0-x1)*(x0-x2)) + y1*(x-x0)*(x-x2)/((x1-x0)*(x1-x2)) + y2*(x-x0)*(x-x1)/((x2-x0)*(x2-x1))  ;
379          dstValue[dstIndex[i]]=y ;
380        }
381
382      }
383    }
384  } 
385
386}
Note: See TracBrowser for help on using the repository browser.