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