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

Last change on this file since 2228 was 2228, checked in by jderouillat, 4 months ago

Fix in exclude masked points from transformation definition (domain_extract)

  • Property svn:eol-style set to native
  • Property svn:executable set to *
File size: 12.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   // Reset geometrical attributes to avoid incompatible (user/domainSource) attributs
44   //   attributs will be defined using domainSource and/or transformation attributs
45   domainDestination->type.reset();
46   domainDestination->ni_glo.reset();
47   domainDestination->nj_glo.reset();
48   
49   domainDestination->i_index.reset();   // defined using domainSource->getLocalElement()
50   domainDestination->j_index.reset();   // "
51   domainDestination->ibegin.reset();    // will be computed in domainDestination->checkDomain() (from checkAttributes())
52   domainDestination->ni.reset();        // "
53   domainDestination->jbegin.reset();    // "
54   domainDestination->nj.reset();        // "
55   
56   domainDestination->mask_1d.reset();   // defined scanning domainSource->getFullView() & domainSource->getWorkflowView() differencies
57   domainDestination->mask_2d.reset();   //   in all case domainDestination->mask_1d used as reference 
58
59   // domainDestination->data_* attributes will be computed in :
60   domainDestination->data_dim.reset();     // domainDestination->checkDomainData() (from checkAttributes())
61   domainDestination->data_ni.reset();
62   domainDestination->data_nj.reset();
63   domainDestination->data_ibegin.reset();
64   domainDestination->data_jbegin.reset();
65   domainDestination->data_i_index.reset(); // domainDestination->checkCompression() (from checkAttributes())
66   domainDestination->data_j_index.reset();
67
68   // Next attributes will be set using domainSource->attributes
69   domainDestination->lonvalue_1d.reset();
70   domainDestination->latvalue_1d.reset();
71   domainDestination->lonvalue_2d.reset();
72   domainDestination->latvalue_2d.reset();
73   domainDestination->nvertex.reset();
74   domainDestination->bounds_lon_1d.reset();
75   domainDestination->bounds_lat_1d.reset();
76   domainDestination->bounds_lon_2d.reset();
77   domainDestination->bounds_lat_2d.reset();
78   domainDestination->area.reset();
79   domainDestination->radius.reset();
80   
81 
82  extractDomain->checkValid(domainSource);
83  extractIBegin_ = extractDomain->ibegin.getValue();
84  extractJBegin_ = extractDomain->jbegin.getValue();
85
86  extractNi_  = extractDomain->ni.getValue();
87  extractNj_  = extractDomain->nj.getValue();
88
89  extractIEnd_ = extractIBegin_ + extractNi_ - 1;
90  extractJEnd_ = extractJBegin_ + extractNj_ - 1;
91
92  if (extractNi_ > domainSource->ni_glo.getValue())
93  {
94    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
95           << "Extract size is greater than size of domain source"
96           << "Size ni_glo of domain source " <<domainSource->getId() << " is " << domainSource->ni_glo.getValue()  << std::endl
97           << "Extract size is " << extractNi_ );
98  }
99
100  if (extractNj_ > domainSource->nj_glo.getValue())
101  {
102    ERROR("CDomainAlgorithmExtract::CDomainAlgorithmExtract(CDomain* domainDestination, CDomain* domainSource, CExtractDomain* extractDomain)",
103           << "Extract size is greater than size of domain source"
104           << "Size nj_glo of domain source " <<domainSource->getId() << " is " << domainSource->nj_glo.getValue()  << std::endl
105           << "Extract size is " << extractNj_ );
106  }
107
108  // Calculate the size of local domain
109  int ind, indLocSrc, indLocDest, iIdxSrc, jIdxSrc, destIBegin = -1, destJBegin = -1, niDest = 0, njDest = 0 ;
110  int indGloDest, indGloSrc, niGloSrc = domainSrc_->ni_glo, iSrc, jSrc;
111  for (int j = 0; j < domainSrc_->nj.getValue(); j++)
112  {
113    for (int i = 0; i < domainSrc_->ni.getValue(); i++)
114    {
115      ind = j*domainSrc_->ni + i;
116      iIdxSrc = domainSrc_->i_index(ind);
117      if ((iIdxSrc >= extractIBegin_) && (iIdxSrc <= extractIEnd_))
118      {
119        jIdxSrc = domainSrc_->j_index(ind);
120        if ((jIdxSrc >= extractJBegin_) && (jIdxSrc <= extractJEnd_))
121        {
122          if ((niDest == 0) && (njDest == 0))
123          {
124            destIBegin = i;
125            destJBegin = j;
126          }
127          if (i == destIBegin) ++njDest;
128        }
129        if (j == destJBegin) ++niDest;
130
131      }
132    }
133  }
134
135   
136  // Set attributes for this transformation
137  domainDest_->type = domainSrc_ -> type ;
138  domainDest_->ni_glo.setValue(extractNi_);
139  domainDest_->nj_glo.setValue(extractNj_);
140  domainDest_->i_index.resize(niDest*njDest);
141  domainDest_->j_index.resize(niDest*njDest);
142
143  // Resize lon/lat, bounds, area arrays to local domain dimensions
144  if (!domainSrc_->lonvalue_1d.isEmpty())
145  {
146    if (domainDest_->type == CDomain::type_attr::rectilinear)
147    {
148      domainDest_->lonvalue_1d.resize(niDest);
149      domainDest_->latvalue_1d.resize(njDest);
150    }
151    else if (domainDest_->type == CDomain::type_attr::unstructured)
152    {
153      domainDest_->lonvalue_1d.resize(niDest);
154      domainDest_->latvalue_1d.resize(niDest);
155    }
156  }
157  else if (!domainSrc_->lonvalue_2d.isEmpty())
158  {
159    domainDest_->lonvalue_2d.resize(niDest,njDest);
160    domainDest_->latvalue_2d.resize(niDest,njDest);
161  }
162  if (domainSrc_->hasBounds)
163  {
164    if (!domainSrc_->bounds_lon_2d.isEmpty())
165    {
166      domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest);
167      domainDest_->bounds_lon_2d.resize(domainDest_->nvertex, niDest, njDest);
168    }
169    else if (!domainSrc_->bounds_lon_1d.isEmpty())
170    {
171      domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest);
172      domainDest_->bounds_lon_1d.resize(domainDest_->nvertex, niDest);
173    }
174  }
175  if (domainSrc_->hasArea) domainDest_->area.resize(niDest,njDest);
176
177
178  // Set attributes required to define domainDestination->localElement_ and associated views, full and workflow)
179  CArray<size_t,1> sourceGlobalIdx = domainSource->getLocalElement()->getGlobalIndex();
180  int indexSize = sourceGlobalIdx.numElements();
181  domainDest_->data_i_index.resize(niDest*njDest);
182  domainDestination->data_i_index = -1; 
183  domainDest_->data_j_index.resize(niDest*njDest);
184  domainDestination->data_j_index = 0; 
185
186  CArray<int,1> sourceWorkflowIdx = domainSource->getLocalView(CElementView::WORKFLOW)->getIndex();
187  int srcWorkflowSize = sourceWorkflowIdx.numElements();
188  for (size_t i = 0; i < srcWorkflowSize ; ++i)
189  {
190    {
191      int iIdxSrc = sourceWorkflowIdx(i)%domainSource->ni-destIBegin;
192      int jIdxSrc = sourceWorkflowIdx(i)/domainSource->ni-destJBegin;
193      int extractedSrcWFIdx = jIdxSrc * niDest + iIdxSrc;
194      if ((extractedSrcWFIdx>=0)&&(extractedSrcWFIdx<niDest*njDest))
195      {
196        domainDest_->data_i_index(extractedSrcWFIdx) = extractedSrcWFIdx;
197      }
198    }
199  }
200
201  int iIdxSrcMin = INT_MAX;
202  int jIdxSrcMin = INT_MAX;
203  for (int countSrc = 0; countSrc < indexSize ; ++countSrc)
204  {
205      if ( sourceGlobalIdx(countSrc)%domainSource->ni_glo < iIdxSrcMin )
206          iIdxSrcMin = sourceGlobalIdx(countSrc)%domainSource->ni_glo;
207      if ( sourceGlobalIdx(countSrc)/domainSource->ni_glo < jIdxSrcMin )
208          jIdxSrcMin = sourceGlobalIdx(countSrc)/domainSource->ni_glo;
209  }
210 
211  int countDest(0); // increment of the position in destination domain
212  for (int countSrc = 0; countSrc < indexSize ; ++countSrc)
213  {
214    int iIdxSrc = sourceGlobalIdx(countSrc)%domainSource->ni_glo;
215    int jIdxSrc = sourceGlobalIdx(countSrc)/domainSource->ni_glo;
216    // check that point countSrc concerned by extract
217    if ( (iIdxSrc >= extractIBegin_) && (iIdxSrc <= extractIEnd_)
218         && (jIdxSrc >= extractJBegin_) && (jIdxSrc <= extractJEnd_) )
219    {
220      // if concerned, convert source the global indexation in the extracted frame
221      domainDest_->i_index(countDest) = iIdxSrc-extractIBegin_;
222      domainDest_->j_index(countDest) = jIdxSrc-extractJBegin_;
223
224      // ------------------ define transformation only if in the WF ------------------
225      int convert_locally_global_idx = (jIdxSrc-jIdxSrcMin)*domainSource->ni + (iIdxSrc-iIdxSrcMin) ;
226      bool concerned_by_WF(false);
227      for ( int i = 0 ; i<sourceWorkflowIdx.numElements() ; ++i )
228      {
229        if (sourceWorkflowIdx(i)==convert_locally_global_idx)
230        {     
231          concerned_by_WF = true;
232          break;
233        }
234      }
235      if (concerned_by_WF)
236      {
237        transformationMapping_[extractNi_*(jIdxSrc-extractJBegin_)+iIdxSrc-extractIBegin_]=sourceGlobalIdx(countSrc);
238      }
239      // -----------------------------------------------------------------------------
240
241      int iIdxDestLocal = countDest%niDest;
242      int jIdxDestLocal = countDest/niDest;
243      int iIdxSrcLocal  = countSrc%domainSource->ni;
244      int jIdxSrcLocal  = countSrc/domainSource->ni;
245
246      // area
247      if (!domainSrc_->area.isEmpty())
248      {
249        domainDest_->area(iIdxDestLocal,jIdxDestLocal) = domainSrc_->area(iIdxSrcLocal,jIdxSrcLocal);
250      }
251
252      // bounds
253      if (!domainDest_->bounds_lon_1d.isEmpty())
254      {
255        for (int n = 0; n < domainSrc_->nvertex; ++n)
256        {
257          domainDest_->bounds_lon_1d(n, countDest) = domainSrc_->bounds_lon_1d(n,countSrc);
258          domainDest_->bounds_lat_1d(n, countDest) = domainSrc_->bounds_lat_1d(n,countSrc);
259        }
260      }
261      else if (!domainDest_->bounds_lon_2d.isEmpty())
262      {
263        for (int n = 0; n < domainSrc_->nvertex; ++n)
264        {
265          domainDest_->bounds_lon_2d(n, iIdxDestLocal, jIdxDestLocal) = domainSrc_->bounds_lon_2d(n, iIdxSrcLocal, jIdxSrcLocal);
266          domainDest_->bounds_lat_2d(n, iIdxDestLocal, jIdxDestLocal) = domainSrc_->bounds_lat_2d(n, iIdxSrcLocal, jIdxSrcLocal);
267        }
268      }
269
270      // lon/lat
271      if (!domainDest_->lonvalue_1d.isEmpty())
272      {
273        if (domainDest_->type == CDomain::type_attr::rectilinear)
274        {
275          // i : scan nbr of points in src
276          domainDest_->lonvalue_1d(iIdxDestLocal) = domainSrc_->lonvalue_1d(iIdxSrcLocal);
277          domainDest_->latvalue_1d(jIdxDestLocal) = domainSrc_->latvalue_1d(jIdxSrcLocal);
278        }
279        else if (domainDest_->type == CDomain::type_attr::unstructured)
280        {
281          domainDest_->lonvalue_1d(countDest) = domainSrc_->lonvalue_1d(countSrc);
282          domainDest_->latvalue_1d(countDest) = domainSrc_->latvalue_1d(countSrc);
283        }
284      }
285      else if (!domainDest_->lonvalue_2d.isEmpty())
286      {
287        if (domainDest_->type == CDomain::type_attr::curvilinear)
288        {
289          domainDest_->lonvalue_2d(iIdxDestLocal, jIdxDestLocal) = domainSrc_->lonvalue_2d(iIdxSrcLocal,jIdxSrcLocal);
290          domainDest_->latvalue_2d(iIdxDestLocal, jIdxDestLocal) = domainSrc_->latvalue_2d(iIdxSrcLocal,jIdxSrcLocal);
291        }
292      }
293     
294      // if point i has been identified as extracted, increment position in destination domain for the next point
295      countDest++;
296    }
297
298  }
299 
300  domainDestination->checkAttributes() ;
301  this->computeAlgorithm(domainSource->getLocalView(CElementView::WORKFLOW), domainDestination->getLocalView(CElementView::WORKFLOW)) ;
302}
303CATCH
304
305
306}
Note: See TracBrowser for help on using the repository browser.