Salome HOME
Refinement per axis + Spread method
[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 void MEDCouplingIMesh::setNodeStruct(const int *nodeStrctStart, const int *nodeStrctStop)
79 {
80   checkSpaceDimension();
81   int sz((int)std::distance(nodeStrctStart,nodeStrctStop));
82   if(sz!=_space_dim)
83     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::setNodeStruct : input vector of node structure has not the right size ! Or change space dimension before calling it !");
84   std::copy(nodeStrctStart,nodeStrctStop,_structure);
85   declareAsNew();
86 }
87
88 std::vector<int> MEDCouplingIMesh::getNodeStruct() const
89 {
90   checkSpaceDimension();
91   return std::vector<int>(_structure,_structure+_space_dim);
92 }
93
94 void MEDCouplingIMesh::setOrigin(const double *originStart, const double *originStop)
95 {
96   checkSpaceDimension();
97   int sz((int)std::distance(originStart,originStop));
98   if(sz!=_space_dim)
99     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::setOrigin : input vector of origin vector has not the right size ! Or change space dimension before calling it !");
100   std::copy(originStart,originStop,_origin);
101   declareAsNew();
102 }
103
104 std::vector<double> MEDCouplingIMesh::getOrigin() const
105 {
106   checkSpaceDimension();
107   return std::vector<double>(_origin,_origin+_space_dim);
108 }
109
110 void MEDCouplingIMesh::setDXYZ(const double *dxyzStart, const double *dxyzStop)
111 {
112   checkSpaceDimension();
113   int sz((int)std::distance(dxyzStart,dxyzStop));
114   if(sz!=_space_dim)
115     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::setDXYZ : input vector of dxyz vector has not the right size ! Or change space dimension before calling it !");
116   std::copy(dxyzStart,dxyzStop,_dxyz);
117   declareAsNew();
118 }
119
120 std::vector<double> MEDCouplingIMesh::getDXYZ() const
121 {
122   checkSpaceDimension();
123   return std::vector<double>(_dxyz,_dxyz+_space_dim);
124 }
125
126 void MEDCouplingIMesh::setAxisUnit(const std::string& unitName)
127 {
128   _axis_unit=unitName;
129   declareAsNew();
130 }
131
132 std::string MEDCouplingIMesh::getAxisUnit() const
133 {
134   return _axis_unit;
135 }
136
137 /*!
138  * This method returns the measure of any cell in \a this.
139  * This specific method of image grid mesh utilizes the fact that any cell in \a this have the same measure.
140  * The value returned by this method is those used to feed the returned field in the MEDCouplingIMesh::getMeasureField.
141  *
142  * \sa getMeasureField
143  */
144 double MEDCouplingIMesh::getMeasureOfAnyCell() const
145 {
146   checkCoherency();
147   int dim(getSpaceDimension());
148   double ret(1.);
149   for(int i=0;i<dim;i++)
150     ret*=fabs(_dxyz[i]);
151   return ret;
152 }
153
154 /*!
155  * This method is allows to convert \a this into MEDCouplingCMesh instance.
156  * This method is the middle level between MEDCouplingIMesh and the most general MEDCouplingUMesh.
157  * This method is useful for MED writers that do not have still the image grid support.
158  *
159  * \sa MEDCouplingMesh::buildUnstructured
160  */
161 MEDCouplingCMesh *MEDCouplingIMesh::convertToCartesian() const
162 {
163   checkCoherency();
164   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCMesh> ret(MEDCouplingCMesh::New());
165   try
166   { ret->copyTinyInfoFrom(this); }
167   catch(INTERP_KERNEL::Exception& ) { }
168   int spaceDim(getSpaceDimension());
169   std::vector<std::string> infos(buildInfoOnComponents());
170   for(int i=0;i<spaceDim;i++)
171     {
172       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr(DataArrayDouble::New()); arr->alloc(_structure[i],1); arr->setInfoOnComponent(0,infos[i]);
173       arr->iota(); arr->applyLin(_dxyz[i],_origin[i]);
174       ret->setCoordsAt(i,arr);
175     }
176   return ret.retn();
177 }
178
179 /*!
180  * This method refines \a this uniformaly along all of its dimensions. In case of success the space covered by \a this will remain
181  * the same before the invocation except that the number of cells will be multiplied by \a factor ^ this->getMeshDimension().
182  * The origin of \a this will be not touched only spacing and node structure will be changed.
183  * This method can be useful for AMR users.
184  */
185 void MEDCouplingIMesh::refineWithFactor(int factor)
186 {
187   if(factor==0)
188     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::refineWithFactor : refinement factor must be != 0 !");
189   checkCoherency();
190   int factAbs(std::abs(factor));
191   double fact2(1./(double)factor);
192   std::transform(_structure,_structure+_space_dim,_structure,std::bind2nd(std::plus<int>(),-1));
193   std::transform(_structure,_structure+_space_dim,_structure,std::bind2nd(std::multiplies<int>(),factAbs));
194   std::transform(_structure,_structure+_space_dim,_structure,std::bind2nd(std::plus<int>(),1));
195   std::transform(_dxyz,_dxyz+_space_dim,_dxyz,std::bind2nd(std::multiplies<double>(),fact2));
196   declareAsNew();
197 }
198
199 /*!
200  * This static method is useful to condense field on cells of a MEDCouplingIMesh instance coming from a refinement ( MEDCouplingIMesh::refineWithFactor for example)
201  * to a coarse MEDCouplingIMesh instance. So this method can be seen as a specialization in P0P0 conservative interpolation non overlaping from fine image mesh
202  * 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.
203  *
204  * \param [in,out] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
205  * \param [in] coarseSt The cell structure of coarse mesh.
206  * \param [in] fineDA The DataArray containing the cell field on uniformly refined mesh
207  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
208  * \param [in] facts The refinement coefficient per axis.
209  * \sa SpreadCoarseToFine
210  */
211 void MEDCouplingIMesh::CondenseFineToCoarse(DataArrayDouble *coarseDA, const std::vector<int>& coarseSt, const DataArrayDouble *fineDA, const std::vector< std::pair<int,int> >& fineLocInCoarse, const std::vector<int>& facts)
212 {
213   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
214     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the parameters 1 or 3 are NULL or not allocated !");
215   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseSt)),nbOfTuplesInFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(fineLocInCoarse));
216   int nbCompo(fineDA->getNumberOfComponents());
217   if(coarseDA->getNumberOfComponents()!=nbCompo)
218     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the number of components of fine DA and coarse one mismatches !");
219   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
220     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) !");
221   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
222     {
223       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarse : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
224       throw INTERP_KERNEL::Exception(oss.str().c_str());
225     }
226   int nbTuplesFine(fineDA->getNumberOfTuples());
227   if(nbTuplesFine%nbOfTuplesInFineExp!=0)
228     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : Invalid nb of tuples in fine DataArray regarding its structure !");
229   int fact(std::accumulate(facts.begin(),facts.end(),1,std::multiplies<int>()));
230   if(nbTuplesFine!=fact*nbOfTuplesInFineExp)
231     {
232       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarse : Invalid number of tuples ("  << nbTuplesFine << ") of fine dataarray is invalid ! Must be " << fact*nbOfTuplesInFineExp << "!";
233       throw INTERP_KERNEL::Exception(oss.str().c_str());
234     }
235   // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
236   double *outPtr(coarseDA->getPointer());
237   const double *inPtr(fineDA->begin());
238   //
239   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
240   switch(meshDim)
241   {
242     case 1:
243       {
244         int offset(fineLocInCoarse[0].first),fact0(facts[0]);
245         for(int i=0;i<dims[0];i++)
246           {
247             double *loc(outPtr+(offset+i)*nbCompo);
248             for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
249               {
250                 if(ifact!=0)
251                   std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
252                 else
253                   std::copy(inPtr,inPtr+nbCompo,loc);
254               }
255           }
256         break;
257       }
258     case 2:
259       {
260         int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first),fact1(facts[1]),fact0(facts[0]);
261         for(int j=0;j<dims[1];j++)
262           {
263             for(int jfact=0;jfact<fact1;jfact++)
264               {
265                 for(int i=0;i<dims[0];i++)
266                   {
267                     double *loc(outPtr+(kk+i)*nbCompo);
268                     for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
269                       {
270                         if(jfact!=0 || ifact!=0)
271                           std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
272                         else
273                           std::copy(inPtr,inPtr+nbCompo,loc);
274                       }
275                   }
276               }
277             kk+=coarseSt[0];
278           }
279         break;
280       }
281     case 3:
282       {
283         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]);
284         for(int k=0;k<dims[2];k++)
285           {
286             for(int kfact=0;kfact<fact2;kfact++)
287               {
288                 for(int j=0;j<dims[1];j++)
289                   {
290                     for(int jfact=0;jfact<fact1;jfact++)
291                       {
292                         for(int i=0;i<dims[0];i++)
293                           {
294                             double *loc(outPtr+(kk+i+j*coarseSt[0])*nbCompo);
295                             for(int ifact=0;ifact<fact0;ifact++,inPtr+=nbCompo)
296                               {
297                                 if(kfact!=0 || jfact!=0 || ifact!=0)
298                                   std::transform(inPtr,inPtr+nbCompo,loc,loc,std::plus<double>());
299                                 else
300                                   std::copy(inPtr,inPtr+nbCompo,loc);
301                               }
302                           }
303                       }
304                   }
305               }
306             kk+=coarseSt[0]*coarseSt[1];
307           }
308         break;
309       }
310     default:
311       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : only dimensions 1, 2 and 3 supported !");
312   }
313 }
314
315 /*!
316  * This method spreads the values of coarse data \a coarseDA into \a fineDA.
317  *
318  * \param [in] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
319  * \param [in] coarseSt The cell structure of coarse mesh.
320  * \param [in,out] fineDA The DataArray containing the cell field on uniformly refined mesh
321  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
322  * \param [in] facts The refinement coefficient per axis.
323  * \sa CondenseFineToCoarse
324  */
325 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)
326 {
327   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
328     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the parameters 1 or 3 are NULL or not allocated !");
329   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseSt)),nbOfTuplesInFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(fineLocInCoarse));
330   int nbCompo(fineDA->getNumberOfComponents());
331   if(coarseDA->getNumberOfComponents()!=nbCompo)
332     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : the number of components of fine DA and coarse one mismatches !");
333   if(meshDim!=(int)fineLocInCoarse.size() || meshDim!=(int)facts.size())
334     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) !");
335   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
336     {
337       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFine : Expecting " << nbOfTuplesInCoarseExp << " tuples having " << coarseDA->getNumberOfTuples() << " !";
338       throw INTERP_KERNEL::Exception(oss.str().c_str());
339     }
340   int nbTuplesFine(fineDA->getNumberOfTuples());
341   if(nbTuplesFine%nbOfTuplesInFineExp!=0)
342     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : Invalid nb of tuples in fine DataArray regarding its structure !");
343   int fact(std::accumulate(facts.begin(),facts.end(),1,std::multiplies<int>()));
344   if(nbTuplesFine!=fact*nbOfTuplesInFineExp)
345     {
346       std::ostringstream oss; oss << "MEDCouplingIMesh::SpreadCoarseToFine : Invalid number of tuples ("  << nbTuplesFine << ") of fine dataarray is invalid ! Must be " << fact*nbOfTuplesInFineExp << "!";
347       throw INTERP_KERNEL::Exception(oss.str().c_str());
348     }
349   // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
350   double *outPtr(fineDA->getPointer());
351   const double *inPtr(coarseDA->begin());
352   //
353   std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
354   switch(meshDim)
355   {
356     case 1:
357       {
358         int offset(fineLocInCoarse[0].first),fact0(facts[0]);
359         for(int i=0;i<dims[0];i++)
360           {
361             const double *loc(inPtr+(offset+i)*nbCompo);
362             for(int ifact=0;ifact<fact0;ifact++)
363               outPtr=std::copy(loc,loc+nbCompo,outPtr);
364           }
365         break;
366       }
367     case 2:
368       {
369         int kk(fineLocInCoarse[0].first+coarseSt[0]*fineLocInCoarse[1].first),fact0(facts[0]),fact1(facts[1]);
370         for(int j=0;j<dims[1];j++)
371           {
372             for(int jfact=0;jfact<fact1;jfact++)
373               {
374                 for(int i=0;i<dims[0];i++)
375                   {
376                     const double *loc(inPtr+(kk+i)*nbCompo);
377                     for(int ifact=0;ifact<fact0;ifact++)
378                       outPtr=std::copy(loc,loc+nbCompo,outPtr);
379                   }
380               }
381             kk+=coarseSt[0];
382           }
383         break;
384       }
385     case 3:
386       {
387         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]);
388         for(int k=0;k<dims[2];k++)
389           {
390             for(int kfact=0;kfact<fact2;kfact++)
391               {
392                 for(int j=0;j<dims[1];j++)
393                   {
394                     for(int jfact=0;jfact<fact1;jfact++)
395                       {
396                         for(int i=0;i<dims[0];i++)
397                           {
398                             const double *loc(inPtr+(kk+i+j*coarseSt[0])*nbCompo);
399                             for(int ifact=0;ifact<fact0;ifact++)
400                               outPtr=std::copy(loc,loc+nbCompo,outPtr);
401                           }
402                       }
403                   }
404               }
405             kk+=coarseSt[0]*coarseSt[1];
406           }
407         break;
408       }
409     default:
410       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::SpreadCoarseToFine : only dimensions 1, 2 and 3 supported !");
411   }
412 }
413
414 void MEDCouplingIMesh::setSpaceDimension(int spaceDim)
415 {
416   if(spaceDim==_space_dim)
417     return ;
418   CheckSpaceDimension(spaceDim);
419   _space_dim=spaceDim;
420   declareAsNew();
421 }
422
423 void MEDCouplingIMesh::updateTime() const
424 {
425 }
426
427 std::size_t MEDCouplingIMesh::getHeapMemorySizeWithoutChildren() const
428 {
429   return MEDCouplingStructuredMesh::getHeapMemorySizeWithoutChildren();
430 }
431
432 std::vector<const BigMemoryObject *> MEDCouplingIMesh::getDirectChildren() const
433 {
434   return std::vector<const BigMemoryObject *>();
435 }
436
437 /*!
438  * This method copyies all tiny strings from other (name and components name).
439  * @throw if other and this have not same mesh type.
440  */
441 void MEDCouplingIMesh::copyTinyStringsFrom(const MEDCouplingMesh *other)
442
443   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
444   if(!otherC)
445     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::copyTinyStringsFrom : meshes have not same type !");
446   MEDCouplingStructuredMesh::copyTinyStringsFrom(other);
447   declareAsNew();
448 }
449
450 bool MEDCouplingIMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
451 {
452   if(!other)
453     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::isEqualIfNotWhy : input other pointer is null !");
454   const MEDCouplingIMesh *otherC(dynamic_cast<const MEDCouplingIMesh *>(other));
455   if(!otherC)
456     {
457       reason="mesh given in input is not castable in MEDCouplingIMesh !";
458       return false;
459     }
460   if(!MEDCouplingStructuredMesh::isEqualIfNotWhy(other,prec,reason))
461     return false;
462   if(!isEqualWithoutConsideringStrInternal(otherC,prec,reason))
463     return false;
464   if(_axis_unit!=otherC->_axis_unit)
465     {
466       reason="The units of axis are not the same !";
467       return false;
468     }
469   return true;
470 }
471
472 bool MEDCouplingIMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
473 {
474   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
475   if(!otherC)
476     return false;
477   std::string tmp;
478   return isEqualWithoutConsideringStrInternal(other,prec,tmp);
479 }
480
481 bool MEDCouplingIMesh::isEqualWithoutConsideringStrInternal(const MEDCouplingMesh *other, double prec, std::string& reason) const
482 {
483   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
484   if(!otherC)
485     return false;
486   if(_space_dim!=otherC->_space_dim)
487     {
488       std::ostringstream oss;
489       oss << "The spaceDimension of this (" << _space_dim << ") is not equal to those of other (" << otherC->_space_dim << ") !";
490       return false;
491     }
492   checkSpaceDimension();
493   for(int i=0;i<_space_dim;i++)
494     {
495       if(fabs(_origin[i]-otherC->_origin[i])>prec)
496         {
497           std::ostringstream oss;
498           oss << "The origin of this and other differs at " << i << " !";
499           reason=oss.str();
500           return false;
501         }
502     }
503   for(int i=0;i<_space_dim;i++)
504     {
505       if(fabs(_dxyz[i]-otherC->_dxyz[i])>prec)
506         {
507           std::ostringstream oss;
508           oss << "The delta of this and other differs at " << i << " !";
509           reason=oss.str();
510           return false;
511         }
512     }
513   for(int i=0;i<_space_dim;i++)
514     {
515       if(_structure[i]!=otherC->_structure[i])
516         {
517           std::ostringstream oss;
518           oss << "The structure of this and other differs at " << i << " !";
519           reason=oss.str();
520           return false;
521         }
522     }
523   return true;
524 }
525
526 void MEDCouplingIMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
527                                             DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const
528 {
529   if(!isEqualWithoutConsideringStr(other,prec))
530     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::checkDeepEquivalWith : Meshes are not the same !");
531 }
532
533 /*!
534  * Nothing is done here (except to check that the other is a ParaMEDMEM::MEDCouplingIMesh instance too).
535  * The user intend that the nodes are the same, so by construction of ParaMEDMEM::MEDCouplingIMesh, \a this and \a other are the same !
536  */
537 void MEDCouplingIMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
538                                                        DataArrayInt *&cellCor) const
539 {
540   if(!isEqualWithoutConsideringStr(other,prec))
541     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::checkDeepEquivalOnSameNodesWith : Meshes are not the same !");
542 }
543
544 void MEDCouplingIMesh::checkCoherency() const
545 {
546   checkSpaceDimension();
547   for(int i=0;i<_space_dim;i++)
548     if(_structure[i]<1)
549       {
550         std::ostringstream oss; oss << "MEDCouplingIMesh::checkCoherency : On axis " << i << "/" << _space_dim << ", number of nodes is equal to " << _structure[i] << " ! must be >=1 !";
551         throw INTERP_KERNEL::Exception(oss.str().c_str());
552       }
553 }
554
555 void MEDCouplingIMesh::checkCoherency1(double eps) const
556 {
557   checkCoherency();
558 }
559
560 void MEDCouplingIMesh::checkCoherency2(double eps) const
561 {
562   checkCoherency1(eps);
563 }
564
565 void MEDCouplingIMesh::getNodeGridStructure(int *res) const
566 {
567   checkSpaceDimension();
568   std::copy(_structure,_structure+_space_dim,res);
569 }
570
571 std::vector<int> MEDCouplingIMesh::getNodeGridStructure() const
572 {
573   checkSpaceDimension();
574   std::vector<int> ret(_structure,_structure+_space_dim);
575   return ret;
576 }
577
578 MEDCouplingStructuredMesh *MEDCouplingIMesh::buildStructuredSubPart(const std::vector< std::pair<int,int> >& cellPart) const
579 {
580   checkCoherency();
581   int dim(getSpaceDimension());
582   if(dim!=(int)cellPart.size())
583     {
584       std::ostringstream oss; oss << "MEDCouplingIMesh::buildStructuredSubPart : the space dimension is " << dim << " and cell part size is " << cellPart.size() << " !";
585       throw INTERP_KERNEL::Exception(oss.str().c_str());
586     }
587   double retOrigin[3]={0.,0.,0.};
588   int retStruct[3]={0,0,0};
589   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> ret(dynamic_cast<MEDCouplingIMesh *>(deepCpy()));
590   for(int i=0;i<dim;i++)
591     {
592       int startNode(cellPart[i].first),endNode(cellPart[i].second+1);
593       int myDelta(endNode-startNode);
594       if(startNode<0 || startNode>=_structure[i])
595         {
596           std::ostringstream oss; oss << "MEDCouplingIMesh::buildStructuredSubPart : At dimension #" << i << " the start node id is " << startNode << " it should be in [0," << _structure[i] << ") !";
597           throw INTERP_KERNEL::Exception(oss.str().c_str());
598         }
599       if(myDelta<0 || myDelta>_structure[i])
600         {
601           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;
602           throw INTERP_KERNEL::Exception(oss.str().c_str());
603         }
604       retOrigin[i]=_origin[i]+startNode*_dxyz[i];
605       retStruct[i]=myDelta;
606     }
607   ret->setNodeStruct(retStruct,retStruct+dim);
608   ret->setOrigin(retOrigin,retOrigin+dim);
609   ret->checkCoherency();
610   return ret.retn();
611 }
612
613 /*!
614  * Return the space dimension of \a this.
615  */
616 int MEDCouplingIMesh::getSpaceDimension() const
617 {
618   return _space_dim;
619 }
620
621 void MEDCouplingIMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const
622 {
623   int tmp[3];
624   int spaceDim(getSpaceDimension());
625   getSplitNodeValues(tmp);
626   int tmp2[3];
627   GetPosFromId(nodeId,spaceDim,tmp,tmp2);
628   for(int j=0;j<spaceDim;j++)
629     coo.push_back(_origin[j]+_dxyz[j]*tmp2[j]);
630 }
631
632 std::string MEDCouplingIMesh::simpleRepr() const
633 {
634   std::ostringstream ret;
635   ret << "Image grid with name : \"" << getName() << "\"\n";
636   ret << "Description of mesh : \"" << getDescription() << "\"\n";
637   int tmpp1,tmpp2;
638   double tt(getTime(tmpp1,tmpp2));
639   int spaceDim(_space_dim);
640   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
641   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
642   ret << "Space dimension : " << spaceDim << "\n";
643   if(spaceDim<0 || spaceDim>3)
644     return ret.str();
645   ret << "The nodal structure is : "; std::copy(_structure,_structure+spaceDim,std::ostream_iterator<int>(ret," ")); ret << "\n";
646   ret << "The origin position is [" << _axis_unit << "]: ";
647   std::copy(_origin,_origin+spaceDim,std::ostream_iterator<double>(ret," ")); ret << "\n";
648   ret << "The intervals along axis are : ";
649   std::copy(_dxyz,_dxyz+spaceDim,std::ostream_iterator<double>(ret," ")); ret << "\n";
650   return ret.str();
651 }
652
653 std::string MEDCouplingIMesh::advancedRepr() const
654 {
655   return simpleRepr();
656 }
657
658 void MEDCouplingIMesh::getBoundingBox(double *bbox) const
659 {
660   checkCoherency();
661   int dim(getSpaceDimension());
662   for(int idim=0; idim<dim; idim++)
663     {
664       bbox[2*idim]=_origin[idim];
665       bbox[2*idim+1]=_origin[idim]+_dxyz[idim]*_structure[idim];
666     }
667 }
668
669 /*!
670  * Returns a new MEDCouplingFieldDouble containing volumes of cells constituting \a this
671  * mesh.<br>
672  * For 1D cells, the returned field contains lengths.<br>
673  * For 2D cells, the returned field contains areas.<br>
674  * For 3D cells, the returned field contains volumes.
675  *  \param [in] isAbs - a not used parameter.
676  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on cells
677  *         and one time . The caller is to delete this field using decrRef() as it is no
678  *         more needed.
679  */
680 MEDCouplingFieldDouble *MEDCouplingIMesh::getMeasureField(bool isAbs) const
681 {
682   checkCoherency();
683   std::string name="MeasureOfMesh_";
684   name+=getName();
685   int nbelem(getNumberOfCells());
686   MEDCouplingFieldDouble *field(MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME));
687   field->setName(name);
688   DataArrayDouble* array(DataArrayDouble::New());
689   array->alloc(nbelem,1);
690   array->fillWithValue(getMeasureOfAnyCell());
691   field->setArray(array) ;
692   array->decrRef();
693   field->setMesh(const_cast<MEDCouplingIMesh *>(this));
694   field->synchronizeTimeWithMesh();
695   return field;
696 }
697
698 /*!
699  * not implemented yet !
700  */
701 MEDCouplingFieldDouble *MEDCouplingIMesh::getMeasureFieldOnNode(bool isAbs) const
702 {
703   throw INTERP_KERNEL::Exception("MEDCouplingIMesh::getMeasureFieldOnNode : not implemented yet !");
704   //return 0;
705 }
706
707 int MEDCouplingIMesh::getCellContainingPoint(const double *pos, double eps) const
708 {
709   int dim(getSpaceDimension()),ret(0),coeff(1);
710   for(int i=0;i<dim;i++)
711     {
712       int nbOfCells(_structure[i]-1);
713       double ref(pos[i]);
714       int tmp((int)((ref-_origin[i])/_dxyz[i]));
715       if(tmp>=0 && tmp<nbOfCells)
716         {
717           ret+=coeff*tmp;
718           coeff*=nbOfCells;
719         }
720       else
721         return -1;
722     }
723   return ret;
724 }
725
726 void MEDCouplingIMesh::rotate(const double *center, const double *vector, double angle)
727 {
728   throw INTERP_KERNEL::Exception("No rotation available on IMesh : Traduce it to unstructured mesh to apply it !");
729 }
730
731 /*!
732  * Translates all nodes of \a this mesh by a given vector. Actually, it adds each
733  * component of the \a vector to all node coordinates of a corresponding axis.
734  *  \param [in] vector - the translation vector whose size must be not less than \a
735  *         this->getSpaceDimension().
736  */
737 void MEDCouplingIMesh::translate(const double *vector)
738 {
739   checkSpaceDimension();
740   int dim(getSpaceDimension());
741   std::transform(_origin,_origin+dim,vector,_origin,std::plus<double>());
742   declareAsNew();
743 }
744
745 /*!
746  * Applies scaling transformation to all nodes of \a this mesh.
747  *  \param [in] point - coordinates of a scaling center. This array is to be of
748  *         size \a this->getSpaceDimension() at least.
749  *  \param [in] factor - a scale factor.
750  */
751 void MEDCouplingIMesh::scale(const double *point, double factor)
752 {
753   checkSpaceDimension();
754   int dim(getSpaceDimension());
755   std::transform(_origin,_origin+dim,point,_origin,std::minus<double>());
756   std::transform(_origin,_origin+dim,_origin,std::bind2nd(std::multiplies<double>(),factor));
757   std::transform(_dxyz,_dxyz+dim,_dxyz,std::bind2nd(std::multiplies<double>(),factor));
758   std::transform(_origin,_origin+dim,point,_origin,std::plus<double>());
759   declareAsNew();
760 }
761
762 MEDCouplingMesh *MEDCouplingIMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
763 {
764   //not implemented yet !
765   return 0;
766 }
767
768 /*!
769  * Returns a new DataArrayDouble holding coordinates of all nodes of \a this mesh.
770  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
771  *          this->getNumberOfNodes() tuples per \a this->getSpaceDimension()
772  *          components. The caller is to delete this array using decrRef() as it is
773  *          no more needed.
774  */
775 DataArrayDouble *MEDCouplingIMesh::getCoordinatesAndOwner() const
776 {
777   checkCoherency();
778   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
779   int spaceDim(getSpaceDimension()),nbNodes(getNumberOfNodes());
780   ret->alloc(nbNodes,spaceDim);
781   double *pt(ret->getPointer());
782   ret->setInfoOnComponents(buildInfoOnComponents());
783   int tmp2[3],tmp[3];
784   getSplitNodeValues(tmp);
785   for(int i=0;i<nbNodes;i++)
786     {
787       GetPosFromId(i,spaceDim,tmp,tmp2);
788       for(int j=0;j<spaceDim;j++)
789         pt[i*spaceDim+j]=_dxyz[j]*tmp2[j]+_origin[j];
790     }
791   return ret.retn();
792 }
793
794 /*!
795  * Returns a new DataArrayDouble holding barycenters of all cells. The barycenter is
796  * computed by averaging coordinates of cell nodes.
797  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
798  *          this->getNumberOfCells() tuples per \a this->getSpaceDimension()
799  *          components. The caller is to delete this array using decrRef() as it is
800  *          no more needed.
801  */
802 DataArrayDouble *MEDCouplingIMesh::getBarycenterAndOwner() const
803 {
804   checkCoherency();
805   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
806   int spaceDim(getSpaceDimension()),nbCells(getNumberOfCells()),tmp[3],tmp2[3];
807   ret->alloc(nbCells,spaceDim);
808   double *pt(ret->getPointer()),shiftOrigin[3];
809   std::transform(_dxyz,_dxyz+spaceDim,shiftOrigin,std::bind2nd(std::multiplies<double>(),0.5));
810   std::transform(_origin,_origin+spaceDim,shiftOrigin,shiftOrigin,std::plus<double>());
811   getSplitCellValues(tmp);
812   ret->setInfoOnComponents(buildInfoOnComponents());
813   for(int i=0;i<nbCells;i++)
814     {
815       GetPosFromId(i,spaceDim,tmp,tmp2);
816       for(int j=0;j<spaceDim;j++)
817         pt[i*spaceDim+j]=_dxyz[j]*tmp2[j]+shiftOrigin[j];
818     }
819   return ret.retn();
820 }
821
822 DataArrayDouble *MEDCouplingIMesh::computeIsoBarycenterOfNodesPerCell() const
823 {
824   return MEDCouplingIMesh::getBarycenterAndOwner();
825 }
826
827 void MEDCouplingIMesh::renumberCells(const int *old2NewBg, bool check)
828 {
829   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for IMesh !");
830 }
831
832 void MEDCouplingIMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
833 {
834   int it,order;
835   double time(getTime(it,order));
836   tinyInfo.clear();
837   tinyInfoD.clear();
838   littleStrings.clear();
839   littleStrings.push_back(getName());
840   littleStrings.push_back(getDescription());
841   littleStrings.push_back(getTimeUnit());
842   littleStrings.push_back(getAxisUnit());
843   tinyInfo.push_back(it);
844   tinyInfo.push_back(order);
845   tinyInfo.push_back(_space_dim);
846   tinyInfo.insert(tinyInfo.end(),_structure,_structure+3);
847   tinyInfoD.push_back(time);
848   tinyInfoD.insert(tinyInfoD.end(),_dxyz,_dxyz+3);
849   tinyInfoD.insert(tinyInfoD.end(),_origin,_origin+3);
850 }
851
852 void MEDCouplingIMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
853 {
854   a1->alloc(0,1);
855   a2->alloc(0,1);
856 }
857
858 void MEDCouplingIMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
859 {
860   a1=DataArrayInt::New();
861   a1->alloc(0,1);
862   a2=DataArrayDouble::New();
863   a2->alloc(0,1);
864 }
865
866 void MEDCouplingIMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
867                                        const std::vector<std::string>& littleStrings)
868 {
869   setName(littleStrings[0]);
870   setDescription(littleStrings[1]);
871   setTimeUnit(littleStrings[2]);
872   setAxisUnit(littleStrings[3]);
873   setTime(tinyInfoD[0],tinyInfo[0],tinyInfo[1]);
874   _space_dim=tinyInfo[2];
875   _structure[0]=tinyInfo[3]; _structure[1]=tinyInfo[4]; _structure[2]=tinyInfo[5];
876   _dxyz[0]=tinyInfoD[1]; _dxyz[1]=tinyInfoD[2]; _dxyz[2]=tinyInfoD[3];
877   _origin[0]=tinyInfoD[4]; _origin[1]=tinyInfoD[5]; _origin[2]=tinyInfoD[6];
878   declareAsNew();
879 }
880
881 void MEDCouplingIMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
882 {
883   checkCoherency();
884   std::ostringstream extent,origin,spacing;
885   for(int i=0;i<3;i++)
886     {
887       if(i<_space_dim)
888         { extent << "0 " <<  _structure[i]-1 << " "; origin << _origin[i] << " "; spacing << _dxyz[i] << " "; }
889       else
890         { extent << "0 0 "; origin << "0 "; spacing << "0 "; }
891     }
892   ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\" Origin=\"" << origin.str() << "\" Spacing=\"" << spacing.str() << "\">\n";
893   ofs << "    <Piece Extent=\"" << extent.str() << "\">\n";
894   ofs << "      <PointData>\n" << pointData << std::endl;
895   ofs << "      </PointData>\n";
896   ofs << "      <CellData>\n" << cellData << std::endl;
897   ofs << "      </CellData>\n";
898   ofs << "      <Coordinates>\n";
899   ofs << "      </Coordinates>\n";
900   ofs << "    </Piece>\n";
901   ofs << "  </" << getVTKDataSetType() << ">\n";
902 }
903
904 void MEDCouplingIMesh::reprQuickOverview(std::ostream& stream) const
905 {
906   stream << "MEDCouplingIMesh C++ instance at " << this << ". Name : \"" << getName() << "\". Space dimension : " << _space_dim << ".";
907   if(_space_dim<0 || _space_dim>3)
908     return ;
909   stream << "\n";
910   std::ostringstream stream0,stream1;
911   int nbNodes(1),nbCells(0);
912   bool isPb(false);
913   for(int i=0;i<_space_dim;i++)
914     {
915       char tmp('X'+i);
916       int tmpNodes(_structure[i]);
917       stream1 << "- Axis " << tmp << " : " << tmpNodes << " nodes (orig=" << _origin[i] << ", inter=" << _dxyz[i] << ").";
918       if(i!=_space_dim-1)
919         stream1 << std::endl;
920       if(tmpNodes>=1)
921         nbNodes*=tmpNodes;
922       else
923         isPb=true;
924       if(tmpNodes>=2)
925         nbCells=nbCells==0?tmpNodes-1:nbCells*(tmpNodes-1);
926     }
927   if(!isPb)
928     {
929       stream0 << "Number of cells : " << nbCells << ", Number of nodes : " << nbNodes;
930       stream << stream0.str();
931       if(_space_dim>0)
932         stream << std::endl;
933     }
934   stream << stream1.str();
935 }
936
937 std::string MEDCouplingIMesh::getVTKDataSetType() const
938 {
939   return std::string("ImageData");
940 }
941
942 std::vector<std::string> MEDCouplingIMesh::buildInfoOnComponents() const
943 {
944   checkSpaceDimension();
945   int dim(getSpaceDimension());
946   std::vector<std::string> ret(dim);
947   for(int i=0;i<dim;i++)
948     {
949       std::ostringstream oss;
950       char tmp('X'+i); oss << tmp;
951       ret[i]=DataArray::BuildInfoFromVarAndUnit(oss.str(),_axis_unit);
952     }
953   return ret;
954 }
955
956 void MEDCouplingIMesh::checkSpaceDimension() const
957 {
958   CheckSpaceDimension(_space_dim);
959 }
960
961 void MEDCouplingIMesh::CheckSpaceDimension(int spaceDim)
962 {
963   if(spaceDim<0 || spaceDim>3)
964     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CheckSpaceDimension : input spaceDim must be in [0,1,2,3] !");
965 }
966
967 int MEDCouplingIMesh::FindIntRoot(int val, int order)
968 {
969   if(order==0)
970     return 1;
971   if(val<0)
972     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : input val is < 0 ! Not possible to compute a root !");
973   if(order==1)
974     return val;
975   if(order!=2 && order!=3)
976     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the order available are 0,1,2 or 3 !");
977   double valf((double)val);
978   if(order==2)
979     {
980       double retf(sqrt(valf));
981       int ret((int)retf);
982       if(ret*ret!=val)
983         throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect square root !");
984       return ret;
985     }
986   else//order==3
987     {
988       double retf(std::pow(val,0.3333333333333333));
989       int ret((int)retf),ret2(ret+1);
990       if(ret*ret*ret!=val && ret2*ret2*ret2!=val)
991         throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect cublic root !");
992       if(ret*ret*ret==val)
993         return ret;
994       else
995         return ret2;
996     }
997 }