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

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

Add missing return statements detected using -fsanitize=return. Return errors at runtime if reached.
  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 12.7 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    if (!interpAxis->coordinate.isEmpty())
85    {
86      coordinateSrc_ = interpAxis->coordinate.getValue();
87      hasCoordinate_=true ;
88      hasCoordinateSrc_=true ;
89    }
90
91    if (!interpAxis->coordinate_src.isEmpty())
92    {
93      coordinateSrc_ = interpAxis->coordinate_src.getValue();
94      hasCoordinate_=true ;
95      hasCoordinateSrc_=true ;
96    }
97
98    if (!interpAxis->coordinate_dst.isEmpty())
99    {
100      coordinateDest_ = interpAxis->coordinate_dst.getValue();
101      hasCoordinate_=true ;
102      hasCoordinateDest_=true ;
103    }
104
105    ngloSrc_=axisSource->n_glo ; 
106    nDest_ =  axisDest_-> getLocalView(CElementView::WORKFLOW)->getSize() ;
107   
108    if (!hasCoordinateDest_)
109    {
110      CArray<double,1> coord ;
111      auto destConnector = make_shared<CLocalConnector>(axisDest_->getLocalView(CElementView::FULL), axisDest_->getLocalView(CElementView::WORKFLOW)) ;
112      destConnector->computeConnector() ;
113      destConnector->transfer(axisDest_->value, coord) ;
114      destCoordinate_ = vector<double>(coord.dataFirst(), coord.dataFirst()+nDest_) ;
115    }
116
117    CArray<size_t,1> globalIndex(ngloSrc_) ;
118    for(int i=0;i<ngloSrc_;i++) 
119    {
120      transformationMapping_[i] = i ;
121      globalIndex(i) = i ;
122    }
123   
124
125    shared_ptr<CLocalElement> axisSourceGlo = make_shared<CLocalElement>(CContext::getCurrent()->getIntraCommRank(), ngloSrc_, globalIndex) ;
126    axisSourceGlo->addFullView() ; 
127
128    this->computeAlgorithm(axisSource->getLocalView(CElementView::WORKFLOW), axisSourceGlo->getView(CElementView::FULL)) ;
129
130    if (!hasCoordinateSrc_)
131    {
132      CArray<double,1> coord ;
133      CArray<double,1> coordGlo ;
134      auto srcConnector = make_shared<CLocalConnector>(axisSrc_->getLocalView(CElementView::FULL), axisSrc_->getLocalView(CElementView::WORKFLOW)) ;
135      srcConnector->computeConnector() ;
136      srcConnector->transfer(axisSrc_->value, coord) ; // full view value -> workflow value
137      transferTransformConnector_ -> transfer(coord, coordGlo) ; // workflow view -> full global view
138      srcCoordinate_ = vector<double>(coordGlo.dataFirst(), coordGlo.dataFirst()+ngloSrc_) ;
139    }
140  }
141  CATCH
142
143  CTransformFilter* CAxisAlgorithmInterpolateCoordinate::createTransformFilter(CGarbageCollector& gc, shared_ptr<CGridAlgorithm> algo, bool detectMissingValues, double defaultValue)
144  {
145    if (hasCoordinateSrc_ && hasCoordinateDest_) return new CTransformFilter(gc, 3, algo, detectMissingValues, defaultValue) ;
146    else return new CTransformFilter(gc, 2, algo, detectMissingValues, defaultValue) ; 
147  }
148
149  void CAxisAlgorithmInterpolateCoordinate::apply(int dimBefore, int dimAfter, const CArray<double,1>& dataIn, 
150                                                   const vector<CArray<double,1>>& auxDataIn, CArray<double,1>& dataOut)
151  {
152    CArray<double,1> dataInTmp;
153    CArray<double,1> auxDataInSrc ;
154    transferTransformConnector_ -> transfer(dimBefore, dimAfter, dataIn, dataInTmp) ;
155    if (hasCoordinateSrc_)  transferTransformConnector_ -> transfer(dimBefore, dimAfter, auxDataIn[0], auxDataInSrc) ;
156
157   
158    dataOut.resize(dimBefore*dimAfter*nDest_) ;
159    const double* pressureSrc ;
160    const double* pressureDest ;
161    if (hasCoordinateSrc_)  pressureSrc = auxDataInSrc.dataFirst() ;
162    if (hasCoordinateDest_ && hasCoordinateSrc_) pressureDest = auxDataIn[1].dataFirst() ;
163    else if (hasCoordinateDest_ && !hasCoordinateSrc_ ) pressureDest = auxDataIn[0].dataFirst() ;
164   
165    const double* in = dataInTmp.dataFirst() ;
166    double* out = dataOut.dataFirst() ;
167
168    size_t sliceSrc  = dimAfter*ngloSrc_ ;
169    size_t sliceDest = dimAfter*nDest_ ;
170    vector<double> srcCoordinate(ngloSrc_) ;
171    vector<double> destCoordinate(nDest_) ;
172    std::vector<int> srcIndex(ngloSrc_);
173    vector<double> srcValue(ngloSrc_) ;
174    vector<double> destValue(nDest_) ;
175    std::vector<int> destIndex(nDest_);
176
177    size_t nsrc, nDest ;
178
179    for(size_t j=0, posJsrc=0, posJdest=0 ;  j<dimBefore ; j++, posJsrc+=sliceSrc, posJdest+=sliceDest )
180      for(size_t k=0, posKsrc=posJsrc, posKdest=posJdest ; k<dimAfter ; k++, posKsrc++,posKdest++)
181      {
182        if (hasCoordinateSrc_)
183        {
184          nsrc=0 ;
185          for(size_t i=0, posIsrc=posKsrc, posIdest=posKdest ; i<ngloSrc_ ; i++, posIsrc+=dimAfter, posIdest+=dimAfter)
186          {
187            if ( !( std::isnan(pressureSrc[posIsrc]) || std::isnan(in[posIsrc]) ) )
188            {
189              srcCoordinate[nsrc]=pressureSrc[posIsrc] ;
190              srcValue[nsrc] = in[posIsrc] ;
191              nsrc++ ;
192            }
193          }
194        }
195        else
196        {
197          nsrc=0 ;
198          for(size_t i=0, posIsrc=posKsrc ; i<ngloSrc_ ; i++, posIsrc+=dimAfter)
199          {
200            if ( !( std::isnan(srcCoordinate_[i]) || std::isnan(in[posIsrc]) ) )
201            {
202              srcCoordinate[nsrc]=srcCoordinate_[i] ;
203              srcValue[nsrc] = in[posIsrc] ;
204              nsrc++ ;
205            }
206          }
207        } 
208
209        if (hasCoordinateDest_)
210        {
211          nDest=0 ;
212          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)
213          {
214            if ( !( std::isnan(pressureDest[posIdest]) ) )
215            {
216              destCoordinate[nDest]=pressureDest[posIdest] ;
217              nDest++ ;
218            }
219          }
220        }
221        else
222        {
223          nDest=0 ;
224          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter)
225          {
226            if ( !( std::isnan(destCoordinate[i]) ) )
227            {
228              destCoordinate[nDest]=destCoordinate_[i] ;
229              nDest++ ;
230            }
231          }
232        }
233 
234        computeInterp(nsrc, srcCoordinate, srcValue, srcIndex, nDest, destCoordinate, destValue, destIndex) ;
235       
236        if (hasCoordinateDest_)
237        {
238          nDest=0 ;
239          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter) 
240          {
241            if ( !( std::isnan(pressureDest[posIdest]) ) ) 
242            {
243              out[posIdest] = destValue[nDest] ;
244              nDest++ ;
245            }
246            else out[posIdest] = std::numeric_limits<double>::quiet_NaN() ;
247          }
248        }
249        else
250        {
251          nDest=0 ;
252          for(size_t i=0, posIdest=posKdest ; i<nDest_ ; i++, posIdest+=dimAfter) 
253          {
254            if ( !( std::isnan(destCoordinate[i]) ) ) 
255            {
256              out[posIdest] = destValue[nDest] ;
257              nDest++ ;
258            }
259            else out[posIdest] = std::numeric_limits<double>::quiet_NaN() ;
260          }
261        }
262      }
263  }
264
265  void CAxisAlgorithmInterpolateCoordinate::computeInterp(int nsrc, vector<double>& srcCoordinate, vector<double>& srcValue, vector<int>& srcIndex,
266                                                          int ndst, vector<double>& dstCoordinate, vector<double>& dstValue, vector<int>& dstIndex)
267  {
268    double x,y ;
269    double d ;
270 
271    iota(srcIndex.data(), srcIndex.data()+nsrc, 0); // sort array and retrive sorted index
272    stable_sort(srcIndex.data(), srcIndex.data()+nsrc, [&srcCoordinate](size_t i1, size_t i2) {return srcCoordinate[i1] < srcCoordinate[i2];});
273
274    iota(dstIndex.data(), dstIndex.data()+ndst, 0);
275    stable_sort(dstIndex.data(), dstIndex.data()+ndst, [&dstCoordinate](size_t i1, size_t i2) {return dstCoordinate[i1] < dstCoordinate[i2];});
276
277    if (order_==1 || nsrc<=2)
278    {
279      if (nsrc<=1) dstValue.assign(ndst,std::numeric_limits<double>::quiet_NaN());
280      else
281      {
282
283        double x0,x1 ;
284        double y0,y1 ;
285        int lastj=0;
286     
287        for(int i=0; i < ndst;i++)
288        {
289          x=dstCoordinate[dstIndex[i]] ;
290          if ( x<=srcCoordinate[srcIndex[0]]) lastj=0 ;
291          else if (x>=srcCoordinate[srcIndex[nsrc-1]]) lastj=nsrc-2 ;
292          else
293          {
294            for(int j=lastj; j<nsrc; j++)
295            { 
296              lastj=j ;
297              if (x >= srcCoordinate[srcIndex[j]] && x<srcCoordinate[srcIndex[j+1]]) break ;
298            } 
299          }
300          x0=srcCoordinate[srcIndex[lastj]] ;
301          x1=srcCoordinate[srcIndex[lastj+1]] ;
302          y0=srcValue[srcIndex[lastj]] ;
303          y1=srcValue[srcIndex[lastj+1]] ;
304          y=((x-x1)*y0-(x-x0)*y1)/(x0-x1) ;
305          dstValue[dstIndex[i]]=y ;
306        }
307      }
308    }
309    else if (order_==2)
310    {
311      double x0,x1,x2 ;
312      double y0,y1,y2 ;
313      int lastj=0, cj ;
314     
315      for(int i=0; i < ndst;i++)
316      {
317        x=dstCoordinate[dstIndex[i]] ;
318        if ( x<=srcCoordinate[srcIndex[0]]) lastj=0 ;
319        else if (x>=srcCoordinate[srcIndex[nsrc-1]]) lastj=nsrc-2 ;
320        else
321        {
322          for(int j=lastj; j<nsrc; j++)
323          { 
324            lastj=j ;
325            if (x >= srcCoordinate[srcIndex[j]] && x<srcCoordinate[srcIndex[j+1]])  break ;
326          } 
327        }
328       
329        if (lastj==0) cj=1 ;
330        else if (lastj>=nsrc-2) cj=nsrc-2 ;
331        else
332        {
333          if ( (x-srcCoordinate[srcIndex[lastj-1]]) > (srcCoordinate[srcIndex[lastj+2]]-x) ) cj=lastj ;
334          else cj=lastj+1 ;
335        } 
336        x0=srcCoordinate[srcIndex[cj-1]] ;
337        x1=srcCoordinate[srcIndex[cj]] ;
338        x2=srcCoordinate[srcIndex[cj+1]] ;
339        y0=srcValue[srcIndex[cj-1]] ;
340        y1=srcValue[srcIndex[cj]] ;
341        y2=srcValue[srcIndex[cj+1]] ;
342           
343        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))  ;
344        dstValue[dstIndex[i]]=y ;
345      }
346    }
347  } 
348
349}
Note: See TracBrowser for help on using the repository browser.