source: XIOS3/trunk/src/transformation/domain_algorithm/domain_algorithm_zoom.cpp @ 2507

Last change on this file since 2507 was 2507, checked in by ymipsl, 13 months ago

Merging XIOS3_ATTACHED branch into XIOS3 trunk.

YM
  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 10.9 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 {
9shared_ptr<CGenericAlgorithmTransformation> CDomainAlgorithmZoom::create(bool isSource, 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 make_shared<CDomainAlgorithmZoom>(isSource, domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], zoomDomain);
28}
29CATCH
30
31bool CDomainAlgorithmZoom::dummyRegistered_ = CDomainAlgorithmZoom::registerTrans();
32bool CDomainAlgorithmZoom::registerTrans()
33TRY
34{
35  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_ZOOM_DOMAIN, create);
36}
37CATCH
38
39CDomainAlgorithmZoom::CDomainAlgorithmZoom(bool isSource, CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)
40: CAlgorithmTransformationTransfer(isSource), domainSrc_(domainSource), domainDest_(domainDestination)
41TRY
42{
43  zoomDomain->checkValid(domainSource);
44  zoomIBegin_ = zoomDomain->ibegin.getValue();
45  zoomJBegin_ = zoomDomain->jbegin.getValue();
46
47  zoomNi_  = zoomDomain->ni.getValue();
48  zoomNj_  = zoomDomain->nj.getValue();
49
50  zoomIEnd_ = zoomIBegin_ + zoomNi_ - 1;
51  zoomJEnd_ = zoomJBegin_ + zoomNj_ - 1;
52
53  if (zoomNi_ > domainSource->ni_glo.getValue())
54  {
55    ERROR("CDomainAlgorithmZoom::CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)",
56           << "Zoom size is greater than size of domain source"
57           << "Size ni_glo of domain source " <<domainSource->getId() << " is " << domainSource->ni_glo.getValue()  << std::endl
58           << "Zoom size is " << zoomNi_ );
59  }
60
61  if (zoomNj_ > domainSource->nj_glo.getValue())
62  {
63    ERROR("CDomainAlgorithmZoom::CDomainAlgorithmZoom(CDomain* domainDestination, CDomain* domainSource, CZoomDomain* zoomDomain)",
64           << "Zoom size is greater than size of domain source"
65           << "Size nj_glo of domain source " <<domainSource->getId() << " is " << domainSource->nj_glo.getValue()  << std::endl
66           << "Zoom size is " << zoomNj_ );
67  }
68
69  // Calculate the size of local domain
70  int ind, indLocSrc, indLocDest, iIdxSrc, jIdxSrc, destIBegin = -1, destJBegin = -1, niDest = 0, njDest = 0, ibeginDest, jbeginDest ;
71  int indGloDest, indGloSrc, niGloSrc = domainSrc_->ni_glo, iSrc, jSrc, nvertex = 0;
72  for (int j = 0; j < domainSrc_->nj.getValue(); j++)
73  {
74    for (int i = 0; i < domainSrc_->ni.getValue(); i++)
75    {
76      ind = j*domainSrc_->ni + i;
77      iIdxSrc = domainSrc_->i_index(ind);
78      if ((iIdxSrc >= zoomIBegin_) && (iIdxSrc <= zoomIEnd_))
79      {
80        jIdxSrc = domainSrc_->j_index(ind);
81        if ((jIdxSrc >= zoomJBegin_) && (jIdxSrc <= zoomJEnd_))
82        {
83          if ((niDest == 0) && (njDest == 0))
84          {
85            destIBegin = i;
86            destJBegin = j;
87          }
88          if (i == destIBegin) ++njDest;
89        }
90        if (j == destJBegin) ++niDest;
91
92      }
93    }
94  }
95  ibeginDest = destIBegin + domainSrc_->ibegin - zoomIBegin_;
96  jbeginDest = destJBegin + domainSrc_->jbegin - zoomJBegin_;
97  domainDest_->ni_glo.setValue(zoomNi_);
98  domainDest_->nj_glo.setValue(zoomNj_);
99  domainDest_->ni.setValue(niDest);
100  domainDest_->nj.setValue(njDest);
101  if ( (niDest==0) || (njDest==0))
102  {
103    domainDest_->ibegin.setValue(0);
104    domainDest_->jbegin.setValue(0);
105  }
106  else
107  {
108    domainDest_->ibegin.setValue(ibeginDest);
109    domainDest_->jbegin.setValue(jbeginDest);
110  }
111  domainDest_->i_index.resize(niDest*njDest);
112  domainDest_->j_index.resize(niDest*njDest);
113
114  domainDest_->data_ni.setValue(niDest);
115  domainDest_->data_nj.setValue(njDest);
116  domainDest_->data_ibegin.setValue(0);
117  domainDest_->data_jbegin.setValue(0);
118  domainDest_->data_i_index.resize(niDest*njDest);
119  domainDest_->data_j_index.resize(niDest*njDest);
120
121  domainDest_->domainMask.resize(niDest*njDest);
122
123  if (!domainSrc_->lonvalue_1d.isEmpty())
124  {
125    if (domainDest_->type == CDomain::type_attr::rectilinear)
126    {
127      domainDest_->lonvalue_1d.resize(niDest);
128      domainDest_->latvalue_1d.resize(njDest);
129    }
130    else if (domainDest_->type == CDomain::type_attr::unstructured)
131    {
132      domainDest_->lonvalue_1d.resize(niDest);
133      domainDest_->latvalue_1d.resize(niDest);
134    }
135    else if (domainDest_->type == CDomain::type_attr::curvilinear)
136    {
137      domainDest_->lonvalue_1d.resize(niDest*njDest);
138      domainDest_->latvalue_1d.resize(niDest*njDest);
139    }
140  }
141  else if (!domainSrc_->lonvalue_2d.isEmpty())
142  {
143    domainDest_->lonvalue_2d.resize(niDest,njDest);
144    domainDest_->latvalue_2d.resize(niDest,njDest);
145  }
146
147  if (domainSrc_->hasBounds)
148  {
149    nvertex = domainSrc_->nvertex;
150    domainDest_->nvertex.setValue(nvertex);
151    if (!domainSrc_->bounds_lon_1d.isEmpty())
152    {
153      if (domainDest_->type == CDomain::type_attr::rectilinear)
154      {
155        domainDest_->bounds_lon_1d.resize(nvertex, niDest);
156        domainDest_->bounds_lat_1d.resize(nvertex, njDest);
157      }
158      else if (domainDest_->type == CDomain::type_attr::unstructured)
159      {
160        domainDest_->bounds_lon_1d.resize(nvertex, niDest);
161        domainDest_->bounds_lat_1d.resize(nvertex, niDest);
162      }
163      else if (domainDest_->type == CDomain::type_attr::curvilinear)
164      {
165        domainDest_->bounds_lon_1d.resize(nvertex, niDest*njDest);
166        domainDest_->bounds_lat_1d.resize(nvertex, niDest*njDest);
167      }
168    }
169    else if (!domainSrc_->bounds_lon_2d.isEmpty())
170    {
171      domainDest_->bounds_lon_2d.resize(nvertex, niDest, njDest);
172      domainDest_->bounds_lat_2d.resize(nvertex, niDest, njDest);
173    }
174  }
175  if (domainSrc_->hasArea)
176  {
177    if (!domainSrc_->area_1d.isEmpty()) domainDest_->area_1d.resize(niDest*njDest);
178    else if (!domainSrc_->area_2d.isEmpty()) domainDest_->area_2d.resize(niDest,njDest);
179  }
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      { 
202         if (!domainSrc_->area_1d.isEmpty())      domainDest_->area_1d(indLocDest)  = domainSrc_->area_1d(ind);
203         else if (!domainSrc_->area_2d.isEmpty()) domainDest_->area_2d(iDest,jDest) = domainSrc_->area_2d(iSrc,jSrc);
204      }
205
206      if (domainSrc_->hasLonLat)
207      {
208        if (!domainSrc_->latvalue_1d.isEmpty())
209        {
210          if (domainDest_->type == CDomain::type_attr::rectilinear)
211          {
212            domainDest_->latvalue_1d(jDest) = domainSrc_->latvalue_1d(jSrc);
213          }
214          else
215          {
216            domainDest_->lonvalue_1d(indLocDest) = domainSrc_->lonvalue_1d(ind);
217            domainDest_->latvalue_1d(indLocDest) = domainSrc_->latvalue_1d(ind);
218          }
219        }
220        else if (!domainSrc_->latvalue_2d.isEmpty())
221        {
222          domainDest_->lonvalue_2d(iDest,jDest) = domainSrc_->lonvalue_2d(iSrc,jSrc);
223          domainDest_->latvalue_2d(iDest,jDest) = domainSrc_->latvalue_2d(iSrc,jSrc);
224        }
225      }
226
227      if (domainSrc_->hasBounds)
228      {
229        if (!domainSrc_->bounds_lon_1d.isEmpty())
230        {
231          if (domainDest_->type == CDomain::type_attr::rectilinear)
232          {
233            for (int n = 0; n < nvertex; ++n)
234              domainDest_->bounds_lat_1d(n,jDest) = domainSrc_->bounds_lat_1d(n,jSrc);
235          }
236          else
237          {
238            for (int n = 0; n < nvertex; ++n)
239            {
240              domainDest_->bounds_lon_1d(n,indLocDest) = domainSrc_->bounds_lon_1d(n,ind);
241              domainDest_->bounds_lat_1d(n,indLocDest) = domainSrc_->bounds_lat_1d(n,ind);
242            }
243          }
244        }
245        else if (!domainSrc_->bounds_lon_2d.isEmpty())
246        {
247          for (int n = 0; n < nvertex; ++n)
248          {
249            domainDest_->bounds_lon_2d(n,iDest,jDest) = domainSrc_->bounds_lon_2d(n,iSrc,jSrc);
250            domainDest_->bounds_lat_2d(n,iDest,jDest) = domainSrc_->bounds_lat_2d(n,iSrc,jSrc);
251          }
252        }
253
254      }
255
256      transformationMapping_[indGloDest]=indGloSrc;
257
258    }
259
260    if (domainSrc_->hasLonLat && !domainSrc_->latvalue_1d.isEmpty())
261    {
262      if (domainDest_->type == CDomain::type_attr::rectilinear)
263      {
264        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
265      }
266    }
267
268    if (domainSrc_->hasBounds && !domainSrc_->bounds_lon_1d.isEmpty())
269    {
270      if (domainDest_->type == CDomain::type_attr::rectilinear)
271      {
272        for (int n = 0; n < nvertex; ++n)
273          domainDest_->bounds_lon_1d(n,iDest) = domainSrc_->bounds_lon_1d(n,iSrc);
274      }
275    }
276  }
277  domainDest_->computeLocalMask();
278
279  domainDestination->checkAttributes() ;
280
281  this->computeAlgorithm(domainSource->getLocalView(CElementView::WORKFLOW), domainDestination->getLocalView(CElementView::WORKFLOW)) ;
282 
283}
284CATCH
285
286}
Note: See TracBrowser for help on using the repository browser.