source: XIOS/trunk/src/transformation/domain_algorithm_extract.cpp @ 1852

Last change on this file since 1852 was 1852, checked in by ymipsl, 4 years ago

Compiler fix : solve the problem of crash occured with recent version of GCC, only in optimised mode > O1
It seems to be due to non return value from a non void function in case of early initialization (static initialization).
Thanks to A. Durocher who find the tip.

YM

File size: 9.4 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(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(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], extractDomain));
28}
29CATCH
30
31bool CDomainAlgorithmExtract::registerTrans()
32TRY
33{
34  return CGridTransformationFactory<CDomain>::registerTransformation(TRANS_EXTRACT_DOMAIN, create);
35}
36CATCH
37
38CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)
39: CDomainAlgorithmTransformation(domainDestination, domainSource)
40TRY
41{
42  extractDomain->checkValid(domainSource);
43  extractIBegin_ = extractDomain->ibegin.getValue();
44  extractJBegin_ = extractDomain->jbegin.getValue();
45
46  extractNi_  = extractDomain->ni.getValue();
47  extractNj_  = extractDomain->nj.getValue();
48
49  extractIEnd_ = extractIBegin_ + extractNi_ - 1;
50  extractJEnd_ = extractJBegin_ + extractNj_ - 1;
51
52  if (extractNi_ > domainSource->ni_glo.getValue())
53  {
54    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
55           << "Extract 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           << "Extract size is " << extractNi_ );
58  }
59
60  if (extractNj_ > domainSource->nj_glo.getValue())
61  {
62    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
63           << "Extract 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           << "Extract size is " << extractNj_ );
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;
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 >= extractIBegin_) && (iIdxSrc <= extractIEnd_))
78      {
79        jIdxSrc = domainSrc_->j_index(ind);
80        if ((jIdxSrc >= extractJBegin_) && (jIdxSrc <= extractJEnd_))
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 - extractIBegin_;
95  jbeginDest = destJBegin + domainSrc_->jbegin - extractJBegin_;
96  domainDest_->ni_glo.setValue(extractNi_);
97  domainDest_->nj_glo.setValue(extractNj_);
98  domainDest_->ni.setValue(niDest);
99  domainDest_->nj.setValue(njDest);
100  domainDest_->ibegin.setValue(ibeginDest);
101  domainDest_->jbegin.setValue(jbeginDest);
102  domainDest_->i_index.resize(niDest*njDest);
103  domainDest_->j_index.resize(niDest*njDest);
104
105  domainDest_->data_ni.setValue(niDest);
106  domainDest_->data_nj.setValue(njDest);
107  domainDest_->data_ibegin.setValue(0);  // local position
108  domainDest_->data_jbegin.setValue(0);  // local position
109  domainDest_->data_i_index.resize(niDest*njDest); // local position
110  domainDest_->data_j_index.resize(niDest*njDest); // local position
111
112  domainDest_->domainMask.resize(niDest*njDest);
113
114  if (!domainSrc_->lonvalue_1d.isEmpty())
115  {
116    if (domainDest_->type == CDomain::type_attr::rectilinear)
117    {
118      domainDest_->lonvalue_1d.resize(niDest);
119      domainDest_->latvalue_1d.resize(njDest);
120    }
121    else if (domainDest_->type == CDomain::type_attr::unstructured)
122    {
123      domainDest_->lonvalue_1d.resize(niDest);
124      domainDest_->latvalue_1d.resize(niDest);
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    if (!domainSrc_->bounds_lon_2d.isEmpty())
136    {
137      domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest);
138      domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest);
139    }
140    else if (!domainSrc_->bounds_lon_1d.isEmpty())
141    {
142      domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest);
143      domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest);
144    }
145  }
146  if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest);
147
148  this->transformationMapping_.resize(1);
149  this->transformationWeight_.resize(1);
150  TransformationIndexMap& transMap = this->transformationMapping_[0];
151  TransformationWeightMap& transWeight = this->transformationWeight_[0];
152
153  for (int iDest = 0; iDest < niDest; iDest++)
154  {
155    iSrc = iDest + destIBegin;
156    for (int jDest = 0; jDest < njDest; jDest++)
157    {
158      jSrc = jDest + destJBegin;
159      ind = jSrc * domainSrc_->ni + iSrc;
160      iIdxSrc = domainSrc_->i_index(ind);
161      jIdxSrc = domainSrc_->j_index(ind);
162      indLocDest = jDest*niDest + iDest;
163      indGloDest = (jDest + jbeginDest)*extractNi_ + (iDest + ibeginDest);
164      indLocSrc = (jDest+destJBegin)*domainSrc_->ni + (iDest+destIBegin);
165      indGloSrc = (jIdxSrc )* niGloSrc + iIdxSrc;
166      domainDest_->i_index(indLocDest) = iDest + ibeginDest;                                             // i_index contains global positions
167      domainDest_->j_index(indLocDest) = jDest + jbeginDest;                                             // i_index contains global positions
168      domainDest_->data_i_index(indLocDest) = (domainSrc_->data_dim == 1) ? indLocDest : iDest;          // data_i_index contains local positions
169      domainDest_->data_j_index(indLocDest) = (domainSrc_->data_dim == 1) ? 0 :jDest;                    // data_i_index contains local positions
170      domainDest_->domainMask(indLocDest) = domainSrc_->domainMask(indLocSrc);
171
172      if (domainSrc_->hasArea)
173        domainDest_->area(iDest,jDest) = domainSrc_->area(iSrc,jSrc);
174
175      if (domainSrc_->hasBounds)
176      {
177        if (!domainSrc_->bounds_lon_2d.isEmpty())
178        {
179          for (int n = 0; n < domainSrc_->nvertex; ++n)
180          {
181            domainDest_->bounds_lon_2d(n,iDest,jDest) = domainSrc_->bounds_lon_2d(n,iSrc,jSrc);
182            domainDest_->bounds_lat_2d(n,iDest,jDest) = domainSrc_->bounds_lat_2d(n,iSrc,jSrc);
183          }
184        }
185        else if (!domainSrc_->bounds_lon_1d.isEmpty())
186        {
187          for (int n = 0; n < domainSrc_->nvertex; ++n)
188          {
189            domainDest_->bounds_lon_1d(n,iDest) = domainSrc_->bounds_lon_1d(n,iSrc);
190            domainDest_->bounds_lat_1d(n,iDest) = domainSrc_->bounds_lat_1d(n,iSrc);
191          }
192        }
193      }
194
195      if (domainSrc_->hasLonLat)
196      {
197        if (domainDest_->type == CDomain::type_attr::rectilinear)
198        {
199          domainDest_->latvalue_1d(jDest) = domainSrc_->latvalue_1d(jSrc);
200        }
201        else if (domainDest_->type == CDomain::type_attr::curvilinear)
202        {
203          domainDest_->lonvalue_2d(iDest,jDest) = domainSrc_->lonvalue_2d(iSrc,jSrc);
204          domainDest_->latvalue_2d(iDest,jDest) = domainSrc_->latvalue_2d(iSrc,jSrc);
205        }
206      }
207
208      transMap[indGloDest].push_back(indGloSrc);
209      transWeight[indGloDest].push_back(1.0);
210
211    }
212    if (domainSrc_->hasLonLat)
213    {
214      if (domainDest_->type == CDomain::type_attr::unstructured)
215      {
216        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
217        domainDest_->latvalue_1d(iDest) = domainSrc_->latvalue_1d(iSrc);
218      }
219      else if (domainDest_->type == CDomain::type_attr::rectilinear)
220      {
221        domainDest_->lonvalue_1d(iDest) = domainSrc_->lonvalue_1d(iSrc);
222      }
223    }
224  }
225
226}
227CATCH
228
229/*!
230  Compute the index mapping between domain on grid source and one on grid destination
231*/
232void CDomainAlgorithmExtract::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs)
233{
234}
235
236
237}
Note: See TracBrowser for help on using the repository browser.