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

Last change on this file since 1612 was 1612, checked in by oabramkina, 10 months ago

Dev: adding exception handling.

To activate it, compilation flag -DXIOS_EXCEPTION should be added.

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