src/MEDCoupling/MEDCouplingIMesh.cxx

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