Salome HOME
Attributes for AMR cartesian mesh developped. Ready to test.
[modules/med.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(nbTuplesFine%nbOfTuplesInFineExp!=0)
298     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : Invalid nb of tuples in fine DataArray regarding its structure !");
299   int fact(std::accumulate(facts.begin(),facts.end(),1,std::multiplies<int>()));
300   if(nbTuplesFine!=fact*nbOfTuplesInFineExp)
301     {
302       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarse : Invalid number of tuples ("  << nbTuplesFine << ") of fine dataarray is invalid ! Must be " << fact*nbOfTuplesInFineExp << "!";
303       throw INTERP_KERNEL::Exception(oss.str().c_str());
304     }
305   // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
306   double *outPtr(coarseDA->getPointer());
307   const double *inPtr(fineDA->begin());
308   //
309   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
310   switch(meshDim)
311   {
312     case 1:
313       {
314         int offset(fineLocInCoarse[0].first),fact0(facts[0]);
315         for(int i=0;i<dims[0];i++)
316           {
317             double *loc(outPtr+(offset+i)*nbCompo);
318             for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
319               {
320                 if(ifact!=0)
321                   std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
322                 else
323                   std::copy(inPtr,inPtr+nbCompo,loc);
324               }
325           }
326         break;
327       }
328     case 2:
329       {
330         int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first),fact1(facts[1]),fact0(facts[0]);
331         for(int j=0;j<dims[1];j++)
332           {
333             for(int jfact=0;jfact<fact1;jfact++)
334               {
335                 for(int i=0;i<dims[0];i++)
336                   {
337                     double *loc(outPtr+(kk+i)*nbCompo);
338                     for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
339                       {
340                         if(jfact!=0 || ifact!=0)
341                           std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
342                         else
343                           std::copy(inPtr,inPtr+nbCompo,loc);
344                       }
345                   }
346               }
347             kk+=coarseSt[0];
348           }
349         break;
350       }
351     case 3:
352       {
353         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]);
354         for(int k=0;k<dims[2];k++)
355           {
356             for(int kfact=0;kfact<fact2;kfact++)
357               {
358                 for(int j=0;j<dims[1];j++)
359                   {
360                     for(int jfact=0;jfact<fact1;jfact++)
361                       {
362                         for(int i=0;i<dims[0];i++)
363                           {
364                             double *loc(outPtr+(kk+i+j*coarseSt[0])*nbCompo);
365                             for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
366                               {
367                                 if(kfact!=0 || jfact!=0 || ifact!=0)
368                                   std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
369                                 else
370                                   std::copy(inPtr,inPtr+nbCompo,loc);
371                               }
372                           }
373                       }
374                   }
375               }
376             kk+=coarseSt[0]*coarseSt[1];
377           }
378         break;
379       }
380     default:
381       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : only dimensions 1, 2 and 3 supported !");
382   }
383 }
384
385 /*!
386  * This static method is useful to condense field on cells of a MEDCouplingIMesh instance coming from a refinement ( MEDCouplingIMesh::refineWithFactor for example)
387  * to a coarse MEDCouplingIMesh instance. So this method can be seen as a specialization in P0P0 conservative interpolation non overlaping from fine image mesh
388  * 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.
389  *
390  * \param [in] coarseSt The cell structure of coarse mesh.
391  * \param [in] fineDA The DataArray containing the cell field on uniformly refined mesh
392  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
393  * \param [in] facts The refinement coefficient per axis.
394  * \param [in,out] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
395  * \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.
396  *
397  * \sa CondenseFineToCoarse,SpreadCoarseToFineGhost
398  */
399 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)
400 {
401   if(ghostSize<0)
402     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : ghost level has to be >= 0 !");
403   if(coarseSt.size()!=fineLocInCoarse.size() || coarseSt.size()!=facts.size())
404     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : All input vectors (dimension) must have the same size !");
405   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
406     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : the parameters 1 or 3 are NULL or not allocated !");
407   std::vector<int> coarseStG(coarseSt.size()); std::transform(coarseSt.begin(),coarseSt.end(),coarseStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
408   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseStG));
409   int nbCompo(fineDA->getNumberOfComponents());
410   if(coarseDA->getNumberOfComponents()!=nbCompo)
411     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : the number of components of fine DA and coarse one mismatches !");
412   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
413     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) !");
414   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
415     {
416       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarseGhost : Expecting " << nbOfTuplesInCoarseExp << " tuples in coarse DataArray having " << coarseDA->getNumberOfTuples() << " !";
417       throw INTERP_KERNEL::Exception(oss.str().c_str());
418     }
419   //
420   std::vector<int> fineStG(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
421   std::transform(fineStG.begin(),fineStG.end(),facts.begin(),fineStG.begin(),std::multiplies<int>());
422   std::transform(fineStG.begin(),fineStG.end(),fineStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
423   int nbTuplesFine(fineDA->getNumberOfTuples()),nbTuplesFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(fineStG));
424   if(fineDA->getNumberOfTuples()!=nbTuplesFineExp)
425     {
426       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarseGhost : Expecting " << nbTuplesFineExp << " tuples in fine DataArray having " << nbTuplesFine << " !";
427       throw INTERP_KERNEL::Exception(oss.str().c_str());
428     }
429   //
430   double *outPtr(coarseDA->getPointer());
431   const double *inPtr(fineDA->begin());
432   //
433   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
434   switch(meshDim)
435   {
436     case 1:
437       {
438         int offset(fineLocInCoarse[0].first+ghostSize),fact0(facts[0]);
439         inPtr+=ghostSize*nbCompo;
440         for(int i=0;i<dims[0];i++)
441           {
442             double *loc(outPtr+(offset+i)*nbCompo);
443             for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
444               {
445                 if(ifact!=0)
446                   std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
447                 else
448                   std::copy(inPtr,inPtr+nbCompo,loc);
449               }
450           }
451         break;
452       }
453     case 2:
454       {
455         int nxwg(coarseSt[0]+2*ghostSize);
456         int kk(fineLocInCoarse[0].first+ghostSize+nxwg*(fineLocInCoarse[1].first+ghostSize)),fact1(facts[1]),fact0(facts[0]);
457         inPtr+=(dims[0]*fact0+2*ghostSize)*nbCompo;
458         for(int j=0;j<dims[1];j++)
459           {
460              for(int jfact=0;jfact<fact1;jfact++)
461               {
462                 inPtr+=ghostSize*nbCompo;
463                 for(int i=0;i<dims[0];i++)
464                   {
465                     double *loc(outPtr+(kk+i)*nbCompo);
466                     for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
467                       {
468                         if(jfact!=0 || ifact!=0)
469                           std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
470                         else
471                           std::copy(inPtr,inPtr+nbCompo,loc);
472                       }
473                   }
474                 inPtr+=ghostSize*nbCompo;
475               }
476             kk+=nxwg;
477           }
478         break;
479       }
480     default:
481       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarseGhost : only dimensions 1, 2 supported !");
482   }
483 }
484
485 /*!
486  * This method spreads the values of coarse data \a coarseDA into \a fineDA.
487  *
488  * \param [in] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
489  * \param [in] coarseSt The cell structure of coarse mesh.
490  * \param [in,out] fineDA The DataArray containing the cell field on uniformly refined mesh
491  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
492  * \param [in] facts The refinement coefficient per axis.
493  * \sa SpreadCoarseToFineGhost, CondenseFineToCoarse
494  */
495 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)
496 {
497   if(coarseSt.size()!=fineLocInCoarse.size() || coarseSt.size()!=facts.size())
498       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : All input vectors (dimension) must have the same size !");
499   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
500     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the parameters 1 or 3 are NULL or not allocated !");
501   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseSt)),nbOfTuplesInFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(fineLocInCoarse));
502   int nbCompo(fineDA->getNumberOfComponents());
503   if(coarseDA->getNumberOfComponents()!=nbCompo)
504     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the number of components of fine DA and coarse one mismatches !");
505   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
506     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) !");
507   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
508     {
509       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFine : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
510       throw INTERP_KERNEL::Exception(oss.str().c_str());
511     }
512   int nbTuplesFine(fineDA->getNumberOfTuples());
513   if(nbTuplesFine%nbOfTuplesInFineExp!=0)
514     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : Invalid nb of tuples in fine DataArray regarding its structure !");
515   int fact(std::accumulate(facts.begin(),facts.end(),1,std::multiplies<int>()));
516   if(nbTuplesFine!=fact*nbOfTuplesInFineExp)
517     {
518       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFine : Invalid number of tuples ("  << nbTuplesFine << ") of fine dataarray is invalid ! Must be " << fact*nbOfTuplesInFineExp << "!";
519       throw INTERP_KERNEL::Exception(oss.str().c_str());
520     }
521   // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
522   double *outPtr(fineDA->getPointer());
523   const double *inPtr(coarseDA->begin());
524   //
525   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
526   switch(meshDim)
527   {
528     case 1:
529       {
530         int offset(fineLocInCoarse[0].first),fact0(facts[0]);
531         for(int i=0;i<dims[0];i++)
532           {
533             const double *loc(inPtr+(offset+i)*nbCompo);
534             for(int ifact=0;ifact<fact0;ifact++)
535               outPtr=std::copy(loc,loc+nbCompo,outPtr);
536           }
537         break;
538       }
539     case 2:
540       {
541         int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first),fact0(facts[0]),fact1(facts[1]);
542         for(int j=0;j<dims[1];j++)
543           {
544             for(int jfact=0;jfact<fact1;jfact++)
545               {
546                 for(int i=0;i<dims[0];i++)
547                   {
548                     const double *loc(inPtr+(kk+i)*nbCompo);
549                     for(int ifact=0;ifact<fact0;ifact++)
550                       outPtr=std::copy(loc,loc+nbCompo,outPtr);
551                   }
552               }
553             kk+=coarseSt[0];
554           }
555         break;
556       }
557     case 3:
558       {
559         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]);
560         for(int k=0;k<dims[2];k++)
561           {
562             for(int kfact=0;kfact<fact2;kfact++)
563               {
564                 for(int j=0;j<dims[1];j++)
565                   {
566                     for(int jfact=0;jfact<fact1;jfact++)
567                       {
568                         for(int i=0;i<dims[0];i++)
569                           {
570                             const double *loc(inPtr+(kk+i+j*coarseSt[0])*nbCompo);
571                             for(int ifact=0;ifact<fact0;ifact++)
572                               outPtr=std::copy(loc,loc+nbCompo,outPtr);
573                           }
574                       }
575                   }
576               }
577             kk+=coarseSt[0]*coarseSt[1];
578           }
579         break;
580       }
581     default:
582       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : only dimensions 1, 2 and 3 supported !");
583   }
584 }
585
586 /*!
587  * This method spreads the values of coarse data \a coarseDA into \a fineDA.
588  *
589  * \param [in] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
590  * \param [in] coarseSt The cell structure of coarse mesh.
591  * \param [in,out] fineDA The DataArray containing the cell field on uniformly refined mesh
592  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
593  * \param [in] facts The refinement coefficient per axis.
594  * \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.
595  *
596  * \sa CondenseFineToCoarse,SpreadCoarseToFineGhostZone
597  */
598 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)
599 {
600   if(ghostSize<0)
601     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : ghost level has to be >= 0 !");
602   if(coarseSt.size()!=fineLocInCoarse.size() || coarseSt.size()!=facts.size())
603     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : All input vectors (dimension) must have the same size !");
604   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
605     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : the parameters 1 or 3 are NULL or not allocated !");
606   std::vector<int> coarseStG(coarseSt.size()); std::transform(coarseSt.begin(),coarseSt.end(),coarseStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
607   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseStG));
608   int nbCompo(fineDA->getNumberOfComponents());
609   if(coarseDA->getNumberOfComponents()!=nbCompo)
610     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : the number of components of fine DA and coarse one mismatches !");
611   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
612     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) !");
613   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
614     {
615       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFineGhost : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
616       throw INTERP_KERNEL::Exception(oss.str().c_str());
617     }
618   //
619   std::vector<int> fineStG(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
620   std::transform(fineStG.begin(),fineStG.end(),facts.begin(),fineStG.begin(),std::multiplies<int>());
621   std::transform(fineStG.begin(),fineStG.end(),fineStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
622   int nbTuplesFine(fineDA->getNumberOfTuples()),nbTuplesFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(fineStG));
623   if(fineDA->getNumberOfTuples()!=nbTuplesFineExp)
624     {
625       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFineGhost : Expecting " << nbTuplesFineExp << " tuples in fine DataArray having " << nbTuplesFine << " !";
626       throw INTERP_KERNEL::Exception(oss.str().c_str());
627     }
628   //
629   double *outPtr(fineDA->getPointer());
630   const double *inPtr(coarseDA->begin());
631   //
632   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
633   switch(meshDim)
634   {
635     case 1:
636       {
637         int offset(fineLocInCoarse[0].first+ghostSize-1),fact0(facts[0]);//offset is always >=0 thanks to the fact that ghostSize>=1 !
638         for(int i=0;i<ghostSize;i++)
639           outPtr=std::copy(inPtr+offset*nbCompo,inPtr+(offset+1)*nbCompo,outPtr);
640         offset=fineLocInCoarse[0].first+ghostSize;
641         for(int i=0;i<dims[0];i++)
642           {
643             const double *loc(inPtr+(offset+i)*nbCompo);
644             for(int ifact=0;ifact<fact0;ifact++)
645               outPtr=std::copy(loc,loc+nbCompo,outPtr);
646           }
647         offset=fineLocInCoarse[0].second+ghostSize;
648         for(int i=0;i<ghostSize;i++)
649           outPtr=std::copy(inPtr+offset*nbCompo,inPtr+(offset+1)*nbCompo,outPtr);
650         break;
651       }
652     case 2:
653       {
654         int nxwg(coarseSt[0]+2*ghostSize),fact0(facts[0]),fact1(facts[1]);
655         int kk(fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].first+ghostSize-1));//kk is always >=0 thanks to the fact that ghostSize>=1 !
656         for(int jg=0;jg<ghostSize;jg++)
657           {
658             for(int ig=0;ig<ghostSize;ig++)
659               outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
660             int kk0(kk+1);
661             for(int ig=0;ig<dims[0];ig++,kk0++)
662               for(int ifact=0;ifact<fact0;ifact++)
663                 outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
664             for(int ik=0;ik<ghostSize;ik++)
665               outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
666           }
667         for(int j=0;j<dims[1];j++)
668           {
669             kk=fineLocInCoarse[0].first-1+ghostSize+nxwg*(fineLocInCoarse[1].first+ghostSize+j);
670             for(int jfact=0;jfact<fact1;jfact++)
671               {
672                 for(int ig=0;ig<ghostSize;ig++)
673                   outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
674                 int kk0(kk+1);//1 not ghost. We make the hypothesis that factors is >= ghostlev
675                 for(int i=0;i<dims[0];i++,kk0++)
676                   {
677                     const double *loc(inPtr+kk0*nbCompo);
678                     for(int ifact=0;ifact<fact0;ifact++)
679                       outPtr=std::copy(loc,loc+nbCompo,outPtr);
680                   }
681                 for(int ig=0;ig<ghostSize;ig++)
682                   outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
683               }
684           }
685         kk=fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].second+ghostSize);
686         for(int jg=0;jg<ghostSize;jg++)
687           {
688             for(int ig=0;ig<ghostSize;ig++)
689               outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
690             int kk0(kk+1);
691             for(int ig=0;ig<dims[0];ig++,kk0++)
692               for(int ifact=0;ifact<fact0;ifact++)
693                 outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
694             for(int ik=0;ik<ghostSize;ik++)
695               outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
696           }
697         break;
698       }
699     default:
700       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhost : only dimensions 1, 2 supported !");
701   }
702 }
703
704 /*!
705  * 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).
706  *
707  * \param [in] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
708  * \param [in] coarseSt The cell structure of coarse mesh.
709  * \param [in,out] fineDA The DataArray containing the cell field on uniformly refined mesh
710  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
711  * \param [in] facts The refinement coefficient per axis.
712  * \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.
713  *
714  * \sa SpreadCoarseToFineGhost
715  */
716 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)
717 {
718   if(ghostSize<0)
719     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : ghost level has to be >= 0 !");
720   if(coarseSt.size()!=fineLocInCoarse.size() || coarseSt.size()!=facts.size())
721     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : All input vectors (dimension) must have the same size !");
722   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
723     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : the parameters 1 or 3 are NULL or not allocated !");
724   std::vector<int> coarseStG(coarseSt.size()); std::transform(coarseSt.begin(),coarseSt.end(),coarseStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
725   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseStG));
726   int nbCompo(fineDA->getNumberOfComponents());
727   if(coarseDA->getNumberOfComponents()!=nbCompo)
728     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : the number of components of fine DA and coarse one mismatches !");
729   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
730     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) !");
731   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
732     {
733       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFineGhostZone : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
734       throw INTERP_KERNEL::Exception(oss.str().c_str());
735     }
736   //
737   std::vector<int> fineStG(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
738   std::transform(fineStG.begin(),fineStG.end(),facts.begin(),fineStG.begin(),std::multiplies<int>());
739   std::transform(fineStG.begin(),fineStG.end(),fineStG.begin(),std::bind2nd(std::plus<int>(),2*ghostSize));
740   int nbTuplesFine(fineDA->getNumberOfTuples()),nbTuplesFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(fineStG));
741   if(fineDA->getNumberOfTuples()!=nbTuplesFineExp)
742     {
743       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFineGhostZone : Expecting " << nbTuplesFineExp << " tuples in fine DataArray having " << nbTuplesFine << " !";
744       throw INTERP_KERNEL::Exception(oss.str().c_str());
745     }
746   //
747   double *outPtr(fineDA->getPointer());
748   const double *inPtr(coarseDA->begin());
749   //
750   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
751   switch(meshDim)
752   {
753     case 1:
754       {
755         int offset(fineLocInCoarse[0].first+ghostSize-1),fact0(facts[0]);//offset is always >=0 thanks to the fact that ghostSize>=1 !
756         for(int i=0;i<ghostSize;i++)
757           outPtr=std::copy(inPtr+offset*nbCompo,inPtr+(offset+1)*nbCompo,outPtr);
758         outPtr+=nbCompo*fact0*dims[0];
759         offset=fineLocInCoarse[0].second+ghostSize;
760         for(int i=0;i<ghostSize;i++)
761           outPtr=std::copy(inPtr+offset*nbCompo,inPtr+(offset+1)*nbCompo,outPtr);
762         break;
763       }
764     case 2:
765       {
766         int nxwg(coarseSt[0]+2*ghostSize),fact0(facts[0]),fact1(facts[1]);
767         int kk(fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].first+ghostSize-1));//kk is always >=0 thanks to the fact that ghostSize>=1 !
768         for(int jg=0;jg<ghostSize;jg++)
769           {
770             for(int ig=0;ig<ghostSize;ig++)
771               outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
772             int kk0(kk+1);
773             for(int ig=0;ig<dims[0];ig++,kk0++)
774               for(int ifact=0;ifact<fact0;ifact++)
775                 outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
776             for(int ik=0;ik<ghostSize;ik++)
777               outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
778           }
779         for(int j=0;j<dims[1];j++)
780           {
781             kk=fineLocInCoarse[0].first-1+ghostSize+nxwg*(fineLocInCoarse[1].first+ghostSize+j);
782             for(int jfact=0;jfact<fact1;jfact++)
783               {
784                 for(int ig=0;ig<ghostSize;ig++)
785                   outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
786                 int kk0(kk+1+dims[0]);//1 not ghost. We make the hypothesis that factors is >= ghostlev
787                 outPtr+=fact0*nbCompo*dims[0];
788                 for(int ig=0;ig<ghostSize;ig++)
789                   outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
790               }
791           }
792         kk=fineLocInCoarse[0].first+ghostSize-1+nxwg*(fineLocInCoarse[1].second+ghostSize);
793         for(int jg=0;jg<ghostSize;jg++)
794           {
795             for(int ig=0;ig<ghostSize;ig++)
796               outPtr=std::copy(inPtr+kk*nbCompo,inPtr+(kk+1)*nbCompo,outPtr);
797             int kk0(kk+1);
798             for(int ig=0;ig<dims[0];ig++,kk0++)
799               for(int ifact=0;ifact<fact0;ifact++)
800                 outPtr=std::copy(inPtr+(kk0)*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
801             for(int ik=0;ik<ghostSize;ik++)
802               outPtr=std::copy(inPtr+kk0*nbCompo,inPtr+(kk0+1)*nbCompo,outPtr);
803           }
804         break;
805       }
806     default:
807       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFineGhostZone : only dimensions 1, 2 supported !");
808   }
809 }
810
811 void MEDCouplingIMesh::setSpaceDimension(int spaceDim)
812 {
813   if(spaceDim==_space_dim)
814     return ;
815   CheckSpaceDimension(spaceDim);
816   _space_dim=spaceDim;
817   declareAsNew();
818 }
819
820 void MEDCouplingIMesh::updateTime() const
821 {
822 }
823
824 std::size_t MEDCouplingIMesh::getHeapMemorySizeWithoutChildren() const
825 {
826   return MEDCouplingStructuredMesh::getHeapMemorySizeWithoutChildren();
827 }
828
829 std::vector<const BigMemoryObject *> MEDCouplingIMesh::getDirectChildren() const
830 {
831   return std::vector<const BigMemoryObject *>();
832 }
833
834 /*!
835  * This method copyies all tiny strings from other (name and components name).
836  * @throw if other and this have not same mesh type.
837  */
838 void MEDCouplingIMesh::copyTinyStringsFrom(const MEDCouplingMesh *other)
839
840   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
841   if(!otherC)
842     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::copyTinyStringsFrom : meshes have not same type !");
843   MEDCouplingStructuredMesh::copyTinyStringsFrom(other);
844   declareAsNew();
845 }
846
847 bool MEDCouplingIMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
848 {
849   if(!other)
850     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::isEqualIfNotWhy : input other pointer is null !");
851   const MEDCouplingIMesh *otherC(dynamic_cast<const MEDCouplingIMesh *>(other));
852   if(!otherC)
853     {
854       reason="mesh given in input is not castable in MEDCouplingIMesh !";
855       return false;
856     }
857   if(!MEDCouplingStructuredMesh::isEqualIfNotWhy(other,prec,reason))
858     return false;
859   if(!isEqualWithoutConsideringStrInternal(otherC,prec,reason))
860     return false;
861   if(_axis_unit!=otherC->_axis_unit)
862     {
863       reason="The units of axis are not the same !";
864       return false;
865     }
866   return true;
867 }
868
869 bool MEDCouplingIMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
870 {
871   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
872   if(!otherC)
873     return false;
874   std::string tmp;
875   return isEqualWithoutConsideringStrInternal(other,prec,tmp);
876 }
877
878 bool MEDCouplingIMesh::isEqualWithoutConsideringStrInternal(const MEDCouplingMesh *other, double prec, std::string& reason) const
879 {
880   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
881   if(!otherC)
882     return false;
883   if(_space_dim!=otherC->_space_dim)
884     {
885       std::ostringstream oss;
886       oss << "The spaceDimension of this (" << _space_dim << ") is not equal to those of other (" << otherC->_space_dim << ") !";
887       return false;
888     }
889   checkSpaceDimension();
890   for(int i=0;i<_space_dim;i++)
891     {
892       if(fabs(_origin[i]-otherC->_origin[i])>prec)
893         {
894           std::ostringstream oss;
895           oss << "The origin of this and other differs at " << i << " !";
896           reason=oss.str();
897           return false;
898         }
899     }
900   for(int i=0;i<_space_dim;i++)
901     {
902       if(fabs(_dxyz[i]-otherC->_dxyz[i])>prec)
903         {
904           std::ostringstream oss;
905           oss << "The delta of this and other differs at " << i << " !";
906           reason=oss.str();
907           return false;
908         }
909     }
910   for(int i=0;i<_space_dim;i++)
911     {
912       if(_structure[i]!=otherC->_structure[i])
913         {
914           std::ostringstream oss;
915           oss << "The structure of this and other differs at " << i << " !";
916           reason=oss.str();
917           return false;
918         }
919     }
920   return true;
921 }
922
923 void MEDCouplingIMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
924                                             DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const
925 {
926   if(!isEqualWithoutConsideringStr(other,prec))
927     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::checkDeepEquivalWith : Meshes are not the same !");
928 }
929
930 /*!
931  * Nothing is done here (except to check that the other is a ParaMEDMEM::MEDCouplingIMesh instance too).
932  * The user intend that the nodes are the same, so by construction of ParaMEDMEM::MEDCouplingIMesh, \a this and \a other are the same !
933  */
934 void MEDCouplingIMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
935                                                        DataArrayInt *&cellCor) const
936 {
937   if(!isEqualWithoutConsideringStr(other,prec))
938     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::checkDeepEquivalOnSameNodesWith : Meshes are not the same !");
939 }
940
941 void MEDCouplingIMesh::checkCoherency() const
942 {
943   checkSpaceDimension();
944   for(int i=0;i<_space_dim;i++)
945     if(_structure[i]<1)
946       {
947         std::ostringstream oss; oss << "MEDCouplingIMesh::checkCoherency : On axis " << i << "/" << _space_dim << ", number of nodes is equal to " << _structure[i] << " ! must be >=1 !";
948         throw INTERP_KERNEL::Exception(oss.str().c_str());
949       }
950 }
951
952 void MEDCouplingIMesh::checkCoherency1(double eps) const
953 {
954   checkCoherency();
955 }
956
957 void MEDCouplingIMesh::checkCoherency2(double eps) const
958 {
959   checkCoherency1(eps);
960 }
961
962 void MEDCouplingIMesh::getNodeGridStructure(int *res) const
963 {
964   checkSpaceDimension();
965   std::copy(_structure,_structure+_space_dim,res);
966 }
967
968 std::vector<int> MEDCouplingIMesh::getNodeGridStructure() const
969 {
970   checkSpaceDimension();
971   std::vector<int> ret(_structure,_structure+_space_dim);
972   return ret;
973 }
974
975 MEDCouplingStructuredMesh *MEDCouplingIMesh::buildStructuredSubPart(const std::vector< std::pair<int,int> >& cellPart) const
976 {
977   checkCoherency();
978   int dim(getSpaceDimension());
979   if(dim!=(int)cellPart.size())
980     {
981       std::ostringstream oss; oss << "MEDCouplingIMesh::buildStructuredSubPart : the space dimension is " << dim << " and cell part size is " << cellPart.size() << " !";
982       throw INTERP_KERNEL::Exception(oss.str().c_str());
983     }
984   double retOrigin[3]={0.,0.,0.};
985   int retStruct[3]={0,0,0};
986   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> ret(dynamic_cast<MEDCouplingIMesh *>(deepCpy()));
987   for(int i=0;i<dim;i++)
988     {
989       int startNode(cellPart[i].first),endNode(cellPart[i].second+1);
990       int myDelta(endNode-startNode);
991       if(startNode<0 || startNode>=_structure[i])
992         {
993           std::ostringstream oss; oss << "MEDCouplingIMesh::buildStructuredSubPart : At dimension #" << i << " the start node id is " << startNode << " it should be in [0," << _structure[i] << ") !";
994           throw INTERP_KERNEL::Exception(oss.str().c_str());
995         }
996       if(myDelta<0 || myDelta>_structure[i])
997         {
998           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;
999           throw INTERP_KERNEL::Exception(oss.str().c_str());
1000         }
1001       retOrigin[i]=_origin[i]+startNode*_dxyz[i];
1002       retStruct[i]=myDelta;
1003     }
1004   ret->setNodeStruct(retStruct,retStruct+dim);
1005   ret->setOrigin(retOrigin,retOrigin+dim);
1006   ret->checkCoherency();
1007   return ret.retn();
1008 }
1009
1010 /*!
1011  * Return the space dimension of \a this.
1012  */
1013 int MEDCouplingIMesh::getSpaceDimension() const
1014 {
1015   return _space_dim;
1016 }
1017
1018 void MEDCouplingIMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const
1019 {
1020   int tmp[3];
1021   int spaceDim(getSpaceDimension());
1022   getSplitNodeValues(tmp);
1023   int tmp2[3];
1024   GetPosFromId(nodeId,spaceDim,tmp,tmp2);
1025   for(int j=0;j<spaceDim;j++)
1026     coo.push_back(_origin[j]+_dxyz[j]*tmp2[j]);
1027 }
1028
1029 std::string MEDCouplingIMesh::simpleRepr() const
1030 {
1031   std::ostringstream ret;
1032   ret << "Image grid with name : \"" << getName() << "\"\n";
1033   ret << "Description of mesh : \"" << getDescription() << "\"\n";
1034   int tmpp1,tmpp2;
1035   double tt(getTime(tmpp1,tmpp2));
1036   int spaceDim(_space_dim);
1037   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
1038   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
1039   ret << "Space dimension : " << spaceDim << "\n";
1040   if(spaceDim<0 || spaceDim>3)
1041     return ret.str();
1042   ret << "The nodal structure is : "; std::copy(_structure,_structure+spaceDim,std::ostream_iterator<int>(ret," ")); ret << "\n";
1043   ret << "The origin position is [" << _axis_unit << "]: ";
1044   std::copy(_origin,_origin+spaceDim,std::ostream_iterator<double>(ret," ")); ret << "\n";
1045   ret << "The intervals along axis are : ";
1046   std::copy(_dxyz,_dxyz+spaceDim,std::ostream_iterator<double>(ret," ")); ret << "\n";
1047   return ret.str();
1048 }
1049
1050 std::string MEDCouplingIMesh::advancedRepr() const
1051 {
1052   return simpleRepr();
1053 }
1054
1055 void MEDCouplingIMesh::getBoundingBox(double *bbox) const
1056 {
1057   checkCoherency();
1058   int dim(getSpaceDimension());
1059   for(int idim=0; idim<dim; idim++)
1060     {
1061       bbox[2*idim]=_origin[idim];
1062       int coeff(_structure[idim]);
1063       if(_structure[idim]<0)
1064         {
1065           std::ostringstream oss; oss << "MEDCouplingIMesh::getBoundingBox : on axis #" << idim << " number of nodes in structure is < 0 !";
1066           throw INTERP_KERNEL::Exception(oss.str().c_str());
1067         }
1068       if(_structure[idim]>1)
1069         coeff=_structure[idim]-1;
1070       bbox[2*idim+1]=_origin[idim]+_dxyz[idim]*coeff;
1071     }
1072 }
1073
1074 /*!
1075  * Returns a new MEDCouplingFieldDouble containing volumes of cells constituting \a this
1076  * mesh.<br>
1077  * For 1D cells, the returned field contains lengths.<br>
1078  * For 2D cells, the returned field contains areas.<br>
1079  * For 3D cells, the returned field contains volumes.
1080  *  \param [in] isAbs - a not used parameter.
1081  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on cells
1082  *         and one time . The caller is to delete this field using decrRef() as it is no
1083  *         more needed.
1084  */
1085 MEDCouplingFieldDouble *MEDCouplingIMesh::getMeasureField(bool isAbs) const
1086 {
1087   checkCoherency();
1088   std::string name="MeasureOfMesh_";
1089   name+=getName();
1090   int nbelem(getNumberOfCells());
1091   MEDCouplingFieldDouble *field(MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME));
1092   field->setName(name);
1093   DataArrayDouble* array(DataArrayDouble::New());
1094   array->alloc(nbelem,1);
1095   array->fillWithValue(getMeasureOfAnyCell());
1096   field->setArray(array) ;
1097   array->decrRef();
1098   field->setMesh(const_cast<MEDCouplingIMesh *>(this));
1099   field->synchronizeTimeWithMesh();
1100   return field;
1101 }
1102
1103 /*!
1104  * not implemented yet !
1105  */
1106 MEDCouplingFieldDouble *MEDCouplingIMesh::getMeasureFieldOnNode(bool isAbs) const
1107 {
1108   throw INTERP_KERNEL::Exception("MEDCouplingIMesh::getMeasureFieldOnNode : not implemented yet !");
1109   //return 0;
1110 }
1111
1112 int MEDCouplingIMesh::getCellContainingPoint(const double *pos, double eps) const
1113 {
1114   int dim(getSpaceDimension()),ret(0),coeff(1);
1115   for(int i=0;i<dim;i++)
1116     {
1117       int nbOfCells(_structure[i]-1);
1118       double ref(pos[i]);
1119       int tmp((int)((ref-_origin[i])/_dxyz[i]));
1120       if(tmp>=0 && tmp<nbOfCells)
1121         {
1122           ret+=coeff*tmp;
1123           coeff*=nbOfCells;
1124         }
1125       else
1126         return -1;
1127     }
1128   return ret;
1129 }
1130
1131 void MEDCouplingIMesh::rotate(const double *center, const double *vector, double angle)
1132 {
1133   throw INTERP_KERNEL::Exception("No rotation available on IMesh : Traduce it to unstructured mesh to apply it !");
1134 }
1135
1136 /*!
1137  * Translates all nodes of \a this mesh by a given vector. Actually, it adds each
1138  * component of the \a vector to all node coordinates of a corresponding axis.
1139  *  \param [in] vector - the translation vector whose size must be not less than \a
1140  *         this->getSpaceDimension().
1141  */
1142 void MEDCouplingIMesh::translate(const double *vector)
1143 {
1144   checkSpaceDimension();
1145   int dim(getSpaceDimension());
1146   std::transform(_origin,_origin+dim,vector,_origin,std::plus<double>());
1147   declareAsNew();
1148 }
1149
1150 /*!
1151  * Applies scaling transformation to all nodes of \a this mesh.
1152  *  \param [in] point - coordinates of a scaling center. This array is to be of
1153  *         size \a this->getSpaceDimension() at least.
1154  *  \param [in] factor - a scale factor.
1155  */
1156 void MEDCouplingIMesh::scale(const double *point, double factor)
1157 {
1158   checkSpaceDimension();
1159   int dim(getSpaceDimension());
1160   std::transform(_origin,_origin+dim,point,_origin,std::minus<double>());
1161   std::transform(_origin,_origin+dim,_origin,std::bind2nd(std::multiplies<double>(),factor));
1162   std::transform(_dxyz,_dxyz+dim,_dxyz,std::bind2nd(std::multiplies<double>(),factor));
1163   std::transform(_origin,_origin+dim,point,_origin,std::plus<double>());
1164   declareAsNew();
1165 }
1166
1167 MEDCouplingMesh *MEDCouplingIMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
1168 {
1169   //not implemented yet !
1170   return 0;
1171 }
1172
1173 /*!
1174  * Returns a new DataArrayDouble holding coordinates of all nodes of \a this mesh.
1175  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
1176  *          this->getNumberOfNodes() tuples per \a this->getSpaceDimension()
1177  *          components. The caller is to delete this array using decrRef() as it is
1178  *          no more needed.
1179  */
1180 DataArrayDouble *MEDCouplingIMesh::getCoordinatesAndOwner() const
1181 {
1182   checkCoherency();
1183   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
1184   int spaceDim(getSpaceDimension()),nbNodes(getNumberOfNodes());
1185   ret->alloc(nbNodes,spaceDim);
1186   double *pt(ret->getPointer());
1187   ret->setInfoOnComponents(buildInfoOnComponents());
1188   int tmp2[3],tmp[3];
1189   getSplitNodeValues(tmp);
1190   for(int i=0;i<nbNodes;i++)
1191     {
1192       GetPosFromId(i,spaceDim,tmp,tmp2);
1193       for(int j=0;j<spaceDim;j++)
1194         pt[i*spaceDim+j]=_dxyz[j]*tmp2[j]+_origin[j];
1195     }
1196   return ret.retn();
1197 }
1198
1199 /*!
1200  * Returns a new DataArrayDouble holding barycenters of all cells. The barycenter is
1201  * computed by averaging coordinates of cell nodes.
1202  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
1203  *          this->getNumberOfCells() tuples per \a this->getSpaceDimension()
1204  *          components. The caller is to delete this array using decrRef() as it is
1205  *          no more needed.
1206  */
1207 DataArrayDouble *MEDCouplingIMesh::getBarycenterAndOwner() const
1208 {
1209   checkCoherency();
1210   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
1211   int spaceDim(getSpaceDimension()),nbCells(getNumberOfCells()),tmp[3],tmp2[3];
1212   ret->alloc(nbCells,spaceDim);
1213   double *pt(ret->getPointer()),shiftOrigin[3];
1214   std::transform(_dxyz,_dxyz+spaceDim,shiftOrigin,std::bind2nd(std::multiplies<double>(),0.5));
1215   std::transform(_origin,_origin+spaceDim,shiftOrigin,shiftOrigin,std::plus<double>());
1216   getSplitCellValues(tmp);
1217   ret->setInfoOnComponents(buildInfoOnComponents());
1218   for(int i=0;i<nbCells;i++)
1219     {
1220       GetPosFromId(i,spaceDim,tmp,tmp2);
1221       for(int j=0;j<spaceDim;j++)
1222         pt[i*spaceDim+j]=_dxyz[j]*tmp2[j]+shiftOrigin[j];
1223     }
1224   return ret.retn();
1225 }
1226
1227 DataArrayDouble *MEDCouplingIMesh::computeIsoBarycenterOfNodesPerCell() const
1228 {
1229   return MEDCouplingIMesh::getBarycenterAndOwner();
1230 }
1231
1232 void MEDCouplingIMesh::renumberCells(const int *old2NewBg, bool check)
1233 {
1234   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for IMesh !");
1235 }
1236
1237 void MEDCouplingIMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
1238 {
1239   int it,order;
1240   double time(getTime(it,order));
1241   tinyInfo.clear();
1242   tinyInfoD.clear();
1243   littleStrings.clear();
1244   littleStrings.push_back(getName());
1245   littleStrings.push_back(getDescription());
1246   littleStrings.push_back(getTimeUnit());
1247   littleStrings.push_back(getAxisUnit());
1248   tinyInfo.push_back(it);
1249   tinyInfo.push_back(order);
1250   tinyInfo.push_back(_space_dim);
1251   tinyInfo.insert(tinyInfo.end(),_structure,_structure+3);
1252   tinyInfoD.push_back(time);
1253   tinyInfoD.insert(tinyInfoD.end(),_dxyz,_dxyz+3);
1254   tinyInfoD.insert(tinyInfoD.end(),_origin,_origin+3);
1255 }
1256
1257 void MEDCouplingIMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
1258 {
1259   a1->alloc(0,1);
1260   a2->alloc(0,1);
1261 }
1262
1263 void MEDCouplingIMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
1264 {
1265   a1=DataArrayInt::New();
1266   a1->alloc(0,1);
1267   a2=DataArrayDouble::New();
1268   a2->alloc(0,1);
1269 }
1270
1271 void MEDCouplingIMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
1272                                        const std::vector<std::string>& littleStrings)
1273 {
1274   setName(littleStrings[0]);
1275   setDescription(littleStrings[1]);
1276   setTimeUnit(littleStrings[2]);
1277   setAxisUnit(littleStrings[3]);
1278   setTime(tinyInfoD[0],tinyInfo[0],tinyInfo[1]);
1279   _space_dim=tinyInfo[2];
1280   _structure[0]=tinyInfo[3]; _structure[1]=tinyInfo[4]; _structure[2]=tinyInfo[5];
1281   _dxyz[0]=tinyInfoD[1]; _dxyz[1]=tinyInfoD[2]; _dxyz[2]=tinyInfoD[3];
1282   _origin[0]=tinyInfoD[4]; _origin[1]=tinyInfoD[5]; _origin[2]=tinyInfoD[6];
1283   declareAsNew();
1284 }
1285
1286 void MEDCouplingIMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
1287 {
1288   checkCoherency();
1289   std::ostringstream extent,origin,spacing;
1290   for(int i=0;i<3;i++)
1291     {
1292       if(i<_space_dim)
1293         { extent << "0 " <<  _structure[i]-1 << " "; origin << _origin[i] << " "; spacing << _dxyz[i] << " "; }
1294       else
1295         { extent << "0 0 "; origin << "0 "; spacing << "0 "; }
1296     }
1297   ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\" Origin=\"" << origin.str() << "\" Spacing=\"" << spacing.str() << "\">\n";
1298   ofs << "    <Piece Extent=\"" << extent.str() << "\">\n";
1299   ofs << "      <PointData>\n" << pointData << std::endl;
1300   ofs << "      </PointData>\n";
1301   ofs << "      <CellData>\n" << cellData << std::endl;
1302   ofs << "      </CellData>\n";
1303   ofs << "      <Coordinates>\n";
1304   ofs << "      </Coordinates>\n";
1305   ofs << "    </Piece>\n";
1306   ofs << "  </" << getVTKDataSetType() << ">\n";
1307 }
1308
1309 void MEDCouplingIMesh::reprQuickOverview(std::ostream& stream) const
1310 {
1311   stream << "MEDCouplingIMesh C++ instance at " << this << ". Name : \"" << getName() << "\". Space dimension : " << _space_dim << ".";
1312   if(_space_dim<0 || _space_dim>3)
1313     return ;
1314   stream << "\n";
1315   std::ostringstream stream0,stream1;
1316   int nbNodes(1),nbCells(0);
1317   bool isPb(false);
1318   for(int i=0;i<_space_dim;i++)
1319     {
1320       char tmp('X'+i);
1321       int tmpNodes(_structure[i]);
1322       stream1 << "- Axis " << tmp << " : " << tmpNodes << " nodes (orig=" << _origin[i] << ", inter=" << _dxyz[i] << ").";
1323       if(i!=_space_dim-1)
1324         stream1 << std::endl;
1325       if(tmpNodes>=1)
1326         nbNodes*=tmpNodes;
1327       else
1328         isPb=true;
1329       if(tmpNodes>=2)
1330         nbCells=nbCells==0?tmpNodes-1:nbCells*(tmpNodes-1);
1331     }
1332   if(!isPb)
1333     {
1334       stream0 << "Number of cells : " << nbCells << ", Number of nodes : " << nbNodes;
1335       stream << stream0.str();
1336       if(_space_dim>0)
1337         stream << std::endl;
1338     }
1339   stream << stream1.str();
1340 }
1341
1342 std::string MEDCouplingIMesh::getVTKDataSetType() const
1343 {
1344   return std::string("ImageData");
1345 }
1346
1347 std::vector<std::string> MEDCouplingIMesh::buildInfoOnComponents() const
1348 {
1349   checkSpaceDimension();
1350   int dim(getSpaceDimension());
1351   std::vector<std::string> ret(dim);
1352   for(int i=0;i<dim;i++)
1353     {
1354       std::ostringstream oss;
1355       char tmp('X'+i); oss << tmp;
1356       ret[i]=DataArray::BuildInfoFromVarAndUnit(oss.str(),_axis_unit);
1357     }
1358   return ret;
1359 }
1360
1361 void MEDCouplingIMesh::checkSpaceDimension() const
1362 {
1363   CheckSpaceDimension(_space_dim);
1364 }
1365
1366 void MEDCouplingIMesh::CheckSpaceDimension(int spaceDim)
1367 {
1368   if(spaceDim<0 || spaceDim>3)
1369     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CheckSpaceDimension : input spaceDim must be in [0,1,2,3] !");
1370 }
1371
1372 int MEDCouplingIMesh::FindIntRoot(int val, int order)
1373 {
1374   if(order==0)
1375     return 1;
1376   if(val<0)
1377     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : input val is < 0 ! Not possible to compute a root !");
1378   if(order==1)
1379     return val;
1380   if(order!=2 && order!=3)
1381     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the order available are 0,1,2 or 3 !");
1382   double valf((double)val);
1383   if(order==2)
1384     {
1385       double retf(sqrt(valf));
1386       int ret((int)retf);
1387       if(ret*ret!=val)
1388         throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect square root !");
1389       return ret;
1390     }
1391   else//order==3
1392     {
1393       double retf(std::pow(val,0.3333333333333333));
1394       int ret((int)retf),ret2(ret+1);
1395       if(ret*ret*ret!=val && ret2*ret2*ret2!=val)
1396         throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect cublic root !");
1397       if(ret*ret*ret==val)
1398         return ret;
1399       else
1400         return ret2;
1401     }
1402 }