Changeset 978 for XIOS/trunk/src

Timestamp:

10/21/16 13:40:33 (8 years ago)

Author:

mhnguyen

Message:

Correcting various bugs relating to transformation

+) Fix the order of transformation selection
+) Correct domain transformation selection
+) Reorganize test_remap

Test
+) On Curie
+) All tests pass

Location:

XIOS/trunk/src

Files:

• test/test_remap.f90 (modified) (7 diffs)
• transformation/domain_algorithm_compute_connectivity.cpp (modified) (1 diff)
• transformation/domain_algorithm_expand.cpp (modified) (1 diff)
• transformation/domain_algorithm_interpolate.cpp (modified) (1 diff)
• transformation/grid_generate.cpp (modified) (1 diff)
• transformation/grid_transformation.cpp (modified) (9 diffs)
• transformation/grid_transformation.hpp (modified) (1 diff)
• transformation/grid_transformation_selector.cpp (modified) (9 diffs)
• transformation/grid_transformation_selector.hpp (modified) (3 diffs)

XIOS/trunk/src/test/test_remap.f90

@@ -20 +20 @@
   DOUBLE PRECISION,ALLOCATABLE :: src_boundslon(:,:), dst_boundslon(:,:)
   DOUBLE PRECISION,ALLOCATABLE :: src_boundslat(:,:), dst_boundslat(:,:)
-  DOUBLE PRECISION,ALLOCATABLE :: src_field(:)
-  DOUBLE PRECISION,ALLOCATABLE :: tmp_field(:), tmp_field_1(:), tmp_field_2(:)
-  DOUBLE PRECISION,ALLOCATABLE :: src_field_3D(:,:), lval(:), lval1(:), src_field_pression(:,:)
+  DOUBLE PRECISION,ALLOCATABLE :: src_field_2D(:)
+  DOUBLE PRECISION,ALLOCATABLE :: tmp_field_0(:), tmp_field_1(:), tmp_field_2(:)
+  DOUBLE PRECISION,ALLOCATABLE :: src_field_3D(:,:), src_field_4D(:,:,:), lval(:), lval1(:), lval2(:), src_field_pression(:,:)
   LOGICAL,ALLOCATABLE :: src_mask_2D(:)
   INTEGER :: src_ni_glo, dst_ni_glo;
@@ -37 +37 @@
   INTEGER :: ts
   INTEGER :: i
-  INTEGER,PARAMETER :: llm=5, interpolatedLlm = 4
+  INTEGER,PARAMETER :: llm=5, interpolatedLlm = 4, llm2 = 6
 
   CALL MPI_INIT(ierr)
@@ -68 +68 @@
   ALLOCATE(src_boundslon(src_nvertex,src_ni))
   ALLOCATE(src_boundslat(src_nvertex,src_ni))
-  ALLOCATE(src_field(src_ni))
+  ALLOCATE(src_field_2D(src_ni))
   ALLOCATE(src_field_3D(src_ni,llm))
+  ALLOCATE(src_field_4D(src_ni,llm,llm2))
   ALLOCATE(src_mask_2D(src_ni))
   ALLOCATE(src_field_pression(src_ni,llm))
   ALLOCATE(lval(llm))
   ALLOCATE(lval1(interpolatedLlm))
+  ALLOCATE(lval2(llm2))
 
   ierr=NF90_INQ_VARID(ncid,"lon",varid)
@@ -84 +86 @@
   ierr=NF90_GET_VAR(ncid,varid, src_boundslat, start=(/1,src_ibegin+1/),count=(/src_nvertex,src_ni/))
   ierr=NF90_INQ_VARID(ncid,"val",varid)
-  ierr=NF90_GET_VAR(ncid,varid, src_field, start=(/src_ibegin+1/),count=(/src_ni/))
+  ierr=NF90_GET_VAR(ncid,varid, src_field_2D, start=(/src_ibegin+1/),count=(/src_ni/))
   DO i=1,src_ni
-    src_field_3D(i,:) = src_field(i)
+    src_field_3D(i,:) = src_field_2D(i)
     IF (MOD(i,10)==0) THEN
       src_mask_2D(i)=.FALSE.
@@ -97 +99 @@
     lval(i) = i*100
     src_field_pression(:,i) = i * 100
-    src_field_3D(:,i) = src_field_3D(:,i) + i * 10
+    src_field_3D(:,i) = src_field_3D(:,i) + i * 10    
+  ENDDO
+
+  DO i=1,llm2
+    src_field_4D(:,:,i) = src_field_3D(:,:)
   ENDDO
 
   DO i=1,interpolatedLlm
+    lval2(i) = i*10
     lval1(i) = i*10 + 2
   ENDDO
@@ -149 +156 @@
 
   CALL xios_set_axis_attr("src_axis", n_glo=llm, value=lval)
+  CALL xios_set_axis_attr("nq", n_glo=llm2, value=lval2)
 
   CALL xios_set_domain_attr("dst_domain", ni_glo=dst_ni_glo, ibegin=dst_ibegin, ni=dst_ni, type="unstructured")
@@ -158 +166 @@
 
   CALL xios_close_context_definition()
-!  CALL xios_get_domain_attr("src_domain_regular_tmp", ni=src_tmp_ni, nj=src_tmp_nj)
-!  ALLOCATE(tmp_field(src_tmp_ni*src_tmp_nj))
-!
-!  CALL xios_get_axis_attr("src_axis_curvilinear", n=src_tmp_n)
-!  CALL xios_get_domain_attr("src_domain_curvilinear", ni=src_tmp_ni, nj=src_tmp_nj)
-!  ALLOCATE(tmp_field_1(src_tmp_ni*src_tmp_nj*src_tmp_n))
-
-!  CALL xios_get_domain_attr("src_domain_unstructured", ni=src_tmp_ni, nj=src_tmp_nj)
-!  ALLOCATE(tmp_field_2(src_tmp_ni*src_tmp_nj))
+  ! CALL xios_get_domain_attr("src_domain_regular_read", ni=src_tmp_ni, nj=src_tmp_nj)
+  ! ALLOCATE(tmp_field_0(src_tmp_ni*src_tmp_nj))
+
+  ! CALL xios_get_axis_attr("src_axis_curvilinear_read", n=src_tmp_n)
+  ! CALL xios_get_domain_attr("src_domain_curvilinear_read", ni=src_tmp_ni, nj=src_tmp_nj)
+  ! ALLOCATE(tmp_field_1(src_tmp_ni*src_tmp_nj*src_tmp_n))
+
+  ! CALL xios_get_domain_attr("src_domain_unstructured_read", ni=src_tmp_ni, nj=src_tmp_nj)
+  ! ALLOCATE(tmp_field_2(src_tmp_ni*src_tmp_nj))
 
   DO ts=1,5
-
-!    CALL xios_recv_field("src_field_regular_tmp", tmp_field)
-!    CALL xios_recv_field("src_field_curvilinear", tmp_field_1)
-!    CALL xios_recv_field("field_src_unstructred", tmp_field_2)
+    ! CALL xios_recv_field("src_field_regular", tmp_field_0)
+    ! CALL xios_recv_field("src_field_curvilinear", tmp_field_1)
+    ! CALL xios_recv_field("src_field_unstructured", tmp_field_2)
     CALL xios_update_calendar(ts)
-    CALL xios_send_field("src_field",src_field)
+    CALL xios_send_field("src_field_2D",src_field_2D)
+    CALL xios_send_field("src_field_2D_clone",src_field_2D)
     CALL xios_send_field("src_field_3D",src_field_3D)
+    CALL xios_send_field("src_field_4D",src_field_4D)
     CALL xios_send_field("src_field_3D_pression",src_field_pression)
-!    CALL xios_send_field("tmp_field",tmp_field)
-!    CALL xios_send_field("tmp_field_1",tmp_field_1)
-!    CALL xios_send_field("tmp_field_2",tmp_field_2)
+    ! CALL xios_send_field("tmp_field_0",tmp_field_0)
+    ! CALL xios_send_field("tmp_field_1",tmp_field_1)
+    ! CALL xios_send_field("tmp_field_2",tmp_field_2)
     CALL wait_us(5000) ;
-  ENDDO
+   ENDDO
 
   CALL xios_context_finalize()
 
-  DEALLOCATE(src_lon, src_lat, src_boundslon,src_boundslat, src_field)
+  DEALLOCATE(src_lon, src_lat, src_boundslon,src_boundslat, src_field_2D)
   DEALLOCATE(dst_lon, dst_lat, dst_boundslon,dst_boundslat)
-!  DEALLOCATE(tmp_field, tmp_field_1, tmp_field_2)
+  ! DEALLOCATE(tmp_field_0, tmp_field_1, tmp_field_2)
 
   CALL MPI_COMM_FREE(comm, ierr)

XIOS/trunk/src/transformation/domain_algorithm_compute_connectivity.cpp

@@ -29 +29 @@
 
   CComputeConnectivityDomain* compute_connectivityDomain = dynamic_cast<CComputeConnectivityDomain*> (transformation);
-  int domainDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
-  int domainSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
 
   return (new CDomainAlgorithmComputeConnectivity(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], compute_connectivityDomain));

XIOS/trunk/src/transformation/domain_algorithm_expand.cpp

@@ -31 +31 @@
 
   CExpandDomain* expandDomain = dynamic_cast<CExpandDomain*> (transformation);
-  int domainDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
-  int domainSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
 
   return (new CDomainAlgorithmExpand(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], expandDomain));

XIOS/trunk/src/transformation/domain_algorithm_interpolate.cpp

@@ -36 +36 @@
 
   CInterpolateDomain* interpolateDomain = dynamic_cast<CInterpolateDomain*> (transformation);
-  int domainDstIndex = elementPositionInGridDst2AxisPosition[elementPositionInGrid];
-  int domainSrcIndex = elementPositionInGridSrc2AxisPosition[elementPositionInGrid];
+  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
+  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
 
   return (new CDomainAlgorithmInterpolate(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], interpolateDomain));

XIOS/trunk/src/transformation/grid_generate.cpp

@@ -103 +103 @@
     int elementPositionInGrid = it->first;
     ETranformationType transType = (it->second).first;
-    int transformationOrder = (it->second).second;
+    int transformationOrder = (it->second).second.first;
+    int algoType = (it->second).second.second;
 
     // First of all, select an algorithm
-    selectAlgo(elementPositionInGrid, transType, transformationOrder, algoTypes_[std::distance(itb, it)]);
+    selectAlgo(elementPositionInGrid, transType, transformationOrder, algoType);
   }
 }

XIOS/trunk/src/transformation/grid_transformation.cpp

@@ -100 +100 @@
   CDomain::TransMapTypes trans = domainListDestP[domainIndex]->getAllTransformations();
   CDomain::TransMapTypes::const_iterator it = trans.begin();
-  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation
+  for (int i = 0; i < transformationOrder; ++i, ++it) {}  // Find the correct transformation  
 
   CGenericAlgorithmTransformation* algo = 0;
@@ -118 +118 @@
 
 /*!
+  Find position of element in a grid as well as its type (domain, axis, scalar) and position in its own element list
+  \return element position: map<int,<int,int> > corresponds to <element position in grid, <element type, element position in element list> >
+*/
+std::map<int,std::pair<int,int> > CGridTransformation::getElementPosition(CGrid* grid)
+{
+  std::vector<CScalar*> scalarListP = grid->getScalars(); 
+  std::vector<CAxis*> axisListP = grid->getAxis();
+  std::vector<CDomain*> domListP = grid->getDomains();  
+  CArray<int,1> axisDomainOrder = grid->axis_domain_order;
+  int scalarIndex = 0, axisIndex = 0, domainIndex = 0;
+  int nbElement = axisDomainOrder.numElements(), elementDim;
+  std::map<int,std::pair<int,int> > elementPosition;
+  for (int idx = 0; idx < nbElement; ++idx)
+  {
+    elementDim = axisDomainOrder(idx);
+    switch (elementDim)
+    {
+      case 2:
+        elementPosition[idx] = std::make_pair(elementDim, domainIndex);
+        ++domainIndex;
+        break;
+      case 1:
+        elementPosition[idx] = std::make_pair(elementDim, axisIndex);
+        ++axisIndex;
+        break;
+      case 0:
+        elementPosition[idx] = std::make_pair(elementDim, scalarIndex);
+        ++scalarIndex;
+        break;
+      default:
+        break;        
+    }
+  }
+
+  return elementPosition;  
+}
+
+/*!
   If there are more than one transformation, a new temporary grid will be created and it will play the role of grid destination.
 This new created one keeps a pointer to the real transformed element of grid destination and generate new copies of other elements from grid source.
@@ -123 +161 @@
   \param [in] transType transformation type
 */
-void CGridTransformation::setUpGridDestination(int elementPositionInGrid, ETranformationType transType, AlgoType algoType)
+void CGridTransformation::setUpGridDestination(int elementPositionInGrid, ETranformationType transType)
 {
   if (isSpecialTransformation(transType)) return;
@@ -130 +168 @@
   {
     tempGridDests_.resize(0);
+  }
+
+  if (1 == getNbAlgo()) 
+  {
+    tmpGridDestination_ = gridDestination_;
+    return;
   }
 
@@ -144 +188 @@
   CArray<int,1> axisDomainOrderDst = gridDestination_->axis_domain_order;
 
-  int scalarIndex = -1, axisIndex = -1, domainIndex = -1;
-  switch (algoType)
-  {
-    case domainType:
-      domainIndex = elementPositionInGridDst2DomainPosition_[elementPositionInGrid];
-      break;
-    case axisType:
-      axisIndex =  elementPositionInGridDst2AxisPosition_[elementPositionInGrid];
-      break;
-    case scalarType:
-      scalarIndex = elementPositionInGridDst2ScalarPosition_[elementPositionInGrid];
-      break;
-    default:
-      break;
-  }
-
+  std::map<int,std::pair<int,int> > elementPositionSrc = getElementPosition(gridSource_);
+  std::map<int,std::pair<int,int> > elementPositionDst = getElementPosition(gridDestination_);
+
+  CArray<int,1> elementOrder(axisDomainOrderDst.numElements());
   for (int idx = 0; idx < axisDomainOrderDst.numElements(); ++idx)
   {
-    int dimElementDst = axisDomainOrderDst(idx);
-    if (2 == dimElementDst)
-    {
-      if (elementPositionInGrid == idx)
-        domainDst.push_back(domListDestP[domainIndex]);
-      else
-        domainDst.push_back(domListSrcP[elementPositionInGridSrc2DomainPosition_[elementPositionInGrid]]);
-    }
-    else if (1 == dimElementDst)
-    {
-      if (elementPositionInGrid == idx)
-        axisDst.push_back(axisListDestP[axisIndex]);
-      else
-        axisDst.push_back(axisListSrcP[elementPositionInGridSrc2AxisPosition_[elementPositionInGrid]]);
+    if (elementPositionInGrid == idx)
+    {
+      int dimElementDst = elementPositionDst[idx].first;
+      int elementIndex  = elementPositionDst[idx].second;
+      switch (dimElementDst) 
+      {
+        case 2:
+          domainDst.push_back(domListDestP[elementIndex]);
+          break;
+        case 1:
+          axisDst.push_back(axisListDestP[elementIndex]);
+          break;
+        case 0:
+          scalarDst.push_back(scalarListDestP[elementIndex]);
+          break;
+        default:
+          break;
+      }
+      elementOrder(idx) = dimElementDst;
     }
     else
-    {
-      if (elementPositionInGrid == idx)
-        scalarDst.push_back(scalarListDestP[scalarIndex]);
-      else
-        scalarDst.push_back(scalarListSrcP[elementPositionInGridSrc2ScalarPosition_[elementPositionInGrid]]);
-    }
-  }
-
-  tmpGridDestination_ = CGrid::createGrid(domainDst, axisDst, scalarDst, gridDestination_->axis_domain_order);
-  tmpGridDestination_->computeGridGlobalDimension(domainDst, axisDst, scalarDst, gridDestination_->axis_domain_order);
+    {      
+      int dimElementSrc = elementPositionSrc[idx].first;
+      int elementIndex  = elementPositionSrc[idx].second;
+      switch (dimElementSrc)
+      {
+        case 2:
+          domainDst.push_back(domListSrcP[elementIndex]);
+          break;
+        case 1:
+          axisDst.push_back(axisListSrcP[elementIndex]);
+          break;
+        case 0:
+          scalarDst.push_back(scalarListSrcP[elementIndex]);
+          break;
+        default:
+          break;
+      }
+      elementOrder(idx) = dimElementSrc;
+    }
+  }
+
+  tmpGridDestination_ = CGrid::createGrid(domainDst, axisDst, scalarDst, elementOrder);
+  tmpGridDestination_->computeGridGlobalDimension(domainDst, axisDst, scalarDst, elementOrder);
   tempGridDests_.push_back(tmpGridDestination_);
 }
@@ -198 +248 @@
   \param [in] transType transformation type
 */
-void CGridTransformation::setUpGridSource(int elementPositionInGrid, AlgoType algoType)
+void CGridTransformation::setUpGridSource(int elementPositionInGrid)
 {
   if (!tempGridSrcs_.empty() && (getNbAlgo()-1) == tempGridSrcs_.size())
@@ -214 +264 @@
   std::vector<CDomain*> domListSrcP = gridSource_->getDomains(), domainSrc;
 
-  CArray<int,1> axisDomainOrderDst = gridDestination_->axis_domain_order;
-
-  int axisIndex = -1, domainIndex = -1, scalarIndex = -1;
-  int axisListIndex = 0, domainListIndex = 0, scalarListIndex = 0;
-  switch (algoType)
-  {
-    case domainType:
-      domainIndex = elementPositionInGridDst2DomainPosition_[elementPositionInGrid];
-      break;
-    case axisType:
-      axisIndex =  elementPositionInGridDst2AxisPosition_[elementPositionInGrid];
-      break;
-    case scalarType:
-      scalarIndex = elementPositionInGridDst2ScalarPosition_[elementPositionInGrid];
-      break;
-    default:
-      break;
-  }
+  CArray<int,1> axisDomainOrderSrc = gridSource_->axis_domain_order;
+  CArray<int,1> axisDomainOrderDst = tmpGridDestination_->axis_domain_order;
+
+  std::map<int,std::pair<int,int> > elementPositionSrc = getElementPosition(gridSource_);
+  std::map<int,std::pair<int,int> > elementPositionDst = getElementPosition(tmpGridDestination_);
 
   for (int idx = 0; idx < axisDomainOrderDst.numElements(); ++idx)
-  {
-    int dimElementDst = axisDomainOrderDst(idx);
-    if (2 == dimElementDst)
-    {
-      if (elementPositionInGrid == idx)
-        domainSrc.push_back(domListDestP[domainIndex]);
-      else
+  {   
+    if (elementPositionInGrid == idx)
+    {
+      int dimElementDst = elementPositionDst[idx].first;
+      int elementIndex  = elementPositionDst[idx].second;
+      if (2 == dimElementDst)
       {
         CDomain* domain = CDomain::createDomain();
-        domain->domain_ref.setValue(domListDestP[domainListIndex]->getId());
+        domain->domain_ref.setValue(domListDestP[elementIndex]->getId());
         domain->solveRefInheritance(true);
         domain->checkAttributesOnClient();
         domainSrc.push_back(domain);
       }
-      ++domainListIndex;
-    }
-    else if (1 == dimElementDst)
-    {
-      if (elementPositionInGrid == idx)
-        axisSrc.push_back(axisListDestP[axisIndex]);
-      else
+      else if (1 == dimElementDst)
       {
         CAxis* axis = CAxis::createAxis();
-        axis->axis_ref.setValue(axisListDestP[axisListIndex]->getId());
+        axis->axis_ref.setValue(axisListDestP[elementIndex]->getId());
         axis->solveRefInheritance(true);
         axis->checkAttributesOnClient();
-        axisSrc.push_back(axis);
-      }
-      ++axisListIndex;
-    }
-    else
-    {
-      if (elementPositionInGrid == idx)
-        scalarSrc.push_back(scalarListDestP[scalarIndex]);
+        axisSrc.push_back(axis); 
+      }
       else
       {
         CScalar* scalar = CScalar::createScalar();
-        scalar->scalar_ref.setValue(scalarListDestP[scalarListIndex]->getId());
+        scalar->scalar_ref.setValue(scalarListDestP[elementIndex]->getId());
         scalar->solveRefInheritance(true);
         scalar->checkAttributesOnClient();
         scalarSrc.push_back(scalar);
       }
-      ++scalarListIndex;
+    }
+    else
+    {      
+      int dimElementDst = elementPositionDst[idx].first;
+      int elementIndex  = elementPositionDst[idx].second;
+      switch (dimElementDst)
+      {
+        case 2:
+          domainSrc.push_back(domListDestP[elementIndex]);
+          break;
+        case 1:
+          axisSrc.push_back(axisListDestP[elementIndex]);
+          break;
+        case 0:
+          scalarSrc.push_back(scalarListDestP[elementIndex]);
+          break;
+        default:
+          break;
+      }
     }
   }
@@ -323 +365 @@
     int elementPositionInGrid = it->first;
     ETranformationType transType = (it->second).first;
-    int transformationOrder = (it->second).second;
+    int transformationOrder = (it->second).second.first;
+    int algoType = ((it->second).second.second); //algoTypes_[std::distance(itb, it)];
     SourceDestinationIndexMap globaIndexWeightFromSrcToDst;
-    AlgoType algoType = algoTypes_[std::distance(itb, it)];
+    
 
     // Create a temporary grid destination which contains transformed element of grid destination and
     // non-transformed elements to grid source
-    setUpGridDestination(elementPositionInGrid, transType, algoType);
+    setUpGridDestination(elementPositionInGrid, transType);
 
     // First of all, select an algorithm
@@ -359 +402 @@
       {
         // Now grid destination becomes grid source in a new transformation
-        if (nbAgloTransformation != (nbNormalAlgos_-1)) setUpGridSource(elementPositionInGrid, algoType);
+        if (nbAgloTransformation != (nbNormalAlgos_-1)) setUpGridSource(elementPositionInGrid);
       }
       ++nbAgloTransformation;

XIOS/trunk/src/transformation/grid_transformation.hpp

@@ -61 +61 @@
   virtual void selectDomainAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder);
 
-  void setUpGridSource(int elementPositionInGrid, AlgoType);
-  void setUpGridDestination(int elementPositionInGrid, ETranformationType, AlgoType);
+  void setUpGridSource(int elementPositionInGrid);
+  void setUpGridDestination(int elementPositionInGrid, ETranformationType);
   void computeTransformationMapping(const SourceDestinationIndexMap& globalIndexWeightFromSrcToDest);
+  std::map<int,std::pair<int,int> > getElementPosition(CGrid* grid);
 
 protected:

XIOS/trunk/src/transformation/grid_transformation_selector.cpp

@@ -102 +102 @@
 
 /*!
-  Initialize the algorithms (transformations)
-*/
-void CGridTransformationSelector::initializeAlgorithms()
+  Update position of elements in grid source and grid destination as well as their positions in element list
+*/
+void CGridTransformationSelector::updateElementPosition()
 {
   int idxScalar = 0, idxAxis = 0, idxDomain = 0;
   CArray<int,1> axisDomainOrderDst = gridDestination_->axis_domain_order;
+  std::map<int, int>().swap(elementPositionInGridDst2DomainPosition_);
+  std::map<int, int>().swap(elementPositionInGridDst2AxisPosition_);
+  std::map<int, int>().swap(elementPositionInGridDst2ScalarPosition_);
   for (int i = 0; i < axisDomainOrderDst.numElements(); ++i)
   {
@@ -130 +133 @@
   idxScalar = idxAxis = idxDomain = 0;
   CArray<int,1> axisDomainOrderSrc = gridSource_->axis_domain_order;
+  std::map<int, int>().swap(elementPositionInGridSrc2DomainPosition_);
+  std::map<int, int>().swap(elementPositionInGridSrc2AxisPosition_);
+  std::map<int, int>().swap(elementPositionInGridSrc2ScalarPosition_);
   for (int i = 0; i < axisDomainOrderSrc.numElements(); ++i)
   {
@@ -149 +155 @@
     }
   }
-
+}
+
+/*!
+  Initialize the algorithms (transformations)
+*/
+void CGridTransformationSelector::initializeAlgorithms()
+{
+  updateElementPosition();
+  CArray<int,1> axisDomainOrderDst = gridDestination_->axis_domain_order;
   for (int i = 0; i < axisDomainOrderDst.numElements(); ++i)
   {
@@ -183 +197 @@
     int scalarDstPos = -1, scalarSrcPos = -1;
     if (0 < elementPositionInGridDst2ScalarPosition_.count(scalarPositionInGrid))
-      scalarDstPos = elementPositionInGridDst2AxisPosition_[scalarPositionInGrid];
+      scalarDstPos = elementPositionInGridDst2ScalarPosition_[scalarPositionInGrid];
     if (0 < elementPositionInGridSrc2ScalarPosition_.count(scalarPositionInGrid))
-      scalarSrcPos = elementPositionInGridSrc2AxisPosition_[scalarPositionInGrid];
+      scalarSrcPos = elementPositionInGridSrc2ScalarPosition_[scalarPositionInGrid];
 
     // If source and destination grid share the same scalar
@@ -199 +213 @@
       for (it = itb; it != ite; ++it)
       {
-        listAlgos_.push_back(std::make_pair(scalarPositionInGrid, std::make_pair(it->first, transformationOrder)));
-        algoTypes_.push_back(scalarType);
+        listAlgos_.push_back(std::make_pair(scalarPositionInGrid, std::make_pair(it->first, std::make_pair(transformationOrder,0))));        
         ++transformationOrder;
         std::vector<StdString> auxInput = (it->second)->checkAuxInputs();
@@ -240 +253 @@
       for (it = itb; it != ite; ++it)
       {
-        listAlgos_.push_back(std::make_pair(axisPositionInGrid, std::make_pair(it->first, transformationOrder)));
-        algoTypes_.push_back(axisType);
+        listAlgos_.push_back(std::make_pair(axisPositionInGrid, std::make_pair(it->first, std::make_pair(transformationOrder,1))));        
         ++transformationOrder;
         std::vector<StdString> auxInput = (it->second)->checkAuxInputs();
@@ -264 +276 @@
     int domDstPos = -1, domSrcPos = -1;
     if (0 < elementPositionInGridDst2DomainPosition_.count(domPositionInGrid))
-      domDstPos = elementPositionInGridDst2AxisPosition_[domPositionInGrid];
+      domDstPos = elementPositionInGridDst2DomainPosition_[domPositionInGrid];
     if (0 < elementPositionInGridSrc2DomainPosition_.count(domPositionInGrid))
-      domSrcPos = elementPositionInGridSrc2AxisPosition_[domPositionInGrid];
+      domSrcPos = elementPositionInGridSrc2DomainPosition_[domPositionInGrid];
 
     // If source and destination grid share the same domain
@@ -280 +292 @@
       for (it = itb; it != ite; ++it)
       {
-        listAlgos_.push_back(std::make_pair(domPositionInGrid, std::make_pair(it->first, transformationOrder)));
-        algoTypes_.push_back(domainType);
+        listAlgos_.push_back(std::make_pair(domPositionInGrid, std::make_pair(it->first, std::make_pair(transformationOrder,2))));        
         ++transformationOrder;
         std::vector<StdString> auxInput = (it->second)->checkAuxInputs();
@@ -297 +308 @@
   \param [in] transType transformation type, for now we have Zoom_axis, inverse_axis
   \param [in] transformationOrder position of the transformation in an element (an element can have several transformation)
-  \param [in] isDomainAlgo flag to specify type of algorithm (for domain or axis)
-*/
-void CGridTransformationSelector::selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, AlgoType algoType)
-{
+  \param [in] algoType flag to specify type of algorithm (2 for domain, 1 for axis and 0 for scalar) 
+*/
+void CGridTransformationSelector::selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, int algoType)
+{
+  updateElementPosition();
   switch (algoType)
   {
-  case scalarType:
+  case 0:
     selectScalarAlgo(elementPositionInGrid, transType, transformationOrder);
     break;
-  case axisType:
+  case 1:
     selectAxisAlgo(elementPositionInGrid, transType, transformationOrder);
     break;
-  case domainType:
+  case 2:
     selectDomainAlgo(elementPositionInGrid, transType, transformationOrder);
     break;

XIOS/trunk/src/transformation/grid_transformation_selector.hpp

@@ -28 +28 @@
 {
 public:
-  typedef std::list<std::pair<int,std::pair<ETranformationType,int> > > ListAlgoType;
+  // Stupid C++98, should be replaced with tuple.
+  // List of algorithm to process. Order of element in this list
+  // <Position of element in grid, <transformationType, <transformation order, alogrithm type> > >
+  typedef std::list<std::pair<int,std::pair<ETranformationType,std::pair<int,int> > > > ListAlgoType;
 protected:
   enum AlgoType {
@@ -50 +53 @@
 
 protected:
+  void updateElementPosition();
   void initializeAlgorithms();
   void initializeDomainAlgorithms(int domPositionInGrid);
@@ -55 +59 @@
   void initializeScalarAlgorithms(int scalarPositionInGrid);
   void initializeTransformations(TransformationType type);
-  void selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, AlgoType algo);
+  void selectAlgo(int elementPositionInGrid, ETranformationType transType, int transformationOrder, int algo);
   bool isSpecialTransformation(ETranformationType transType);
   void registerTransformations();