source: XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp @ 2297

Last change on this file since 2297 was 1937, checked in by ymipsl, 4 years ago

Fix problem when using area for domain interpolation.

YM

File size: 37.4 KB
RevLine 
[657]1/*!
2   \file domain_algorithm_interpolate_from_file.cpp
3   \author Ha NGUYEN
4   \since 09 Jul 2015
[689]5   \date 15 Sep 2015
[657]6
7   \brief Algorithm for interpolation on a domain.
8 */
[689]9#include "domain_algorithm_interpolate.hpp"
[1542]10#include <unordered_map>
[657]11#include "context.hpp"
12#include "context_client.hpp"
13#include "distribution_client.hpp"
14#include "client_server_mapping_distributed.hpp"
[660]15#include "netcdf.hpp"
[688]16#include "mapper.hpp"
[821]17#include "mpi_tag.hpp"
[933]18#include "domain.hpp"
19#include "grid_transformation_factory_impl.hpp"
20#include "interpolate_domain.hpp"
21#include "grid.hpp"
[657]22
23namespace xios {
[933]24CGenericAlgorithmTransformation* CDomainAlgorithmInterpolate::create(CGrid* gridDst, CGrid* gridSrc,
25                                                                     CTransformation<CDomain>* transformation,
26                                                                     int elementPositionInGrid,
27                                                                     std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
28                                                                     std::map<int, int>& elementPositionInGridSrc2AxisPosition,
29                                                                     std::map<int, int>& elementPositionInGridSrc2DomainPosition,
30                                                                     std::map<int, int>& elementPositionInGridDst2ScalarPosition,
31                                                                     std::map<int, int>& elementPositionInGridDst2AxisPosition,
32                                                                     std::map<int, int>& elementPositionInGridDst2DomainPosition)
[1622]33TRY
[933]34{
35  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
36  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
[657]37
[933]38  CInterpolateDomain* interpolateDomain = dynamic_cast<CInterpolateDomain*> (transformation);
[978]39  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
40  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
[933]41
42  return (new CDomainAlgorithmInterpolate(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], interpolateDomain));
43}
[1622]44CATCH
[933]45
46bool CDomainAlgorithmInterpolate::registerTrans()
[1622]47TRY
[933]48{
[1852]49  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_INTERPOLATE_DOMAIN, create);
[933]50}
[1622]51CATCH
[933]52
[689]53CDomainAlgorithmInterpolate::CDomainAlgorithmInterpolate(CDomain* domainDestination, CDomain* domainSource, CInterpolateDomain* interpDomain)
[1021]54: CDomainAlgorithmTransformation(domainDestination, domainSource), interpDomain_(interpDomain), writeToFile_(false), readFromFile_(false)
[1622]55TRY
[657]56{
[1021]57  CContext* context = CContext::getCurrent();
[657]58  interpDomain_->checkValid(domainSource);
[1264]59
60  detectMissingValue = interpDomain_->detect_missing_value ;
61  renormalize = interpDomain_->renormalize ;
62  quantity = interpDomain_->quantity ;
[1269]63
64  if (interpDomain_->read_write_convention == CInterpolateDomain::read_write_convention_attr::fortran) fortranConvention=true ;
65  else fortranConvention=false ;
[1264]66 
[1021]67  fileToReadWrite_ = "xios_interpolation_weights_";
68
69  if (interpDomain_->weight_filename.isEmpty())
70  {
71    fileToReadWrite_ += context->getId() + "_" + 
72                    domainSource->getDomainOutputName() + "_" + 
73                    domainDestination->getDomainOutputName() + ".nc";   
74  }
75  else 
76    fileToReadWrite_ = interpDomain_->weight_filename;
77
78  ifstream f(fileToReadWrite_.c_str()); 
79  switch (interpDomain_->mode)
80  {
81    case CInterpolateDomain::mode_attr::read:
82      readFromFile_ = true;     
83      break;
84    case CInterpolateDomain::mode_attr::compute:
85      readFromFile_ = false;
86      break;
87    case CInterpolateDomain::mode_attr::read_or_compute:     
88      if (!f.good())
89        readFromFile_ = false;
90      else
91        readFromFile_ = true;
92      break;
93    default:
94      break;
95  } 
96
97  writeToFile_ = interpDomain_->write_weight; 
98   
[657]99}
[1622]100CATCH
[657]101
[689]102/*!
103  Compute remap with integrated remap calculation module
104*/
105void CDomainAlgorithmInterpolate::computeRemap()
[1622]106TRY
[688]107{
108  using namespace sphereRemap;
109
110  CContext* context = CContext::getCurrent();
111  CContextClient* client=context->client;
112  int clientRank = client->clientRank;
113  int i, j, k, idx;
[689]114  std::vector<double> srcPole(3,0), dstPole(3,0);
[915]115  int orderInterp = interpDomain_->order.getValue();
[688]116
[846]117
[743]118  const double poleValue = 90.0;
119  const int constNVertex = 4; // Value by default number of vertex for rectangular domain
[688]120  int nVertexSrc, nVertexDest;
121  nVertexSrc = nVertexDest = constNVertex;
122
123  // First of all, try to retrieve the boundary values of domain source and domain destination
124  int localDomainSrcSize = domainSrc_->i_index.numElements();
125  int niSrc = domainSrc_->ni.getValue(), njSrc = domainSrc_->nj.getValue();
126  bool hasBoundSrc = domainSrc_->hasBounds;
127  if (hasBoundSrc) nVertexSrc = domainSrc_->nvertex.getValue();
128  CArray<double,2> boundsLonSrc(nVertexSrc,localDomainSrcSize);
129  CArray<double,2> boundsLatSrc(nVertexSrc,localDomainSrcSize);
[1937]130  CArray<double,1> areaSrc;
[688]131
[953]132  if (domainSrc_->hasPole) srcPole[2] = 1;
[688]133  if (hasBoundSrc)  // Suppose that domain source is curvilinear or unstructured
134  {
135    if (!domainSrc_->bounds_lon_2d.isEmpty())
136    {
137      for (j = 0; j < njSrc; ++j)
138        for (i = 0; i < niSrc; ++i)
139        {
140          k=j*niSrc+i;
141          for(int n=0;n<nVertexSrc;++n)
142          {
143            boundsLonSrc(n,k) = domainSrc_->bounds_lon_2d(n,i,j);
144            boundsLatSrc(n,k) = domainSrc_->bounds_lat_2d(n,i,j);
145          }
146        }
147    }
148    else
149    {
150      boundsLonSrc = domainSrc_->bounds_lon_1d;
151      boundsLatSrc = domainSrc_->bounds_lat_1d;
152    }
153  }
154  else // if domain source is rectilinear, not do anything now
155  {
[809]156    CArray<double,1> lon_g ;
157    CArray<double,1> lat_g ;
[753]158
[809]159    if (!domainSrc_->lonvalue_1d.isEmpty() && !domainSrc_->latvalue_1d.isEmpty())
[808]160    {
[915]161      domainSrc_->AllgatherRectilinearLonLat(domainSrc_->lonvalue_1d,domainSrc_->latvalue_1d, lon_g,lat_g) ;
162    }
163    else if (! domainSrc_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainSrc_->lonvalue_rectilinear_read_from_file.isEmpty() )
[808]164    {
[915]165      lat_g=domainSrc_->latvalue_rectilinear_read_from_file ;
166      lon_g=domainSrc_->lonvalue_rectilinear_read_from_file ;
167    }
168    else if (!domainSrc_->lon_start.isEmpty() && !domainSrc_->lon_end.isEmpty() &&
169             !domainSrc_->lat_start.isEmpty() && !domainSrc_->lat_end.isEmpty())
170    {
171      double step=(domainSrc_->lon_end-domainSrc_->lon_start)/domainSrc_->ni_glo ;
172      for (int i=0; i<domainSrc_->ni_glo; ++i) lon_g(i)=domainSrc_->lon_start+i*step ;
173      step=(domainSrc_->lat_end-domainSrc_->lat_start)/domainSrc_->nj_glo ;
174      for (int i=0; i<domainSrc_->ni_glo; ++i) lat_g(i)=domainSrc_->lat_start+i*step ;
175    }
176    else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
[808]177
[688]178    nVertexSrc = constNVertex;
[809]179    domainSrc_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonSrc, boundsLatSrc);
[688]180  }
[1937]181 
182  if (domainSrc_->hasArea) 
183  {
184    areaSrc.resize(niSrc*njSrc);
185    for (int j = 0; j < njSrc; ++j)
186      for (int i = 0; i < niSrc; ++i)
187      {
188        int k = j * niSrc + i;
189        areaSrc(k) = domainSrc_->area(i,j);
190      }
191  }
192 
[743]193  std::map<int,std::vector<std::pair<int,double> > > interpMapValueNorthPole;
194  std::map<int,std::vector<std::pair<int,double> > > interpMapValueSouthPole;
195
[688]196  int localDomainDestSize = domainDest_->i_index.numElements();
197  int niDest = domainDest_->ni.getValue(), njDest = domainDest_->nj.getValue();
198  bool hasBoundDest = domainDest_->hasBounds;
199  if (hasBoundDest) nVertexDest = domainDest_->nvertex.getValue();
200  CArray<double,2> boundsLonDest(nVertexDest,localDomainDestSize);
201  CArray<double,2> boundsLatDest(nVertexDest,localDomainDestSize);
[1937]202  CArray<double,1> areaDest;
[688]203
[953]204  if (domainDest_->hasPole) dstPole[2] = 1;
[688]205  if (hasBoundDest)
206  {
207    if (!domainDest_->bounds_lon_2d.isEmpty())
208    {
209      for (j = 0; j < njDest; ++j)
210        for (i = 0; i < niDest; ++i)
211        {
212          k=j*niDest+i;
213          for(int n=0;n<nVertexDest;++n)
214          {
215            boundsLonDest(n,k) = domainDest_->bounds_lon_2d(n,i,j);
216            boundsLatDest(n,k) = domainDest_->bounds_lat_2d(n,i,j);
217          }
218        }
219    }
220    else
221    {
222      boundsLonDest = domainDest_->bounds_lon_1d;
223      boundsLatDest = domainDest_->bounds_lat_1d;
224    }
225  }
226  else
227  {
[753]228    bool isNorthPole = false;
229    bool isSouthPole = false;
[809]230
231    CArray<double,1> lon_g ;
232    CArray<double,1> lat_g ;
233
234    if (!domainDest_->lonvalue_1d.isEmpty() && !domainDest_->latvalue_1d.isEmpty())
235    {
[949]236      domainDest_->AllgatherRectilinearLonLat(domainDest_->lonvalue_1d,domainDest_->latvalue_1d, lon_g,lat_g) ;
237    }
238    else if (! domainDest_->latvalue_rectilinear_read_from_file.isEmpty() && ! domainDest_->lonvalue_rectilinear_read_from_file.isEmpty() )
[809]239    {
[949]240      lat_g=domainDest_->latvalue_rectilinear_read_from_file ;
241      lon_g=domainDest_->lonvalue_rectilinear_read_from_file ;
242    }
243    else if (!domainDest_->lon_start.isEmpty() && !domainDest_->lon_end.isEmpty() &&
244             !domainDest_->lat_start.isEmpty() && !domainDest_->lat_end.isEmpty())
245    {
246      double step=(domainDest_->lon_end-domainDest_->lon_start)/domainDest_->ni_glo ;
247      for(int i=0; i<domainDest_->ni_glo; ++i) lon_g(i)=domainDest_->lon_start+i*step ;
248      step=(domainDest_->lat_end-domainDest_->lat_start)/domainDest_->nj_glo ;
249      for(int i=0; i<domainDest_->ni_glo; ++i) lat_g(i)=domainDest_->lat_start+i*step ;
250    }
251    else ERROR("void CDomainAlgorithmInterpolate::computeRemap()",<<"Cannot compute bounds for rectilinear domain") ;
252   
[809]253    if (std::abs(poleValue - std::abs(lat_g(0))) < NumTraits<double>::epsilon()) isNorthPole = true;
254    if (std::abs(poleValue - std::abs(lat_g(domainDest_->nj_glo-1))) < NumTraits<double>::epsilon()) isSouthPole = true;
255
[821]256
257
258
[743]259    if (isNorthPole && (0 == domainDest_->jbegin.getValue()))
260    {
261      int ibegin = domainDest_->ibegin.getValue();
262      for (i = 0; i < niDest; ++i)
263      {
264        interpMapValueNorthPole[i+ibegin];
265      }
266    }
267
268    if (isSouthPole && (domainDest_->nj_glo.getValue() == (domainDest_->jbegin.getValue() + njDest)))
269    {
270      int ibegin = domainDest_->ibegin.getValue();
271      int njGlo = domainDest_->nj_glo.getValue();
272      int niGlo = domainDest_->ni_glo.getValue();
273      for (i = 0; i < niDest; ++i)
274      {
275        k = (njGlo - 1)*niGlo + i + ibegin;
276        interpMapValueSouthPole[k];
277      }
278    }
279
[688]280    // Ok, fill in boundary values for rectangular domain
[689]281    nVertexDest = constNVertex;
[809]282    domainDest_->fillInRectilinearBoundLonLat(lon_g,lat_g, boundsLonDest, boundsLatDest);
[688]283  }
[1937]284 
285  if (domainDest_->hasArea) 
286  {
287    areaDest.resize(niDest*njDest);
288    for (int j = 0; j < njDest; ++j)
289      for (int i = 0; i < niDest; ++i)
290      {
291        int k = j * niDest + i;
292        areaDest(k) = domainDest_->area(i,j);
293      }
294  }
295 
[688]296
[709]297
[688]298  // Ok, now use mapper to calculate
299  int nSrcLocal = domainSrc_->i_index.numElements();
300  int nDstLocal = domainDest_->i_index.numElements();
[709]301  long int * globalSrc = new long int [nSrcLocal];
302  long int * globalDst = new long int [nDstLocal];
303
304  long int globalIndex;
305  int i_ind, j_ind;
306  for (int idx = 0; idx < nSrcLocal; ++idx)
307  {
308    i_ind=domainSrc_->i_index(idx) ;
309    j_ind=domainSrc_->j_index(idx) ;
310
311    globalIndex = i_ind + j_ind * domainSrc_->ni_glo;
312    globalSrc[idx] = globalIndex;
313  }
314
315  for (int idx = 0; idx < nDstLocal; ++idx)
316  {
317    i_ind=domainDest_->i_index(idx) ;
318    j_ind=domainDest_->j_index(idx) ;
319
320    globalIndex = i_ind + j_ind * domainDest_->ni_glo;
321    globalDst[idx] = globalIndex;
322  }
323
324
325  // Calculate weight index
[688]326  Mapper mapper(client->intraComm);
327  mapper.setVerbosity(PROGRESS) ;
328
[880]329     
[846]330  // supress masked data for the source
331  int nSrcLocalUnmasked = 0 ;
[893]332  for (int idx=0 ; idx < nSrcLocal; idx++) if (domainSrc_->localMask(idx)) ++nSrcLocalUnmasked ;
[846]333
[880]334
[846]335  CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
336  CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked);
[1615]337  CArray<double,1> areaSrcUnmasked(nSrcLocalUnmasked);
338 
[846]339  long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked];
340
341  nSrcLocalUnmasked=0 ;
[1617]342  bool hasSrcArea=domainSrc_->hasArea && !domainSrc_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true  ;
[1615]343  double srcAreaFactor ;
344  if (hasSrcArea) srcAreaFactor=1./(domainSrc_->radius*domainSrc_->radius) ;
345 
[846]346  for (int idx=0 ; idx < nSrcLocal; idx++)
347  {
[893]348    if (domainSrc_->localMask(idx))
[846]349    {
350      for(int n=0;n<nVertexSrc;++n)
351      {
352        boundsLonSrcUnmasked(n,nSrcLocalUnmasked) = boundsLonSrc(n,idx) ;
353        boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ;
354      }
[1937]355      if (hasSrcArea) areaSrcUnmasked(nSrcLocalUnmasked) = areaSrc(idx)*srcAreaFactor ;
[846]356      globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ;
357      ++nSrcLocalUnmasked ;
358    }
359  }
[1615]360 
[846]361
362  int nDstLocalUnmasked = 0 ;
[893]363  for (int idx=0 ; idx < nDstLocal; idx++) if (domainDest_->localMask(idx)) ++nDstLocalUnmasked ;
[846]364
365  CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked);
366  CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked);
[1617]367  CArray<double,1>   areaDstUnmasked(nDstLocalUnmasked);
[1615]368
[846]369  long int * globalDstUnmasked = new long int [nDstLocalUnmasked];
370
371  nDstLocalUnmasked=0 ;
[1617]372  bool hasDstArea=domainDest_->hasArea && !domainDest_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true ;
[1615]373  double dstAreaFactor ;
374  if (hasDstArea) dstAreaFactor=1./(domainDest_->radius*domainDest_->radius) ;
[846]375  for (int idx=0 ; idx < nDstLocal; idx++)
376  {
[893]377    if (domainDest_->localMask(idx))
[846]378    {
379      for(int n=0;n<nVertexDest;++n)
380      {
381        boundsLonDestUnmasked(n,nDstLocalUnmasked) = boundsLonDest(n,idx) ;
382        boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ;
383      }
[1937]384      if (hasDstArea) areaDstUnmasked(nDstLocalUnmasked) = areaDest(idx)*dstAreaFactor ;
[846]385      globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ;
386      ++nDstLocalUnmasked ;
387    }
388  }
[856]389
[1615]390  double* ptAreaSrcUnmasked = NULL ;
391  if (hasSrcArea) ptAreaSrcUnmasked=areaSrcUnmasked.dataFirst() ;
[846]392
[1615]393  double* ptAreaDstUnmasked = NULL ;
394  if (hasDstArea) ptAreaDstUnmasked=areaDstUnmasked.dataFirst() ;
395
396  mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), ptAreaSrcUnmasked, nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked);
397  mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), ptAreaDstUnmasked, nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked);
398
[1158]399  std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity);
[846]400
[709]401  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
[743]402  std::map<int,std::vector<std::pair<int,double> > >::const_iterator iteNorthPole = interpMapValueNorthPole.end(),
403                                                                     iteSouthPole = interpMapValueSouthPole.end();
[688]404  for (int idx = 0;  idx < mapper.nWeights; ++idx)
405  {
[709]406    interpMapValue[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
[743]407    if (iteNorthPole != interpMapValueNorthPole.find(mapper.targetWeightId[idx]))
408    {
409      interpMapValueNorthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
410    }
411
412    if (iteSouthPole != interpMapValueSouthPole.find(mapper.targetWeightId[idx]))
413    {
414      interpMapValueSouthPole[mapper.targetWeightId[idx]].push_back(make_pair(mapper.sourceWeightId[idx],mapper.remapMatrix[idx]));
415    }
[688]416  }
[743]417  int niGloDst = domainDest_->ni_glo.getValue();
418  processPole(interpMapValueNorthPole, niGloDst);
419  processPole(interpMapValueSouthPole, niGloDst);
420
421  if (!interpMapValueNorthPole.empty())
422  {
423     std::map<int,std::vector<std::pair<int,double> > >::iterator itNorthPole = interpMapValueNorthPole.begin();
424     for (; itNorthPole != iteNorthPole; ++itNorthPole)
425     {
426       if (!(itNorthPole->second.empty()))
427        itNorthPole->second.swap(interpMapValue[itNorthPole->first]);
428     }
429  }
430
431  if (!interpMapValueSouthPole.empty())
432  {
433     std::map<int,std::vector<std::pair<int,double> > >::iterator itSouthPole = interpMapValueSouthPole.begin();
434     for (; itSouthPole != iteSouthPole; ++itSouthPole)
435     {
436       if (!(itSouthPole->second.empty()))
437        itSouthPole->second.swap(interpMapValue[itSouthPole->first]);
438     }
439  }
440
[1173]441  if (writeToFile_ && !readFromFile_) writeRemapInfo(interpMapValue);
442//  exchangeRemapInfo(interpMapValue);
443  convertRemapInfo(interpMapValue) ;
[709]444
445  delete [] globalSrc;
[846]446  delete [] globalSrcUnmasked;
[709]447  delete [] globalDst;
[846]448  delete [] globalDstUnmasked;
449
[688]450}
[1622]451CATCH
[688]452
[743]453void CDomainAlgorithmInterpolate::processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole,
454                                              int nbGlobalPointOnPole)
[1622]455TRY
[743]456{
457  CContext* context = CContext::getCurrent();
458  CContextClient* client=context->client;
459
[1639]460  MPI_Comm poleComme(MPI_COMM_NULL);
[743]461  MPI_Comm_split(client->intraComm, interMapValuePole.empty() ? MPI_UNDEFINED : 1, 0, &poleComme);
[1639]462  if (MPI_COMM_NULL != poleComme)
[743]463  {
464    int nbClientPole;
[1639]465    MPI_Comm_size(poleComme, &nbClientPole);
[743]466
467    std::map<int,std::vector<std::pair<int,double> > >::iterator itePole = interMapValuePole.end(), itPole,
468                                                                 itbPole = interMapValuePole.begin();
469
470    int nbWeight = 0;
471    for (itPole = itbPole; itPole != itePole; ++itPole)
472       nbWeight += itPole->second.size();
473
474    std::vector<int> recvCount(nbClientPole,0);
475    std::vector<int> displ(nbClientPole,0);
[1639]476    MPI_Allgather(&nbWeight,1,MPI_INT,&recvCount[0],1,MPI_INT,poleComme) ;
[743]477
478    displ[0]=0;
479    for(int n=1;n<nbClientPole;++n) displ[n]=displ[n-1]+recvCount[n-1] ;
480    int recvSize=displ[nbClientPole-1]+recvCount[nbClientPole-1] ;
481
482    std::vector<int> sendSourceIndexBuff(nbWeight);
483    std::vector<double> sendSourceWeightBuff(nbWeight);
484    int k = 0;
485    for (itPole = itbPole; itPole != itePole; ++itPole)
486    {
487      for (int idx = 0; idx < itPole->second.size(); ++idx)
488      {
489        sendSourceIndexBuff[k] = (itPole->second)[idx].first;
490        sendSourceWeightBuff[k] = (itPole->second)[idx].second;
491        ++k;
492      }
493    }
494
495    std::vector<int> recvSourceIndexBuff(recvSize);
496    std::vector<double> recvSourceWeightBuff(recvSize);
497
498    // Gather all index and weight for pole
[1639]499    MPI_Allgatherv(&sendSourceIndexBuff[0],nbWeight,MPI_INT,&recvSourceIndexBuff[0],&recvCount[0],&displ[0],MPI_INT,poleComme);
500    MPI_Allgatherv(&sendSourceWeightBuff[0],nbWeight,MPI_DOUBLE,&recvSourceWeightBuff[0],&recvCount[0],&displ[0],MPI_DOUBLE,poleComme);
[743]501
502    std::map<int,double> recvTemp;
503    for (int idx = 0; idx < recvSize; ++idx)
504    {
505      if (recvTemp.end() != recvTemp.find(recvSourceIndexBuff[idx]))
506        recvTemp[recvSourceIndexBuff[idx]] += recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
507      else
[1317]508        recvTemp[recvSourceIndexBuff[idx]] = recvSourceWeightBuff[idx]/nbGlobalPointOnPole;
[743]509    }
510
511    std::map<int,double>::const_iterator itRecvTemp, itbRecvTemp = recvTemp.begin(), iteRecvTemp = recvTemp.end();
512
513    for (itPole = itbPole; itPole != itePole; ++itPole)
514    {
515      itPole->second.clear();
516      for (itRecvTemp = itbRecvTemp; itRecvTemp != iteRecvTemp; ++itRecvTemp)
517          itPole->second.push_back(make_pair(itRecvTemp->first, itRecvTemp->second));
518    }
519  }
520
521}
[1622]522CATCH
[743]523
[689]524/*!
525  Compute the index mapping between domain on grid source and one on grid destination
526*/
[827]527void CDomainAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
[1622]528TRY
[688]529{
[1267]530  if (readFromFile_) 
[689]531    readRemapInfo();
532  else
[1021]533  {
534    computeRemap(); 
535  }
[688]536}
[1622]537CATCH
[688]538
[1021]539void CDomainAlgorithmInterpolate::writeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
[1622]540TRY
[1021]541{ 
542  writeInterpolationInfo(fileToReadWrite_, interpMapValue);
543}
[1622]544CATCH
[1021]545
[689]546void CDomainAlgorithmInterpolate::readRemapInfo()
[1622]547TRY
[1021]548{ 
[657]549  std::map<int,std::vector<std::pair<int,double> > > interpMapValue;
[1021]550  readInterpolationInfo(fileToReadWrite_, interpMapValue);
[657]551
[709]552  exchangeRemapInfo(interpMapValue);
553}
[1622]554CATCH
[709]555
[1173]556void CDomainAlgorithmInterpolate::convertRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
[1622]557TRY
[1173]558{
559  CContext* context = CContext::getCurrent();
560  CContextClient* client=context->client;
561  int clientRank = client->clientRank;
[709]562
[1173]563  this->transformationMapping_.resize(1);
564  this->transformationWeight_.resize(1);
565
566  TransformationIndexMap& transMap = this->transformationMapping_[0];
567  TransformationWeightMap& transWeight = this->transformationWeight_[0];
568
569  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
570                                                                     ite = interpMapValue.end();
571 
572  for (it = itb; it != ite; ++it)
573  {   
574    const std::vector<std::pair<int,double> >& tmp = it->second;
575    for (int i = 0; i < tmp.size(); ++i)
576    {
577      transMap[it->first].push_back(tmp[i].first);
578      transWeight[it->first].push_back(tmp[i].second);
579    }     
580  }     
581}
[1622]582CATCH
[1173]583
[709]584/*!
585  Read remap information from file then distribute it among clients
586*/
[856]587void CDomainAlgorithmInterpolate::exchangeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
[1622]588TRY
[709]589{
590  CContext* context = CContext::getCurrent();
591  CContextClient* client=context->client;
592  int clientRank = client->clientRank;
593
[827]594  this->transformationMapping_.resize(1);
595  this->transformationWeight_.resize(1);
596
[833]597  TransformationIndexMap& transMap = this->transformationMapping_[0];
598  TransformationWeightMap& transWeight = this->transformationWeight_[0];
[827]599
[1542]600  std::unordered_map<size_t,int> globalIndexOfDomainDest;
[657]601  int ni = domainDest_->ni.getValue();
602  int nj = domainDest_->nj.getValue();
603  int ni_glo = domainDest_->ni_glo.getValue();
604  size_t globalIndex;
605  int nIndexSize = domainDest_->i_index.numElements(), i_ind, j_ind;
606  for (int idx = 0; idx < nIndexSize; ++idx)
607  {
608    i_ind=domainDest_->i_index(idx) ;
609    j_ind=domainDest_->j_index(idx) ;
610
611    globalIndex = i_ind + j_ind * ni_glo;
612    globalIndexOfDomainDest[globalIndex] = clientRank;
613  }
614
615  CClientServerMappingDistributed domainIndexClientClientMapping(globalIndexOfDomainDest,
616                                                                 client->intraComm,
617                                                                 true);
618  CArray<size_t,1> globalIndexInterp(interpMapValue.size());
619  std::map<int,std::vector<std::pair<int,double> > >::const_iterator itb = interpMapValue.begin(), it,
620                                                                     ite = interpMapValue.end();
621  size_t globalIndexCount = 0;
622  for (it = itb; it != ite; ++it)
623  {
624    globalIndexInterp(globalIndexCount) = it->first;
625    ++globalIndexCount;
626  }
627
[1336]628  domainIndexClientClientMapping.computeServerIndexMapping(globalIndexInterp, client->serverSize);
[829]629  const CClientServerMapping::GlobalIndexMap& globalIndexInterpSendToClient = domainIndexClientClientMapping.getGlobalIndexOnServer();
[657]630
631  //Inform each client number of index they will receive
632  int nbClient = client->clientSize;
633  int* sendBuff = new int[nbClient];
634  int* recvBuff = new int[nbClient];
[1856]635
636  int* sendParticipants = new int[nbClient];
637  int* recvParticipants = new int[nbClient];
638 
[709]639  for (int i = 0; i < nbClient; ++i)
640  {
641    sendBuff[i] = 0;
642    recvBuff[i] = 0;
[1856]643    sendParticipants[i]=0 ;
644    recvParticipants[i]=0 ;
[709]645  }
[657]646  int sendBuffSize = 0;
[829]647  CClientServerMapping::GlobalIndexMap::const_iterator itbMap = globalIndexInterpSendToClient.begin(), itMap,
648                                                       iteMap = globalIndexInterpSendToClient.end();
[657]649  for (itMap = itbMap; itMap != iteMap; ++itMap)
650  {
[709]651    const std::vector<size_t>& tmp = itMap->second;
[657]652    int sizeIndex = 0, mapSize = (itMap->second).size();
653    for (int idx = 0; idx < mapSize; ++idx)
654    {
[856]655//      sizeIndex += interpMapValue.at((itMap->second)[idx]).size();
656      sizeIndex += (interpMapValue[(int)(itMap->second)[idx]]).size();
[657]657    }
658    sendBuff[itMap->first] = sizeIndex;
[1856]659    sendParticipants[itMap->first] = 1 ;
[657]660    sendBuffSize += sizeIndex;
661  }
662
[1639]663  MPI_Allreduce(sendBuff, recvBuff, nbClient, MPI_INT, MPI_SUM, client->intraComm);
[1856]664  MPI_Allreduce(sendParticipants, recvParticipants, nbClient, MPI_INT, MPI_SUM, client->intraComm);
[657]665
666  int* sendIndexDestBuff = new int [sendBuffSize];
667  int* sendIndexSrcBuff  = new int [sendBuffSize];
668  double* sendWeightBuff = new double [sendBuffSize];
669
[1639]670  std::vector<MPI_Request> sendRequest;
[709]671
672  int sendOffSet = 0, l = 0;
[657]673  for (itMap = itbMap; itMap != iteMap; ++itMap)
674  {
[709]675    const std::vector<size_t>& indexToSend = itMap->second;
676    int mapSize = indexToSend.size();
[657]677    int k = 0;
678    for (int idx = 0; idx < mapSize; ++idx)
679    {
[856]680      std::vector<std::pair<int,double> >& interpMap = interpMapValue[(int)indexToSend[idx]]; //interpMapValue.at(indexToSend[idx]);
[657]681      for (int i = 0; i < interpMap.size(); ++i)
682      {
[709]683        sendIndexDestBuff[l] = indexToSend[idx];
684        sendIndexSrcBuff[l]  = interpMap[i].first;
685        sendWeightBuff[l]    = interpMap[i].second;
[657]686        ++k;
[709]687        ++l;
[657]688      }
689    }
[1856]690   
691    sendRequest.push_back(MPI_Request());
692    MPI_Isend(&clientRank,
693             1,
694             MPI_INT,
695             itMap->first,
696             MPI_DOMAIN_INTERPOLATION_SOURCE_RANK,
697             client->intraComm,
698             &sendRequest.back());
[657]699
[1639]700    sendRequest.push_back(MPI_Request());
701    MPI_Isend(sendIndexDestBuff + sendOffSet,
[657]702             k,
[1639]703             MPI_INT,
[657]704             itMap->first,
[821]705             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
[657]706             client->intraComm,
707             &sendRequest.back());
[1639]708    sendRequest.push_back(MPI_Request());
709    MPI_Isend(sendIndexSrcBuff + sendOffSet,
[657]710             k,
[1639]711             MPI_INT,
[657]712             itMap->first,
[821]713             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
[657]714             client->intraComm,
715             &sendRequest.back());
[1639]716    sendRequest.push_back(MPI_Request());
717    MPI_Isend(sendWeightBuff + sendOffSet,
[657]718             k,
[1639]719             MPI_DOUBLE,
[657]720             itMap->first,
[821]721             MPI_DOMAIN_INTERPOLATION_WEIGHT,
[657]722             client->intraComm,
723             &sendRequest.back());
724    sendOffSet += k;
725  }
726
727  int recvBuffSize = recvBuff[clientRank];
[1856]728  int numberOfParticipants = recvParticipants[clientRank] ;
729
[657]730  int* recvIndexDestBuff = new int [recvBuffSize];
731  int* recvIndexSrcBuff  = new int [recvBuffSize];
732  double* recvWeightBuff = new double [recvBuffSize];
733  int receivedSize = 0;
734  int clientSrcRank;
[1856]735 
736 
737 // this part is done to impose a specific order for reception to retrive reproducibility
738  set<int> rankOrder ; 
739  for (int np=0 ; np < numberOfParticipants; ++np)
[657]740  {
[1639]741    MPI_Status recvStatus;
[1856]742    int rank ;
743    MPI_Recv(&rank,
744             1,
745             MPI_INT,
746             MPI_ANY_SOURCE,
747             MPI_DOMAIN_INTERPOLATION_SOURCE_RANK,
748             client->intraComm,
749             &recvStatus);
750    rankOrder.insert(rank) ;
751  }
752 
753  for (auto  clientSrcRank : rankOrder)
754  {
755    MPI_Status recvStatus;
[657]756    MPI_Recv((recvIndexDestBuff + receivedSize),
757             recvBuffSize,
[1639]758             MPI_INT,
[1856]759             clientSrcRank,
[821]760             MPI_DOMAIN_INTERPOLATION_DEST_INDEX,
[657]761             client->intraComm,
762             &recvStatus);
763
764    int countBuff = 0;
[1639]765    MPI_Get_count(&recvStatus, MPI_INT, &countBuff);
[1856]766//    clientSrcRank = recvStatus.MPI_SOURCE;
[657]767
[1639]768    MPI_Recv((recvIndexSrcBuff + receivedSize),
[657]769             recvBuffSize,
[1639]770             MPI_INT,
[657]771             clientSrcRank,
[821]772             MPI_DOMAIN_INTERPOLATION_SRC_INDEX,
[657]773             client->intraComm,
774             &recvStatus);
775
[1639]776    MPI_Recv((recvWeightBuff + receivedSize),
[657]777             recvBuffSize,
[1639]778             MPI_DOUBLE,
[657]779             clientSrcRank,
[821]780             MPI_DOMAIN_INTERPOLATION_WEIGHT,
[657]781             client->intraComm,
782             &recvStatus);
783
784    for (int idx = 0; idx < countBuff; ++idx)
785    {
[827]786      transMap[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvIndexSrcBuff + receivedSize + idx));
787      transWeight[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvWeightBuff + receivedSize + idx));
[657]788    }
789    receivedSize += countBuff;
790  }
[1856]791 
[709]792  MPI_Waitall(sendRequest.size(), &sendRequest[0], MPI_STATUS_IGNORE);
[657]793
[1856]794  delete [] sendParticipants ;
795  delete [] recvParticipants ;
[657]796  delete [] sendIndexDestBuff;
797  delete [] sendIndexSrcBuff;
798  delete [] sendWeightBuff;
799  delete [] recvIndexDestBuff;
800  delete [] recvIndexSrcBuff;
801  delete [] recvWeightBuff;
802  delete [] sendBuff;
803  delete [] recvBuff;
804}
[1622]805CATCH
[1021]806 
807/*! Redefined some functions of CONetCDF4 to make use of them */
[1639]808CDomainAlgorithmInterpolate::WriteNetCdf::WriteNetCdf(const StdString& filename, const MPI_Comm comm)
[1158]809  : CNc4DataOutput(NULL, filename, false, false, true, comm, false, true) {}
[1021]810int CDomainAlgorithmInterpolate::WriteNetCdf::addDimensionWrite(const StdString& name, 
811                                                                const StdSize size)
[1622]812TRY
[1021]813{
[1158]814  return CONetCDF4::addDimension(name, size); 
[1021]815}
[1622]816CATCH
[657]817
[1021]818int CDomainAlgorithmInterpolate::WriteNetCdf::addVariableWrite(const StdString& name, nc_type type,
819                                                               const std::vector<StdString>& dim)
[1622]820TRY
[1021]821{
[1158]822  return CONetCDF4::addVariable(name, type, dim);
[1021]823}
[1622]824CATCH
[1021]825
[1158]826void CDomainAlgorithmInterpolate::WriteNetCdf::endDefinition()
[1622]827TRY
[1158]828{
829  CONetCDF4::definition_end();
830}
[1622]831CATCH
[1158]832
[1021]833void CDomainAlgorithmInterpolate::WriteNetCdf::writeDataIndex(const CArray<int,1>& data, const StdString& name,
834                                                              bool collective, StdSize record,
835                                                              const std::vector<StdSize>* start,
836                                                              const std::vector<StdSize>* count)
[1622]837TRY
[1021]838{
839  CONetCDF4::writeData<int,1>(data, name, collective, record, start, count);
840}
[1622]841CATCH
[1021]842
843void CDomainAlgorithmInterpolate::WriteNetCdf::writeDataIndex(const CArray<double,1>& data, const StdString& name,
844                                                              bool collective, StdSize record,
845                                                              const std::vector<StdSize>* start,
846                                                              const std::vector<StdSize>* count)
[1622]847TRY
[1021]848{
849  CONetCDF4::writeData<double,1>(data, name, collective, record, start, count);
850}
[1622]851CATCH
[1021]852
853/*
854   Write interpolation weights into a file
855   \param [in] filename name of output file
856   \param interpMapValue mapping of global index of domain destination and domain source as well as the corresponding weight
857*/
858void CDomainAlgorithmInterpolate::writeInterpolationInfo(std::string& filename,
859                                                         std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
[1622]860TRY
[1021]861{
862  CContext* context = CContext::getCurrent();
863  CContextClient* client=context->client;
864
865  size_t n_src = domainSrc_->ni_glo * domainSrc_->nj_glo;
866  size_t n_dst = domainDest_->ni_glo * domainDest_->nj_glo;
867
868  long localNbWeight = 0;
869  long globalNbWeight;
870  long startIndex;
871  typedef std::map<int,std::vector<std::pair<int,double> > > IndexRemap;
872  IndexRemap::iterator itb = interpMapValue.begin(), it,
873                       ite = interpMapValue.end();
874  for (it = itb; it!=ite; ++it)
875  {
876    localNbWeight += (it->second).size();
877  }
878
879  CArray<int,1> src_idx(localNbWeight);
880  CArray<int,1> dst_idx(localNbWeight);
881  CArray<double,1> weights(localNbWeight);
882
883  int index = 0;
[1269]884  int indexOffset=0 ;
885  if (fortranConvention) indexOffset=1 ;
[1021]886  for (it = itb; it !=ite; ++it)
887  {
888    std::vector<std::pair<int,double> >& tmp = it->second;
889    for (int idx = 0; idx < tmp.size(); ++idx)
890    {
[1269]891      dst_idx(index) = it->first + indexOffset;
892      src_idx(index) = tmp[idx].first + indexOffset;
[1021]893      weights(index) = tmp[idx].second;
894      ++index;
895    }   
896  }
897
[1639]898  MPI_Allreduce(&localNbWeight, &globalNbWeight, 1, MPI_LONG, MPI_SUM, client->intraComm);
899  MPI_Scan(&localNbWeight, &startIndex, 1, MPI_LONG, MPI_SUM, client->intraComm);
[1021]900 
[1158]901  if (0 == globalNbWeight)
902  {
903    info << "There is no interpolation weights calculated between "
904         << "domain source: " << domainSrc_->getDomainOutputName()
905         << " and domain destination: " << domainDest_->getDomainOutputName()
906         << std::endl;
907    return;
908  }
909
[1021]910  std::vector<StdSize> start(1, startIndex - localNbWeight);
911  std::vector<StdSize> count(1, localNbWeight);
[1158]912 
913  WriteNetCdf netCdfWriter(filename, client->intraComm); 
[1021]914
915  // Define some dimensions
916  netCdfWriter.addDimensionWrite("n_src", n_src);
917  netCdfWriter.addDimensionWrite("n_dst", n_dst);
918  netCdfWriter.addDimensionWrite("n_weight", globalNbWeight);
919 
920  std::vector<StdString> dims(1,"n_weight");
921
922  // Add some variables
923  netCdfWriter.addVariableWrite("src_idx", NC_INT, dims);
924  netCdfWriter.addVariableWrite("dst_idx", NC_INT, dims);
925  netCdfWriter.addVariableWrite("weight", NC_DOUBLE, dims);
926
[1158]927  // End of definition
928  netCdfWriter.endDefinition();
929
[1021]930  // // Write variables
[1158]931  if (0 != localNbWeight)
932  {
933    netCdfWriter.writeDataIndex(src_idx, "src_idx", false, 0, &start, &count);
934    netCdfWriter.writeDataIndex(dst_idx, "dst_idx", false, 0, &start, &count);
935    netCdfWriter.writeDataIndex(weights, "weight", false, 0, &start, &count);
936  }
[1021]937
938  netCdfWriter.closeFile();
939}
[1622]940CATCH
[1021]941
[657]942/*!
943  Read interpolation information from a file
944  \param [in] filename interpolation file
945  \param [in/out] interpMapValue Mapping between (global) index of domain on grid destination and
946         corresponding global index of domain and associated weight value on grid source
947*/
[689]948void CDomainAlgorithmInterpolate::readInterpolationInfo(std::string& filename,
[709]949                                                        std::map<int,std::vector<std::pair<int,double> > >& interpMapValue)
[1622]950TRY
[657]951{
[660]952  int ncid ;
953  int weightDimId ;
954  size_t nbWeightGlo ;
[657]955
[1268]956
[660]957  CContext* context = CContext::getCurrent();
958  CContextClient* client=context->client;
959  int clientRank = client->clientRank;
960  int clientSize = client->clientSize;
[663]961
[1268]962
963  {
964    ifstream f(filename.c_str());
965    if (!f.good()) ERROR("void CDomainAlgorithmInterpolate::readInterpolationInfo",
966                      << "Attempt to read file weight :"  << filename << " which doesn't seem to exist." << std::endl
967                      << "Please check this file ");
968  }
969                 
[660]970  nc_open(filename.c_str(),NC_NOWRITE, &ncid) ;
971  nc_inq_dimid(ncid,"n_weight",&weightDimId) ;
972  nc_inq_dimlen(ncid,weightDimId,&nbWeightGlo) ;
973
974  size_t nbWeight ;
975  size_t start ;
976  size_t div = nbWeightGlo/clientSize ;
977  size_t mod = nbWeightGlo%clientSize ;
978  if (clientRank < mod)
979  {
980    nbWeight=div+1 ;
981    start=clientRank*(div+1) ;
982  }
983  else
984  {
985    nbWeight=div ;
986    start= mod * (div+1) + (clientRank-mod) * div ;
987  }
988
989  double* weight=new double[nbWeight] ;
990  int weightId ;
991  nc_inq_varid (ncid, "weight", &weightId) ;
992  nc_get_vara_double(ncid, weightId, &start, &nbWeight, weight) ;
993
[663]994  long* srcIndex=new long[nbWeight] ;
[660]995  int srcIndexId ;
996  nc_inq_varid (ncid, "src_idx", &srcIndexId) ;
997  nc_get_vara_long(ncid, srcIndexId, &start, &nbWeight, srcIndex) ;
998
[663]999  long* dstIndex=new long[nbWeight] ;
[660]1000  int dstIndexId ;
1001  nc_inq_varid (ncid, "dst_idx", &dstIndexId) ;
1002  nc_get_vara_long(ncid, dstIndexId, &start, &nbWeight, dstIndex) ;
[663]1003
[1269]1004  int indexOffset=0 ;
1005  if (fortranConvention) indexOffset=1 ;
1006    for(size_t ind=0; ind<nbWeight;++ind)
1007      interpMapValue[dstIndex[ind]-indexOffset].push_back(make_pair(srcIndex[ind]-indexOffset,weight[ind]));
1008 }
[1622]1009CATCH
[657]1010
[1264]1011void CDomainAlgorithmInterpolate::apply(const std::vector<std::pair<int,double> >& localIndex,
1012                                            const double* dataInput,
1013                                            CArray<double,1>& dataOut,
1014                                            std::vector<bool>& flagInitial,
1015                                            bool ignoreMissingValue, bool firstPass  )
[1622]1016TRY
[1264]1017{
1018  int nbLocalIndex = localIndex.size();   
1019  double defaultValue = std::numeric_limits<double>::quiet_NaN();
1020   
1021  if (detectMissingValue)
1022  {
1023     if (firstPass && renormalize)
1024     {
1025       renormalizationFactor.resize(dataOut.numElements()) ;
1026       renormalizationFactor=1 ;
1027     }
1028
[1480]1029    if (firstPass)
1030    {
1031      allMissing.resize(dataOut.numElements()) ;
1032      allMissing=true ;
1033    }
1034
[1264]1035    for (int idx = 0; idx < nbLocalIndex; ++idx)
1036    {
[1474]1037      if (NumTraits<double>::isNan(*(dataInput + idx)))
[1264]1038      {
[1480]1039        allMissing(localIndex[idx].first) = allMissing(localIndex[idx].first) && true;
[1264]1040        if (renormalize) renormalizationFactor(localIndex[idx].first)-=localIndex[idx].second ;
1041      }
1042      else
1043      {
1044        dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
[1480]1045        allMissing(localIndex[idx].first) = allMissing(localIndex[idx].first) && false; // Reset flag to indicate not all data source are nan
[1264]1046      }
1047    }
1048
1049  }
1050  else
1051  {
1052    for (int idx = 0; idx < nbLocalIndex; ++idx)
1053    {
1054      dataOut(localIndex[idx].first) += *(dataInput + idx) * localIndex[idx].second;
1055    }
1056  }
[657]1057}
[1622]1058CATCH
[1264]1059
1060void CDomainAlgorithmInterpolate::updateData(CArray<double,1>& dataOut)
[1622]1061TRY
[1264]1062{
[1480]1063  if (detectMissingValue)
[1273]1064  {
[1480]1065    double defaultValue = std::numeric_limits<double>::quiet_NaN();
1066    size_t nbIndex=dataOut.numElements() ; 
1067
[1615]1068    if (allMissing.numElements()>0)
[1480]1069    {
[1615]1070      for (int idx = 0; idx < nbIndex; ++idx)
1071      {
1072        if (allMissing(idx)) dataOut(idx) = defaultValue; // If all data source are nan then data destination must be nan
1073      }
[1480]1074    }
[1767]1075    else
1076    {
1077      for (int idx = 0; idx < nbIndex; ++idx)
1078      {
1079        dataOut(idx) = defaultValue; // If allMissing has no element means you have no source element so target is NaN
1080      }
1081    }
[1480]1082    if (renormalize)
1083    {
1084      if (renormalizationFactor.numElements()>0) dataOut/=renormalizationFactor ; // In some case, process doesn't received any data for interpolation (mask)
1085                                                                                 // so renormalizationFactor is not initialized
1086    }
[1273]1087  }
[1264]1088}
[1622]1089CATCH
[1264]1090
1091}
Note: See TracBrowser for help on using the repository browser.