src/MEDCoupling/MEDCouplingIMesh.cxx

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