source: XIOS/dev/dev_olga/src/transformation/domain_algorithm_zoom.cpp @ 1609

Last change on this file since 1609 was 1609, checked in by oabramkina, 17 months ago

Bugfix for domain zoom.

Values of ibegin and jbegin were not set properly for the destination domain. This was causing an error later on in case of multiple transformations when the destination grid of a transformation becomes the source grid for the next transformation.

File size: 10.7 KB
Line 
1#include "domain_algorithm_zoom.hpp"
2#include "zoom_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* CDomainAlgorithmZoom::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)
18{
19  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
20  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
21
22  CZoomDomain* zoomDomain = dynamic_cast<CZoomDomain*> (transformation);
23  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
24  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
25
26  return (new CDomainAlgorithmZoom(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], zoomDomain));
27}
28
29bool CDomainAlgorithmZoom::registerTrans()
30{
31  CGridTransformationFactory<CDomain>::registerTransformation(TRANS_ZOOM_DOMAIN, create);
32}
33
34CDomainAlgorithmZoom::CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)
35: CDomainAlgorithmTransformation(domainDestination, domainSource)
36{
37  zoomDomain->checkValid(domainSource);
38  zoomIBegin_ = zoomDomain->ibegin.getValue();
39  zoomJBegin_ = zoomDomain->jbegin.getValue();
40
41  zoomNi_  = zoomDomain->ni.getValue();
42  zoomNj_  = zoomDomain->nj.getValue();
43
44  zoomIEnd_ = zoomIBegin_ + zoomNi_ - 1;
45  zoomJEnd_ = zoomJBegin_ + zoomNj_ - 1;
46
47  if (zoomNi_ > domainSource->ni_glo.getValue())
48  {
49    ERROR("CDomainAlgorithmZoom::CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)",
50           << "Zoom size is greater than size of domain source"
51           << "Size ni_glo of domain source " <<domainSource->getId() << " is " << domainSource->ni_glo.getValue()  << std::endl
52           << "Zoom size is " << zoomNi_ );
53  }
54
55  if (zoomNj_ > domainSource->nj_glo.getValue())
56  {
57    ERROR("CDomainAlgorithmZoom::CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)",
58           << "Zoom size is greater than size of domain source"
59           << "Size nj_glo of domain source " <<domainSource->getId() << " is " << domainSource->nj_glo.getValue()  << std::endl
60           << "Zoom size is " << zoomNj_ );
61  }
62
63  // Calculate the size of local domain
64  int ind, indLocSrc, indLocDest, iIdxSrc, jIdxSrc, destIBegin = -1, destJBegin = -1, niDest = 0, njDest = 0, ibeginDest, jbeginDest ;
65  int indGloDest, indGloSrc, niGloSrc = domainSrc_->ni_glo, iSrc, jSrc, nvertex = 0;
66  for (int j = 0; j < domainSrc_->nj.getValue(); j++)
67  {
68    for (int i = 0; i < domainSrc_->ni.getValue(); i++)
69    {
70      ind = j*domainSrc_->ni + i;
71      iIdxSrc = domainSrc_->i_index(ind);
72      if ((iIdxSrc >= zoomIBegin_) && (iIdxSrc <= zoomIEnd_))
73      {
74        jIdxSrc = domainSrc_->j_index(ind);
75        if ((jIdxSrc >= zoomJBegin_) && (jIdxSrc <= zoomJEnd_))
76        {
77          if ((niDest == 0) && (njDest == 0))
78          {
79            destIBegin = i;
80            destJBegin = j;
81          }
82          if (i == destIBegin) ++njDest;
83        }
84        if (j == destJBegin) ++niDest;
85
86      }
87    }
88  }
89  ibeginDest = destIBegin + domainSrc_->ibegin - zoomIBegin_;
90  jbeginDest = destJBegin + domainSrc_->jbegin - zoomJBegin_;
91  domainDest_->ni_glo.setValue(zoomNi_);
92  domainDest_->nj_glo.setValue(zoomNj_);
93  domainDest_->ni.setValue(niDest);
94  domainDest_->nj.setValue(njDest);
95  if ( (niDest==0) || (njDest==0))
96  {
97    domainDest_->ibegin.setValue(0);
98    domainDest_->jbegin.setValue(0);
99  }
100  else
101  {
102    domainDest_->ibegin.setValue(ibeginDest);
103    domainDest_->jbegin.setValue(jbeginDest);
104  }
105  domainDest_->i_index.resize(niDest*njDest);
106  domainDest_->j_index.resize(niDest*njDest);
107
108  domainDest_->data_ni.setValue(niDest);
109  domainDest_->data_nj.setValue(njDest);
110  domainDest_->data_ibegin.setValue(0);  // local position
111  domainDest_->data_jbegin.setValue(0);  // local position
112  domainDest_->data_i_index.resize(niDest*njDest); // local position
113  domainDest_->data_j_index.resize(niDest*njDest); // local position
114
115  domainDest_->domainMask.resize(niDest*njDest);
116
117  if (!domainSrc_->lonvalue_1d.isEmpty())
118  {
119    if (domainDest_->type == CDomain::type_attr::rectilinear)
120    {
121      domainDest_->lonvalue_1d.resize(niDest);
122      domainDest_->latvalue_1d.resize(njDest);
123    }
124    else if (domainDest_->type == CDomain::type_attr::unstructured)
125    {
126      domainDest_->lonvalue_1d.resize(niDest);
127      domainDest_->latvalue_1d.resize(niDest);
128    }
129    else if (domainDest_->type == CDomain::type_attr::curvilinear)
130    {
131      domainDest_->lonvalue_1d.resize(niDest*njDest);
132      domainDest_->latvalue_1d.resize(niDest*njDest);
133    }
134  }
135  else if (!domainSrc_->lonvalue_2d.isEmpty())
136  {
137    domainDest_->lonvalue_2d.resize(niDest,njDest);
138    domainDest_->latvalue_2d.resize(niDest,njDest);
139  }
140
141  if (domainSrc_->hasBounds)
142  {
143    nvertex = domainSrc_->nvertex;
144    domainDest_->nvertex.setValue(nvertex);
145    if (!domainSrc_->bounds_lon_1d.isEmpty())
146    {
147      if (domainDest_->type == CDomain::type_attr::rectilinear)
148      {
149        domainDest_->bounds_lon_1d.resize(nvertex, niDest);
150        domainDest_->bounds_lat_1d.resize(nvertex, njDest);
151      }
152      else if (domainDest_->type == CDomain::type_attr::unstructured)
153      {
154        domainDest_->bounds_lon_1d.resize(nvertex, niDest);
155        domainDest_->bounds_lat_1d.resize(nvertex, niDest);
156      }
157      else if (domainDest_->type == CDomain::type_attr::curvilinear)
158      {
159        domainDest_->bounds_lon_1d.resize(nvertex, niDest*njDest);
160        domainDest_->bounds_lat_1d.resize(nvertex, niDest*njDest);
161      }
162    }
163    else if (!domainSrc_->bounds_lon_2d.isEmpty())
164    {
165      domainDest_->bounds_lon_2d.resize(nvertex, niDest, njDest);
166      domainDest_->bounds_lat_2d.resize(nvertex, niDest, njDest);
167    }
168  }
169  if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest);
170
171  this->transformationMapping_.resize(1);
172  this->transformationWeight_.resize(1);
173  TransformationIndexMap& transMap = this->transformationMapping_[0];
174  TransformationWeightMap& transWeight = this->transformationWeight_[0];
175
176  for (int iDest = 0; iDest < niDest; iDest++)
177  {
178    iSrc = iDest + destIBegin;
179    for (int jDest = 0; jDest < njDest; jDest++)
180    {
181      jSrc = jDest + destJBegin;
182      ind = jSrc * domainSrc_->ni + iSrc;
183      iIdxSrc = domainSrc_->i_index(ind);
184      jIdxSrc = domainSrc_->j_index(ind);
185      indLocDest = jDest*niDest + iDest;
186      indGloDest = (jDest + jbeginDest)*zoomNi_ + (iDest + ibeginDest);
187      indLocSrc = (jDest+destJBegin)*domainSrc_->ni + (iDest+destIBegin);
188      indGloSrc = (jIdxSrc )* niGloSrc + iIdxSrc;
189      domainDest_->i_index(indLocDest) = iDest + ibeginDest;                                             // i_index contains global positions
190      domainDest_->j_index(indLocDest) = jDest + jbeginDest;                                             // i_index contains global positions
191      domainDest_->data_i_index(indLocDest) = (domainSrc_->data_dim == 1) ? indLocDest : iDest;          // data_i_index contains local positions
192      domainDest_->data_j_index(indLocDest) = (domainSrc_->data_dim == 1) ? 0 :jDest;                    // data_i_index contains local positions
193      domainDest_->domainMask(indLocDest) = domainSrc_->domainMask(indLocSrc);
194
195      if (domainSrc_->hasArea)
196        domainDest_->area(iDest,jDest) = domainSrc_->area(iSrc,jSrc);
197
198      if (domainSrc_->hasLonLat)
199      {
200        if (!domainSrc_->latvalue_1d.isEmpty())
201        {
202          if (domainDest_->type == CDomain::type_attr::rectilinear)
203          {
204            domainDest_->latvalue_1d(jDest) = domainSrc_->latvalue_1d(jSrc);
205          }
206          else
207          {
208            domainDest_->lonvalue_1d(indLocDest) = domainSrc_->lonvalue_1d(ind);
209            domainDest_->latvalue_1d(indLocDest) = domainSrc_->latvalue_1d(ind);
210          }
211        }
212        else if (!domainSrc_->latvalue_2d.isEmpty())
213        {
214          domainDest_->lonvalue_2d(iDest,jDest) = domainSrc_->lonvalue_2d(iSrc,jSrc);
215          domainDest_->latvalue_2d(iDest,jDest) = domainSrc_->latvalue_2d(iSrc,jSrc);
216        }
217      }
218
219      if (domainSrc_->hasBounds)
220      {
221        if (!domainSrc_->bounds_lon_1d.isEmpty())
222        {
223          if (domainDest_->type == CDomain::type_attr::rectilinear)
224          {
225            for (int n = 0; n < nvertex; ++n)
226              domainDest_->bounds_lat_1d(n,jDest) = domainSrc_->bounds_lat_1d(n,jSrc);
227          }
228          else
229          {
230            for (int n = 0; n < nvertex; ++n)
231            {
232              domainDest_->bounds_lon_1d(n,indLocDest) = domainSrc_->bounds_lon_1d(n,ind);
233              domainDest_->bounds_lat_1d(n,indLocDest) = domainSrc_->bounds_lat_1d(n,ind);
234            }
235          }
236        }
237        else if (!domainSrc_->bounds_lon_2d.isEmpty())
238        {
239          for (int n = 0; n < nvertex; ++n)
240          {
241            domainDest_->bounds_lon_2d(n,iDest,jDest) = domainSrc_->bounds_lon_2d(n,iSrc,jSrc);
242            domainDest_->bounds_lat_2d(n,iDest,jDest) = domainSrc_->bounds_lat_2d(n,iSrc,jSrc);
243          }
244        }
245
246      }
247
248      transMap[indGloDest].push_back(indGloSrc);
249      transWeight[indGloDest].push_back(1.0);
250    }
251
252    if (domainSrc_->hasLonLat && !domainSrc_->latvalue_1d.isEmpty())
253    {
254      if (domainDest_->type == CDomain::type_attr::rectilinear)
255      {
256        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
257      }
258    }
259
260    if (domainSrc_->hasBounds && !domainSrc_->bounds_lon_1d.isEmpty())
261    {
262      if (domainDest_->type == CDomain::type_attr::rectilinear)
263      {
264        for (int n = 0; n < nvertex; ++n)
265          domainDest_->bounds_lon_1d(n,iDest) = domainSrc_->bounds_lon_1d(n,iSrc);
266      }
267    }
268  }
269
270}
271
272/*!
273  Compute the index mapping between domain on grid source and one on grid destination
274*/
275void CDomainAlgorithmZoom::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
276{
277}
278
279
280}
Note: See TracBrowser for help on using the repository browser.