source: XIOS/trunk/src/transformation/domain_algorithm_zoom.cpp @ 1583

Last change on this file since 1583 was 1583, checked in by oabramkina, 3 years ago

Bugfix on domain zoom for curvilinear domains with 1D lon/lat.

File size: 10.5 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  domainDest_->ibegin.setValue(ibeginDest);
96  domainDest_->jbegin.setValue(jbeginDest);
97  domainDest_->i_index.resize(niDest*njDest);
98  domainDest_->j_index.resize(niDest*njDest);
99
100  domainDest_->data_ni.setValue(niDest);
101  domainDest_->data_nj.setValue(njDest);
102  domainDest_->data_ibegin.setValue(0);  // local position
103  domainDest_->data_jbegin.setValue(0);  // local position
104  domainDest_->data_i_index.resize(niDest*njDest); // local position
105  domainDest_->data_j_index.resize(niDest*njDest); // local position
106
107  domainDest_->domainMask.resize(niDest*njDest);
108
109  if (!domainSrc_->lonvalue_1d.isEmpty())
110  {
111    if (domainDest_->type == CDomain::type_attr::rectilinear)
112    {
113      domainDest_->lonvalue_1d.resize(niDest);
114      domainDest_->latvalue_1d.resize(njDest);
115    }
116    else if (domainDest_->type == CDomain::type_attr::unstructured)
117    {
118      domainDest_->lonvalue_1d.resize(niDest);
119      domainDest_->latvalue_1d.resize(niDest);
120    }
121    else if (domainDest_->type == CDomain::type_attr::curvilinear)
122    {
123      domainDest_->lonvalue_1d.resize(niDest*njDest);
124      domainDest_->latvalue_1d.resize(niDest*njDest);
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    nvertex = domainSrc_->nvertex;
136    domainDest_->nvertex.setValue(nvertex);
137    if (!domainSrc_->bounds_lon_1d.isEmpty())
138    {
139      if (domainDest_->type == CDomain::type_attr::rectilinear)
140      {
141        domainDest_->bounds_lon_1d.resize(nvertex, niDest);
142        domainDest_->bounds_lat_1d.resize(nvertex, njDest);
143      }
144      else if (domainDest_->type == CDomain::type_attr::unstructured)
145      {
146        domainDest_->bounds_lon_1d.resize(nvertex, niDest);
147        domainDest_->bounds_lat_1d.resize(nvertex, niDest);
148      }
149      else if (domainDest_->type == CDomain::type_attr::curvilinear)
150      {
151        domainDest_->bounds_lon_1d.resize(nvertex, niDest*njDest);
152        domainDest_->bounds_lat_1d.resize(nvertex, niDest*njDest);
153      }
154    }
155    else if (!domainSrc_->bounds_lon_2d.isEmpty())
156    {
157      domainDest_->bounds_lon_2d.resize(nvertex, niDest, njDest);
158      domainDest_->bounds_lat_2d.resize(nvertex, niDest, njDest);
159    }
160  }
161  if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest);
162
163  this->transformationMapping_.resize(1);
164  this->transformationWeight_.resize(1);
165  TransformationIndexMap& transMap = this->transformationMapping_[0];
166  TransformationWeightMap& transWeight = this->transformationWeight_[0];
167
168  for (int iDest = 0; iDest < niDest; iDest++)
169  {
170    iSrc = iDest + destIBegin;
171    for (int jDest = 0; jDest < njDest; jDest++)
172    {
173      jSrc = jDest + destJBegin;
174      ind = jSrc * domainSrc_->ni + iSrc;
175      iIdxSrc = domainSrc_->i_index(ind);
176      jIdxSrc = domainSrc_->j_index(ind);
177      indLocDest = jDest*niDest + iDest;
178      indGloDest = (jDest + jbeginDest)*zoomNi_ + (iDest + ibeginDest);
179      indLocSrc = (jDest+destJBegin)*domainSrc_->ni + (iDest+destIBegin);
180      indGloSrc = (jIdxSrc )* niGloSrc + iIdxSrc;
181      domainDest_->i_index(indLocDest) = iDest + ibeginDest;                                             // i_index contains global positions
182      domainDest_->j_index(indLocDest) = jDest + jbeginDest;                                             // i_index contains global positions
183      domainDest_->data_i_index(indLocDest) = (domainSrc_->data_dim == 1) ? indLocDest : iDest;          // data_i_index contains local positions
184      domainDest_->data_j_index(indLocDest) = (domainSrc_->data_dim == 1) ? 0 :jDest;                    // data_i_index contains local positions
185      domainDest_->domainMask(indLocDest) = domainSrc_->domainMask(indLocSrc);
186
187      if (domainSrc_->hasArea)
188        domainDest_->area(iDest,jDest) = domainSrc_->area(iSrc,jSrc);
189
190      if (domainSrc_->hasLonLat)
191      {
192        if (!domainSrc_->latvalue_1d.isEmpty())
193        {
194          if (domainDest_->type == CDomain::type_attr::rectilinear)
195          {
196            domainDest_->latvalue_1d(jDest) = domainSrc_->latvalue_1d(jSrc);
197          }
198          else
199          {
200            domainDest_->lonvalue_1d(indLocDest) = domainSrc_->lonvalue_1d(ind);
201            domainDest_->latvalue_1d(indLocDest) = domainSrc_->latvalue_1d(ind);
202          }
203        }
204        else if (!domainSrc_->latvalue_2d.isEmpty())
205        {
206          domainDest_->lonvalue_2d(iDest,jDest) = domainSrc_->lonvalue_2d(iSrc,jSrc);
207          domainDest_->latvalue_2d(iDest,jDest) = domainSrc_->latvalue_2d(iSrc,jSrc);
208        }
209      }
210
211      if (domainSrc_->hasBounds)
212      {
213        if (!domainSrc_->bounds_lon_1d.isEmpty())
214        {
215          if (domainDest_->type == CDomain::type_attr::rectilinear)
216          {
217            for (int n = 0; n < nvertex; ++n)
218              domainDest_->bounds_lat_1d(n,jDest) = domainSrc_->bounds_lat_1d(n,jSrc);
219          }
220          else
221          {
222            for (int n = 0; n < nvertex; ++n)
223            {
224              domainDest_->bounds_lon_1d(n,indLocDest) = domainSrc_->bounds_lon_1d(n,ind);
225              domainDest_->bounds_lat_1d(n,indLocDest) = domainSrc_->bounds_lat_1d(n,ind);
226            }
227          }
228        }
229        else if (!domainSrc_->bounds_lon_2d.isEmpty())
230        {
231          for (int n = 0; n < nvertex; ++n)
232          {
233            domainDest_->bounds_lon_2d(n,iDest,jDest) = domainSrc_->bounds_lon_2d(n,iSrc,jSrc);
234            domainDest_->bounds_lat_2d(n,iDest,jDest) = domainSrc_->bounds_lat_2d(n,iSrc,jSrc);
235          }
236        }
237
238      }
239
240      transMap[indGloDest].push_back(indGloSrc);
241      transWeight[indGloDest].push_back(1.0);
242    }
243
244    if (domainSrc_->hasLonLat && !domainSrc_->latvalue_1d.isEmpty())
245    {
246      if (domainDest_->type == CDomain::type_attr::rectilinear)
247      {
248        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
249      }
250    }
251
252    if (domainSrc_->hasBounds && !domainSrc_->bounds_lon_1d.isEmpty())
253    {
254      if (domainDest_->type == CDomain::type_attr::rectilinear)
255      {
256        for (int n = 0; n < nvertex; ++n)
257          domainDest_->bounds_lon_1d(n,iDest) = domainSrc_->bounds_lon_1d(n,iSrc);
258      }
259    }
260  }
261
262}
263
264/*!
265  Compute the index mapping between domain on grid source and one on grid destination
266*/
267void CDomainAlgorithmZoom::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
268{
269}
270
271
272}
Note: See TracBrowser for help on using the repository browser.