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