source: XIOS/dev/dev_ym/XIOS_COUPLING/src/transformation/domain_algorithm/domain_algorithm_extract.cpp @ 2006

Last change on this file since 2006 was 2006, checked in by ymipsl, 3 years ago

Some fix domain_algorithm_extract. Now the domain destination not need to be a reference of the source grid.
YM
  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 9.8 KB
Line 
1#include "domain_algorithm_extract.hpp"
2#include "extract_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* CDomainAlgorithmExtract::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  CExtractDomain* extractDomain = dynamic_cast<CExtractDomain*> (transformation);
24  int domainDstIndex = elementPositionInGridDst2DomainPosition[elementPositionInGrid];
25  int domainSrcIndex = elementPositionInGridSrc2DomainPosition[elementPositionInGrid];
26
27  return (new CDomainAlgorithmExtract(isSource, domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], extractDomain));
28}
29CATCH
30
31bool CDomainAlgorithmExtract::dummyRegistered_ = CDomainAlgorithmExtract::registerTrans();
32bool CDomainAlgorithmExtract::registerTrans()
33TRY
34{
35  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_EXTRACT_DOMAIN, create);
36}
37CATCH
38
39CDomainAlgorithmExtract::CDomainAlgorithmExtract(bool isSource, CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)
40: CAlgorithmTransformationTransfer(isSource), domainSrc_(domainSource), domainDest_(domainDestination)
41TRY
42{
43  extractDomain->checkValid(domainSource);
44  extractIBegin_ = extractDomain->ibegin.getValue();
45  extractJBegin_ = extractDomain->jbegin.getValue();
46
47  extractNi_  = extractDomain->ni.getValue();
48  extractNj_  = extractDomain->nj.getValue();
49
50  extractIEnd_ = extractIBegin_ + extractNi_ - 1;
51  extractJEnd_ = extractJBegin_ + extractNj_ - 1;
52
53  if (extractNi_ > domainSource->ni_glo.getValue())
54  {
55    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
56           << "Extract 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           << "Extract size is " << extractNi_ );
59  }
60
61  if (extractNj_ > domainSource->nj_glo.getValue())
62  {
63    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
64           << "Extract 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           << "Extract size is " << extractNj_ );
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;
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 >= extractIBegin_) && (iIdxSrc <= extractIEnd_))
79      {
80        jIdxSrc = domainSrc_->j_index(ind);
81        if ((jIdxSrc >= extractJBegin_) && (jIdxSrc <= extractJEnd_))
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 - extractIBegin_;
96  jbeginDest = destJBegin + domainSrc_->jbegin - extractJBegin_;
97  if (niDest==0) ibeginDest=0 ;
98  if (njDest==0) jbeginDest=0 ;
99 
100  domainDest_->type = domainSrc_ -> type ;
101  domainDest_->data_dim = domainSrc_->data_dim ;
102  domainDest_->ni_glo.setValue(extractNi_);
103  domainDest_->nj_glo.setValue(extractNj_);
104  domainDest_->ni.setValue(niDest);
105  domainDest_->nj.setValue(njDest);
106  domainDest_->ibegin.setValue(ibeginDest);
107  domainDest_->jbegin.setValue(jbeginDest);
108  domainDest_->i_index.resize(niDest*njDest);
109  domainDest_->j_index.resize(niDest*njDest);
110
111  domainDest_->data_ni.setValue(niDest);
112  domainDest_->data_nj.setValue(njDest);
113  domainDest_->data_ibegin.setValue(0);  // local position
114  domainDest_->data_jbegin.setValue(0);  // local position
115  domainDest_->data_i_index.resize(niDest*njDest); // local position
116  domainDest_->data_j_index.resize(niDest*njDest); // local position
117
118  //domainDest_->domainMask.resize(niDest*njDest);
119
120  if (!domainSrc_->mask_2d.isEmpty()) domainDest_->mask_2d.resize(niDest,njDest);
121  if (!domainSrc_->mask_1d.isEmpty()) domainDest_->mask_1d.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  }
136  else if (!domainSrc_->lonvalue_2d.isEmpty())
137  {
138    domainDest_->lonvalue_2d.resize(niDest,njDest);
139    domainDest_->latvalue_2d.resize(niDest,njDest);
140  }
141
142  if (domainSrc_->hasBounds)
143  {
144    if (!domainSrc_->bounds_lon_2d.isEmpty())
145    {
146      domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest);
147      domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest);
148    }
149    else if (!domainSrc_->bounds_lon_1d.isEmpty())
150    {
151      domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest);
152      domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest);
153    }
154  }
155  if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest);
156 
157  for (int iDest = 0; iDest < niDest; iDest++)
158  {
159    iSrc = iDest + destIBegin;
160    for (int jDest = 0; jDest < njDest; jDest++)
161    {
162      jSrc = jDest + destJBegin;
163      ind = jSrc * domainSrc_->ni + iSrc;
164      iIdxSrc = domainSrc_->i_index(ind);
165      jIdxSrc = domainSrc_->j_index(ind);
166      indLocDest = jDest*niDest + iDest;
167      indGloDest = (jDest + jbeginDest)*extractNi_ + (iDest + ibeginDest);
168      indLocSrc = (jDest+destJBegin)*domainSrc_->ni + (iDest+destIBegin);
169      indGloSrc = (jIdxSrc )* niGloSrc + iIdxSrc;
170      domainDest_->i_index(indLocDest) = iDest + ibeginDest;                                             // i_index contains global positions
171      domainDest_->j_index(indLocDest) = jDest + jbeginDest;                                             // i_index contains global positions
172      domainDest_->data_i_index(indLocDest) = (domainSrc_->data_dim == 1) ? indLocDest : iDest;          // data_i_index contains local positions
173      domainDest_->data_j_index(indLocDest) = (domainSrc_->data_dim == 1) ? 0 :jDest;                    // data_i_index contains local positions
174      //domainDest_->domainMask(indLocDest) = domainSrc_->domainMask(indLocSrc);
175
176      if (!domainSrc_->mask_2d.isEmpty())
177        domainDest_->mask_2d(iDest,jDest) = domainSrc_->mask_2d(iSrc,jSrc);
178
179      if (!domainSrc_->mask_1d.isEmpty())
180        domainDest_->mask_1d(indLocDest) = domainSrc_->mask_1d(indLocSrc);
181
182      if (domainSrc_->hasArea)
183        domainDest_->area(iDest,jDest) = domainSrc_->area(iSrc,jSrc);
184
185
186      if (domainSrc_->hasBounds)
187      {
188        if (!domainSrc_->bounds_lon_2d.isEmpty())
189        {
190          for (int n = 0; n < domainSrc_->nvertex; ++n)
191          {
192            domainDest_->bounds_lon_2d(n,iDest,jDest) = domainSrc_->bounds_lon_2d(n,iSrc,jSrc);
193            domainDest_->bounds_lat_2d(n,iDest,jDest) = domainSrc_->bounds_lat_2d(n,iSrc,jSrc);
194          }
195        }
196        else if (!domainSrc_->bounds_lon_1d.isEmpty())
197        {
198          for (int n = 0; n < domainSrc_->nvertex; ++n)
199          {
200            domainDest_->bounds_lon_1d(n,iDest) = domainSrc_->bounds_lon_1d(n,iSrc);
201            domainDest_->bounds_lat_1d(n,iDest) = domainSrc_->bounds_lat_1d(n,iSrc);
202          }
203        }
204      }
205
206      if (domainSrc_->hasLonLat)
207      {
208        if (domainDest_->type == CDomain::type_attr::rectilinear)
209        {
210          domainDest_->latvalue_1d(jDest) = domainSrc_->latvalue_1d(jSrc);
211        }
212        else if (domainDest_->type == CDomain::type_attr::curvilinear)
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      transformationMapping_[indGloDest]=indGloSrc;
220
221    }
222    if (domainSrc_->hasLonLat)
223    {
224      if (domainDest_->type == CDomain::type_attr::unstructured)
225      {
226        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
227        domainDest_->latvalue_1d(iDest) = domainSrc_->latvalue_1d(iSrc);
228      }
229      else if (domainDest_->type == CDomain::type_attr::rectilinear)
230      {
231        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
232      }
233    }
234  }
235 
236  domainDestination->checkAttributes() ;
237  this->computeAlgorithm(domainSource->getLocalView(CElementView::WORKFLOW), domainDestination->getLocalView(CElementView::WORKFLOW)) ;
238}
239CATCH
240
241
242}
Note: See TracBrowser for help on using the repository browser.