Changeset 827 for XIOS


Ignore:
Timestamp:
03/23/16 16:10:45 (8 years ago)
Author:
mhnguyen
Message:

Implementing dynamic interpolation on axis

+) Change grid transformation to make it more flexible
+) Make some small improvements

Test
+) On Curie
+) All test pass

Location:
XIOS/trunk
Files:
30 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/inputs/Version2/iodef.xml

    r821 r827  
    1515     <field id="field_Domain_transformed_Interpolated"  operation="average" freq_op="3600s" field_ref="field_A" grid_ref="grid_Domain_tranformed_Interpolated" /> 
    1616     <field id="field_Scalar" operation="average" freq_op="3600s" grid_ref="ScalarGrid" /> 
     17     <field id="field_Value"  operation="average" freq_op="3600s" grid_ref="grid_A" /> 
    1718   </field_definition> 
    1819 
     
    5354     </axis> 
    5455     <axis id="axis_F" axis_ref="axis_C"> 
    55         <inverse_axis /> 
    56         <zoom_axis begin="0" n="4" /> 
    57         <inverse_axis /> 
     56       <interpolate_axis type="polynomial" order="3" coordinate="field_Value"/> 
     57<!--       <inverse_axis />--> 
     58<!--        <zoom_axis begin="0" n="4" />--> 
     59<!--        <inverse_axis />--> 
    5860     </axis> 
    5961   </axis_definition> 
  • XIOS/trunk/src/config/interpolate_axis_attribute.conf

    r630 r827  
    22DECLARE_ATTRIBUTE(StdString, type) 
    33DECLARE_ATTRIBUTE(int, order) 
     4DECLARE_ATTRIBUTE(StdString, coordinate) 
  • XIOS/trunk/src/filter/filter.cpp

    r639 r827  
    77    , COutputPin() 
    88    , engine(engine) 
     9    , inputSlotCount(inputSlotCount) 
    910  { /* Nothing to do */ } 
    1011 
  • XIOS/trunk/src/filter/filter.hpp

    r639 r827  
    2828    protected: 
    2929      IFilterEngine* engine; //!< The filter engine, might be the filter itself 
     30      size_t inputSlotCount; //!< Number of slot on filter 
    3031 
    3132      /*! 
  • XIOS/trunk/src/filter/spatial_transform_filter.cpp

    r790 r827  
    66namespace xios 
    77{ 
    8   CSpatialTransformFilter::CSpatialTransformFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine) 
    9     : CFilter(gc, 1, engine) 
     8  CSpatialTransformFilter::CSpatialTransformFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, size_t inputSlotsCount) 
     9    : CFilter(gc, inputSlotsCount, engine) 
    1010  { /* Nothing to do */ } 
    1111 
     
    2424      CGridTransformation* gridTransformation = destGrid->getTransformations(); 
    2525      CSpatialTransformFilterEngine* engine = CSpatialTransformFilterEngine::get(destGrid->getTransformations()); 
    26       boost::shared_ptr<CSpatialTransformFilter> filter(new CSpatialTransformFilter(gc, engine)); 
     26      const std::vector<StdString>& auxInputs = gridTransformation->getAuxInputs(); 
     27      size_t inputCount = 1 + (auxInputs.empty() ? 0 : auxInputs.size()); 
     28      boost::shared_ptr<CSpatialTransformFilter> filter(new CSpatialTransformFilter(gc, engine, inputCount)); 
    2729 
    2830      if (!lastFilter) 
     
    3234 
    3335      firstFilter = filter; 
     36      for (size_t idx = 0; idx < auxInputs.size(); ++idx) 
     37      { 
     38        CField* fieldAuxInput = CField::get(auxInputs[idx]); 
     39        fieldAuxInput->buildFilterGraph(gc, false); 
     40        fieldAuxInput->getInstantDataFilter()->connectOutput(firstFilter,idx+1); 
     41      } 
     42 
    3443      destGrid = gridTransformation->getGridSource(); 
    3544    } 
     
    7483    if (packet->status == CDataPacket::NO_ERROR) 
    7584    { 
     85      if (1 < data.size())  // Dynamical transformations 
     86      { 
     87        std::vector<CArray<double,1>* > dataAuxInputs(data.size()-1); 
     88        for (size_t idx = 0; idx < dataAuxInputs.size(); ++idx) dataAuxInputs[idx] = &(data[idx+1]->data); 
     89        gridTransformation->computeAll(dataAuxInputs); 
     90      } 
    7691      packet->data.resize(gridTransformation->getGridDestination()->storeIndex_client.numElements()); 
    7792      apply(data[0]->data, packet->data); 
  • XIOS/trunk/src/filter/spatial_transform_filter.hpp

    r644 r827  
    1111 
    1212  /*! 
    13    * A generic filter with one input slot wrapping any type of spatial transformations. 
     13   * A generic filter with multiple input slots wrapping any type of spatial transformations. 
    1414   */ 
    1515  class CSpatialTransformFilter : public CFilter 
     
    2121       * \param gc the associated garbage collector 
    2222       * \param engine the engine defining the spatial transformation 
     23       * \param [in] inputSlotsCount number of input, by default there is only one for field src 
    2324       */ 
    24       CSpatialTransformFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine); 
     25      CSpatialTransformFilter(CGarbageCollector& gc, CSpatialTransformFilterEngine* engine, size_t inputSlotsCount = 1); 
    2526 
    2627      /*! 
  • XIOS/trunk/src/node/interpolate_axis.cpp

    r630 r827  
    11#include "interpolate_axis.hpp" 
    22#include "type.hpp" 
     3#include "field.hpp" 
    34 
    45namespace xios { 
     
    2324  ENodeType CInterpolateAxis::GetType(void)    { return eInterpolateAxis; } 
    2425 
    25   void CInterpolateAxis::checkValid(CAxis* axisDest) 
     26  void CInterpolateAxis::checkValid(CAxis* axisSrc) 
    2627  { 
     28    if (this->order.isEmpty()) this->order.setValue(2); 
     29    int order = this->order.getValue(); 
     30    if (order >= axisSrc->n_glo.getValue()) 
     31    { 
     32      ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     33             << "Order of interpolation is greater than global size of axis source" 
     34             << "Size of axis source " <<axisSrc->getId() << " is " << axisSrc->n_glo.getValue()  << std::endl 
     35             << "Order of interpolation is " << order ); 
     36    } 
     37 
     38 
     39    if (!this->coordinate.isEmpty()) 
     40    { 
     41      StdString coordinate = this->coordinate.getValue(); 
     42      if (!CField::has(coordinate)) 
     43        ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     44               << "Coordinate field whose id " << coordinate << "does not exist " 
     45               << "Please define one"); 
     46    } 
    2747  } 
    2848 
     49  std::vector<StdString> CInterpolateAxis::checkAuxInputs_() 
     50  { 
     51    std::vector<StdString> auxInputs; 
     52    if (!this->coordinate.isEmpty()) 
     53    { 
     54      StdString coordinate = this->coordinate.getValue(); 
     55      if (!CField::has(coordinate)) 
     56        ERROR("CInterpolateAxis::checkValid(CAxis* axisSrc)", 
     57               << "Coordinate field whose id " << coordinate << "does not exist " 
     58               << "Please define one"); 
     59      auxInputs.push_back(coordinate); 
     60    } 
     61 
     62    return auxInputs; 
     63  } 
    2964} 
  • XIOS/trunk/src/node/interpolate_axis.hpp

    r630 r827  
    5252      virtual void checkValid(CAxis* axisDest); 
    5353 
     54      std::vector<StdString> checkAuxInputs_(); 
     55 
    5456      /// Accesseurs statiques /// 
    5557      static StdString GetName(void); 
  • XIOS/trunk/src/node/transformation.hpp

    r625 r827  
    2424      virtual void checkValid(T* dest) = 0; 
    2525 
     26      std::vector<StdString> checkAuxInputs() { return checkAuxInputs_(); } 
     27 
    2628      /// Destructeur /// 
    2729      virtual ~CTransformation(void) {} 
     30 
     31    protected: 
     32      virtual std::vector<StdString> checkAuxInputs_() { return std::vector<StdString>(); } 
    2833  }; // class CTransformation 
    2934 
  • XIOS/trunk/src/test/test_new_features.f90

    r821 r827  
    2828  DOUBLE PRECISION,DIMENSION(ni_glo,nj_glo) :: lon_glo,lat_glo 
    2929  DOUBLE PRECISION,DIMENSION(4,ni_glo,nj_glo) :: bnds_lon_glo, bnds_lat_glo 
    30   DOUBLE PRECISION :: field_A_glo(ni_glo,nj_glo,llm), lval_ni_glo(ni_glo), lval_nj_glo(nj_glo) 
     30  DOUBLE PRECISION :: field_A_glo(ni_glo,nj_glo,llm), lval_ni_glo(ni_glo), lval_nj_glo(nj_glo), field_Value_glo(ni_glo,nj_glo,llm) 
    3131  DOUBLE PRECISION,ALLOCATABLE :: lon(:,:),lat(:,:),field_A(:,:,:), field_All_Axis(:,:,:), lonvalue(:) , & 
    3232                                  field_Axis(:), lvaln(:), lval_ni(:), lval_nj(:), field_Two_Axis(:,:), lvalnInterp(:), & 
    3333                                  lontransformed(:,:), lattransformed(:,:), lon_glotransformed(:,:), lat_glotransformed(:,:), & 
    34                                   bnds_lon(:,:,:), bnds_lat(:,:,:) 
     34                                  bnds_lon(:,:,:), bnds_lat(:,:,:), field_value(:,:,:) 
    3535  INTEGER :: ni,ibegin,iend,nj,jbegin,jend, nAxis, axisBegin, axisEnd 
    3636  INTEGER :: axisterpBegin, nAxisinterp, axisinterpEnd 
     
    6868      DO l=1,llm 
    6969        field_A_glo(i,j,l)=(i-1)+(j-1)*ni_glo+10000*l 
     70        field_Value_glo(i,j,l)=l*100 
    7071      ENDDO 
    7172    ENDDO 
     
    8081 
    8182  DO j=1,llm 
    82     lval(j) = j *10 
     83    lval(j) = j *100 
    8384  ENDDO 
    8485  axisBegin = 0 
     
    109110          lvaln(nAxis), lval_ni(ni), lval_nj(nj), lvalnInterp(nAxisinterp), & 
    110111          lontransformed(niDomInterp, njDomInterp), lattransformed(niDomInterp, njDomInterp), & 
    111           bnds_lon(4,ni,nj), bnds_lat(4,ni,nj)) 
     112          bnds_lon(4,ni,nj), bnds_lat(4,ni,nj), field_value(0:ni+1,-1:nj+2,llm)) 
    112113 
    113114  ALLOCATE(mask(nj)) 
     
    126127  lattransformed(:,:) = lat_glotransformed(ibeginDomInterp+1:iendDomInterp+1,jbeginDomInterp+1:jendDomInterp+1) 
    127128  field_A(1:ni,1:nj,:) = field_A_glo(ibegin+1:iend+1,jbegin+1:jend+1,:) 
     129  field_value(1:ni,1:nj,:) = field_Value_glo(ibegin+1:iend+1,jbegin+1:jend+1,:) 
    128130  field_Axis(1:nAxis)  = field_A_glo(1,1,axisBegin+1:axisEnd+1) 
    129131  field_Two_Axis(:,1:nj)  = field_A_glo(:,jbegin+1:jend+1,1) 
     
    201203 
    202204  PRINT*,"field field_A is active ? ",xios_field_is_active("field_A") 
    203   DO ts=1,24*10 
     205  DO ts=1,24*1 
    204206    CALL xios_update_calendar(ts) 
    205207    CALL xios_send_field("field_A",field_A) 
     208    CALL xios_send_field("field_Value",field_value) 
    206209    CALL xios_send_field("field_Axis",field_Axis) 
    207210    CALL xios_send_field("field_Two_Axis",field_Two_Axis) 
  • XIOS/trunk/src/transformation/axis_algorithm_interpolate.cpp

    r668 r827  
    1212#include "context_client.hpp" 
    1313#include "utils.hpp" 
     14#include "grid.hpp" 
     15#include "distribution_client.hpp" 
    1416 
    1517namespace xios { 
    1618 
    1719CAxisAlgorithmInterpolate::CAxisAlgorithmInterpolate(CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis) 
    18 : CAxisAlgorithmTransformation(axisDestination, axisSource) 
     20: CAxisAlgorithmTransformation(axisDestination, axisSource), coordinate_(), transPosition_() 
    1921{ 
    2022  interpAxis->checkValid(axisSource); 
    2123  order_ = interpAxis->order.getValue(); 
    22   if (order_ >= axisSource->n_glo.getValue()) 
    23   { 
    24     ERROR("CAxisAlgorithmInterpolate::CAxisAlgorithmInterpolate(CAxis* axisDestination, CAxis* axisSource, CInterpolateAxis* interpAxis)", 
    25            << "Order of interpolation is greater than global size of axis source" 
    26            << "Size of axis source " <<axisSource->getId() << " is " << axisSource->n_glo.getValue()  << std::endl 
    27            << "Order of interpolation is " << order_ ); 
    28   } 
    29  
    30   computeIndexSourceMapping(); 
     24  if (!interpAxis->coordinate.isEmpty()) 
     25  { 
     26    coordinate_ = interpAxis->coordinate.getValue(); 
     27    this->idAuxInputs_.resize(1); 
     28    this->idAuxInputs_[0] = coordinate_; 
     29  } 
     30 
     31//  computeIndexSourceMapping(); 
    3132} 
    3233 
     
    3435  Compute the index mapping between axis on grid source and one on grid destination 
    3536*/ 
    36 void CAxisAlgorithmInterpolate::computeIndexSourceMapping() 
    37 { 
    38   CArray<double,1>& axisValue = axisSrc_->value; 
    39   CArray<bool,1>& axisMask = axisSrc_->mask; 
    40  
     37void CAxisAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 
     38{ 
    4139  CContext* context = CContext::getCurrent(); 
    4240  CContextClient* client=context->client; 
    4341  int nbClient = client->clientSize; 
    44  
     42  CArray<bool,1>& axisMask = axisSrc_->mask; 
    4543  int srcSize  = axisSrc_->n_glo.getValue(); 
    46   int numValue = axisValue.numElements(); 
    47  
    48   std::vector<double> recvBuff(srcSize); 
    49   std::vector<int> indexVec(srcSize); 
    50  
    51   retrieveAllAxisValue(recvBuff, indexVec); 
    52   XIOSAlgorithms::sortWithIndex<double, CVectorStorage>(recvBuff, indexVec); 
    53   computeInterpolantPoint(recvBuff, indexVec); 
     44  std::vector<CArray<double,1> > vecAxisValue; 
     45 
     46  // Fill in axis value from coordinate 
     47  fillInAxisValue(vecAxisValue, dataAuxInputs); 
     48 
     49  for (int idx = 0; idx < vecAxisValue.size(); ++idx) 
     50  { 
     51    CArray<double,1>& axisValue = vecAxisValue[idx]; 
     52    std::vector<double> recvBuff(srcSize); 
     53    std::vector<int> indexVec(srcSize); 
     54    retrieveAllAxisValue(axisValue, axisMask, recvBuff, indexVec); 
     55    XIOSAlgorithms::sortWithIndex<double, CVectorStorage>(recvBuff, indexVec); 
     56    computeInterpolantPoint(recvBuff, indexVec, idx); 
     57  } 
    5458} 
    5559 
     
    6064  \param [in] indexVec permutation index of axisValue 
    6165*/ 
    62 void CAxisAlgorithmInterpolate::computeInterpolantPoint(const std::vector<double>& axisValue, const std::vector<int>& indexVec) 
     66void CAxisAlgorithmInterpolate::computeInterpolantPoint(const std::vector<double>& axisValue, const std::vector<int>& indexVec, int transPos) 
    6367{ 
    6468  std::vector<double>::const_iterator itb = axisValue.begin(), ite = axisValue.end(); 
     
    7882    itUpperBound = std::upper_bound(itb, ite, destValue); 
    7983    if ((ite != itUpperBound) && (sfmax == *itUpperBound)) itUpperBound = ite; 
    80  
    8184 
    8285    // If the value is not in the range, that means we'll do extra-polation 
     
    118121    } 
    119122  } 
    120   computeWeightedValueAndMapping(interpolatingIndexValues); 
     123  computeWeightedValueAndMapping(interpolatingIndexValues, transPos); 
    121124} 
    122125 
     
    125128  \param [in] interpolatingIndexValues the necessary axis value to calculate the coeffs 
    126129*/ 
    127 void CAxisAlgorithmInterpolate::computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues) 
    128 { 
    129   std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
    130   std::map<int, std::vector<double> >& transWeight = this->transformationWeight_; 
     130void CAxisAlgorithmInterpolate::computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos) 
     131{ 
     132  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[transPos]; 
     133  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[transPos]; 
    131134  std::map<int, std::vector<std::pair<int,double> > >::const_iterator itb = interpolatingIndexValues.begin(), it, 
    132135                                                                      ite = interpolatingIndexValues.end(); 
     
    138141    const std::vector<std::pair<int,double> >& interpVal = it->second; 
    139142    int interpSize = interpVal.size(); 
     143    transMap[globalIndexDest].resize(interpSize); 
     144    transWeight[globalIndexDest].resize(interpSize); 
    140145    for (int idx = 0; idx < interpSize; ++idx) 
    141146    { 
     
    149154        weight /= (interpVal[idx].second - interpVal[k].second); 
    150155      } 
    151       transMap[globalIndexDest].push_back(index); 
    152       transWeight[globalIndexDest].push_back(weight); 
     156      transMap[globalIndexDest][idx] = index; 
     157      transWeight[globalIndexDest][idx] = weight; 
     158      if (!transPosition_.empty()) 
     159      { 
     160        (this->transformationPosition_[transPos])[globalIndexDest] = transPosition_[transPos]; 
     161      } 
    153162    } 
    154163  } 
     
    160169  \param [in/out] indexVec mapping between values and global index of axis 
    161170*/ 
    162 void CAxisAlgorithmInterpolate::retrieveAllAxisValue(std::vector<double>& recvBuff, std::vector<int>& indexVec) 
    163 { 
    164   CArray<double,1>& axisValue = axisSrc_->value; 
    165   CArray<bool,1>& axisMask = axisSrc_->mask; 
    166  
     171void CAxisAlgorithmInterpolate::retrieveAllAxisValue(const CArray<double,1>& axisValue, const CArray<bool,1>& axisMask, 
     172                                                     std::vector<double>& recvBuff, std::vector<int>& indexVec) 
     173{ 
    167174  CContext* context = CContext::getCurrent(); 
    168175  CContextClient* client=context->client; 
     
    230237} 
    231238 
    232 } 
     239/*! 
     240  Fill in axis value dynamically from a field whose grid is composed of a domain and an axis 
     241  \param [in/out] vecAxisValue vector axis value filled in from input field 
     242*/ 
     243void CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue, 
     244                                                const std::vector<CArray<double,1>* >& dataAuxInputs) 
     245{ 
     246  if (coordinate_.empty()) 
     247  { 
     248    vecAxisValue.resize(1); 
     249    vecAxisValue[0].resize(axisSrc_->value.numElements()); 
     250    vecAxisValue[0] = axisSrc_->value; 
     251    this->transformationMapping_.resize(1); 
     252    this->transformationWeight_.resize(1); 
     253  } 
     254  else 
     255  { 
     256    CField* field = CField::get(coordinate_); 
     257    CGrid* grid = field->grid; 
     258 
     259    std::vector<CDomain*> domListP = grid->getDomains(); 
     260    std::vector<CAxis*> axisListP = grid->getAxis(); 
     261    if (domListP.empty() || axisListP.empty() || (1 < domListP.size()) || (1 < axisListP.size())) 
     262      ERROR("CAxisAlgorithmInterpolate::fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue)", 
     263             << "XIOS only supports dynamic interpolation with coordinate (field) associated with grid composed of a domain and an axis" 
     264             << "Coordinate (field) id = " <<field->getId() << std::endl 
     265             << "Associated grid id = " << grid->getId()); 
     266 
     267    CDomain* dom = domListP[0]; 
     268    size_t vecAxisValueSize = dom->i_index.numElements(); 
     269    vecAxisValue.resize(vecAxisValueSize); 
     270    if (transPosition_.empty()) 
     271    { 
     272      transPosition_.resize(vecAxisValueSize); 
     273      for (size_t idx = 0; idx < vecAxisValueSize; ++idx) 
     274      { 
     275        transPosition_[idx].resize(2); 
     276        transPosition_[idx][0] = (dom->i_index)(idx); 
     277        transPosition_[idx][1] = (dom->j_index)(idx); 
     278      } 
     279    } 
     280    this->transformationMapping_.resize(vecAxisValueSize); 
     281    this->transformationWeight_.resize(vecAxisValueSize); 
     282    this->transformationPosition_.resize(vecAxisValueSize); 
     283 
     284    const std::vector<size_t>& globalIndexSendToServer = grid->getDistributionClient()->getGlobalDataIndexSendToServer(); 
     285    const std::vector<int>& localDataSendToServer = grid->getDistributionClient()->getLocalDataIndexSendToServer(); 
     286    std::vector<int> permutIndex(globalIndexSendToServer.size()); 
     287    std::vector<int>::iterator itVec; 
     288    XIOSAlgorithms::fillInIndex(globalIndexSendToServer.size(), permutIndex); 
     289    XIOSAlgorithms::sortWithIndex<size_t, CVectorStorage>(globalIndexSendToServer, permutIndex); 
     290    size_t axisSrcSize = axisSrc_->index.numElements(); 
     291    std::vector<int> globalDimension = grid->getGlobalDimension(); 
     292    XIOSBinarySearchWithIndex<size_t> binSearch(globalIndexSendToServer); 
     293    for (size_t idx = 0; idx < vecAxisValueSize; ++idx) 
     294    { 
     295      size_t axisValueSize = 0; 
     296      for (size_t jdx = 0; jdx < axisSrcSize; ++jdx) 
     297      { 
     298        size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1]; 
     299        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndex, itVec)) 
     300        { 
     301          ++axisValueSize; 
     302        } 
     303      } 
     304 
     305      vecAxisValue[idx].resize(axisValueSize); 
     306      axisValueSize = 0; 
     307      for (size_t jdx = 0; jdx < axisSrcSize; ++jdx) 
     308      { 
     309        size_t globalIndex = ((dom->i_index)(idx) + (dom->j_index)(idx)*globalDimension[0]) + (axisSrc_->index)(jdx)*globalDimension[0]*globalDimension[1]; 
     310        if (binSearch.search(permutIndex.begin(), permutIndex.end(), globalIndex, itVec)) 
     311        { 
     312          vecAxisValue[idx](axisValueSize) = (*dataAuxInputs[0])(localDataSendToServer[*itVec]); 
     313          ++axisValueSize; 
     314        } 
     315      } 
     316    } 
     317  } 
     318} 
     319 
     320} 
  • XIOS/trunk/src/transformation/axis_algorithm_interpolate.hpp

    r630 r827  
    2727  virtual ~CAxisAlgorithmInterpolate() {} 
    2828 
    29   virtual void computeIndexSourceMapping(); 
     29protected: 
     30  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); 
    3031 
    3132private: 
    32   void retrieveAllAxisValue(std::vector<double>& recvBuff, std::vector<int>& indexVec); 
    33   void computeInterpolantPoint(const std::vector<double>& recvBuff, const std::vector<int>& indexVec); 
    34   void computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues); 
     33  void retrieveAllAxisValue(const CArray<double,1>& axisValue, const CArray<bool,1>& axisMask, 
     34                            std::vector<double>& recvBuff, std::vector<int>& indexVec); 
     35  void computeInterpolantPoint(const std::vector<double>& recvBuff, const std::vector<int>& indexVec, int transPos = 0); 
     36  void computeWeightedValueAndMapping(const std::map<int, std::vector<std::pair<int,double> > >& interpolatingIndexValues, int transPos = 0); 
     37  void fillInAxisValue(std::vector<CArray<double,1> >& vecAxisValue, 
     38                       const std::vector<CArray<double,1>* >& dataAuxInputs); 
    3539 
    3640private: 
    3741  // Interpolation order 
    3842  int order_; 
     43  StdString coordinate_; 
     44  std::vector<std::vector<int> > transPosition_; 
    3945}; 
    4046 
  • XIOS/trunk/src/transformation/axis_algorithm_inverse.cpp

    r666 r827  
    2525  } 
    2626 
    27   this->computeIndexSourceMapping(); 
     27//  this->computeIndexSourceMapping(); 
     28} 
     29 
     30void CAxisAlgorithmInverse::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 
     31{ 
     32  this->transformationMapping_.resize(1); 
     33  this->transformationWeight_.resize(1); 
     34 
     35  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[0]; 
     36  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[0]; 
     37 
     38  int globalIndexSize = axisDestGlobalIndex_.size(); 
     39  for (int idx = 0; idx < globalIndexSize; ++idx) 
     40  { 
     41    transMap[axisDestGlobalIndex_[idx]].push_back(axisDestGlobalSize_-axisDestGlobalIndex_[idx]-1); 
     42    transWeight[axisDestGlobalIndex_[idx]].push_back(1.0); 
     43  } 
     44 
    2845  int niSrc   = axisSrc_->n.getValue(); 
    2946  int sizeSrc = axisSrc_->n_glo.getValue(); 
     
    3552      axisDest_->value(idx) = axisSrc_->value(sizeSrc-idx-1); 
    3653    } 
    37   } 
    38 } 
    39  
    40 void CAxisAlgorithmInverse::computeIndexSourceMapping() 
    41 { 
    42   std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
    43   std::map<int, std::vector<double> >& transWeight = this->transformationWeight_; 
    44  
    45   int globalIndexSize = axisDestGlobalIndex_.size(); 
    46   for (int idx = 0; idx < globalIndexSize; ++idx) 
    47   { 
    48     transMap[axisDestGlobalIndex_[idx]].push_back(axisDestGlobalSize_-axisDestGlobalIndex_[idx]-1); 
    49     transWeight[axisDestGlobalIndex_[idx]].push_back(1.0); 
    5054  } 
    5155} 
     
    6569  CTransformationMapping transformationMap(axisDest_, axisSrc_); 
    6670 
    67   std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
    68   std::map<int, std::vector<double> >& transWeight = this->transformationWeight_; 
     71  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[0]; 
     72  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[0]; 
    6973 
    7074  std::map<size_t, std::vector<std::pair<size_t,double> > > globaIndexMapFromDestToSource; 
  • XIOS/trunk/src/transformation/axis_algorithm_inverse.hpp

    r630 r827  
    2525  virtual ~CAxisAlgorithmInverse() {} 
    2626 
    27   virtual void computeIndexSourceMapping(); 
     27protected: 
     28  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); 
    2829 
    2930private: 
  • XIOS/trunk/src/transformation/axis_algorithm_transformation.cpp

    r668 r827  
    2929} 
    3030 
    31 void CAxisAlgorithmTransformation::computeIndexSourceMapping() 
     31void CAxisAlgorithmTransformation::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 
    3232{ 
    3333} 
     
    4545void CAxisAlgorithmTransformation::computeGlobalGridIndexFromGlobalIndexElement(int axisDestGlobalIndex, 
    4646                                                                          const std::vector<int>& axisSrcGlobalIndex, 
     47                                                                          const std::vector<int>& destGlobalIndexPositionInGrid, 
    4748                                                                          int axisPositionInGrid, 
    4849                                                                          const std::vector<int>& gridDestGlobalDim, 
     
    5253                                                                          std::vector<std::vector<size_t> >& globalIndexSrcGrid) 
    5354{ 
     55  bool hasDestGlobalIndexPos = !destGlobalIndexPositionInGrid.empty(); 
    5456  int globalDim = gridDestGlobalDim.size(); 
    5557 
     
    6264  { 
    6365    if (axisPositionInGrid == i) gridAxisGlobalDim[i] = 1; 
    64     else gridAxisGlobalDim[i] = gridDestGlobalDim[i]; 
     66    else 
     67    { 
     68      if (!hasDestGlobalIndexPos) gridAxisGlobalDim[i] = gridDestGlobalDim[i]; 
     69      else gridAxisGlobalDim[i] = 1; 
     70    } 
    6571    ssize *= gridAxisGlobalDim[i]; 
    6672  } 
     
    8086    } 
    8187 
    82     for (int i = 0; i < globalDim; ++i) currentIndex[i] = idxLoop[i]; 
     88    int j = 0; 
     89    for (int i = 0; i < globalDim; ++i) 
     90    { 
     91      if (!hasDestGlobalIndexPos) currentIndex[i] = idxLoop[i]; 
     92      else 
     93      { 
     94        if (axisPositionInGrid == i) currentIndex[i] = axisDestGlobalIndex; 
     95        else 
     96        { 
     97          currentIndex[i] = destGlobalIndexPositionInGrid[j]; 
     98          ++j; 
     99        } 
     100      } 
     101    } 
     102 
    83103    currentIndex[axisPositionInGrid] = axisDestGlobalIndex; 
    84104 
     
    118138    } 
    119139 
    120     for (int i = 0; i < globalDim; ++i) currentIndex[i] = idxLoop[i]; 
     140    int j = 0; 
     141    for (int i = 0; i < globalDim; ++i) 
     142    { 
     143      if (!hasDestGlobalIndexPos) currentIndex[i] = idxLoop[i]; 
     144      else 
     145      { 
     146        if (axisPositionInGrid == i) currentIndex[i] = axisDestGlobalIndex; 
     147        else 
     148        { 
     149          currentIndex[i] = destGlobalIndexPositionInGrid[j]; 
     150          ++j; 
     151        } 
     152      } 
     153    } 
    121154    currentIndex[axisPositionInGrid] = axisDestGlobalIndex; 
    122155 
     
    150183    ++idx; 
    151184  } 
    152 } 
    153 } 
     185 
     186 
     187 
     188//  int globalDim = gridDestGlobalDim.size(); 
     189// 
     190//  std::vector<size_t> currentIndex(globalDim); 
     191//  std::vector<int> gridAxisGlobalDim(globalDim); 
     192//  std::vector<int> idxLoop(globalDim, 0); 
     193// 
     194//  size_t ssize = 1, idx = 0, realGlobalIndexSize = 0; 
     195//  for (int i = 0; i< globalDim; ++i) 
     196//  { 
     197//    if (axisPositionInGrid == i) gridAxisGlobalDim[i] = 1; 
     198//    else gridAxisGlobalDim[i] = gridDestGlobalDim[i]; 
     199//    ssize *= gridAxisGlobalDim[i]; 
     200//  } 
     201// 
     202//  std::vector<size_t>::const_iterator itbArr = globalIndexGridDestSendToServer.begin(), itArr, 
     203//                                      iteArr = globalIndexGridDestSendToServer.end(); 
     204// 
     205//  while (idx < ssize) 
     206//  { 
     207//    for (int i = 0; i < globalDim-1; ++i) 
     208//    { 
     209//      if (idxLoop[i] == gridAxisGlobalDim[i]) 
     210//      { 
     211//        idxLoop[i] = 0; 
     212//        ++idxLoop[i+1]; 
     213//      } 
     214//    } 
     215// 
     216//    for (int i = 0; i < globalDim; ++i) currentIndex[i] = idxLoop[i]; 
     217//    currentIndex[axisPositionInGrid] = axisDestGlobalIndex; 
     218// 
     219//    size_t globIndex = currentIndex[0]; 
     220//    size_t mulDim = 1; 
     221//    for (int k = 1; k < globalDim; ++k) 
     222//    { 
     223//      mulDim *= gridDestGlobalDim[k-1]; 
     224//      globIndex += (currentIndex[k])*mulDim; 
     225//    } 
     226// 
     227//    if (std::binary_search(itbArr, iteArr, globIndex)) ++realGlobalIndexSize; 
     228//    ++idxLoop[0]; 
     229//    ++idx; 
     230//  } 
     231// 
     232//  if (globalIndexDestGrid.numElements() != realGlobalIndexSize) 
     233//    globalIndexDestGrid.resize(realGlobalIndexSize); 
     234// 
     235//  if (realGlobalIndexSize != globalIndexSrcGrid.size()) globalIndexSrcGrid.resize(realGlobalIndexSize); 
     236//  for (int i = 0; i < globalIndexSrcGrid.size(); ++i) 
     237//    if (globalIndexSrcGrid[i].size() != axisSrcGlobalIndex.size()) 
     238//      globalIndexSrcGrid[i].resize(axisSrcGlobalIndex.size()); 
     239// 
     240//  size_t realGlobalIndex = 0; 
     241//  idx = 0; 
     242//  idxLoop.assign(globalDim, 0); 
     243//  while (idx < ssize) 
     244//  { 
     245//    for (int i = 0; i < globalDim-1; ++i) 
     246//    { 
     247//      if (idxLoop[i] == gridAxisGlobalDim[i]) 
     248//      { 
     249//        idxLoop[i] = 0; 
     250//        ++idxLoop[i+1]; 
     251//      } 
     252//    } 
     253// 
     254//    for (int i = 0; i < globalDim; ++i) currentIndex[i] = idxLoop[i]; 
     255//    currentIndex[axisPositionInGrid] = axisDestGlobalIndex; 
     256// 
     257//    size_t globIndex = currentIndex[0]; 
     258//    size_t mulDim = 1; 
     259//    for (int k = 1; k < globalDim; ++k) 
     260//    { 
     261//      mulDim *= gridDestGlobalDim[k-1]; 
     262//      globIndex += (currentIndex[k])*mulDim; 
     263//    } 
     264// 
     265//    if (std::binary_search(itbArr, iteArr, globIndex)) 
     266//    { 
     267//      globalIndexDestGrid(realGlobalIndex) = globIndex; 
     268//      for (int i = 0; i < globalIndexSrcGrid[realGlobalIndex].size(); ++i) 
     269//      { 
     270//        currentIndex[axisPositionInGrid] = axisSrcGlobalIndex[i]; 
     271//        globIndex = currentIndex[0]; 
     272//        mulDim = 1; 
     273//        for (int k = 1; k < globalDim; ++k) 
     274//        { 
     275//          mulDim *= gridDestGlobalDim[k-1]; 
     276//          globIndex += (currentIndex[k])*mulDim; 
     277//        } 
     278//        (globalIndexSrcGrid[realGlobalIndex])[i] = globIndex; 
     279//      } 
     280//      ++realGlobalIndex; 
     281//    } 
     282// 
     283//    ++idxLoop[0]; 
     284//    ++idx; 
     285//  } 
     286} 
     287} 
  • XIOS/trunk/src/transformation/axis_algorithm_transformation.hpp

    r668 r827  
    2929  virtual void computeGlobalGridIndexFromGlobalIndexElement(int axisDestGlobalIndex, 
    3030                                                        const std::vector<int>& axisSrcGlobalIndex, 
     31                                                        const std::vector<int>& destGlobalIndexPositionInGrid, 
    3132                                                        int axisPositionInGrid, 
    3233                                                        const std::vector<int>& gridDestGlobalDim, 
     
    3536                                                        CArray<size_t,1>& globalIndexDestGrid, 
    3637                                                        std::vector<std::vector<size_t> >& globalIndexSrcGrid); 
    37   void computeIndexSourceMapping(); 
     38 
     39  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); 
    3840 
    3941protected: 
  • XIOS/trunk/src/transformation/axis_algorithm_zoom.cpp

    r821 r827  
    2727  } 
    2828 
    29   computeIndexSourceMapping(); 
     29//  computeIndexSourceMapping(); 
    3030} 
    3131 
     
    3333  Compute the index mapping between axis on grid source and one on grid destination 
    3434*/ 
    35 void CAxisAlgorithmZoom::computeIndexSourceMapping() 
     35void CAxisAlgorithmZoom::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 
    3636{ 
    3737  StdSize niSource = axisSrc_->n.getValue(); 
     
    4444  if (iend < ibegin) ni = 0; 
    4545 
    46   std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
    47   std::map<int, std::vector<double> >& transWeight = this->transformationWeight_; 
     46  this->transformationMapping_.resize(1); 
     47  this->transformationWeight_.resize(1); 
     48 
     49  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[0]; 
     50  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[0]; 
     51 
    4852  for (StdSize idx = 0; idx < ni; ++idx) 
    4953  { 
     
    7579  StdSize iBeginMask = axisDest_->begin.getValue(); 
    7680  StdSize globalIndexMask = 0; 
    77   std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
     81  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[0]; 
    7882  std::map<int, std::vector<int> >::const_iterator ite = (transMap).end(); 
    7983  for (StdSize idx = 0; idx < niMask; ++idx) 
  • XIOS/trunk/src/transformation/axis_algorithm_zoom.hpp

    r630 r827  
    2828  virtual ~CAxisAlgorithmZoom() {} 
    2929 
    30   virtual void computeIndexSourceMapping(); 
     30protected: 
     31  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); 
    3132 
    3233private: 
  • XIOS/trunk/src/transformation/domain_algorithm_generate_rectilinear.cpp

    r775 r827  
    3333  Compute the index mapping between domain on grid source and one on grid destination 
    3434*/ 
    35 void CDomainAlgorithmGenerateRectilinear::computeIndexSourceMapping() 
     35void CDomainAlgorithmGenerateRectilinear::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 
    3636{ 
    3737 /* Nothing to do */ 
  • XIOS/trunk/src/transformation/domain_algorithm_generate_rectilinear.hpp

    r775 r827  
    3333  virtual ~CDomainAlgorithmGenerateRectilinear() {} 
    3434 
    35   virtual void computeIndexSourceMapping(); 
     35protected: 
     36  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); 
    3637 
    3738private: 
  • XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

    r821 r827  
    2323{ 
    2424  interpDomain_->checkValid(domainSource); 
    25   computeIndexSourceMapping(); 
     25//  computeIndexSourceMapping(); 
    2626} 
    2727 
     
    351351  Compute the index mapping between domain on grid source and one on grid destination 
    352352*/ 
    353 void CDomainAlgorithmInterpolate::computeIndexSourceMapping() 
     353void CDomainAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 
    354354{ 
    355355  if (!interpDomain_->file.isEmpty()) 
     
    381381  CContextClient* client=context->client; 
    382382  int clientRank = client->clientRank; 
     383 
     384  this->transformationMapping_.resize(1); 
     385  this->transformationWeight_.resize(1); 
     386 
     387  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[0]; 
     388  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[0]; 
    383389 
    384390  boost::unordered_map<size_t,int> globalIndexOfDomainDest; 
     
    531537    for (int idx = 0; idx < countBuff; ++idx) 
    532538    { 
    533       transformationMapping_[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvIndexSrcBuff + receivedSize + idx)); 
    534       transformationWeight_[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvWeightBuff + receivedSize + idx)); 
     539      transMap[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvIndexSrcBuff + receivedSize + idx)); 
     540      transWeight[*(recvIndexDestBuff + receivedSize + idx)].push_back(*(recvWeightBuff + receivedSize + idx)); 
    535541    } 
    536542    receivedSize += countBuff; 
  • XIOS/trunk/src/transformation/domain_algorithm_interpolate.hpp

    r743 r827  
    2626  virtual ~CDomainAlgorithmInterpolate() {} 
    2727 
    28   virtual void computeIndexSourceMapping(); 
     28protected: 
     29  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); 
    2930 
    3031private: 
  • XIOS/trunk/src/transformation/domain_algorithm_transformation.cpp

    r668 r827  
    2121} 
    2222 
    23 void CDomainAlgorithmTransformation::computeIndexSourceMapping() 
     23void CDomainAlgorithmTransformation::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 
    2424{ 
    2525} 
     
    3838void CDomainAlgorithmTransformation::computeGlobalGridIndexFromGlobalIndexElement(int domainDestGlobalIndex, 
    3939                                                                          const std::vector<int>& domainSrcGlobalIndex, 
     40                                                                          const std::vector<int>& destGlobalIndexPositionInGrid, 
    4041                                                                          int domainPositionInGrid, 
    4142                                                                          const std::vector<int>& gridDestGlobalDim, 
  • XIOS/trunk/src/transformation/domain_algorithm_transformation.hpp

    r668 r827  
    2929  virtual void computeGlobalGridIndexFromGlobalIndexElement(int domainDestGlobalIndex, 
    3030                                                        const std::vector<int>& domainSrcGlobalIndex, 
     31                                                        const std::vector<int>& destGlobalIndexPositionInGrid, 
    3132                                                        int domainPositionInGrid, 
    3233                                                        const std::vector<int>& gridDestGlobalDim, 
     
    3536                                                        CArray<size_t,1>& globalIndexDestGrid, 
    3637                                                        std::vector<std::vector<size_t> >& globalIndexSrcGrid); 
    37   void computeIndexSourceMapping(); 
     38 
     39  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >&); 
    3840 
    3941protected: 
  • XIOS/trunk/src/transformation/domain_algorithm_zoom.cpp

    r787 r827  
    4040  } 
    4141 
    42   computeIndexSourceMapping(); 
     42//  computeIndexSourceMapping(); 
    4343} 
    4444 
     
    4646  Compute the index mapping between domain on grid source and one on grid destination 
    4747*/ 
    48 void CDomainAlgorithmZoom::computeIndexSourceMapping() 
     48void CDomainAlgorithmZoom::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 
    4949{ 
    5050  int niSource = domainSrc_->ni.getValue(); 
     
    6868  int niGlob = domainSrc_->ni_glo.getValue(); 
    6969  int njGlob = domainSrc_->nj_glo.getValue(); 
    70   std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
    71   std::map<int, std::vector<double> >& transWeight = this->transformationWeight_; 
     70 
     71  this->transformationMapping_.resize(1); 
     72  this->transformationWeight_.resize(1); 
     73 
     74  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[0]; 
     75  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_[0]; 
     76 
     77//  std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
     78//  std::map<int, std::vector<double> >& transWeight = this->transformationWeight_; 
    7279  int domainGlobalIndex; 
    7380  for (int j = 0; j < nj; ++j) 
     
    110117  int globalIndexMask = 0; 
    111118 
    112   std::map<int, std::vector<int> >& transMap = this->transformationMapping_; 
     119  std::map<int, std::vector<int> >& transMap = this->transformationMapping_[0]; 
    113120  std::map<int, std::vector<int> >::const_iterator ite = (transMap).end(); 
    114121  for (int j = 0; j < njMask; ++j) 
  • XIOS/trunk/src/transformation/domain_algorithm_zoom.hpp

    r631 r827  
    2828  virtual ~CDomainAlgorithmZoom() {} 
    2929 
    30   virtual void computeIndexSourceMapping(); 
     30protected: 
     31  void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs); 
    3132 
    3233private: 
  • XIOS/trunk/src/transformation/generic_algorithm_transformation.cpp

    r709 r827  
    1212 
    1313CGenericAlgorithmTransformation::CGenericAlgorithmTransformation() 
    14  : transformationMapping_(), transformationWeight_() 
     14 : transformationMapping_(), transformationWeight_(), transformationPosition_(), idAuxInputs_() 
    1515{ 
    1616} 
     
    3232                                                             std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexWeightFromDestToSource) 
    3333{ 
    34   std::map<int, std::vector<int> >::const_iterator itbTransMap = transformationMapping_.begin(), itTransMap, 
    35                                                    iteTransMap = transformationMapping_.end(); 
    36   std::map<int, std::vector<double> >::const_iterator itTransWeight = transformationWeight_.begin(); 
    37   std::vector<std::vector<size_t> > globalIndexSrcGrid; 
    38   CArray<size_t,1> globalIndexDestGrid; 
     34  bool isTransPosEmpty = transformationPosition_.empty(); 
    3935  std::vector<size_t> vecGlobalIndexGridSendToServer(globalIndexGridDestSendToServer.begin(), globalIndexGridDestSendToServer.end()); 
    4036  std::sort(vecGlobalIndexGridSendToServer.begin(), vecGlobalIndexGridSendToServer.end()); 
    41   for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight) 
     37  for (size_t idxTrans = 0; idxTrans < transformationMapping_.size(); ++idxTrans) 
    4238  { 
    43     this->computeGlobalGridIndexFromGlobalIndexElement(itTransMap->first, 
    44                                                    itTransMap->second, 
    45                                                    elementPositionInGrid, 
    46                                                    gridDestGlobalDim, 
    47                                                    gridSrcGlobalDim, 
    48                                                    vecGlobalIndexGridSendToServer, 
    49                                                    globalIndexDestGrid, 
    50                                                    globalIndexSrcGrid); 
    51     size_t globalIndexSize = globalIndexDestGrid.numElements(); 
    52     const std::vector<double>& currentVecWeight = itTransWeight->second; 
    53     for (size_t idx = 0; idx < globalIndexSize; ++idx) 
     39    std::map<int, std::vector<int> >::const_iterator itbTransMap = transformationMapping_[idxTrans].begin(), itTransMap, 
     40                                                     iteTransMap = transformationMapping_[idxTrans].end(); 
     41    std::map<int, std::vector<double> >::const_iterator itTransWeight = transformationWeight_[idxTrans].begin(); 
     42 
     43    // If transformation position exists 
     44    std::map<int, std::vector<int> >::const_iterator itTransPos, iteTransPos; 
     45    if (!isTransPosEmpty) 
    5446    { 
    55       size_t globalIndex = globalIndexDestGrid(idx); 
    56       for (int i = 0; i < globalIndexSrcGrid[idx].size(); ++i) 
     47      itTransPos  = transformationPosition_[idxTrans].begin(), 
     48      iteTransPos = transformationPosition_[idxTrans].end(); 
     49    } 
     50    std::vector<int> emptyTransPos; 
     51 
     52    std::vector<std::vector<size_t> > globalIndexSrcGrid; 
     53    CArray<size_t,1> globalIndexDestGrid; 
     54//    std::vector<size_t> vecGlobalIndexGridSendToServer(globalIndexGridDestSendToServer.begin(), globalIndexGridDestSendToServer.end()); 
     55//    std::sort(vecGlobalIndexGridSendToServer.begin(), vecGlobalIndexGridSendToServer.end()); 
     56    for (itTransMap = itbTransMap; itTransMap != iteTransMap; ++itTransMap, ++itTransWeight) 
     57    { 
     58      if (!isTransPosEmpty) 
    5759      { 
    58         globaIndexWeightFromDestToSource[globalIndex].push_back(make_pair(globalIndexSrcGrid[idx][i], currentVecWeight[i])); 
     60        this->computeGlobalGridIndexFromGlobalIndexElement(itTransMap->first, 
     61                                                       itTransMap->second, 
     62                                                       itTransPos->second, 
     63                                                       elementPositionInGrid, 
     64                                                       gridDestGlobalDim, 
     65                                                       gridSrcGlobalDim, 
     66                                                       vecGlobalIndexGridSendToServer, 
     67                                                       globalIndexDestGrid, 
     68                                                       globalIndexSrcGrid); 
     69        ++itTransPos; 
     70      } 
     71      else 
     72      { 
     73        this->computeGlobalGridIndexFromGlobalIndexElement(itTransMap->first, 
     74                                                       itTransMap->second, 
     75                                                       emptyTransPos, 
     76                                                       elementPositionInGrid, 
     77                                                       gridDestGlobalDim, 
     78                                                       gridSrcGlobalDim, 
     79                                                       vecGlobalIndexGridSendToServer, 
     80                                                       globalIndexDestGrid, 
     81                                                       globalIndexSrcGrid); 
     82      } 
     83      size_t globalIndexSize = globalIndexDestGrid.numElements(); 
     84      const std::vector<double>& currentVecWeight = itTransWeight->second; 
     85      for (size_t idx = 0; idx < globalIndexSize; ++idx) 
     86      { 
     87        size_t globalIndex = globalIndexDestGrid(idx); 
     88        for (int i = 0; i < globalIndexSrcGrid[idx].size(); ++i) 
     89        { 
     90          globaIndexWeightFromDestToSource[globalIndex].push_back(make_pair(globalIndexSrcGrid[idx][i], currentVecWeight[i])); 
     91        } 
    5992      } 
    6093    } 
     
    6295} 
    6396 
     97void CGenericAlgorithmTransformation::computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs) 
     98{ 
     99  computeIndexSourceMapping_(dataAuxInputs); 
    64100} 
     101 
     102std::vector<StdString> CGenericAlgorithmTransformation::getIdAuxInputs() 
     103{ 
     104  return idAuxInputs_; 
     105} 
     106 
     107} 
  • XIOS/trunk/src/transformation/generic_algorithm_transformation.hpp

    r668 r827  
    3232                                std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexWeightFromDestToSource); 
    3333 
     34  std::vector<StdString> getIdAuxInputs(); 
     35 
    3436  /*! 
    3537  Compute global index mapping from one element of destination grid to the corresponding element of source grid 
    3638  */ 
    37   virtual void computeIndexSourceMapping() = 0; 
     39  void computeIndexSourceMapping(const std::vector<CArray<double,1>* >& dataAuxInputs = std::vector<CArray<double,1>* >()); 
    3840 
    3941protected: 
     
    5052  virtual void computeGlobalGridIndexFromGlobalIndexElement(int destGlobalIndex, 
    5153                                                        const std::vector<int>& srcGlobalIndex, 
     54                                                        const std::vector<int>& destGlobalIndexPositionInGrid, 
    5255                                                        int elementPositionInGrid, 
    5356                                                        const std::vector<int>& gridDestGlobalDim, 
     
    5760                                                        std::vector<std::vector<size_t> >& globalIndexSrcGrid) = 0; 
    5861 
    59  
     62  virtual void computeIndexSourceMapping_(const std::vector<CArray<double,1>* >&) = 0; 
    6063 
    6164protected: 
    62   std::map<int, std::vector<int> > transformationMapping_; 
    63   std::map<int, std::vector<double> > transformationWeight_; 
     65  //! Map between global index of destination element and source element 
     66  std::vector<std::map<int, std::vector<int> > > transformationMapping_; 
     67  //! Weight corresponding of source to destination 
     68  std::vector<std::map<int, std::vector<double> > > transformationWeight_; 
     69  //! Map of global index of destination element and corresponding global index of other elements in the same grid 
     70  //! By default, one index of an element corresponds to all index of remaining element in the grid. So it's empty 
     71  std::vector<std::map<int, std::vector<int> > > transformationPosition_; 
     72 
     73  //! Id of auxillary inputs which help doing transformation dynamically 
     74  std::vector<StdString> idAuxInputs_; 
    6475}; 
    6576 
  • XIOS/trunk/src/transformation/grid_transformation.cpp

    r824 r827  
    2323CGridTransformation::CGridTransformation(CGrid* destination, CGrid* source) 
    2424: gridSource_(source), gridDestination_(destination), originalGridSource_(source), 
    25   algoTypes_(), nbAlgos_(0), currentGridIndexToOriginalGridIndex_() 
     25  algoTypes_(), nbAlgos_(0), currentGridIndexToOriginalGridIndex_(), tempGrids_(), 
     26  auxInputs_(), dynamicalTransformation_(false) 
    2627 
    2728{ 
     
    6869    if (!isSpecialTransformation(transType)) ++nbAlgos_; 
    6970  } 
    70  
    71   if (1<nbAlgos_)  // Only when there are more than 1 algorithm, will we create temporary grid source 
    72   { 
    73     std::vector<CAxis*> axisSrcTmp = gridSource_->getAxis(), axisSrc; 
    74     std::vector<CDomain*> domainSrcTmp = gridSource_->getDomains(), domainSrc; 
    75     for (int idx = 0; idx < axisSrcTmp.size(); ++idx) 
    76     { 
    77       CAxis* axis = CAxis::createAxis(); 
    78       axis->axis_ref.setValue(axisSrcTmp[idx]->getId()); 
    79       axis->solveRefInheritance(true); 
    80       axis->solveInheritanceTransformation(); 
    81       axis->checkAttributesOnClient(); 
    82       axisSrc.push_back(axis); 
    83     } 
    84  
    85     for (int idx = 0; idx < domainSrcTmp.size(); ++idx) 
    86     { 
    87       CDomain* domain = CDomain::createDomain(); 
    88       domain->domain_ref.setValue(domainSrcTmp[idx]->getId()); 
    89       domain->solveRefInheritance(true); 
    90       domain->solveInheritanceTransformation(); 
    91       domain->checkAttributesOnClient(); 
    92       domainSrc.push_back(domain); 
    93     } 
    94  
    95     gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order); 
    96     gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order); 
    97   } 
    9871} 
    9972 
    10073CGridTransformation::~CGridTransformation() 
    10174{ 
    102   std::list<CGenericAlgorithmTransformation*>::const_iterator itb = algoTransformation_.begin(), it, 
    103                                                               ite = algoTransformation_.end(); 
     75  std::vector<CGenericAlgorithmTransformation*>::const_iterator itb = algoTransformation_.begin(), it, 
     76                                                                ite = algoTransformation_.end(); 
    10477  for (it = itb; it != ite; ++it) delete (*it); 
    10578} 
     
    181154        algoTypes_.push_back(false); 
    182155        ++transformationOrder; 
     156        std::vector<StdString> auxInput = (it->second)->checkAuxInputs(); 
     157        for (int idx = 0; idx < auxInput.size(); ++idx) auxInputs_.push_back(auxInput[idx]); 
    183158      } 
    184159    } 
     
    208183        algoTypes_.push_back(true); 
    209184        ++transformationOrder; 
     185        std::vector<StdString> auxInput = (it->second)->checkAuxInputs(); 
     186        for (int idx = 0; idx < auxInput.size(); ++idx) auxInputs_.push_back(auxInput[idx]); 
    210187      } 
    211188    } 
     
    313290  \param [in] transType transformation type 
    314291*/ 
    315 void CGridTransformation::setUpGrid(int elementPositionInGrid, ETranformationType transType) 
    316 { 
     292void CGridTransformation::setUpGrid(int elementPositionInGrid, ETranformationType transType, int nbTransformation) 
     293{ 
     294  if (!tempGrids_.empty() && (getNbAlgo()-1) == tempGrids_.size()) 
     295  { 
     296    gridSource_ = tempGrids_[nbTransformation]; 
     297    return; 
     298  } 
     299 
    317300  std::vector<CAxis*> axisListDestP = gridDestination_->getAxis(); 
    318   std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(); 
     301  std::vector<CAxis*> axisListSrcP = gridSource_->getAxis(), axisSrc; 
    319302 
    320303  std::vector<CDomain*> domListDestP = gridDestination_->getDomains(); 
    321   std::vector<CDomain*> domListSrcP = gridSource_->getDomains(); 
    322  
    323   int axisIndex, domainIndex; 
     304  std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainSrc; 
     305 
     306  int axisIndex = -1, domainIndex = -1; 
    324307  switch (transType) 
    325308  { 
     
    327310    case TRANS_ZOOM_DOMAIN: 
    328311      domainIndex = elementPosition2DomainPositionInGrid_[elementPositionInGrid]; 
    329       domListSrcP[domainIndex]->clearAllAttributes(); 
    330       domListSrcP[domainIndex]->duplicateAttributes(domListDestP[domainIndex]); 
    331312      break; 
    332313 
     
    335316    case TRANS_INVERSE_AXIS: 
    336317      axisIndex =  elementPosition2AxisPositionInGrid_[elementPositionInGrid]; 
    337       axisListSrcP[axisIndex]->clearAllAttributes(); 
    338       axisListSrcP[axisIndex]->duplicateAttributes(axisListDestP[axisIndex]); 
    339318      break; 
    340319    default: 
    341320      break; 
    342321  } 
    343   gridSource_->createMask(); 
    344   gridSource_->computeGridGlobalDimension(domListSrcP, axisListSrcP, gridSource_->axis_domain_order); 
     322 
     323  for (int idx = 0; idx < axisListSrcP.size(); ++idx) 
     324  { 
     325    CAxis* axis = CAxis::createAxis(); 
     326    if (axisIndex != idx) axis->axis_ref.setValue(axisListSrcP[idx]->getId()); 
     327    else axis->axis_ref.setValue(axisListDestP[idx]->getId()); 
     328    axis->solveRefInheritance(true); 
     329    axis->checkAttributesOnClient(); 
     330    axisSrc.push_back(axis); 
     331  } 
     332 
     333  for (int idx = 0; idx < domListSrcP.size(); ++idx) 
     334  { 
     335    CDomain* domain = CDomain::createDomain(); 
     336    if (domainIndex != idx) domain->domain_ref.setValue(domListSrcP[idx]->getId()); 
     337    else domain->domain_ref.setValue(domListDestP[idx]->getId()); 
     338    domain->solveRefInheritance(true); 
     339    domain->checkAttributesOnClient(); 
     340    domainSrc.push_back(domain); 
     341  } 
     342 
     343  gridSource_ = CGrid::createGrid(domainSrc, axisSrc, gridDestination_->axis_domain_order); 
     344  gridSource_->computeGridGlobalDimension(domainSrc, axisSrc, gridDestination_->axis_domain_order); 
     345 
     346  tempGrids_.push_back(gridSource_); 
    345347} 
    346348 
     
    353355  -) Make current grid destination become grid source in the next transformation 
    354356*/ 
    355 void CGridTransformation::computeAll() 
     357void CGridTransformation::computeAll(const std::vector<CArray<double,1>* >& dataAuxInputs) 
    356358{ 
    357359  if (nbAlgos_ < 1) return; 
     360  if (!auxInputs_.empty() && !dynamicalTransformation_) { dynamicalTransformation_ = true; return; } 
     361  if (dynamicalTransformation_) GlobalIndexMap().swap(currentGridIndexToOriginalGridIndex_);  // Reset map 
    358362 
    359363  CContext* context = CContext::getCurrent(); 
     
    373377 
    374378    // First of all, select an algorithm 
    375     selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]); 
    376     algo = algoTransformation_.back(); 
     379    if (!dynamicalTransformation_ || (algoTransformation_.size() < listAlgos_.size())) 
     380    { 
     381      selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]); 
     382      algo = algoTransformation_.back(); 
     383    } 
     384    else 
     385      algo = algoTransformation_[std::distance(itb, it)]; 
    377386 
    378387    if (0 != algo) // Only registered transformation can be executed 
    379388    { 
     389      algo->computeIndexSourceMapping(dataAuxInputs); 
     390 
    380391      // Recalculate the distribution of grid destination 
    381392      CDistributionClient distributionClientDest(client->clientRank, gridDestination_); 
     
    398409 
    399410        // Now grid destination becomes grid source in a new transformation 
    400         if (nbAgloTransformation != (nbAlgos_-1)) setUpGrid(elementPositionInGrid, transType); 
     411        if (nbAgloTransformation != (nbAlgos_-1)) setUpGrid(elementPositionInGrid, transType, nbAgloTransformation); 
    401412      } 
    402413      else 
     
    741752    int sourceRank = itMapRecv->first; 
    742753    int numGlobalIndex = (itMapRecv->second).size(); 
     754    localIndexToReceiveOnGridDest_[sourceRank].resize(numGlobalIndex); 
    743755    for (int i = 0; i < numGlobalIndex; ++i) 
    744756    { 
     
    754766        } 
    755767      } 
    756       localIndexToReceiveOnGridDest_[sourceRank].push_back(tmpVec); 
     768      localIndexToReceiveOnGridDest_[sourceRank][i] = tmpVec; 
    757769    } 
    758770  } 
  • XIOS/trunk/src/transformation/grid_transformation.hpp

    r824 r827  
    4141  ~CGridTransformation(); 
    4242 
    43   void computeAll(); 
     43  void computeAll(const std::vector<CArray<double,1>* >& dataAuxInput=std::vector<CArray<double,1>* >()); 
    4444 
    4545  const std::map<int, CArray<int,1> >& getLocalIndexToSendFromGridSource() const; 
     
    4949  ListAlgoType getAlgoList() const {return listAlgos_; } 
    5050  int getNbAlgo() { return nbAlgos_; } 
     51  const std::vector<StdString>& getAuxInputs() const { return auxInputs_; } 
    5152 
    5253protected: 
     
    6061  void selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder); 
    6162  void selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, bool isDomainAlgo); 
    62   void setUpGrid(int elementPositionInGrid, ETranformationType transType); 
     63  void setUpGrid(int elementPositionInGrid, ETranformationType transType, int nbTransformation); 
    6364  void computeFinalTransformationMapping(); 
    6465  void computeTransformationFromOriginalGridSource(const std::map<size_t, std::vector<std::pair<size_t,double> > >& globaIndexMapFromDestToSource); 
     
    8990 
    9091  // Mapping between position of an element in grid and its transformation (if any) 
    91   std::list<CGenericAlgorithmTransformation*> algoTransformation_; 
     92  std::vector<CGenericAlgorithmTransformation*> algoTransformation_; 
    9293 
    9394  //! Mapping of (grid) global index representing tranformation. 
     
    105106  //! (Grid) Global index of grid source 
    106107  GlobalIndexMap currentGridIndexToOriginalGridIndex_; 
     108 
     109  std::vector<CGrid*> tempGrids_; 
     110  std::vector<StdString> auxInputs_; 
     111  bool dynamicalTransformation_; 
    107112}; 
    108113 
Note: See TracChangeset for help on using the changeset viewer.