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