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

Last change on this file since 1620 was 1620, checked in by oabramkina, 16 months ago

Dev: adding exception handling.

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

File size: 10.7 KB
RevLine 
[631]1#include "domain_algorithm_zoom.hpp"
[933]2#include "zoom_domain.hpp"
3#include "domain.hpp"
4#include "grid.hpp"
5#include "grid_transformation_factory_impl.hpp"
[1553]6#include "attribute_template.hpp"
[631]7
8namespace xios {
[933]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)
[1612]18TRY
[933]19{
20  std::vector<CDomain*> domainListDestP = gridDst->getDomains();
21  std::vector<CDomain*> domainListSrcP  = gridSrc->getDomains();
[631]22
[933]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}
[1612]29CATCH
[933]30
31bool CDomainAlgorithmZoom::registerTrans()
[1612]32TRY
[933]33{
34  CGridTransformationFactory<CDomain>::registerTransformation(TRANS_ZOOM_DOMAIN, create);
35}
[1612]36CATCH
[933]37
[631]38CDomainAlgorithmZoom::CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)
39: CDomainAlgorithmTransformation(domainDestination, domainSource)
[1612]40TRY
[631]41{
42  zoomDomain->checkValid(domainSource);
[787]43  zoomIBegin_ = zoomDomain->ibegin.getValue();
44  zoomJBegin_ = zoomDomain->jbegin.getValue();
[631]45
[787]46  zoomNi_  = zoomDomain->ni.getValue();
47  zoomNj_  = zoomDomain->nj.getValue();
[631]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
[1553]68  // Calculate the size of local domain
69  int ind, indLocSrc, indLocDest, iIdxSrc, jIdxSrc, destIBegin = -1, destJBegin = -1, niDest = 0, njDest = 0, ibeginDest, jbeginDest ;
[1583]70  int indGloDest, indGloSrc, niGloSrc = domainSrc_->ni_glo, iSrc, jSrc, nvertex = 0;
[1553]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;
[631]90
[1553]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);
[1609]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  }
[1553]110  domainDest_->i_index.resize(niDest*njDest);
111  domainDest_->j_index.resize(niDest*njDest);
[827]112
[1553]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    }
[1583]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    }
[1553]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  {
[1583]148    nvertex = domainSrc_->nvertex;
149    domainDest_->nvertex.setValue(nvertex);
150    if (!domainSrc_->bounds_lon_1d.isEmpty())
[1553]151    {
[1583]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      }
[1553]167    }
[1583]168    else if (!domainSrc_->bounds_lon_2d.isEmpty())
[1553]169    {
[1583]170      domainDest_->bounds_lon_2d.resize(nvertex, niDest, njDest);
171      domainDest_->bounds_lat_2d.resize(nvertex, niDest, njDest);
[1553]172    }
173  }
174  if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest);
175
[827]176  this->transformationMapping_.resize(1);
177  this->transformationWeight_.resize(1);
[833]178  TransformationIndexMap& transMap = this->transformationMapping_[0];
179  TransformationWeightMap& transWeight = this->transformationWeight_[0];
[827]180
[1553]181  for (int iDest = 0; iDest < niDest; iDest++)
[631]182  {
[1553]183    iSrc = iDest + destIBegin;
184    for (int jDest = 0; jDest < njDest; jDest++)
[631]185    {
[1553]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
[1583]203      if (domainSrc_->hasLonLat)
[1553]204      {
[1583]205        if (!domainSrc_->latvalue_1d.isEmpty())
[1553]206        {
[1583]207          if (domainDest_->type == CDomain::type_attr::rectilinear)
[1553]208          {
[1583]209            domainDest_->latvalue_1d(jDest) = domainSrc_->latvalue_1d(jSrc);
[1553]210          }
[1583]211          else
[1553]212          {
[1583]213            domainDest_->lonvalue_1d(indLocDest) = domainSrc_->lonvalue_1d(ind);
214            domainDest_->latvalue_1d(indLocDest) = domainSrc_->latvalue_1d(ind);
[1553]215          }
216        }
[1583]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        }
[1553]222      }
223
[1583]224      if (domainSrc_->hasBounds)
[1553]225      {
[1583]226        if (!domainSrc_->bounds_lon_1d.isEmpty())
[1553]227        {
[1583]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          }
[1553]241        }
[1583]242        else if (!domainSrc_->bounds_lon_2d.isEmpty())
[1553]243        {
[1583]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          }
[1553]249        }
[1583]250
[1553]251      }
252
253      transMap[indGloDest].push_back(indGloSrc);
254      transWeight[indGloDest].push_back(1.0);
[1583]255    }
[1553]256
[1583]257    if (domainSrc_->hasLonLat && !domainSrc_->latvalue_1d.isEmpty())
[1553]258    {
[1583]259      if (domainDest_->type == CDomain::type_attr::rectilinear)
[1553]260      {
261        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
262      }
[1583]263    }
264
265    if (domainSrc_->hasBounds && !domainSrc_->bounds_lon_1d.isEmpty())
266    {
267      if (domainDest_->type == CDomain::type_attr::rectilinear)
[1553]268      {
[1583]269        for (int n = 0; n < nvertex; ++n)
270          domainDest_->bounds_lon_1d(n,iDest) = domainSrc_->bounds_lon_1d(n,iSrc);
[1553]271      }
272    }
[631]273  }
[1620]274  domainDest_->computeLocalMask();
[631]275}
[1612]276CATCH
[631]277
278/*!
[1553]279  Compute the index mapping between domain on grid source and one on grid destination
[631]280*/
[1553]281void CDomainAlgorithmZoom::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
[631]282{
283}
284
285
286}
Note: See TracBrowser for help on using the repository browser.