src/MEDCoupling/MEDCouplingIMesh.cxx

   1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License, or (at your option) any later version.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "MEDCouplingIMesh.hxx"
  22 #include "MEDCouplingCMesh.hxx"
  23 #include "MEDCouplingMemArray.hxx"
  24 #include "MEDCouplingFieldDouble.hxx"
  25
  26 #include <functional>
  27 #include <algorithm>
  28 #include <sstream>
  29 #include <numeric>
  30
  31 using namespace ParaMEDMEM;
  32
  33 MEDCouplingIMesh::MEDCouplingIMesh():_space_dim(-1)
  34 {
  35   _origin[0]=0.; _origin[1]=0.; _origin[2]=0.;
  36   _dxyz[0]=0.; _dxyz[1]=0.; _dxyz[2]=0.;
  37   _structure[0]=0; _structure[1]=0; _structure[2]=0;
  38 }
  39
  40 MEDCouplingIMesh::MEDCouplingIMesh(const MEDCouplingIMesh& other, bool deepCopy):MEDCouplingStructuredMesh(other,deepCopy),_space_dim(other._space_dim),_axis_unit(other._axis_unit)
  41 {
  42   _origin[0]=other._origin[0]; _origin[1]=other._origin[1]; _origin[2]=other._origin[2];
  43   _dxyz[0]=other._dxyz[0]; _dxyz[1]=other._dxyz[1]; _dxyz[2]=other._dxyz[2];
  44   _structure[0]=other._structure[0]; _structure[1]=other._structure[1]; _structure[2]=other._structure[2];
  45 }
  46
  47 MEDCouplingIMesh::~MEDCouplingIMesh()
  48 {
  49 }
  50
  51 MEDCouplingIMesh *MEDCouplingIMesh::New()
  52 {
  53   return new MEDCouplingIMesh;
  54 }
  55
  56 MEDCouplingIMesh *MEDCouplingIMesh::New(const std::string& meshName, int spaceDim, const int *nodeStrctStart, const int *nodeStrctStop,
  57                                         const double *originStart, const double *originStop, const double *dxyzStart, const double *dxyzStop)
  58 {
  59   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> ret(new MEDCouplingIMesh);
  60   ret->setName(meshName);
  61   ret->setSpaceDimension(spaceDim);
  62   ret->setNodeStruct(nodeStrctStart,nodeStrctStop);
  63   ret->setOrigin(originStart,originStop);
  64   ret->setDXYZ(dxyzStart,dxyzStop);
  65   return ret.retn();
  66 }
  67
  68 MEDCouplingMesh *MEDCouplingIMesh::deepCpy() const
  69 {
  70   return clone(true);
  71 }
  72
  73 MEDCouplingIMesh *MEDCouplingIMesh::clone(bool recDeepCpy) const
  74 {
  75   return new MEDCouplingIMesh(*this,recDeepCpy);
  76 }
  77
  78 void MEDCouplingIMesh::setNodeStruct(const int *nodeStrctStart, const int *nodeStrctStop)
  79 {
  80   checkSpaceDimension();
  81   int sz((int)std::distance(nodeStrctStart,nodeStrctStop));
  82   if(sz!=_space_dim)
  83     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::setNodeStruct : input vector of node structure has not the right size ! Or change space dimension before calling it !");
  84   std::copy(nodeStrctStart,nodeStrctStop,_structure);
  85   declareAsNew();
  86 }
  87
  88 std::vector<int> MEDCouplingIMesh::getNodeStruct() const
  89 {
  90   checkSpaceDimension();
  91   return std::vector<int>(_structure,_structure+_space_dim);
  92 }
  93
  94 void MEDCouplingIMesh::setOrigin(const double *originStart, const double *originStop)
  95 {
  96   checkSpaceDimension();
  97   int sz((int)std::distance(originStart,originStop));
  98   if(sz!=_space_dim)
  99     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::setOrigin : input vector of origin vector has not the right size ! Or change space dimension before calling it !");
 100   std::copy(originStart,originStop,_origin);
 101   declareAsNew();
 102 }
 103
 104 std::vector<double> MEDCouplingIMesh::getOrigin() const
 105 {
 106   checkSpaceDimension();
 107   return std::vector<double>(_origin,_origin+_space_dim);
 108 }
 109
 110 void MEDCouplingIMesh::setDXYZ(const double *dxyzStart, const double *dxyzStop)
 111 {
 112   checkSpaceDimension();
 113   int sz((int)std::distance(dxyzStart,dxyzStop));
 114   if(sz!=_space_dim)
 115     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::setDXYZ : input vector of dxyz vector has not the right size ! Or change space dimension before calling it !");
 116   std::copy(dxyzStart,dxyzStop,_dxyz);
 117   declareAsNew();
 118 }
 119
 120 std::vector<double> MEDCouplingIMesh::getDXYZ() const
 121 {
 122   checkSpaceDimension();
 123   return std::vector<double>(_dxyz,_dxyz+_space_dim);
 124 }
 125
 126 void MEDCouplingIMesh::setAxisUnit(const std::string& unitName)
 127 {
 128   _axis_unit=unitName;
 129   declareAsNew();
 130 }
 131
 132 std::string MEDCouplingIMesh::getAxisUnit() const
 133 {
 134   return _axis_unit;
 135 }
 136
 137 /*!
 138  * This method returns the measure of any cell in \a this.
 139  * This specific method of image grid mesh utilizes the fact that any cell in \a this have the same measure.
 140  * The value returned by this method is those used to feed the returned field in the MEDCouplingIMesh::getMeasureField.
 141  *
 142  * \sa getMeasureField
 143  */
 144 double MEDCouplingIMesh::getMeasureOfAnyCell() const
 145 {
 146   checkCoherency();
 147   int dim(getSpaceDimension());
 148   double ret(1.);
 149   for(int i=0;i<dim;i++)
 150     ret*=fabs(_dxyz[i]);
 151   return ret;
 152 }
 153
 154 /*!
 155  * This method is allows to convert \a this into MEDCouplingCMesh instance.
 156  * This method is the middle level between MEDCouplingIMesh and the most general MEDCouplingUMesh.
 157  * This method is useful for MED writers that do not have still the image grid support.
 158  *
 159  * \sa MEDCouplingMesh::buildUnstructured
 160  */
 161 MEDCouplingCMesh *MEDCouplingIMesh::convertToCartesian() const
 162 {
 163   checkCoherency();
 164   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCMesh> ret(MEDCouplingCMesh::New());
 165   try
 166   { ret->copyTinyInfoFrom(this); }
 167   catch(INTERP_KERNEL::Exception& ) { }
 168   int spaceDim(getSpaceDimension());
 169   std::vector<std::string> infos(buildInfoOnComponents());
 170   for(int i=0;i<spaceDim;i++)
 171     {
 172       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> arr(DataArrayDouble::New()); arr->alloc(_structure[i],1); arr->setInfoOnComponent(0,infos[i]);
 173       arr->iota(); arr->applyLin(_dxyz[i],_origin[i]);
 174       ret->setCoordsAt(i,arr);
 175     }
 176   return ret.retn();
 177 }
 178
 179 /*!
 180  * This method refines \a this uniformaly along all of its dimensions. In case of success the space covered by \a this will remain
 181  * the same before the invocation except that the number of cells will be multiplied by \a factor ^ this->getMeshDimension().
 182  * The origin of \a this will be not touched only spacing and node structure will be changed.
 183  * This method can be useful for AMR users.
 184  */
 185 void MEDCouplingIMesh::refineWithFactor(int factor)
 186 {
 187   if(factor==0)
 188     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::refineWithFactor : refinement factor must be != 0 !");
 189   checkCoherency();
 190   int factAbs(std::abs(factor));
 191   double fact2(1./(double)factor);
 192   std::transform(_structure,_structure+_space_dim,_structure,std::bind2nd(std::plus<int>(),-1));
 193   std::transform(_structure,_structure+_space_dim,_structure,std::bind2nd(std::multiplies<int>(),factAbs));
 194   std::transform(_structure,_structure+_space_dim,_structure,std::bind2nd(std::plus<int>(),1));
 195   std::transform(_dxyz,_dxyz+_space_dim,_dxyz,std::bind2nd(std::multiplies<double>(),fact2));
 196   declareAsNew();
 197 }
 198
 199 /*!
 200  * This static method is useful to condense field on cells of a MEDCouplingIMesh instance coming from a refinement ( MEDCouplingIMesh::refineWithFactor for example)
 201  * to a coarse MEDCouplingIMesh instance. So this method can be seen as a specialization in P0P0 conservative interpolation non overlaping from fine image mesh
 202  * to a coarse image mesh. Only tuples ( deduced from \a fineLocInCoarse ) of \a coarseDA will be modified. Other tuples of \a coarseDA will be let unchanged.
 203  *
 204  * \param [in,out] coarseDA The DataArrayDouble corresponding to the a cell field of a coarse mesh whose cell structure is defined by \a coarseSt.
 205  * \param [in] coarseSt The cell structure of coarse mesh.
 206  * \param [in] fineDA The DataArray containing the cell field on uniformly refined mesh
 207  * \param [in] fineLocInCoarse The cell localization of refined mesh into the coarse one.
 208  */
 209 void MEDCouplingIMesh::CondenseFineToCoarse(DataArrayDouble *coarseDA, const std::vector<int>& coarseSt, const DataArrayDouble *fineDA, const std::vector< std::pair<int,int> >& fineLocInCoarse)
 210 {
 211   if(!coarseDA || !coarseDA->isAllocated() || !fineDA || !fineDA->isAllocated())
 212     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the parameters 1 or 3 are NULL or not allocated !");
 213   int meshDim((int)coarseSt.size()),nbOfTuplesInCoarseExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenStructure(coarseSt)),nbOfTuplesInFineExp(MEDCouplingStructuredMesh::DeduceNumberOfGivenRangeInCompactFrmt(fineLocInCoarse));
 214   int nbCompo(fineDA->getNumberOfComponents());
 215   if(coarseDA->getNumberOfComponents()!=nbCompo)
 216     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the number of components of fine DA and coarse one mismatches !");
 217   if(meshDim!=(int)fineLocInCoarse.size())
 218     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : the size of fineLocInCoarse (4th param) must be equal to the sier of coarseSt (2nd param) !");
 219   if(coarseDA->getNumberOfTuples()!=nbOfTuplesInCoarseExp)
 220     {
 221       std::ostringstream oss; oss << "MEDCouplingIMesh::CondenseFineToCoarse : Expecting " << nbOfTuplesInCoarseExp << " having " << coarseDA->getNumberOfTuples() << " !";
 222       throw INTERP_KERNEL::Exception(oss.str().c_str());
 223     }
 224   int nbTuplesFine(fineDA->getNumberOfTuples());
 225   if(nbTuplesFine%nbOfTuplesInCoarseExp!=0)
 226     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : Invalid nb of tuples in fine DataArray regarding its structure !");
 227   int factN(nbTuplesFine/nbOfTuplesInFineExp);
 228   int fact(FindIntRoot(factN,meshDim));
 229   // to improve use jump-iterator. Factorizes with SwitchOnIdsFrom BuildExplicitIdsFrom
 230   MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids(BuildExplicitIdsFrom(coarseSt,fineLocInCoarse));
 231   const int *idsPtr(ids->begin());
 232   double *outPtr(coarseDA->getPointer());
 233   const double *inPtr(fineDA->begin());
 234   coarseDA->setPartOfValuesSimple3(0.,ids->begin(),ids->end(),0,nbCompo,1);
 235   //
 236   switch(meshDim)
 237   {
 238     case 2:
 239       {
 240         int kk(0);
 241         std::vector<int> dims(MEDCouplingStructuredMesh::GetDimensionsFromCompactFrmt(fineLocInCoarse));
 242         for(int it=0;it<dims[1];it++)
 243           {
 244             for(int i=0;i<fact;i++)
 245               {
 246                 for(int j=0;j<dims[0];j++,inPtr+=fact)
 247                   {
 248                     double *loc(outPtr+idsPtr[kk+j]*nbCompo);
 249                     std::transform(inPtr,inPtr+fact,loc,loc,std::plus<double>());
 250                   }
 251               }
 252             kk+=it;
 253           }
 254         break;
 255       }
 256     default:
 257       throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CondenseFineToCoarse : only dimensions 2 supported !");
 258   }
 259 }
 260
 261 void MEDCouplingIMesh::setSpaceDimension(int spaceDim)
 262 {
 263   if(spaceDim==_space_dim)
 264     return ;
 265   CheckSpaceDimension(spaceDim);
 266   _space_dim=spaceDim;
 267   declareAsNew();
 268 }
 269
 270 void MEDCouplingIMesh::updateTime() const
 271 {
 272 }
 273
 274 std::size_t MEDCouplingIMesh::getHeapMemorySizeWithoutChildren() const
 275 {
 276   return MEDCouplingStructuredMesh::getHeapMemorySizeWithoutChildren();
 277 }
 278
 279 std::vector<const BigMemoryObject *> MEDCouplingIMesh::getDirectChildren() const
 280 {
 281   return std::vector<const BigMemoryObject *>();
 282 }
 283
 284 /*!
 285  * This method copyies all tiny strings from other (name and components name).
 286  * @throw if other and this have not same mesh type.
 287  */
 288 void MEDCouplingIMesh::copyTinyStringsFrom(const MEDCouplingMesh *other)
 289 {
 290   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
 291   if(!otherC)
 292     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::copyTinyStringsFrom : meshes have not same type !");
 293   MEDCouplingStructuredMesh::copyTinyStringsFrom(other);
 294   declareAsNew();
 295 }
 296
 297 bool MEDCouplingIMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
 298 {
 299   if(!other)
 300     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::isEqualIfNotWhy : input other pointer is null !");
 301   const MEDCouplingIMesh *otherC(dynamic_cast<const MEDCouplingIMesh *>(other));
 302   if(!otherC)
 303     {
 304       reason="mesh given in input is not castable in MEDCouplingIMesh !";
 305       return false;
 306     }
 307   if(!MEDCouplingStructuredMesh::isEqualIfNotWhy(other,prec,reason))
 308     return false;
 309   if(!isEqualWithoutConsideringStrInternal(otherC,prec,reason))
 310     return false;
 311   if(_axis_unit!=otherC->_axis_unit)
 312     {
 313       reason="The units of axis are not the same !";
 314       return false;
 315     }
 316   return true;
 317 }
 318
 319 bool MEDCouplingIMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
 320 {
 321   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
 322   if(!otherC)
 323     return false;
 324   std::string tmp;
 325   return isEqualWithoutConsideringStrInternal(other,prec,tmp);
 326 }
 327
 328 bool MEDCouplingIMesh::isEqualWithoutConsideringStrInternal(const MEDCouplingMesh *other, double prec, std::string& reason) const
 329 {
 330   const MEDCouplingIMesh *otherC=dynamic_cast<const MEDCouplingIMesh *>(other);
 331   if(!otherC)
 332     return false;
 333   if(_space_dim!=otherC->_space_dim)
 334     {
 335       std::ostringstream oss;
 336       oss << "The spaceDimension of this (" << _space_dim << ") is not equal to those of other (" << otherC->_space_dim << ") !";
 337       return false;
 338     }
 339   checkSpaceDimension();
 340   for(int i=0;i<_space_dim;i++)
 341     {
 342       if(fabs(_origin[i]-otherC->_origin[i])>prec)
 343         {
 344           std::ostringstream oss;
 345           oss << "The origin of this and other differs at " << i << " !";
 346           reason=oss.str();
 347           return false;
 348         }
 349     }
 350   for(int i=0;i<_space_dim;i++)
 351     {
 352       if(fabs(_dxyz[i]-otherC->_dxyz[i])>prec)
 353         {
 354           std::ostringstream oss;
 355           oss << "The delta of this and other differs at " << i << " !";
 356           reason=oss.str();
 357           return false;
 358         }
 359     }
 360   for(int i=0;i<_space_dim;i++)
 361     {
 362       if(_structure[i]!=otherC->_structure[i])
 363         {
 364           std::ostringstream oss;
 365           oss << "The structure of this and other differs at " << i << " !";
 366           reason=oss.str();
 367           return false;
 368         }
 369     }
 370   return true;
 371 }
 372
 373 void MEDCouplingIMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
 374                                             DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const
 375 {
 376   if(!isEqualWithoutConsideringStr(other,prec))
 377     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::checkDeepEquivalWith : Meshes are not the same !");
 378 }
 379
 380 /*!
 381  * Nothing is done here (except to check that the other is a ParaMEDMEM::MEDCouplingIMesh instance too).
 382  * The user intend that the nodes are the same, so by construction of ParaMEDMEM::MEDCouplingIMesh, \a this and \a other are the same !
 383  */
 384 void MEDCouplingIMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
 385                                                        DataArrayInt *&cellCor) const
 386 {
 387   if(!isEqualWithoutConsideringStr(other,prec))
 388     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::checkDeepEquivalOnSameNodesWith : Meshes are not the same !");
 389 }
 390
 391 void MEDCouplingIMesh::checkCoherency() const
 392 {
 393   checkSpaceDimension();
 394   for(int i=0;i<_space_dim;i++)
 395     if(_structure[i]<1)
 396       {
 397         std::ostringstream oss; oss << "MEDCouplingIMesh::checkCoherency : On axis " << i << "/" << _space_dim << ", number of nodes is equal to " << _structure[i] << " ! must be >=1 !";
 398         throw INTERP_KERNEL::Exception(oss.str().c_str());
 399       }
 400 }
 401
 402 void MEDCouplingIMesh::checkCoherency1(double eps) const
 403 {
 404   checkCoherency();
 405 }
 406
 407 void MEDCouplingIMesh::checkCoherency2(double eps) const
 408 {
 409   checkCoherency1(eps);
 410 }
 411
 412 void MEDCouplingIMesh::getNodeGridStructure(int *res) const
 413 {
 414   checkSpaceDimension();
 415   std::copy(_structure,_structure+_space_dim,res);
 416 }
 417
 418 std::vector<int> MEDCouplingIMesh::getNodeGridStructure() const
 419 {
 420   checkSpaceDimension();
 421   std::vector<int> ret(_structure,_structure+_space_dim);
 422   return ret;
 423 }
 424
 425 MEDCouplingStructuredMesh *MEDCouplingIMesh::buildStructuredSubPart(const std::vector< std::pair<int,int> >& cellPart) const
 426 {
 427   checkCoherency();
 428   int dim(getSpaceDimension());
 429   if(dim!=(int)cellPart.size())
 430     {
 431       std::ostringstream oss; oss << "MEDCouplingIMesh::buildStructuredSubPart : the space dimension is " << dim << " and cell part size is " << cellPart.size() << " !";
 432       throw INTERP_KERNEL::Exception(oss.str().c_str());
 433     }
 434   double retOrigin[3]={0.,0.,0.};
 435   int retStruct[3]={0,0,0};
 436   MEDCouplingAutoRefCountObjectPtr<MEDCouplingIMesh> ret(dynamic_cast<MEDCouplingIMesh *>(deepCpy()));
 437   for(int i=0;i<dim;i++)
 438     {
 439       int startNode(cellPart[i].first),endNode(cellPart[i].second+1);
 440       int myDelta(endNode-startNode);
 441       if(startNode<0 || startNode>=_structure[i])
 442         {
 443           std::ostringstream oss; oss << "MEDCouplingIMesh::buildStructuredSubPart : At dimension #" << i << " the start node id is " << startNode << " it should be in [0," << _structure[i] << ") !";
 444           throw INTERP_KERNEL::Exception(oss.str().c_str());
 445         }
 446       if(myDelta<0 || myDelta>_structure[i])
 447         {
 448           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;
 449           throw INTERP_KERNEL::Exception(oss.str().c_str());
 450         }
 451       retOrigin[i]=_origin[i]+startNode*_dxyz[i];
 452       retStruct[i]=myDelta;
 453     }
 454   ret->setNodeStruct(retStruct,retStruct+dim);
 455   ret->setOrigin(retOrigin,retOrigin+dim);
 456   ret->checkCoherency();
 457   return ret.retn();
 458 }
 459
 460 /*!
 461  * Return the space dimension of \a this.
 462  */
 463 int MEDCouplingIMesh::getSpaceDimension() const
 464 {
 465   return _space_dim;
 466 }
 467
 468 void MEDCouplingIMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const
 469 {
 470   int tmp[3];
 471   int spaceDim(getSpaceDimension());
 472   getSplitNodeValues(tmp);
 473   int tmp2[3];
 474   GetPosFromId(nodeId,spaceDim,tmp,tmp2);
 475   for(int j=0;j<spaceDim;j++)
 476     coo.push_back(_origin[j]+_dxyz[j]*tmp2[j]);
 477 }
 478
 479 std::string MEDCouplingIMesh::simpleRepr() const
 480 {
 481   std::ostringstream ret;
 482   ret << "Image grid with name : \"" << getName() << "\"\n";
 483   ret << "Description of mesh : \"" << getDescription() << "\"\n";
 484   int tmpp1,tmpp2;
 485   double tt(getTime(tmpp1,tmpp2));
 486   int spaceDim(_space_dim);
 487   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
 488   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
 489   ret << "Space dimension : " << spaceDim << "\n";
 490   if(spaceDim<0 || spaceDim>3)
 491     return ret.str();
 492   ret << "The nodal structure is : "; std::copy(_structure,_structure+spaceDim,std::ostream_iterator<int>(ret," ")); ret << "\n";
 493   ret << "The origin position is [" << _axis_unit << "]: ";
 494   std::copy(_origin,_origin+spaceDim,std::ostream_iterator<double>(ret," ")); ret << "\n";
 495   ret << "The intervals along axis are : ";
 496   std::copy(_dxyz,_dxyz+spaceDim,std::ostream_iterator<double>(ret," ")); ret << "\n";
 497   return ret.str();
 498 }
 499
 500 std::string MEDCouplingIMesh::advancedRepr() const
 501 {
 502   return simpleRepr();
 503 }
 504
 505 void MEDCouplingIMesh::getBoundingBox(double *bbox) const
 506 {
 507   checkCoherency();
 508   int dim(getSpaceDimension());
 509   for(int idim=0; idim<dim; idim++)
 510     {
 511       bbox[2*idim]=_origin[idim];
 512       bbox[2*idim+1]=_origin[idim]+_dxyz[idim]*_structure[idim];
 513     }
 514 }
 515
 516 /*!
 517  * Returns a new MEDCouplingFieldDouble containing volumes of cells constituting \a this
 518  * mesh.<br>
 519  * For 1D cells, the returned field contains lengths.<br>
 520  * For 2D cells, the returned field contains areas.<br>
 521  * For 3D cells, the returned field contains volumes.
 522  *  \param [in] isAbs - a not used parameter.
 523  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on cells
 524  *         and one time . The caller is to delete this field using decrRef() as it is no
 525  *         more needed.
 526  */
 527 MEDCouplingFieldDouble *MEDCouplingIMesh::getMeasureField(bool isAbs) const
 528 {
 529   checkCoherency();
 530   std::string name="MeasureOfMesh_";
 531   name+=getName();
 532   int nbelem(getNumberOfCells());
 533   MEDCouplingFieldDouble *field(MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME));
 534   field->setName(name);
 535   DataArrayDouble* array(DataArrayDouble::New());
 536   array->alloc(nbelem,1);
 537   array->fillWithValue(getMeasureOfAnyCell());
 538   field->setArray(array) ;
 539   array->decrRef();
 540   field->setMesh(const_cast<MEDCouplingIMesh *>(this));
 541   field->synchronizeTimeWithMesh();
 542   return field;
 543 }
 544
 545 /*!
 546  * not implemented yet !
 547  */
 548 MEDCouplingFieldDouble *MEDCouplingIMesh::getMeasureFieldOnNode(bool isAbs) const
 549 {
 550   throw INTERP_KERNEL::Exception("MEDCouplingIMesh::getMeasureFieldOnNode : not implemented yet !");
 551   //return 0;
 552 }
 553
 554 int MEDCouplingIMesh::getCellContainingPoint(const double *pos, double eps) const
 555 {
 556   int dim(getSpaceDimension()),ret(0),coeff(1);
 557   for(int i=0;i<dim;i++)
 558     {
 559       int nbOfCells(_structure[i]-1);
 560       double ref(pos[i]);
 561       int tmp((int)((ref-_origin[i])/_dxyz[i]));
 562       if(tmp>=0 && tmp<nbOfCells)
 563         {
 564           ret+=coeff*tmp;
 565           coeff*=nbOfCells;
 566         }
 567       else
 568         return -1;
 569     }
 570   return ret;
 571 }
 572
 573 void MEDCouplingIMesh::rotate(const double *center, const double *vector, double angle)
 574 {
 575   throw INTERP_KERNEL::Exception("No rotation available on IMesh : Traduce it to unstructured mesh to apply it !");
 576 }
 577
 578 /*!
 579  * Translates all nodes of \a this mesh by a given vector. Actually, it adds each
 580  * component of the \a vector to all node coordinates of a corresponding axis.
 581  *  \param [in] vector - the translation vector whose size must be not less than \a
 582  *         this->getSpaceDimension().
 583  */
 584 void MEDCouplingIMesh::translate(const double *vector)
 585 {
 586   checkSpaceDimension();
 587   int dim(getSpaceDimension());
 588   std::transform(_origin,_origin+dim,vector,_origin,std::plus<double>());
 589   declareAsNew();
 590 }
 591
 592 /*!
 593  * Applies scaling transformation to all nodes of \a this mesh.
 594  *  \param [in] point - coordinates of a scaling center. This array is to be of
 595  *         size \a this->getSpaceDimension() at least.
 596  *  \param [in] factor - a scale factor.
 597  */
 598 void MEDCouplingIMesh::scale(const double *point, double factor)
 599 {
 600   checkSpaceDimension();
 601   int dim(getSpaceDimension());
 602   std::transform(_origin,_origin+dim,point,_origin,std::minus<double>());
 603   std::transform(_origin,_origin+dim,_origin,std::bind2nd(std::multiplies<double>(),factor));
 604   std::transform(_dxyz,_dxyz+dim,_dxyz,std::bind2nd(std::multiplies<double>(),factor));
 605   std::transform(_origin,_origin+dim,point,_origin,std::plus<double>());
 606   declareAsNew();
 607 }
 608
 609 MEDCouplingMesh *MEDCouplingIMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
 610 {
 611   //not implemented yet !
 612   return 0;
 613 }
 614
 615 /*!
 616  * Returns a new DataArrayDouble holding coordinates of all nodes of \a this mesh.
 617  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
 618  *          this->getNumberOfNodes() tuples per \a this->getSpaceDimension()
 619  *          components. The caller is to delete this array using decrRef() as it is
 620  *          no more needed.
 621  */
 622 DataArrayDouble *MEDCouplingIMesh::getCoordinatesAndOwner() const
 623 {
 624   checkCoherency();
 625   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
 626   int spaceDim(getSpaceDimension()),nbNodes(getNumberOfNodes());
 627   ret->alloc(nbNodes,spaceDim);
 628   double *pt(ret->getPointer());
 629   ret->setInfoOnComponents(buildInfoOnComponents());
 630   int tmp2[3],tmp[3];
 631   getSplitNodeValues(tmp);
 632   for(int i=0;i<nbNodes;i++)
 633     {
 634       GetPosFromId(i,spaceDim,tmp,tmp2);
 635       for(int j=0;j<spaceDim;j++)
 636         pt[i*spaceDim+j]=_dxyz[j]*tmp2[j]+_origin[j];
 637     }
 638   return ret.retn();
 639 }
 640
 641 /*!
 642  * Returns a new DataArrayDouble holding barycenters of all cells. The barycenter is
 643  * computed by averaging coordinates of cell nodes.
 644  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
 645  *          this->getNumberOfCells() tuples per \a this->getSpaceDimension()
 646  *          components. The caller is to delete this array using decrRef() as it is
 647  *          no more needed.
 648  */
 649 DataArrayDouble *MEDCouplingIMesh::getBarycenterAndOwner() const
 650 {
 651   checkCoherency();
 652   MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> ret(DataArrayDouble::New());
 653   int spaceDim(getSpaceDimension()),nbCells(getNumberOfCells()),tmp[3],tmp2[3];
 654   ret->alloc(nbCells,spaceDim);
 655   double *pt(ret->getPointer()),shiftOrigin[3];
 656   std::transform(_dxyz,_dxyz+spaceDim,shiftOrigin,std::bind2nd(std::multiplies<double>(),0.5));
 657   std::transform(_origin,_origin+spaceDim,shiftOrigin,shiftOrigin,std::plus<double>());
 658   getSplitCellValues(tmp);
 659   ret->setInfoOnComponents(buildInfoOnComponents());
 660   for(int i=0;i<nbCells;i++)
 661     {
 662       GetPosFromId(i,spaceDim,tmp,tmp2);
 663       for(int j=0;j<spaceDim;j++)
 664         pt[i*spaceDim+j]=_dxyz[j]*tmp2[j]+shiftOrigin[j];
 665     }
 666   return ret.retn();
 667 }
 668
 669 DataArrayDouble *MEDCouplingIMesh::computeIsoBarycenterOfNodesPerCell() const
 670 {
 671   return MEDCouplingIMesh::getBarycenterAndOwner();
 672 }
 673
 674 void MEDCouplingIMesh::renumberCells(const int *old2NewBg, bool check)
 675 {
 676   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for IMesh !");
 677 }
 678
 679 void MEDCouplingIMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 680 {
 681   int it,order;
 682   double time(getTime(it,order));
 683   tinyInfo.clear();
 684   tinyInfoD.clear();
 685   littleStrings.clear();
 686   littleStrings.push_back(getName());
 687   littleStrings.push_back(getDescription());
 688   littleStrings.push_back(getTimeUnit());
 689   littleStrings.push_back(getAxisUnit());
 690   tinyInfo.push_back(it);
 691   tinyInfo.push_back(order);
 692   tinyInfo.push_back(_space_dim);
 693   tinyInfo.insert(tinyInfo.end(),_structure,_structure+3);
 694   tinyInfoD.push_back(time);
 695   tinyInfoD.insert(tinyInfoD.end(),_dxyz,_dxyz+3);
 696   tinyInfoD.insert(tinyInfoD.end(),_origin,_origin+3);
 697 }
 698
 699 void MEDCouplingIMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 700 {
 701   a1->alloc(0,1);
 702   a2->alloc(0,1);
 703 }
 704
 705 void MEDCouplingIMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
 706 {
 707   a1=DataArrayInt::New();
 708   a1->alloc(0,1);
 709   a2=DataArrayDouble::New();
 710   a2->alloc(0,1);
 711 }
 712
 713 void MEDCouplingIMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
 714                                        const std::vector<std::string>& littleStrings)
 715 {
 716   setName(littleStrings[0]);
 717   setDescription(littleStrings[1]);
 718   setTimeUnit(littleStrings[2]);
 719   setAxisUnit(littleStrings[3]);
 720   setTime(tinyInfoD[0],tinyInfo[0],tinyInfo[1]);
 721   _space_dim=tinyInfo[2];
 722   _structure[0]=tinyInfo[3]; _structure[1]=tinyInfo[4]; _structure[2]=tinyInfo[5];
 723   _dxyz[0]=tinyInfoD[1]; _dxyz[1]=tinyInfoD[2]; _dxyz[2]=tinyInfoD[3];
 724   _origin[0]=tinyInfoD[4]; _origin[1]=tinyInfoD[5]; _origin[2]=tinyInfoD[6];
 725   declareAsNew();
 726 }
 727
 728 void MEDCouplingIMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
 729 {
 730   checkCoherency();
 731   std::ostringstream extent,origin,spacing;
 732   for(int i=0;i<3;i++)
 733     {
 734       if(i<_space_dim)
 735         { extent << "0 " <<  _structure[i]-1 << " "; origin << _origin[i] << " "; spacing << _dxyz[i] << " "; }
 736       else
 737         { extent << "0 0 "; origin << "0 "; spacing << "0 "; }
 738     }
 739   ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\" Origin=\"" << origin.str() << "\" Spacing=\"" << spacing.str() << "\">\n";
 740   ofs << "    <Piece Extent=\"" << extent.str() << "\">\n";
 741   ofs << "      <PointData>\n" << pointData << std::endl;
 742   ofs << "      </PointData>\n";
 743   ofs << "      <CellData>\n" << cellData << std::endl;
 744   ofs << "      </CellData>\n";
 745   ofs << "      <Coordinates>\n";
 746   ofs << "      </Coordinates>\n";
 747   ofs << "    </Piece>\n";
 748   ofs << "  </" << getVTKDataSetType() << ">\n";
 749 }
 750
 751 void MEDCouplingIMesh::reprQuickOverview(std::ostream& stream) const
 752 {
 753   stream << "MEDCouplingIMesh C++ instance at " << this << ". Name : \"" << getName() << "\". Space dimension : " << _space_dim << ".";
 754   if(_space_dim<0 || _space_dim>3)
 755     return ;
 756   stream << "\n";
 757   std::ostringstream stream0,stream1;
 758   int nbNodes(1),nbCells(0);
 759   bool isPb(false);
 760   for(int i=0;i<_space_dim;i++)
 761     {
 762       char tmp('X'+i);
 763       int tmpNodes(_structure[i]);
 764       stream1 << "- Axis " << tmp << " : " << tmpNodes << " nodes (orig=" << _origin[i] << ", inter=" << _dxyz[i] << ").";
 765       if(i!=_space_dim-1)
 766         stream1 << std::endl;
 767       if(tmpNodes>=1)
 768         nbNodes*=tmpNodes;
 769       else
 770         isPb=true;
 771       if(tmpNodes>=2)
 772         nbCells=nbCells==0?tmpNodes-1:nbCells*(tmpNodes-1);
 773     }
 774   if(!isPb)
 775     {
 776       stream0 << "Number of cells : " << nbCells << ", Number of nodes : " << nbNodes;
 777       stream << stream0.str();
 778       if(_space_dim>0)
 779         stream << std::endl;
 780     }
 781   stream << stream1.str();
 782 }
 783
 784 std::string MEDCouplingIMesh::getVTKDataSetType() const
 785 {
 786   return std::string("ImageData");
 787 }
 788
 789 std::vector<std::string> MEDCouplingIMesh::buildInfoOnComponents() const
 790 {
 791   checkSpaceDimension();
 792   int dim(getSpaceDimension());
 793   std::vector<std::string> ret(dim);
 794   for(int i=0;i<dim;i++)
 795     {
 796       std::ostringstream oss;
 797       char tmp('X'+i); oss << tmp;
 798       ret[i]=DataArray::BuildInfoFromVarAndUnit(oss.str(),_axis_unit);
 799     }
 800   return ret;
 801 }
 802
 803 void MEDCouplingIMesh::checkSpaceDimension() const
 804 {
 805   CheckSpaceDimension(_space_dim);
 806 }
 807
 808 void MEDCouplingIMesh::CheckSpaceDimension(int spaceDim)
 809 {
 810   if(spaceDim<0 || spaceDim>3)
 811     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::CheckSpaceDimension : input spaceDim must be in [0,1,2,3] !");
 812 }
 813
 814 int MEDCouplingIMesh::FindIntRoot(int val, int order)
 815 {
 816   if(order==0)
 817     return 1;
 818   if(val<0)
 819     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : input val is < 0 ! Not possible to compute a root !");
 820   if(order==1)
 821     return val;
 822   if(order!=2 && order!=3)
 823     throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the order available are 0,1,2 or 3 !");
 824   double valf((double)val);
 825   if(order==2)
 826     {
 827       double retf(sqrt(valf));
 828       int ret((int)retf);
 829       if(ret*ret!=val)
 830         throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect square root !");
 831       return ret;
 832     }
 833   else//order==3
 834     {
 835       double retf(std::pow(val,0.3333333333333333));
 836       int ret((int)retf),ret2(ret+1);
 837       if(ret*ret*ret!=val || ret2*ret2*ret2!=val)
 838         throw INTERP_KERNEL::Exception("MEDCouplingIMesh::FindIntRoot : the input val is not a perfect cublic root !");
 839       if(ret*ret*ret==val)
 840         return ret;
 841       else
 842         return ret2;
 843     }
 844 }