source: XIOS/trunk/src/transformation/domain_algorithm_transformation.cpp @ 827

Last change on this file since 827 was 827, checked in by mhnguyen, 8 years ago

Implementing dynamic interpolation on axis

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

Test
+) On Curie
+) All test pass

File size: 7.0 KB
Line 
1/*!
2   \file domain_algorithm_transformation.hpp
3   \author Ha NGUYEN
4   \since 02 Jul 2015
5   \date 02 Jul 2015
6
7   \brief Interface for all domain transformation algorithms.
8 */
9
10#include "domain_algorithm_transformation.hpp"
11
12namespace xios {
13
14CDomainAlgorithmTransformation::CDomainAlgorithmTransformation(CDomain* domainDestination, CDomain* domainSource)
15 : CGenericAlgorithmTransformation(), domainDest_(domainDestination), domainSrc_(domainSource)
16{
17}
18
19CDomainAlgorithmTransformation::~CDomainAlgorithmTransformation()
20{
21}
22
23void CDomainAlgorithmTransformation::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
24{
25}
26
27/*!
28  Compute an array of global index from a global index on a domain
29  \param[in] domainDestGlobalIndex global index on an domain of destination grid
30  \param[in] domainSrcGlobalIndex global index on an domain of source grid (which are needed by one index on domain destination)
31  \param[in] domainPositionInGrid position of the domain in the grid
32  \param[in] gridDestGlobalDim dimension size of destination grid (it should share the same size for all dimension, maybe except the domain on which transformation is performed)
33  \param[in] gridSrcGlobalDim dimension size of source grid (it should share the same size for all dimension, maybe except the domain on which transformation is performed)
34  \param[in] globalIndexGridDestSendToServer global index of destination grid which are to be sent to server(s), this array is already acsending sorted
35  \param[in/out] globalIndexDestGrid array of global index (for 2d grid, this array is a line, for 3d, this array represents a plan). It should be preallocated
36  \param[in/out] globalIndexSrcGrid array of global index of source grid (for 2d grid, this array is a line, for 3d, this array represents a plan). It should be preallocated
37*/
38void CDomainAlgorithmTransformation::computeGlobalGridIndexFromGlobalIndexElement(int domainDestGlobalIndex,
39                                                                          const std::vector<int>& domainSrcGlobalIndex,
40                                                                          const std::vector<int>& destGlobalIndexPositionInGrid,
41                                                                          int domainPositionInGrid,
42                                                                          const std::vector<int>& gridDestGlobalDim,
43                                                                          const std::vector<int>& gridSrcGlobalDim,
44                                                                          const std::vector<size_t>& globalIndexGridDestSendToServer,
45                                                                          CArray<size_t,1>& globalIndexDestGrid,
46                                                                          std::vector<std::vector<size_t> >& globalIndexSrcGrid)
47{
48  int globalDim = gridDestGlobalDim.size();
49  int numElement = globalDim - 1;
50  int domainElementPosition = 0;
51  std::vector<int> currentIndex(globalDim);
52  std::vector<int> gridDomainGlobalDim(numElement), indexMap(numElement);
53  std::vector<int> idxLoop(numElement, 0), domainSrcGlobalDim(2), domainDestGlobalDim(2);
54
55
56  size_t ssize = 1, idx = 0, realGlobalIndexSize = 0;
57  for (int i = 0; i< globalDim; ++i)
58  {
59    if (domainPositionInGrid == i) continue;
60    if (domainPositionInGrid == (i+1))
61    {
62      domainDestGlobalDim[0] = gridDestGlobalDim[i];
63      domainDestGlobalDim[1] = gridDestGlobalDim[i+1];
64      domainSrcGlobalDim[0] = gridSrcGlobalDim[i];
65      domainSrcGlobalDim[1] = gridSrcGlobalDim[i+1];
66      gridDomainGlobalDim[idx] = 1;
67      domainElementPosition = idx;
68    }
69    else
70    {
71      gridDomainGlobalDim[idx] = gridDestGlobalDim[i];
72    }
73    indexMap[idx] = i;
74    ++idx;
75  }
76  int iIndex = domainDestGlobalIndex % domainDestGlobalDim[0];
77  int jIndex = domainDestGlobalIndex / domainDestGlobalDim[0];
78
79  for (int i = 0; i< numElement; ++i) ssize *= gridDomainGlobalDim[i];
80
81  std::vector<size_t>::const_iterator itbArr = globalIndexGridDestSendToServer.begin(), itArr,
82                                      iteArr = globalIndexGridDestSendToServer.end();
83  idx = 0;
84  while (idx < ssize)
85  {
86    for (int i = 0; i < numElement; ++i)
87    {
88      if (idxLoop[i] == gridDomainGlobalDim[i])
89      {
90        idxLoop[i] = 0;
91        ++idxLoop[i+1];
92      }
93    }
94
95    for (int i = 0; i < numElement; ++i)
96    {
97      if (domainElementPosition == i)
98      {
99        currentIndex[indexMap[i]]   = iIndex;
100        currentIndex[indexMap[i]+1] = jIndex;
101      }
102      else
103        currentIndex[indexMap[i]] = idxLoop[i];
104    }
105
106    size_t globIndex = currentIndex[0];
107    size_t mulDim = 1;
108    for (int k = 1; k < globalDim; ++k)
109    {
110      mulDim *= gridDestGlobalDim[k-1];
111      globIndex += (currentIndex[k])*mulDim;
112    }
113
114    if (std::binary_search(itbArr, iteArr, globIndex)) ++realGlobalIndexSize;
115    ++idxLoop[0];
116    ++idx;
117  }
118
119  if (globalIndexDestGrid.numElements() != realGlobalIndexSize)
120    globalIndexDestGrid.resize(realGlobalIndexSize);
121
122  if (realGlobalIndexSize != globalIndexSrcGrid.size()) globalIndexSrcGrid.resize(realGlobalIndexSize);
123  for (int i = 0; i < globalIndexSrcGrid.size(); ++i)
124    if (globalIndexSrcGrid[i].size() != domainSrcGlobalIndex.size())
125      globalIndexSrcGrid[i].resize(domainSrcGlobalIndex.size());
126
127  size_t realGlobalIndex = 0;
128  idx = 0;
129  idxLoop.assign(globalDim, 0);
130  while (idx < ssize)
131  {
132    for (int i = 0; i < numElement; ++i)
133    {
134      if (idxLoop[i] == gridDomainGlobalDim[i])
135      {
136        idxLoop[i] = 0;
137        ++idxLoop[i+1];
138      }
139    }
140
141    for (int i = 0; i < numElement; ++i)
142    {
143      if (domainElementPosition == i)
144      {
145        currentIndex[indexMap[i]]   = iIndex;
146        currentIndex[indexMap[i]+1] = jIndex;
147      }
148      else
149        currentIndex[indexMap[i]] = idxLoop[i];
150    }
151
152    size_t globIndex = currentIndex[0];
153    size_t mulDim = 1;
154    for (int k = 1; k < globalDim; ++k)
155    {
156      mulDim *= gridDestGlobalDim[k-1];
157      globIndex += (currentIndex[k])*mulDim;
158    }
159
160    if (std::binary_search(itbArr, iteArr, globIndex))
161    {
162      globalIndexDestGrid(realGlobalIndex) = globIndex;
163      for (int i = 0; i < globalIndexSrcGrid[realGlobalIndex].size(); ++i)
164      {
165        domainGlobalIndex(domainSrcGlobalIndex[i], domainSrcGlobalDim[0], domainSrcGlobalDim[1],
166                          currentIndex[indexMap[domainElementPosition]],
167                          currentIndex[indexMap[domainElementPosition]+1]);
168
169        globIndex = currentIndex[0];
170        mulDim = 1;
171        for (int k = 1; k < globalDim; ++k)
172        {
173          mulDim *= gridSrcGlobalDim[k-1];
174          globIndex += (currentIndex[k])*mulDim;
175        }
176        (globalIndexSrcGrid[realGlobalIndex])[i] = globIndex;
177      }
178      ++realGlobalIndex;
179    }
180
181    ++idxLoop[0];
182    ++idx;
183  }
184}
185
186void CDomainAlgorithmTransformation::domainGlobalIndex(const int& index, const int& niGlob, const int& njGlob,
187                                                       int& iIndex, int& jIndex)
188{
189   iIndex = index % niGlob;
190   jIndex = index / niGlob;
191}
192}
Note: See TracBrowser for help on using the repository browser.