source: XIOS2/trunk/src/transformation/domain_algorithm_extract.cpp @ 2442

Last change on this file since 2442 was 2442, checked in by jderouillat, 19 months ago

Fix bounds management in transformations. Disabled temporarily bounds settings in the rectilinear case of the generic_testcase
File size: 9.5 KB
Line 
1#include "domain_algorithm_extract.hpp"
2#include "extract_domain.hpp"
3#include "domain.hpp"
4#include "grid.hpp"
5#include "grid_transformation_factory_impl.hpp"
6#include "attribute_template.hpp"
7
8namespace xios {
9CGenericAlgorithmTransformation* CDomainAlgorithmExtract::create(CGrid* gridDst, CGrid* gridSrc,
10                                                             CTransformation<CDomain>* transformation,
11                                                             int elementPositionInGrid,
12                                                             std::map<int, int>& elementPositionInGridSrc2ScalarPosition,
13                                                             std::map<int, int>& elementPositionInGridSrc2AxisPosition,
14                                                             std::map<int, int>& elementPositionInGridSrc2DomainPosition,
15                                                             std::map<int, int>& elementPositionInGridDst2ScalarPosition,
16                                                             std::map<int, int>& elementPositionInGridDst2AxisPosition,
17                                                             std::map<int, int>& elementPositionInGridDst2DomainPosition)
18TRY
19{
20  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
21  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
22
23  CExtractDomain* extractDomain = dynamic_cast<CExtractDomain*> (transformation);
24  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
25  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
26
27  return (new CDomainAlgorithmExtract(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], extractDomain));
28}
29CATCH
30
31bool CDomainAlgorithmExtract::registerTrans()
32TRY
33{
34  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_EXTRACT_DOMAIN, create);
35}
36CATCH
37
38CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)
39: CDomainAlgorithmTransformation(domainDestination, domainSource)
40TRY
41{
42  extractDomain->checkValid(domainSource);
43  extractIBegin_ = extractDomain->ibegin.getValue();
44  extractJBegin_ = extractDomain->jbegin.getValue();
45
46  extractNi_  = extractDomain->ni.getValue();
47  extractNj_  = extractDomain->nj.getValue();
48
49  extractIEnd_ = extractIBegin_ + extractNi_ - 1;
50  extractJEnd_ = extractJBegin_ + extractNj_ - 1;
51
52  if (extractNi_ > domainSource->ni_glo.getValue())
53  {
54    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
55           << "Extract size is greater than size of domain source"
56           << "Size ni_glo of domain source " <<domainSource->getId() << " is " << domainSource->ni_glo.getValue()  << std::endl
57           << "Extract size is " << extractNi_ );
58  }
59
60  if (extractNj_ > domainSource->nj_glo.getValue())
61  {
62    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
63           << "Extract size is greater than size of domain source"
64           << "Size nj_glo of domain source " <<domainSource->getId() << " is " << domainSource->nj_glo.getValue()  << std::endl
65           << "Extract size is " << extractNj_ );
66  }
67
68  // Calculate the size of local domain
69  int ind, indLocSrc, indLocDest, iIdxSrc, jIdxSrc, destIBegin = -1, destJBegin = -1, niDest = 0, njDest = 0, ibeginDest, jbeginDest ;
70  int indGloDest, indGloSrc, niGloSrc = domainSrc_->ni_glo, iSrc, jSrc;
71  for (int j = 0; j < domainSrc_->nj.getValue(); j++)
72  {
73    for (int i = 0; i < domainSrc_->ni.getValue(); i++)
74    {
75      ind = j*domainSrc_->ni + i;
76      iIdxSrc = domainSrc_->i_index(ind);
77      if ((iIdxSrc >= extractIBegin_) && (iIdxSrc <= extractIEnd_))
78      {
79        jIdxSrc = domainSrc_->j_index(ind);
80        if ((jIdxSrc >= extractJBegin_) && (jIdxSrc <= extractJEnd_))
81        {
82          if ((niDest == 0) && (njDest == 0))
83          {
84            destIBegin = i;
85            destJBegin = j;
86          }
87          if (i == destIBegin) ++njDest;
88        }
89        if (j == destJBegin) ++niDest;
90
91      }
92    }
93  }
94  ibeginDest = destIBegin + domainSrc_->ibegin - extractIBegin_;
95  jbeginDest = destJBegin + domainSrc_->jbegin - extractJBegin_;
96  domainDest_->ni_glo.setValue(extractNi_);
97  domainDest_->nj_glo.setValue(extractNj_);
98  domainDest_->ni.setValue(niDest);
99  domainDest_->nj.setValue(njDest);
100  domainDest_->ibegin.setValue(ibeginDest);
101  domainDest_->jbegin.setValue(jbeginDest);
102  domainDest_->i_index.resize(niDest*njDest);
103  domainDest_->j_index.resize(niDest*njDest);
104
105  domainDest_->data_ni.setValue(niDest);
106  domainDest_->data_nj.setValue(njDest);
107  domainDest_->data_ibegin.setValue(0);  // local position
108  domainDest_->data_jbegin.setValue(0);  // local position
109  domainDest_->data_i_index.resize(niDest*njDest); // local position
110  domainDest_->data_j_index.resize(niDest*njDest); // local position
111
112  domainDest_->domainMask.resize(niDest*njDest);
113
114  if (!domainSrc_->lonvalue_1d.isEmpty())
115  {
116    if (domainDest_->type == CDomain::type_attr::rectilinear)
117    {
118      domainDest_->lonvalue_1d.resize(niDest);
119      domainDest_->latvalue_1d.resize(njDest);
120    }
121    else if (domainDest_->type == CDomain::type_attr::unstructured)
122    {
123      domainDest_->lonvalue_1d.resize(niDest);
124      domainDest_->latvalue_1d.resize(niDest);
125    }
126  }
127  else if (!domainSrc_->lonvalue_2d.isEmpty())
128  {
129    domainDest_->lonvalue_2d.resize(niDest,njDest);
130    domainDest_->latvalue_2d.resize(niDest,njDest);
131  }
132
133  if (domainSrc_->hasBounds)
134  {
135    if (!domainSrc_->bounds_lon_2d.isEmpty())
136    {
137      domainDest_->bounds_lon_2d.resize(domainSrc_->bounds_lon_2d.shape()[0], niDest, njDest);
138      domainDest_->bounds_lon_2d.resize(domainSrc_->bounds_lon_2d.shape()[0], niDest, njDest);
139    }
140    else if (!domainSrc_->bounds_lon_1d.isEmpty())
141    {
142      domainDest_->bounds_lon_1d.resize(domainSrc_->bounds_lon_1d.shape()[0], niDest);
143      domainDest_->bounds_lon_1d.resize(domainSrc_->bounds_lon_1d.shape()[0], niDest);
144    }
145  }
146  if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest);
147
148  this->transformationMapping_.resize(1);
149  this->transformationWeight_.resize(1);
150  TransformationIndexMap& transMap = this->transformationMapping_[0];
151  TransformationWeightMap& transWeight = this->transformationWeight_[0];
152
153  for (int iDest = 0; iDest < niDest; iDest++)
154  {
155    iSrc = iDest + destIBegin;
156    for (int jDest = 0; jDest < njDest; jDest++)
157    {
158      jSrc = jDest + destJBegin;
159      ind = jSrc * domainSrc_->ni + iSrc;
160      iIdxSrc = domainSrc_->i_index(ind);
161      jIdxSrc = domainSrc_->j_index(ind);
162      indLocDest = jDest*niDest + iDest;
163      indGloDest = (jDest + jbeginDest)*extractNi_ + (iDest + ibeginDest);
164      indLocSrc = (jDest+destJBegin)*domainSrc_->ni + (iDest+destIBegin);
165      indGloSrc = (jIdxSrc )* niGloSrc + iIdxSrc;
166      domainDest_->i_index(indLocDest) = iDest + ibeginDest;                                             // i_index contains global positions
167      domainDest_->j_index(indLocDest) = jDest + jbeginDest;                                             // i_index contains global positions
168      domainDest_->data_i_index(indLocDest) = (domainSrc_->data_dim == 1) ? indLocDest : iDest;          // data_i_index contains local positions
169      domainDest_->data_j_index(indLocDest) = (domainSrc_->data_dim == 1) ? 0 :jDest;                    // data_i_index contains local positions
170      domainDest_->domainMask(indLocDest) = domainSrc_->domainMask(indLocSrc);
171
172      if (domainSrc_->hasArea)
173        domainDest_->area(iDest,jDest) = domainSrc_->area(iSrc,jSrc);
174
175      if (domainSrc_->hasBounds)
176      {
177        if (!domainSrc_->bounds_lon_2d.isEmpty())
178        {
179          for (int n = 0; n < domainSrc_->bounds_lon_2d.shape()[0]; ++n)
180          {
181            domainDest_->bounds_lon_2d(n,iDest,jDest) = domainSrc_->bounds_lon_2d(n,iSrc,jSrc);
182            domainDest_->bounds_lat_2d(n,iDest,jDest) = domainSrc_->bounds_lat_2d(n,iSrc,jSrc);
183          }
184        }
185        else if (!domainSrc_->bounds_lon_1d.isEmpty())
186        {
187          for (int n = 0; n < domainSrc_->bounds_lon_1d.shape()[0]; ++n)
188          {
189            domainDest_->bounds_lon_1d(n,iDest) = domainSrc_->bounds_lon_1d(n,iSrc);
190            domainDest_->bounds_lat_1d(n,iDest) = domainSrc_->bounds_lat_1d(n,iSrc);
191          }
192        }
193      }
194
195      if (domainSrc_->hasLonLat)
196      {
197        if (domainDest_->type == CDomain::type_attr::rectilinear)
198        {
199          domainDest_->latvalue_1d(jDest) = domainSrc_->latvalue_1d(jSrc);
200        }
201        else if (domainDest_->type == CDomain::type_attr::curvilinear)
202        {
203          domainDest_->lonvalue_2d(iDest,jDest) = domainSrc_->lonvalue_2d(iSrc,jSrc);
204          domainDest_->latvalue_2d(iDest,jDest) = domainSrc_->latvalue_2d(iSrc,jSrc);
205        }
206      }
207
208      transMap[indGloDest].push_back(indGloSrc);
209      transWeight[indGloDest].push_back(1.0);
210
211    }
212    if (domainSrc_->hasLonLat)
213    {
214      if (domainDest_->type == CDomain::type_attr::unstructured)
215      {
216        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
217        domainDest_->latvalue_1d(iDest) = domainSrc_->latvalue_1d(iSrc);
218      }
219      else if (domainDest_->type == CDomain::type_attr::rectilinear)
220      {
221        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
222      }
223    }
224  }
225
226}
227CATCH
228
229/*!
230  Compute the index mapping between domain on grid source and one on grid destination
231*/
232void CDomainAlgorithmExtract::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
233{
234}
235
236
237}
Note: See TracBrowser for help on using the repository browser.