Salome HOME
Debug on limit cases when input part are empty. Correct bug of DeduceNumberOfEntities.
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingIMesh.cxx
1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19 // Author : Anthony Geay (CEA/DEN)
20
21 #include "MEDCouplingIMesh.hxx"
22 #include "MEDCouplingCMesh.hxx"
23 #include "MEDCouplingMemArray.hxx"
24 #include "MEDCouplingFieldDouble.hxx"
25
26 #include <functional>
27 #include <algorithm>
28 #include <sstream>
29 #include <numeric>
30
31 using namespace ParaMEDMEM;
32
33 MEDCouplingIMesh::MEDCouplingIMesh():_space_dim(-1)
34 {
35   _origin[0]=0.; _origin[1]=0.; _origin[2]=0.;
36   _dxyz[0]=0.; _dxyz[1]=0.; _dxyz[2]=0.;
37   _structure[0]=0; _structure[1]=0; _structure[2]=0;
38 }
39
40 MEDCouplingIMesh::MEDCouplingIMesh(const MEDCouplingIMesh& other, bool deepCopy):MEDCouplingStructuredMesh(other,deepCopy),_space_dim(other._space_dim),_axis_unit(other._axis_unit)
41 {
42   _origin[0]=other._origin[0]; _origin[1]=other._origin[1]; _origin[2]=other._origin[2];
43   _dxyz[0]=other._dxyz[0]; _dxyz[1]=other._dxyz[1]; _dxyz[2]=other._dxyz[2];
44   _structure[0]=other._structure[0]; _structure[1]=other._structure[1]; _structure[2]=other._structure[2];
45 }
46
47 MEDCouplingIMesh::~MEDCouplingIMesh()
48 {
49 }
50
51 MEDCouplingIMesh *MEDCouplingIMesh::New()
52 {
53   return new MEDCouplingIMesh;
54 }
55
56 MEDCouplingIMesh *MEDCouplingIMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
57                                         const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
58 {
59   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> ret(new MEDCouplingIMesh);
60   ret->setName(meshName);
61   ret->setSpaceDimension(spaceDim);
62   ret->setNodeStruct(nodeStrctStart,nodeStrctStop);
63   ret->setOrigin(originStart,originStop);
64   ret->setDXYZ(dxyzStart,dxyzStop);
65   return ret.retn();
66 }
67
68 MEDCouplingMesh *MEDCouplingIMesh::deepCpy() const
69 {
70   return clone(true);
71 }
72
73 MEDCouplingIMesh *MEDCouplingIMesh::clone(bool recDeepCpy) const
74 {
75   return new MEDCouplingIMesh(*this,recDeepCpy);
76 }
77
78 /*!
79  * This method creates a copy of \a this enlarged by \a ghostLev cells on each axis.
80  * If \a ghostLev equal to 0 this method behaves as MEDCouplingIMesh::clone.
81  *
82  * \param [in] ghostLev - the ghost level expected
83  * \return MEDCouplingIMesh * - a newly alloacted object to be managed by the caller.
84  * \throw if \a ghostLev < 0.
85  */
86 MEDCouplingIMesh *MEDCouplingIMesh::buildWithGhost(int ghostLev) const
87 {
88   if(ghostLev<0)
89     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::buildWithGhost : the ghostLev must be >= 0 !");
90   checkCoherency();
91   int spaceDim(getSpaceDimension());
92   double origin[3],dxyz[3];
93   int structure[3];
94   for(int i=0;i<spaceDim;i++)
95     {
96       origin[i]=_origin[i]-ghostLev*_dxyz[i];
97       dxyz[i]=_dxyz[i];
98       structure[i]=_structure[i]+2*ghostLev;
99     }
100   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> ret(MEDCouplingIMesh::New(getName(),spaceDim,structure,structure+spaceDim,origin,origin+spaceDim,dxyz,dxyz+spaceDim));
101   ret->copyTinyInfoFrom(this);
102   return ret.retn();
103 }
104
105 void MEDCouplingIMesh::setNodeStruct(const int *nodeStrctStart, const int *nodeStrctStop)
106 {
107   checkSpaceDimension();
108   int sz((int)std::distance(nodeStrctStart,nodeStrctStop));
109   if(sz!=_space_dim)
110     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::setNodeStruct : input vector of node structure has not the right size ! Or change space dimension before calling it !");
111   std::copy(nodeStrctStart,nodeStrctStop,_structure);
112   declareAsNew();
113 }
114
115 std::vector<int> MEDCouplingIMesh::getNodeStruct() const
116 {
117   checkSpaceDimension();
118   return std::vector<int>(_structure,_structure+_space_dim);
119 }
120
121 void MEDCouplingIMesh::setOrigin(const double *originStart, const double *originStop)
122 {
123   checkSpaceDimension();
124   int sz((int)std::distance(originStart,originStop));
125   if(sz!=_space_dim)
126     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::setOrigin : input vector of origin vector has not the right size ! Or change space dimension before calling it !");
127   std::copy(originStart,originStop,_origin);
128   declareAsNew();
129 }
130
131 std::vector<double> MEDCouplingIMesh::getOrigin() const
132 {
133   checkSpaceDimension();
134   return std::vector<double>(_origin,_origin+_space_dim);
135 }
136
137 void MEDCouplingIMesh::setDXYZ(const double *dxyzStart, const double *dxyzStop)
138 {
139   checkSpaceDimension();
140   int sz((int)std::distance(dxyzStart,dxyzStop));
141   if(sz!=_space_dim)
142     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::setDXYZ : input vector of dxyz vector has not the right size ! Or change space dimension before calling it !");
143   std::copy(dxyzStart,dxyzStop,_dxyz);
144   declareAsNew();
145 }
146
147 std::vector<double> MEDCouplingIMesh::getDXYZ() const
148 {
149   checkSpaceDimension();
150   return std::vector<double>(_dxyz,_dxyz+_space_dim);
151 }
152
153 void MEDCouplingIMesh::setAxisUnit(const std::string& unitName)
154 {
155   _axis_unit=unitName;
156   declareAsNew();
157 }
158
159 std::string MEDCouplingIMesh::getAxisUnit() const
160 {
161   return _axis_unit;
162 }
163
164 /*!
165  * This method returns the measure of any cell in \a this.
166  * This specific method of image grid mesh utilizes the fact that any cell in \a this have the same measure.
167  * The value returned by this method is those used to feed the returned field in the MEDCouplingIMesh::getMeasureField.
168  *
169  * \sa getMeasureField
170  */
171 double MEDCouplingIMesh::getMeasureOfAnyCell() const
172 {
173   checkCoherency();
174   int dim(getSpaceDimension());
175   double ret(1.);
176   for(int i=0;i<dim;i++)
177     ret*=fabs(_dxyz[i]);
178   return ret;
179 }
180
181 /*!
182  * This method is allows to convert \a this into MEDCouplingCMesh instance.
183  * This method is the middle level between MEDCouplingIMesh and the most general MEDCouplingUMesh.
184  * This method is useful for MED writers that do not have still the image grid support.
185  *
186  * \sa MEDCouplingMesh::buildUnstructured
187  */
188 MEDCouplingCMesh *MEDCouplingIMesh::convertToCartesian() const
189 {
190   checkCoherency();
191   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCMesh> ret(MEDCouplingCMesh::New());
192   try
193   { ret->copyTinyInfoFrom(this); }
194   catch(INTERP_KERNEL::Exception& ) { }
195   int spaceDim(getSpaceDimension());
196   std::vector<std::string> infos(buildInfoOnComponents());
197   for(int i=0;i<spaceDim;i++)
198     {
199       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr(DataArrayDouble::New()); arr->alloc(_structure[i],1); arr->setInfoOnComponent(0,infos[i]);
200       arr->iota(); arr->applyLin(_dxyz[i],_origin[i]);
201       ret->setCoordsAt(i,arr);
202     }
203   return ret.retn();
204 }
205
206 /*!
207  * This method refines \a this uniformaly along all of its dimensions. In case of success the space covered by \a this will remain
208  * the same before the invocation except that the number of cells will be multiplied by \a factor ^ this->getMeshDimension().
209  * The origin of \a this will be not touched only spacing and node structure will be changed.
210  * This method can be useful for AMR users.
211  */
212 void MEDCouplingIMesh::refineWithFactor(const std::vector<int>& factors)
213 {
214   if((int)factors.size()!=_space_dim)
215     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::refineWithFactor : refinement factors must have size equal to spaceDim !");
216   checkCoherency();
217   std::vector<int> structure(_structure,_structure+3);
218   std::vector<double> dxyz(_dxyz,_dxyz+3);
219   for(int i=0;i<_space_dim;i++)
220     {
221       if(factors[i]<=0)
222         {
223           std::ostringstream oss; oss << "MEDCouplingIMesh::refineWithFactor : factor for axis #" << i << " (" << factors[i] << ")is invalid ! Must be > 0 !";
224           throw INTERP_KERNEL::Exception(oss.str().c_str());
225         }
226       int factAbs(std::abs(factors[i]));
227       double fact2(1./(double)factors[i]);
228       structure[i]=(_structure[i]-1)*factAbs+1;
229       dxyz[i]=fact2*_dxyz[i];
230     }
231   std::copy(structure.begin(),structure.end(),_structure);
232   std::copy(dxyz.begin(),dxyz.end(),_dxyz);
233   declareAsNew();
234 }
235
236 /*!
237  * This method returns a newly created mesh containing a single cell in it. This returned cell covers exactly the space covered by \a this.
238  *
239  * \return MEDCouplingIMesh * - A newly created object (to be managed by the caller with decrRef) containing simply one cell.
240  *
241  * \throw if \a this does not pass the \c checkCoherency test.
242  */
243 MEDCouplingIMesh *MEDCouplingIMesh::asSingleCell() const
244 {
245   checkCoherency();
246   int spaceDim(getSpaceDimension()),nodeSt[3];
247   double dxyz[3];
248   for(int i=0;i<spaceDim;i++)
249     {
250       if(_structure[i]>=2)
251         {
252           nodeSt[i]=2;
253           dxyz[i]=(_structure[i]-1)*_dxyz[i];
254         }
255       else
256         {
257           nodeSt[i]=_structure[i];
258           dxyz[i]=_dxyz[i];
259         }
260     }
261   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> ret(MEDCouplingIMesh::New(getName(),getSpaceDimension(),nodeSt,nodeSt+spaceDim,_origin,_origin+spaceDim,dxyz,dxyz+spaceDim));
262   ret->copyTinyInfoFrom(this);
263   return ret.retn();
264 }
265
266 /*!
267  * This static method is useful to condense field on cells of a MEDCouplingIMesh instance coming from a refinement ( MEDCouplingIMesh::refineWithFactor for example)
268  * to a coarse MEDCouplingIMesh instance. So this method can be seen as a specialization in P0P0 conservative interpolation non overlaping from fine image mesh
269  * to a coarse image mesh. Only tuples ( deduced from \a fineLocInCoarse ) of \a coarseDA will be modified. Other tuples of \a coarseDA will be let unchanged.
270  *
271  * \param [in] coarseSt The cell structure of coarse mesh.
272  * \param [in] fineDA The DataArray containing the cell field on uniformly refined mesh
273  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
274  * \param [in] facts The refinement coefficient per axis.
275  * \param [in,out] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
276  *
277  * \sa CondenseFineToCoarseGhost,SpreadCoarseToFine
278  */
279 void MEDCouplingIMesh::CondenseFineToCoarse(const std::vector<int>& coarseSt, const DataArrayDouble *fineDA, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts, DataArrayDouble *coarseDA)
280 {
281   if(coarseSt.size()!=fineLocInCoarse.size() || coarseSt.size()!=facts.size())
282     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : All input vectors (dimension) must have the same size !");
283   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
284     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the parameters 1 or 3 are NULL or not allocated !");
285   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseSt)),nbOfTuplesInFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(fineLocInCoarse));
286   int nbCompo(fineDA->getNumberOfComponents());
287   if(coarseDA->getNumberOfComponents()!=nbCompo)
288     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the number of components of fine DA and coarse one mismatches !");
289   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
290     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the size of fineLocInCoarse (4th param) and facts (5th param) must be equal to the sier of coarseSt (2nd param) !");
291   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
292     {
293       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarse : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
294       throw INTERP_KERNEL::Exception(oss.str().c_str());
295     }
296   int nbTuplesFine(fineDA->getNumberOfTuples());
297   if(nbOfTuplesInFineExp==0)
298     return ;
299   if(nbTuplesFine%nbOfTuplesInFineExp!=0)
300     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : Invalid nb of tuples in fine DataArray regarding its structure !");
301   int fact(std::accumulate(facts.begin(),facts.end(),1,std::multiplies<int>()));
302   if(nbTuplesFine!=fact*nbOfTuplesInFineExp)
303     {
304       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarse : Invalid number of tuples ("  << nbTuplesFine << ") of fine dataarray is invalid ! Must be " << fact*nbOfTuplesInFineExp << "!";
305       throw INTERP_KERNEL::Exception(oss.str().c_str());
306     }
307   // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
308   double *outPtr(coarseDA->getPointer());
309   const double *inPtr(fineDA->begin());
310   //
311   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
312   switch(meshDim)
313   {
314     case 1:
315       {
316         int offset(fineLocInCoarse[0].first),fact0(facts[0]);
317         for(int i=0;i<dims[0];i++)
318           {
319             double *loc(outPtr+(offset+i)*nbCompo);
320             for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
321               {
322                 if(ifact!=0)
323                   std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
324                 else
325                   std::copy(inPtr,inPtr+nbCompo,loc);
326               }
327           }
328         break;
329       }
330     case 2:
331       {
332         int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first),fact1(facts[1]),fact0(facts[0]);
333         for(int j=0;j<dims[1];j++)
334           {
335             for(int jfact=0;jfact<fact1;jfact++)
336               {
337                 for(int i=0;i<dims[0];i++)
338                   {
339                     double *loc(outPtr+(kk+i)*nbCompo);
340                     for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
341                       {
342                         if(jfact!=0 || ifact!=0)
343                           std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
344                         else
345                           std::copy(inPtr,inPtr+nbCompo,loc);
346                       }
347                   }
348               }
349             kk+=coarseSt[0];
350           }
351         break;
352       }
353     case 3:
354       {
355         int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first+coarseSt[0]*coarseSt[1]*fineLocInCoarse[2].first),fact2(facts[2]),fact1(facts[1]),fact0(facts[0]);
356         for(int k=0;k<dims[2];k++)
357           {
358             for(int kfact=0;kfact<fact2;kfact++)
359               {
360                 for(int j=0;j<dims[1];j++)
361                   {
362                     for(int jfact=0;jfact<fact1;jfact++)
363                       {
364                         for(int i=0;i<dims[0];i++)
365                           {
366                             double *loc(outPtr+(kk+i+j*coarseSt[0])*nbCompo);
367                             for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
368                               {
369                                 if(kfact!=0 || jfact!=0 || ifact!=0)
370                                   std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
371                                 else
372                                   std::copy(inPtr,inPtr+nbCompo,loc);
373                               }
374                           }
375                       }
376                   }
377               }
378             kk+=coarseSt[0]*coarseSt[1];
379           }
380         break;
381       }
382     default:
383       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : only dimensions 1, 2 and 3 supported !");
384   }
385 }
386
387 /*!
388  * This static method is useful to condense field on cells of a MEDCouplingIMesh instance coming from a refinement ( MEDCouplingIMesh::refineWithFactor for example)
389  * to a coarse MEDCouplingIMesh instance. So this method can be seen as a specialization in P0P0 conservative interpolation non overlaping from fine image mesh
390  * to a coarse image mesh. Only tuples ( deduced from \a fineLocInCoarse ) of \a coarseDA will be modified. Other tuples of \a coarseDA will be let unchanged.
391  *
392  * \param [in] coarseSt The cell structure of coarse mesh.
393  * \param [in] fineDA The DataArray containing the cell field on uniformly refined mesh
394  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
395  * \param [in] facts The refinement coefficient per axis.
396  * \param [in,out] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
397  * \param [in] ghostSize - The size of the ghost zone. The ghost zone is expected to be the same for all axis and both for coarse and fine meshes.
398  *
399  * \sa CondenseFineToCoarse,SpreadCoarseToFineGhost
400  */
401 void MEDCouplingIMesh::CondenseFineToCoarseGhost(const std::vector<int>& coarseSt, const DataArrayDouble *fineDA, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts, DataArrayDouble *coarseDA, int ghostSize)
402 {
403   if(ghostSize<0)
404     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : ghost level has to be >= 0 !");
405   if(coarseSt.size()!=fineLocInCoarse.size() || coarseSt.size()!=facts.size())
406     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : All input vectors (dimension) must have the same size !");
407   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
408     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : the parameters 1 or 3 are NULL or not allocated !");
409   std::vector<int> coarseStG(coarseSt.size()); std::transform(coarseSt.begin(),coarseSt.end(),coarseStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
410   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseStG));
411   int nbCompo(fineDA->getNumberOfComponents());
412   if(coarseDA->getNumberOfComponents()!=nbCompo)
413     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : the number of components of fine DA and coarse one mismatches !");
414   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
415     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : the size of fineLocInCoarse (4th param) and facts (5th param) must be equal to the sier of coarseSt (2nd param) !");
416   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
417     {
418       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarseGhost : Expecting " << nbOfTuplesInCoarseExp << " tuples in coarse DataArray having " << coarseDA->getNumberOfTuples() << " !";
419       throw INTERP_KERNEL::Exception(oss.str().c_str());
420     }
421   //
422   std::vector<int> fineStG(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
423   std::transform(fineStG.begin(),fineStG.end(),facts.begin(),fineStG.begin(),std::multiplies<int>());
424   std::transform(fineStG.begin(),fineStG.end(),fineStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
425   int nbTuplesFine(fineDA->getNumberOfTuples()),nbTuplesFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(fineStG));
426   if(fineDA->getNumberOfTuples()!=nbTuplesFineExp)
427     {
428       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarseGhost : Expecting " << nbTuplesFineExp << " tuples in fine DataArray having " << nbTuplesFine << " !";
429       throw INTERP_KERNEL::Exception(oss.str().c_str());
430     }
431   //
432   double *outPtr(coarseDA->getPointer());
433   const double *inPtr(fineDA->begin());
434   //
435   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
436   switch(meshDim)
437   {
438     case 1:
439       {
440         int offset(fineLocInCoarse[0].first+ghostSize),fact0(facts[0]);
441         inPtr+=ghostSize*nbCompo;
442         for(int i=0;i<dims[0];i++)
443           {
444             double *loc(outPtr+(offset+i)*nbCompo);
445             for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
446               {
447                 if(ifact!=0)
448                   std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
449                 else
450                   std::copy(inPtr,inPtr+nbCompo,loc);
451               }
452           }
453         break;
454       }
455     case 2:
456       {
457         int nxwg(coarseSt[0]+2*ghostSize);
458         int kk(fineLocInCoarse[0].first+ghostSize+nxwg*(fineLocInCoarse[1].first+ghostSize)),fact1(facts[1]),fact0(facts[0]);
459         inPtr+=(dims[0]*fact0+2*ghostSize)*nbCompo;
460         for(int j=0;j<dims[1];j++)
461           {
462              for(int jfact=0;jfact<fact1;jfact++)
463               {
464                 inPtr+=ghostSize*nbCompo;
465                 for(int i=0;i<dims[0];i++)
466                   {
467                     double *loc(outPtr+(kk+i)*nbCompo);
468                     for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
469                       {
470                         if(jfact!=0 || ifact!=0)
471                           std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
472                         else
473                           std::copy(inPtr,inPtr+nbCompo,loc);
474                       }
475                   }
476                 inPtr+=ghostSize*nbCompo;
477               }
478             kk+=nxwg;
479           }
480         break;
481       }
482     default:
483       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : only dimensions 1, 2 supported !");
484   }
485 }
486
487 /*!
488  * This method spreads the values of coarse data \a coarseDA into \a fineDA.
489  *
490  * \param [in] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
491  * \param [in] coarseSt The cell structure of coarse mesh.
492  * \param [in,out] fineDA The DataArray containing the cell field on uniformly refined mesh
493  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
494  * \param [in] facts The refinement coefficient per axis.
495  * \sa SpreadCoarseToFineGhost, CondenseFineToCoarse
496  */
497 void MEDCouplingIMesh::SpreadCoarseToFine(const DataArrayDouble *coarseDA, const std::vector<int>& coarseSt, DataArrayDouble *fineDA, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts)
498 {
499   if(coarseSt.size()!=fineLocInCoarse.size() || coarseSt.size()!=facts.size())
500       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : All input vectors (dimension) must have the same size !");
501   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
502     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the parameters 1 or 3 are NULL or not allocated !");
503   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseSt)),nbOfTuplesInFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(fineLocInCoarse));
504   int nbCompo(fineDA->getNumberOfComponents());
505   if(coarseDA->getNumberOfComponents()!=nbCompo)
506     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the number of components of fine DA and coarse one mismatches !");
507   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
508     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the size of fineLocInCoarse (4th param) and facts (5th param) must be equal to the sier of coarseSt (2nd param) !");
509   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
510     {
511       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFine : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
512       throw INTERP_KERNEL::Exception(oss.str().c_str());
513     }
514   int nbTuplesFine(fineDA->getNumberOfTuples());
515   if(nbTuplesFine%nbOfTuplesInFineExp!=0)
516     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : Invalid nb of tuples in fine DataArray regarding its structure !");
517   int fact(std::accumulate(facts.begin(),facts.end(),1,std::multiplies<int>()));
518   if(nbTuplesFine!=fact*nbOfTuplesInFineExp)
519     {
520       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFine : Invalid number of tuples ("  << nbTuplesFine << ") of fine dataarray is invalid ! Must be " << fact*nbOfTuplesInFineExp << "!";
521       throw INTERP_KERNEL::Exception(oss.str().c_str());
522     }
523   // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
524   double *outPtr(fineDA->getPointer());
525   const double *inPtr(coarseDA->begin());
526   //
527   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
528   switch(meshDim)
529   {
530     case 1:
531       {
532         int offset(fineLocInCoarse[0].first),fact0(facts[0]);
533         for(int i=0;i<dims[0];i++)
534           {
535             const double *loc(inPtr+(offset+i)*nbCompo);
536             for(int ifact=0;ifact<fact0;ifact++)
537               outPtr=std::copy(loc,loc+nbCompo,outPtr);
538           }
539         break;
540       }
541     case 2:
542       {
543         int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first),fact0(facts[0]),fact1(facts[1]);
544         for(int j=0;j<dims[1];j++)
545           {
546             for(int jfact=0;jfact<fact1;jfact++)
547               {
548                 for(int i=0;i<dims[0];i++)
549                   {
550                     const double *loc(inPtr+(kk+i)*nbCompo);
551                     for(int ifact=0;ifact<fact0;ifact++)
552                       outPtr=std::copy(loc,loc+nbCompo,outPtr);
553                   }
554               }
555             kk+=coarseSt[0];
556           }
557         break;
558       }
559     case 3:
560       {
561         int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first+coarseSt[0]*coarseSt[1]*fineLocInCoarse[2].first),fact0(facts[0]),fact1(facts[2]),fact2(facts[2]);
562         for(int k=0;k<dims[2];k++)
563           {
564             for(int kfact=0;kfact<fact2;kfact++)
565               {
566                 for(int j=0;j<dims[1];j++)
567                   {
568                     for(int jfact=0;jfact<fact1;jfact++)
569                       {
570                         for(int i=0;i<dims[0];i++)
571                           {
572                             const double *loc(inPtr+(kk+i+j*coarseSt[0])*nbCompo);
573                             for(int ifact=0;ifact<fact0;ifact++)
574                               outPtr=std::copy(loc,loc+nbCompo,outPtr);
575                           }
576                       }
577                   }
578               }
579             kk+=coarseSt[0]*coarseSt[1];
580           }
581         break;
582       }
583     default:
584       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : only dimensions 1, 2 and 3 supported !");
585   }
586 }
587
588 /*!
589  * This method spreads the values of coarse data \a coarseDA into \a fineDA.
590  *
591  * \param [in] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
592  * \param [in] coarseSt The cell structure of coarse mesh.
593  * \param [in,out] fineDA The DataArray containing the cell field on uniformly refined mesh
594  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
595  * \param [in] facts The refinement coefficient per axis.
596  * \param [in] ghostSize - The size of the ghost zone. The ghost zone is expected to be the same for all axis and both for coarse and fine meshes.
597  *
598  * \sa CondenseFineToCoarse, SpreadCoarseToFineGhostZone
599  */
600 void MEDCouplingIMesh::SpreadCoarseToFineGhost(const DataArrayDouble *coarseDA, const std::vector<int>& coarseSt, DataArrayDouble *fineDA, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts, int ghostSize)
601 {
602   if(ghostSize<0)
603     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : ghost level has to be >= 0 !");
604   if(coarseSt.size()!=fineLocInCoarse.size() || coarseSt.size()!=facts.size())
605     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : All input vectors (dimension) must have the same size !");
606   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
607     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : the parameters 1 or 3 are NULL or not allocated !");
608   std::vector<int> coarseStG(coarseSt.size()); std::transform(coarseSt.begin(),coarseSt.end(),coarseStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
609   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseStG));
610   int nbCompo(fineDA->getNumberOfComponents());
611   if(coarseDA->getNumberOfComponents()!=nbCompo)
612     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : the number of components of fine DA and coarse one mismatches !");
613   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
614     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : the size of fineLocInCoarse (4th param) and facts (5th param) must be equal to the sier of coarseSt (2nd param) !");
615   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
616     {
617       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFineGhost : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
618       throw INTERP_KERNEL::Exception(oss.str().c_str());
619     }
620   //
621   std::vector<int> fineStG(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
622   std::transform(fineStG.begin(),fineStG.end(),facts.begin(),fineStG.begin(),std::multiplies<int>());
623   std::transform(fineStG.begin(),fineStG.end(),fineStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
624   int nbTuplesFine(fineDA->getNumberOfTuples()),nbTuplesFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(fineStG));
625   if(fineDA->getNumberOfTuples()!=nbTuplesFineExp)
626     {
627       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFineGhost : Expecting " << nbTuplesFineExp << " tuples in fine DataArray having " << nbTuplesFine << " !";
628       throw INTERP_KERNEL::Exception(oss.str().c_str());
629     }
630   //
631   double *outPtr(fineDA->getPointer());
632   const double *inPtr(coarseDA->begin());
633   //
634   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
635   switch(meshDim)
636   {
637     case 1:
638       {
639         int offset(fineLocInCoarse[0].first+ghostSize-1),fact0(facts[0]);//offset is always >=0 thanks to the fact that ghostSize>=1 !
640         for(int i=0;i<ghostSize;i++)
641           outPtr=std::copy(inPtr+offset*nbCompo,inPtr+(offset+1)*nbCompo,outPtr);
642         offset=fineLocInCoarse[0].first+ghostSize;
643         for(int i=0;i<dims[0];i++)
644           {
645             const double *loc(inPtr+(offset+i)*nbCompo);
646             for(int ifact=0;ifact<fact0;ifact++)
647               outPtr=std::copy(loc,loc+nbCompo,outPtr);
648           }
649         offset=fineLocInCoarse[0].second+ghostSize;
650         for(int i=0;i<ghostSize;i++)
651           outPtr=std::copy(inPtr+offset*nbCompo,inPtr+(offset+1)*nbCompo,outPtr);
652         break;
653       }
654     case 2:
655       {
656         int nxwg(coarseSt[0]+2*ghostSize),fact0(facts[0]),fact1(facts[1]);
657         int kk(fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].first+ghostSize-1));//kk is always >=0 thanks to the fact that ghostSize>=1 !
658         for(int jg=0;jg<ghostSize;jg++)
659           {
660             for(int ig=0;ig<ghostSize;ig++)
661               outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
662             int kk0(kk+1);
663             for(int ig=0;ig<dims[0];ig++,kk0++)
664               for(int ifact=0;ifact<fact0;ifact++)
665                 outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
666             for(int ik=0;ik<ghostSize;ik++)
667               outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
668           }
669         for(int j=0;j<dims[1];j++)
670           {
671             kk=fineLocInCoarse[0].first-1+ghostSize+nxwg*(fineLocInCoarse[1].first+ghostSize+j);
672             for(int jfact=0;jfact<fact1;jfact++)
673               {
674                 for(int ig=0;ig<ghostSize;ig++)
675                   outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
676                 int kk0(kk+1);//1 not ghost. We make the hypothesis that factors is >= ghostlev
677                 for(int i=0;i<dims[0];i++,kk0++)
678                   {
679                     const double *loc(inPtr+kk0*nbCompo);
680                     for(int ifact=0;ifact<fact0;ifact++)
681                       outPtr=std::copy(loc,loc+nbCompo,outPtr);
682                   }
683                 for(int ig=0;ig<ghostSize;ig++)
684                   outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
685               }
686           }
687         kk=fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].second+ghostSize);
688         for(int jg=0;jg<ghostSize;jg++)
689           {
690             for(int ig=0;ig<ghostSize;ig++)
691               outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
692             int kk0(kk+1);
693             for(int ig=0;ig<dims[0];ig++,kk0++)
694               for(int ifact=0;ifact<fact0;ifact++)
695                 outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
696             for(int ik=0;ik<ghostSize;ik++)
697               outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
698           }
699         break;
700       }
701     default:
702       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : only dimensions 1, 2 supported !");
703   }
704 }
705
706 /*!
707  * This method spreads the values of coarse data \a coarseDA into \a fineDA \b ONLY \b in \b the \b ghost \b zone (contrary to SpreadCoarseToFineGhost that spread the values everywhere).
708  *
709  * \param [in] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
710  * \param [in] coarseSt The cell structure of coarse mesh.
711  * \param [in,out] fineDA The DataArray containing the cell field on uniformly refined mesh
712  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
713  * \param [in] facts The refinement coefficient per axis.
714  * \param [in] ghostSize - The size of the ghost zone. The ghost zone is expected to be the same for all axis and both for coarse and fine meshes.
715  *
716  * \sa SpreadCoarseToFineGhost
717  */
718 void MEDCouplingIMesh::SpreadCoarseToFineGhostZone(const DataArrayDouble *coarseDA, const std::vector<int>& coarseSt, DataArrayDouble *fineDA, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts, int ghostSize)
719 {
720   if(ghostSize<0)
721     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : ghost level has to be >= 0 !");
722   if(coarseSt.size()!=fineLocInCoarse.size() || coarseSt.size()!=facts.size())
723     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : All input vectors (dimension) must have the same size !");
724   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
725     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : the parameters 1 or 3 are NULL or not allocated !");
726   std::vector<int> coarseStG(coarseSt.size()); std::transform(coarseSt.begin(),coarseSt.end(),coarseStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
727   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseStG));
728   int nbCompo(fineDA->getNumberOfComponents());
729   if(coarseDA->getNumberOfComponents()!=nbCompo)
730     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : the number of components of fine DA and coarse one mismatches !");
731   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
732     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : the size of fineLocInCoarse (4th param) and facts (5th param) must be equal to the sier of coarseSt (2nd param) !");
733   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
734     {
735       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFineGhostZone : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
736       throw INTERP_KERNEL::Exception(oss.str().c_str());
737     }
738   //
739   std::vector<int> fineStG(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
740   std::transform(fineStG.begin(),fineStG.end(),facts.begin(),fineStG.begin(),std::multiplies<int>());
741   std::transform(fineStG.begin(),fineStG.end(),fineStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
742   int nbTuplesFine(fineDA->getNumberOfTuples()),nbTuplesFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(fineStG));
743   if(fineDA->getNumberOfTuples()!=nbTuplesFineExp)
744     {
745       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFineGhostZone : Expecting " << nbTuplesFineExp << " tuples in fine DataArray having " << nbTuplesFine << " !";
746       throw INTERP_KERNEL::Exception(oss.str().c_str());
747     }
748   //
749   double *outPtr(fineDA->getPointer());
750   const double *inPtr(coarseDA->begin());
751   //
752   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
753   switch(meshDim)
754   {
755     case 1:
756       {
757         int offset(fineLocInCoarse[0].first+ghostSize-1),fact0(facts[0]);//offset is always >=0 thanks to the fact that ghostSize>=1 !
758         for(int i=0;i<ghostSize;i++)
759           outPtr=std::copy(inPtr+offset*nbCompo,inPtr+(offset+1)*nbCompo,outPtr);
760         outPtr+=nbCompo*fact0*dims[0];
761         offset=fineLocInCoarse[0].second+ghostSize;
762         for(int i=0;i<ghostSize;i++)
763           outPtr=std::copy(inPtr+offset*nbCompo,inPtr+(offset+1)*nbCompo,outPtr);
764         break;
765       }
766     case 2:
767       {
768         int nxwg(coarseSt[0]+2*ghostSize),fact0(facts[0]),fact1(facts[1]);
769         int kk(fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].first+ghostSize-1));//kk is always >=0 thanks to the fact that ghostSize>=1 !
770         for(int jg=0;jg<ghostSize;jg++)
771           {
772             for(int ig=0;ig<ghostSize;ig++)
773               outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
774             int kk0(kk+1);
775             for(int ig=0;ig<dims[0];ig++,kk0++)
776               for(int ifact=0;ifact<fact0;ifact++)
777                 outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
778             for(int ik=0;ik<ghostSize;ik++)
779               outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
780           }
781         for(int j=0;j<dims[1];j++)
782           {
783             kk=fineLocInCoarse[0].first-1+ghostSize+nxwg*(fineLocInCoarse[1].first+ghostSize+j);
784             for(int jfact=0;jfact<fact1;jfact++)
785               {
786                 for(int ig=0;ig<ghostSize;ig++)
787                   outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
788                 int kk0(kk+1+dims[0]);//1 not ghost. We make the hypothesis that factors is >= ghostlev
789                 outPtr+=fact0*nbCompo*dims[0];
790                 for(int ig=0;ig<ghostSize;ig++)
791                   outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
792               }
793           }
794         kk=fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].second+ghostSize);
795         for(int jg=0;jg<ghostSize;jg++)
796           {
797             for(int ig=0;ig<ghostSize;ig++)
798               outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
799             int kk0(kk+1);
800             for(int ig=0;ig<dims[0];ig++,kk0++)
801               for(int ifact=0;ifact<fact0;ifact++)
802                 outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
803             for(int ik=0;ik<ghostSize;ik++)
804               outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
805           }
806         break;
807       }
808     default:
809       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : only dimensions 1, 2 supported !");
810   }
811 }
812
813 void MEDCouplingIMesh::setSpaceDimension(int spaceDim)
814 {
815   if(spaceDim==_space_dim)
816     return ;
817   CheckSpaceDimension(spaceDim);
818   _space_dim=spaceDim;
819   declareAsNew();
820 }
821
822 void MEDCouplingIMesh::updateTime() const
823 {
824 }
825
826 std::size_t MEDCouplingIMesh::getHeapMemorySizeWithoutChildren() const
827 {
828   return MEDCouplingStructuredMesh::getHeapMemorySizeWithoutChildren();
829 }
830
831 std::vector<const BigMemoryObject *> MEDCouplingIMesh::getDirectChildren() const
832 {
833   return std::vector<const BigMemoryObject *>();
834 }
835
836 /*!
837  * This method copyies all tiny strings from other (name and components name).
838  * @throw if other and this have not same mesh type.
839  */
840 void MEDCouplingIMesh::copyTinyStringsFrom(const MEDCouplingMesh *other)
841
842   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
843   if(!otherC)
844     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::copyTinyStringsFrom : meshes have not same type !");
845   MEDCouplingStructuredMesh::copyTinyStringsFrom(other);
846   declareAsNew();
847 }
848
849 bool MEDCouplingIMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
850 {
851   if(!other)
852     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::isEqualIfNotWhy : input other pointer is null !");
853   const MEDCouplingIMesh *otherC(dynamic_cast<const MEDCouplingIMesh *>(other));
854   if(!otherC)
855     {
856       reason="mesh given in input is not castable in MEDCouplingIMesh !";
857       return false;
858     }
859   if(!MEDCouplingStructuredMesh::isEqualIfNotWhy(other,prec,reason))
860     return false;
861   if(!isEqualWithoutConsideringStrInternal(otherC,prec,reason))
862     return false;
863   if(_axis_unit!=otherC->_axis_unit)
864     {
865       reason="The units of axis are not the same !";
866       return false;
867     }
868   return true;
869 }
870
871 bool MEDCouplingIMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
872 {
873   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
874   if(!otherC)
875     return false;
876   std::string tmp;
877   return isEqualWithoutConsideringStrInternal(other,prec,tmp);
878 }
879
880 bool MEDCouplingIMesh::isEqualWithoutConsideringStrInternal(const MEDCouplingMesh *other, double prec, std::string& reason) const
881 {
882   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
883   if(!otherC)
884     return false;
885   if(_space_dim!=otherC->_space_dim)
886     {
887       std::ostringstream oss;
888       oss << "The spaceDimension of this (" << _space_dim << ") is not equal to those of other (" << otherC->_space_dim << ") !";
889       return false;
890     }
891   checkSpaceDimension();
892   for(int i=0;i<_space_dim;i++)
893     {
894       if(fabs(_origin[i]-otherC->_origin[i])>prec)
895         {
896           std::ostringstream oss;
897           oss << "The origin of this and other differs at " << i << " !";
898           reason=oss.str();
899           return false;
900         }
901     }
902   for(int i=0;i<_space_dim;i++)
903     {
904       if(fabs(_dxyz[i]-otherC->_dxyz[i])>prec)
905         {
906           std::ostringstream oss;
907           oss << "The delta of this and other differs at " << i << " !";
908           reason=oss.str();
909           return false;
910         }
911     }
912   for(int i=0;i<_space_dim;i++)
913     {
914       if(_structure[i]!=otherC->_structure[i])
915         {
916           std::ostringstream oss;
917           oss << "The structure of this and other differs at " << i << " !";
918           reason=oss.str();
919           return false;
920         }
921     }
922   return true;
923 }
924
925 void MEDCouplingIMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
926                                             DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const
927 {
928   if(!isEqualWithoutConsideringStr(other,prec))
929     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::checkDeepEquivalWith : Meshes are not the same !");
930 }
931
932 /*!
933  * Nothing is done here (except to check that the other is a ParaMEDMEM::MEDCouplingIMesh instance too).
934  * The user intend that the nodes are the same, so by construction of ParaMEDMEM::MEDCouplingIMesh, \a this and \a other are the same !
935  */
936 void MEDCouplingIMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
937                                                        DataArrayInt *&cellCor) const
938 {
939   if(!isEqualWithoutConsideringStr(other,prec))
940     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::checkDeepEquivalOnSameNodesWith : Meshes are not the same !");
941 }
942
943 void MEDCouplingIMesh::checkCoherency() const
944 {
945   checkSpaceDimension();
946   for(int i=0;i<_space_dim;i++)
947     if(_structure[i]<1)
948       {
949         std::ostringstream oss; oss << "MEDCouplingIMesh::checkCoherency : On axis " << i << "/" << _space_dim << ", number of nodes is equal to " << _structure[i] << " ! must be >=1 !";
950         throw INTERP_KERNEL::Exception(oss.str().c_str());
951       }
952 }
953
954 void MEDCouplingIMesh::checkCoherency1(double eps) const
955 {
956   checkCoherency();
957 }
958
959 void MEDCouplingIMesh::checkCoherency2(double eps) const
960 {
961   checkCoherency1(eps);
962 }
963
964 void MEDCouplingIMesh::getNodeGridStructure(int *res) const
965 {
966   checkSpaceDimension();
967   std::copy(_structure,_structure+_space_dim,res);
968 }
969
970 std::vector<int> MEDCouplingIMesh::getNodeGridStructure() const
971 {
972   checkSpaceDimension();
973   std::vector<int> ret(_structure,_structure+_space_dim);
974   return ret;
975 }
976
977 MEDCouplingStructuredMesh *MEDCouplingIMesh::buildStructuredSubPart(const std::vector< std::pair<int,int> >& cellPart) const
978 {
979   checkCoherency();
980   int dim(getSpaceDimension());
981   if(dim!=(int)cellPart.size())
982     {
983       std::ostringstream oss; oss << "MEDCouplingIMesh::buildStructuredSubPart : the space dimension is " << dim << " and cell part size is " << cellPart.size() << " !";
984       throw INTERP_KERNEL::Exception(oss.str().c_str());
985     }
986   double retOrigin[3]={0.,0.,0.};
987   int retStruct[3]={0,0,0};
988   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> ret(dynamic_cast<MEDCouplingIMesh *>(deepCpy()));
989   for(int i=0;i<dim;i++)
990     {
991       int startNode(cellPart[i].first),endNode(cellPart[i].second+1);
992       int myDelta(endNode-startNode);
993       if(startNode<0 || startNode>=_structure[i])
994         {
995           std::ostringstream oss; oss << "MEDCouplingIMesh::buildStructuredSubPart : At dimension #" << i << " the start node id is " << startNode << " it should be in [0," << _structure[i] << ") !";
996           throw INTERP_KERNEL::Exception(oss.str().c_str());
997         }
998       if(myDelta<0 || myDelta>_structure[i])
999         {
1000           std::ostringstream oss; oss << "MEDCouplingIMesh::buildStructuredSubPart : Along dimension #" << i << " the number of nodes is " << _structure[i] << ", and you are requesting for " << myDelta << " nodes wide range !" << std::endl;
1001           throw INTERP_KERNEL::Exception(oss.str().c_str());
1002         }
1003       retOrigin[i]=_origin[i]+startNode*_dxyz[i];
1004       retStruct[i]=myDelta;
1005     }
1006   ret->setNodeStruct(retStruct,retStruct+dim);
1007   ret->setOrigin(retOrigin,retOrigin+dim);
1008   ret->checkCoherency();
1009   return ret.retn();
1010 }
1011
1012 /*!
1013  * Return the space dimension of \a this.
1014  */
1015 int MEDCouplingIMesh::getSpaceDimension() const
1016 {
1017   return _space_dim;
1018 }
1019
1020 void MEDCouplingIMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const
1021 {
1022   int tmp[3];
1023   int spaceDim(getSpaceDimension());
1024   getSplitNodeValues(tmp);
1025   int tmp2[3];
1026   GetPosFromId(nodeId,spaceDim,tmp,tmp2);
1027   for(int j=0;j<spaceDim;j++)
1028     coo.push_back(_origin[j]+_dxyz[j]*tmp2[j]);
1029 }
1030
1031 std::string MEDCouplingIMesh::simpleRepr() const
1032 {
1033   std::ostringstream ret;
1034   ret << "Image grid with name : \"" << getName() << "\"\n";
1035   ret << "Description of mesh : \"" << getDescription() << "\"\n";
1036   int tmpp1,tmpp2;
1037   double tt(getTime(tmpp1,tmpp2));
1038   int spaceDim(_space_dim);
1039   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
1040   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
1041   ret << "Space dimension : " << spaceDim << "\n";
1042   if(spaceDim<0 || spaceDim>3)
1043     return ret.str();
1044   ret << "The nodal structure is : "; std::copy(_structure,_structure+spaceDim,std::ostream_iterator<int>(ret," ")); ret << "\n";
1045   ret << "The origin position is [" << _axis_unit << "]: ";
1046   std::copy(_origin,_origin+spaceDim,std::ostream_iterator<double>(ret," ")); ret << "\n";
1047   ret << "The intervals along axis are : ";
1048   std::copy(_dxyz,_dxyz+spaceDim,std::ostream_iterator<double>(ret," ")); ret << "\n";
1049   return ret.str();
1050 }
1051
1052 std::string MEDCouplingIMesh::advancedRepr() const
1053 {
1054   return simpleRepr();
1055 }
1056
1057 void MEDCouplingIMesh::getBoundingBox(double *bbox) const
1058 {
1059   checkCoherency();
1060   int dim(getSpaceDimension());
1061   for(int idim=0; idim<dim; idim++)
1062     {
1063       bbox[2*idim]=_origin[idim];
1064       int coeff(_structure[idim]);
1065       if(_structure[idim]<0)
1066         {
1067           std::ostringstream oss; oss << "MEDCouplingIMesh::getBoundingBox : on axis #" << idim << " number of nodes in structure is < 0 !";
1068           throw INTERP_KERNEL::Exception(oss.str().c_str());
1069         }
1070       if(_structure[idim]>1)
1071         coeff=_structure[idim]-1;
1072       bbox[2*idim+1]=_origin[idim]+_dxyz[idim]*coeff;
1073     }
1074 }
1075
1076 /*!
1077  * Returns a new MEDCouplingFieldDouble containing volumes of cells constituting \a this
1078  * mesh.<br>
1079  * For 1D cells, the returned field contains lengths.<br>
1080  * For 2D cells, the returned field contains areas.<br>
1081  * For 3D cells, the returned field contains volumes.
1082  *  \param [in] isAbs - a not used parameter.
1083  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on cells
1084  *         and one time . The caller is to delete this field using decrRef() as it is no
1085  *         more needed.
1086  */
1087 MEDCouplingFieldDouble *MEDCouplingIMesh::getMeasureField(bool isAbs) const
1088 {
1089   checkCoherency();
1090   std::string name="MeasureOfMesh_";
1091   name+=getName();
1092   int nbelem(getNumberOfCells());
1093   MEDCouplingFieldDouble *field(MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME));
1094   field->setName(name);
1095   DataArrayDouble* array(DataArrayDouble::New());
1096   array->alloc(nbelem,1);
1097   array->fillWithValue(getMeasureOfAnyCell());
1098   field->setArray(array) ;
1099   array->decrRef();
1100   field->setMesh(const_cast<MEDCouplingIMesh *>(this));
1101   field->synchronizeTimeWithMesh();
1102   return field;
1103 }
1104
1105 /*!
1106  * not implemented yet !
1107  */
1108 MEDCouplingFieldDouble *MEDCouplingIMesh::getMeasureFieldOnNode(bool isAbs) const
1109 {
1110   throw INTERP_KERNEL::Exception("MEDCouplingIMesh::getMeasureFieldOnNode : not implemented yet !");
1111   //return 0;
1112 }
1113
1114 int MEDCouplingIMesh::getCellContainingPoint(const double *pos, double eps) const
1115 {
1116   int dim(getSpaceDimension()),ret(0),coeff(1);
1117   for(int i=0;i<dim;i++)
1118     {
1119       int nbOfCells(_structure[i]-1);
1120       double ref(pos[i]);
1121       int tmp((int)((ref-_origin[i])/_dxyz[i]));
1122       if(tmp>=0 && tmp<nbOfCells)
1123         {
1124           ret+=coeff*tmp;
1125           coeff*=nbOfCells;
1126         }
1127       else
1128         return -1;
1129     }
1130   return ret;
1131 }
1132
1133 void MEDCouplingIMesh::rotate(const double *center, const double *vector, double angle)
1134 {
1135   throw INTERP_KERNEL::Exception("No rotation available on IMesh : Traduce it to unstructured mesh to apply it !");
1136 }
1137
1138 /*!
1139  * Translates all nodes of \a this mesh by a given vector. Actually, it adds each
1140  * component of the \a vector to all node coordinates of a corresponding axis.
1141  *  \param [in] vector - the translation vector whose size must be not less than \a
1142  *         this->getSpaceDimension().
1143  */
1144 void MEDCouplingIMesh::translate(const double *vector)
1145 {
1146   checkSpaceDimension();
1147   int dim(getSpaceDimension());
1148   std::transform(_origin,_origin+dim,vector,_origin,std::plus<double>());
1149   declareAsNew();
1150 }
1151
1152 /*!
1153  * Applies scaling transformation to all nodes of \a this mesh.
1154  *  \param [in] point - coordinates of a scaling center. This array is to be of
1155  *         size \a this->getSpaceDimension() at least.
1156  *  \param [in] factor - a scale factor.
1157  */
1158 void MEDCouplingIMesh::scale(const double *point, double factor)
1159 {
1160   checkSpaceDimension();
1161   int dim(getSpaceDimension());
1162   std::transform(_origin,_origin+dim,point,_origin,std::minus<double>());
1163   std::transform(_origin,_origin+dim,_origin,std::bind2nd(std::multiplies<double>(),factor));
1164   std::transform(_dxyz,_dxyz+dim,_dxyz,std::bind2nd(std::multiplies<double>(),factor));
1165   std::transform(_origin,_origin+dim,point,_origin,std::plus<double>());
1166   declareAsNew();
1167 }
1168
1169 MEDCouplingMesh *MEDCouplingIMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
1170 {
1171   //not implemented yet !
1172   return 0;
1173 }
1174
1175 /*!
1176  * Returns a new DataArrayDouble holding coordinates of all nodes of \a this mesh.
1177  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
1178  *          this->getNumberOfNodes() tuples per \a this->getSpaceDimension()
1179  *          components. The caller is to delete this array using decrRef() as it is
1180  *          no more needed.
1181  */
1182 DataArrayDouble *MEDCouplingIMesh::getCoordinatesAndOwner() const
1183 {
1184   checkCoherency();
1185   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
1186   int spaceDim(getSpaceDimension()),nbNodes(getNumberOfNodes());
1187   ret->alloc(nbNodes,spaceDim);
1188   double *pt(ret->getPointer());
1189   ret->setInfoOnComponents(buildInfoOnComponents());
1190   int tmp2[3],tmp[3];
1191   getSplitNodeValues(tmp);
1192   for(int i=0;i<nbNodes;i++)
1193     {
1194       GetPosFromId(i,spaceDim,tmp,tmp2);
1195       for(int j=0;j<spaceDim;j++)
1196         pt[i*spaceDim+j]=_dxyz[j]*tmp2[j]+_origin[j];
1197     }
1198   return ret.retn();
1199 }
1200
1201 /*!
1202  * Returns a new DataArrayDouble holding barycenters of all cells. The barycenter is
1203  * computed by averaging coordinates of cell nodes.
1204  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
1205  *          this->getNumberOfCells() tuples per \a this->getSpaceDimension()
1206  *          components. The caller is to delete this array using decrRef() as it is
1207  *          no more needed.
1208  */
1209 DataArrayDouble *MEDCouplingIMesh::getBarycenterAndOwner() const
1210 {
1211   checkCoherency();
1212   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
1213   int spaceDim(getSpaceDimension()),nbCells(getNumberOfCells()),tmp[3],tmp2[3];
1214   ret->alloc(nbCells,spaceDim);
1215   double *pt(ret->getPointer()),shiftOrigin[3];
1216   std::transform(_dxyz,_dxyz+spaceDim,shiftOrigin,std::bind2nd(std::multiplies<double>(),0.5));
1217   std::transform(_origin,_origin+spaceDim,shiftOrigin,shiftOrigin,std::plus<double>());
1218   getSplitCellValues(tmp);
1219   ret->setInfoOnComponents(buildInfoOnComponents());
1220   for(int i=0;i<nbCells;i++)
1221     {
1222       GetPosFromId(i,spaceDim,tmp,tmp2);
1223       for(int j=0;j<spaceDim;j++)
1224         pt[i*spaceDim+j]=_dxyz[j]*tmp2[j]+shiftOrigin[j];
1225     }
1226   return ret.retn();
1227 }
1228
1229 DataArrayDouble *MEDCouplingIMesh::computeIsoBarycenterOfNodesPerCell() const
1230 {
1231   return MEDCouplingIMesh::getBarycenterAndOwner();
1232 }
1233
1234 void MEDCouplingIMesh::renumberCells(const int *old2NewBg, bool check)
1235 {
1236   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for IMesh !");
1237 }
1238
1239 void MEDCouplingIMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
1240 {
1241   int it,order;
1242   double time(getTime(it,order));
1243   tinyInfo.clear();
1244   tinyInfoD.clear();
1245   littleStrings.clear();
1246   littleStrings.push_back(getName());
1247   littleStrings.push_back(getDescription());
1248   littleStrings.push_back(getTimeUnit());
1249   littleStrings.push_back(getAxisUnit());
1250   tinyInfo.push_back(it);
1251   tinyInfo.push_back(order);
1252   tinyInfo.push_back(_space_dim);
1253   tinyInfo.insert(tinyInfo.end(),_structure,_structure+3);
1254   tinyInfoD.push_back(time);
1255   tinyInfoD.insert(tinyInfoD.end(),_dxyz,_dxyz+3);
1256   tinyInfoD.insert(tinyInfoD.end(),_origin,_origin+3);
1257 }
1258
1259 void MEDCouplingIMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
1260 {
1261   a1->alloc(0,1);
1262   a2->alloc(0,1);
1263 }
1264
1265 void MEDCouplingIMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
1266 {
1267   a1=DataArrayInt::New();
1268   a1->alloc(0,1);
1269   a2=DataArrayDouble::New();
1270   a2->alloc(0,1);
1271 }
1272
1273 void MEDCouplingIMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
1274                                        const std::vector<std::string>& littleStrings)
1275 {
1276   setName(littleStrings[0]);
1277   setDescription(littleStrings[1]);
1278   setTimeUnit(littleStrings[2]);
1279   setAxisUnit(littleStrings[3]);
1280   setTime(tinyInfoD[0],tinyInfo[0],tinyInfo[1]);
1281   _space_dim=tinyInfo[2];
1282   _structure[0]=tinyInfo[3]; _structure[1]=tinyInfo[4]; _structure[2]=tinyInfo[5];
1283   _dxyz[0]=tinyInfoD[1]; _dxyz[1]=tinyInfoD[2]; _dxyz[2]=tinyInfoD[3];
1284   _origin[0]=tinyInfoD[4]; _origin[1]=tinyInfoD[5]; _origin[2]=tinyInfoD[6];
1285   declareAsNew();
1286 }
1287
1288 void MEDCouplingIMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
1289 {
1290   checkCoherency();
1291   std::ostringstream extent,origin,spacing;
1292   for(int i=0;i<3;i++)
1293     {
1294       if(i<_space_dim)
1295         { extent << "0 " <<  _structure[i]-1 << " "; origin << _origin[i] << " "; spacing << _dxyz[i] << " "; }
1296       else
1297         { extent << "0 0 "; origin << "0 "; spacing << "0 "; }
1298     }
1299   ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\" Origin=\"" << origin.str() << "\" Spacing=\"" << spacing.str() << "\">\n";
1300   ofs << "    <Piece Extent=\"" << extent.str() << "\">\n";
1301   ofs << "      <PointData>\n" << pointData << std::endl;
1302   ofs << "      </PointData>\n";
1303   ofs << "      <CellData>\n" << cellData << std::endl;
1304   ofs << "      </CellData>\n";
1305   ofs << "      <Coordinates>\n";
1306   ofs << "      </Coordinates>\n";
1307   ofs << "    </Piece>\n";
1308   ofs << "  </" << getVTKDataSetType() << ">\n";
1309 }
1310
1311 void MEDCouplingIMesh::reprQuickOverview(std::ostream& stream) const
1312 {
1313   stream << "MEDCouplingIMesh C++ instance at " << this << ". Name : \"" << getName() << "\". Space dimension : " << _space_dim << ".";
1314   if(_space_dim<0 || _space_dim>3)
1315     return ;
1316   stream << "\n";
1317   std::ostringstream stream0,stream1;
1318   int nbNodes(1),nbCells(0);
1319   bool isPb(false);
1320   for(int i=0;i<_space_dim;i++)
1321     {
1322       char tmp('X'+i);
1323       int tmpNodes(_structure[i]);
1324       stream1 << "- Axis " << tmp << " : " << tmpNodes << " nodes (orig=" << _origin[i] << ", inter=" << _dxyz[i] << ").";
1325       if(i!=_space_dim-1)
1326         stream1 << std::endl;
1327       if(tmpNodes>=1)
1328         nbNodes*=tmpNodes;
1329       else
1330         isPb=true;
1331       if(tmpNodes>=2)
1332         nbCells=nbCells==0?tmpNodes-1:nbCells*(tmpNodes-1);
1333     }
1334   if(!isPb)
1335     {
1336       stream0 << "Number of cells : " << nbCells << ", Number of nodes : " << nbNodes;
1337       stream << stream0.str();
1338       if(_space_dim>0)
1339         stream << std::endl;
1340     }
1341   stream << stream1.str();
1342 }
1343
1344 std::string MEDCouplingIMesh::getVTKDataSetType() const
1345 {
1346   return std::string("ImageData");
1347 }
1348
1349 std::vector<std::string> MEDCouplingIMesh::buildInfoOnComponents() const
1350 {
1351   checkSpaceDimension();
1352   int dim(getSpaceDimension());
1353   std::vector<std::string> ret(dim);
1354   for(int i=0;i<dim;i++)
1355     {
1356       std::ostringstream oss;
1357       char tmp('X'+i); oss << tmp;
1358       ret[i]=DataArray::BuildInfoFromVarAndUnit(oss.str(),_axis_unit);
1359     }
1360   return ret;
1361 }
1362
1363 void MEDCouplingIMesh::checkSpaceDimension() const
1364 {
1365   CheckSpaceDimension(_space_dim);
1366 }
1367
1368 void MEDCouplingIMesh::CheckSpaceDimension(int spaceDim)
1369 {
1370   if(spaceDim<0 || spaceDim>3)
1371     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CheckSpaceDimension : input spaceDim must be in [0,1,2,3] !");
1372 }
1373
1374 int MEDCouplingIMesh::FindIntRoot(int val, int order)
1375 {
1376   if(order==0)
1377     return 1;
1378   if(val<0)
1379     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : input val is < 0 ! Not possible to compute a root !");
1380   if(order==1)
1381     return val;
1382   if(order!=2 && order!=3)
1383     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the order available are 0,1,2 or 3 !");
1384   double valf((double)val);
1385   if(order==2)
1386     {
1387       double retf(sqrt(valf));
1388       int ret((int)retf);
1389       if(ret*ret!=val)
1390         throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect square root !");
1391       return ret;
1392     }
1393   else//order==3
1394     {
1395       double retf(std::pow(val,0.3333333333333333));
1396       int ret((int)retf),ret2(ret+1);
1397       if(ret*ret*ret!=val && ret2*ret2*ret2!=val)
1398         throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect cublic root !");
1399       if(ret*ret*ret==val)
1400         return ret;
1401       else
1402         return ret2;
1403     }
1404 }