src/MEDCoupling/MEDCouplingCMesh.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 "MEDCouplingCMesh.hxx"
  22 #include "MEDCouplingMemArray.hxx"
  23 #include "MEDCouplingFieldDouble.hxx"
  24
  25 #include <functional>
  26 #include <algorithm>
  27 #include <sstream>
  28 #include <numeric>
  29
  30 using namespace ParaMEDMEM;
  31
  32 MEDCouplingCMesh::MEDCouplingCMesh():_x_array(0),_y_array(0),_z_array(0)
  33 {
  34 }
  35
  36 MEDCouplingCMesh::MEDCouplingCMesh(const MEDCouplingCMesh& other, bool deepCopy):MEDCouplingStructuredMesh(other,deepCopy)
  37 {
  38   if(deepCopy)
  39     {
  40       if(other._x_array)
  41         _x_array=other._x_array->deepCpy();
  42       else
  43         _x_array=0;
  44       if(other._y_array)
  45         _y_array=other._y_array->deepCpy();
  46       else
  47         _y_array=0;
  48       if(other._z_array)
  49         _z_array=other._z_array->deepCpy();
  50       else
  51         _z_array=0;
  52     }
  53   else
  54     {
  55       _x_array=other._x_array;
  56       if(_x_array)
  57         _x_array->incrRef();
  58       _y_array=other._y_array;
  59       if(_y_array)
  60         _y_array->incrRef();
  61       _z_array=other._z_array;
  62       if(_z_array)
  63         _z_array->incrRef();
  64     }
  65 }
  66
  67 MEDCouplingCMesh::~MEDCouplingCMesh()
  68 {
  69   if(_x_array)
  70     _x_array->decrRef();
  71   if(_y_array)
  72     _y_array->decrRef();
  73   if(_z_array)
  74     _z_array->decrRef();
  75 }
  76
  77 MEDCouplingCMesh *MEDCouplingCMesh::New()
  78 {
  79   return new MEDCouplingCMesh;
  80 }
  81
  82 MEDCouplingCMesh *MEDCouplingCMesh::New(const std::string& meshName)
  83 {
  84   MEDCouplingCMesh *ret=new MEDCouplingCMesh;
  85   ret->setName(meshName);
  86   return ret;
  87 }
  88
  89 MEDCouplingMesh *MEDCouplingCMesh::deepCpy() const
  90 {
  91   return clone(true);
  92 }
  93
  94 MEDCouplingCMesh *MEDCouplingCMesh::clone(bool recDeepCpy) const
  95 {
  96   return new MEDCouplingCMesh(*this,recDeepCpy);
  97 }
  98
  99 void MEDCouplingCMesh::updateTime() const
 100 {
 101   if(_x_array)
 102     updateTimeWith(*_x_array);
 103   if(_y_array)
 104     updateTimeWith(*_y_array);
 105   if(_z_array)
 106     updateTimeWith(*_z_array);
 107 }
 108
 109 std::size_t MEDCouplingCMesh::getHeapMemorySizeWithoutChildren() const
 110 {
 111   return MEDCouplingStructuredMesh::getHeapMemorySizeWithoutChildren();
 112 }
 113
 114 std::vector<const BigMemoryObject *> MEDCouplingCMesh::getDirectChildren() const
 115 {
 116   std::vector<const BigMemoryObject *> ret;
 117   if(_x_array)
 118     ret.push_back(_x_array);
 119   if(_y_array)
 120     ret.push_back(_y_array);
 121   if(_z_array)
 122     ret.push_back(_z_array);
 123   return ret;
 124 }
 125
 126 /*!
 127  * This method copyies all tiny strings from other (name and components name).
 128  * @throw if other and this have not same mesh type.
 129  */
 130 void MEDCouplingCMesh::copyTinyStringsFrom(const MEDCouplingMesh *other)
 131 {
 132   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
 133   if(!otherC)
 134     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::copyTinyStringsFrom : meshes have not same type !");
 135   MEDCouplingStructuredMesh::copyTinyStringsFrom(other);
 136   if(_x_array && otherC->_x_array)
 137     _x_array->copyStringInfoFrom(*otherC->_x_array);
 138   if(_y_array && otherC->_y_array)
 139     _y_array->copyStringInfoFrom(*otherC->_y_array);
 140   if(_z_array && otherC->_z_array)
 141     _z_array->copyStringInfoFrom(*otherC->_z_array);
 142 }
 143
 144 bool MEDCouplingCMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const
 145 {
 146   if(!other)
 147     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::isEqualIfNotWhy : input other pointer is null !");
 148   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
 149   if(!otherC)
 150     {
 151       reason="mesh given in input is not castable in MEDCouplingCMesh !";
 152       return false;
 153     }
 154   if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
 155     return false;
 156   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 157   const DataArrayDouble *otherArr[3]={otherC->_x_array,otherC->_y_array,otherC->_z_array};
 158   std::ostringstream oss; oss.precision(15);
 159   for(int i=0;i<3;i++)
 160     {
 161       if((thisArr[i]!=0 && otherArr[i]==0) || (thisArr[i]==0 && otherArr[i]!=0))
 162         {
 163           oss << "Only one CMesh between the two this and other has its coordinates of rank" << i << " defined !";
 164           reason=oss.str();
 165           return false;
 166         }
 167       if(thisArr[i])
 168         if(!thisArr[i]->isEqualIfNotWhy(*otherArr[i],prec,reason))
 169           {
 170             oss << "Coordinates DataArrayDouble of rank #" << i << " differ :";
 171             reason.insert(0,oss.str());
 172             return false;
 173           }
 174     }
 175   return true;
 176 }
 177
 178 bool MEDCouplingCMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
 179 {
 180   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
 181   if(!otherC)
 182     return false;
 183   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 184   const DataArrayDouble *otherArr[3]={otherC->_x_array,otherC->_y_array,otherC->_z_array};
 185   for(int i=0;i<3;i++)
 186     {
 187       if((thisArr[i]!=0 && otherArr[i]==0) || (thisArr[i]==0 && otherArr[i]!=0))
 188         return false;
 189       if(thisArr[i])
 190         if(!thisArr[i]->isEqualWithoutConsideringStr(*otherArr[i],prec))
 191           return false;
 192     }
 193   return true;
 194 }
 195
 196 void MEDCouplingCMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
 197                                             DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const
 198 {
 199   if(!isEqualWithoutConsideringStr(other,prec))
 200     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkDeepEquivalWith : Meshes are not the same !");
 201 }
 202
 203 /*!
 204  * Nothing is done here (except to check that the other is a ParaMEDMEM::MEDCouplingCMesh instance too).
 205  * The user intend that the nodes are the same, so by construction of ParaMEDMEM::MEDCouplingCMesh, 'this' and 'other' are the same !
 206  */
 207 void MEDCouplingCMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
 208                                                        DataArrayInt *&cellCor) const
 209 {
 210   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
 211   if(!otherC)
 212     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkDeepEquivalOnSameNodesWith : other is NOT a cartesian mesh ! Impossible to check equivalence !");
 213 }
 214
 215 void MEDCouplingCMesh::checkCoherency() const
 216 {
 217   const char msg0[]="Invalid ";
 218   const char msg1[]=" array ! Must contain more than 1 element.";
 219   const char msg2[]=" array ! Must be with only one component.";
 220   getSpaceDimension();// here to check that no holes in arrays !
 221   if(_x_array)
 222     {
 223       if(_x_array->getNbOfElems()<2)
 224         {
 225           std::ostringstream os; os << msg0 << 'X' << msg1;
 226           throw INTERP_KERNEL::Exception(os.str().c_str());
 227         }
 228       if(_x_array->getNumberOfComponents()!=1)
 229         {
 230           std::ostringstream os; os << msg0 << 'X' << msg2;
 231           throw INTERP_KERNEL::Exception(os.str().c_str());
 232         }
 233     }
 234   if(_y_array)
 235     {
 236       if(_y_array->getNbOfElems()<2)
 237         {
 238           std::ostringstream os; os << msg0 << 'Y' << msg1;
 239           throw INTERP_KERNEL::Exception(os.str().c_str());
 240         }
 241       if(_y_array->getNumberOfComponents()!=1)
 242         {
 243           std::ostringstream os; os << msg0 << 'Y' << msg2;
 244           throw INTERP_KERNEL::Exception(os.str().c_str());
 245         }
 246     }
 247   if(_z_array)
 248     {
 249       if(_z_array->getNbOfElems()<2)
 250         {
 251           std::ostringstream os; os << msg0 << 'Z' << msg1;
 252           throw INTERP_KERNEL::Exception(os.str().c_str());
 253         }
 254       if(_z_array->getNumberOfComponents()!=1)
 255         {
 256           std::ostringstream os; os << msg0 << 'Z' << msg2;
 257           throw INTERP_KERNEL::Exception(os.str().c_str());
 258         }
 259     }
 260 }
 261
 262 void MEDCouplingCMesh::checkCoherency1(double eps) const
 263 {
 264   checkCoherency();
 265   if(_x_array)
 266     _x_array->checkMonotonic(true, eps);
 267   if(_y_array)
 268     _y_array->checkMonotonic(true, eps);
 269   if(_z_array)
 270     _z_array->checkMonotonic(true, eps);
 271 }
 272
 273 void MEDCouplingCMesh::checkCoherency2(double eps) const
 274 {
 275   checkCoherency1(eps);
 276 }
 277
 278 int MEDCouplingCMesh::getNumberOfCells() const
 279 {
 280   int ret=1;
 281   if(_x_array)
 282     ret*=_x_array->getNbOfElems()-1;
 283   if(_y_array)
 284     ret*=_y_array->getNbOfElems()-1;
 285   if(_z_array)
 286     ret*=_z_array->getNbOfElems()-1;
 287   return ret;
 288 }
 289
 290 int MEDCouplingCMesh::getNumberOfNodes() const
 291 {
 292   int ret=1;
 293   if(_x_array)
 294     ret*=_x_array->getNbOfElems();
 295   if(_y_array)
 296     ret*=_y_array->getNbOfElems();
 297   if(_z_array)
 298     ret*=_z_array->getNbOfElems();
 299   return ret;
 300 }
 301
 302 void MEDCouplingCMesh::getSplitCellValues(int *res) const
 303 {
 304   int meshDim(getMeshDimension());
 305   for(int l=0;l<meshDim;l++)
 306     {
 307       int val=1;
 308       for(int p=0;p<meshDim-l-1;p++)
 309         val*=getCoordsAt(p)->getNbOfElems()-1;
 310       res[meshDim-l-1]=val;
 311     }
 312 }
 313
 314 void MEDCouplingCMesh::getSplitNodeValues(int *res) const
 315 {
 316   int spaceDim(getSpaceDimension());
 317   for(int l=0;l<spaceDim;l++)
 318     {
 319       int val=1;
 320       for(int p=0;p<spaceDim-l-1;p++)
 321         val*=getCoordsAt(p)->getNbOfElems();
 322       res[spaceDim-l-1]=val;
 323     }
 324 }
 325
 326 void MEDCouplingCMesh::getNodeGridStructure(int *res) const
 327 {
 328   int spaceDim(getSpaceDimension());
 329   for(int i=0;i<spaceDim;i++)
 330     res[i]=getCoordsAt(i)->getNbOfElems();
 331 }
 332
 333 std::vector<int> MEDCouplingCMesh::getNodeGridStructure() const
 334 {
 335   std::vector<int> ret(getSpaceDimension());
 336   getNodeGridStructure(&ret[0]);
 337   return ret;
 338 }
 339
 340 MEDCouplingStructuredMesh *MEDCouplingCMesh::buildStructuredSubPart(const std::vector< std::pair<int,int> >& cellPart) const
 341 {
 342   checkCoherency();
 343   int dim(getMeshDimension());
 344   if(dim!=(int)cellPart.size())
 345     {
 346       std::ostringstream oss; oss << "MEDCouplingCMesh::buildStructuredSubPart : the mesh dimension is " << dim << " and cell part size is " << cellPart.size() << " !";
 347       throw INTERP_KERNEL::Exception(oss.str().c_str());
 348     }
 349   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCMesh> ret(dynamic_cast<MEDCouplingCMesh *>(deepCpy()));
 350   for(int i=0;i<dim;i++)
 351     {
 352       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> tmp(ret->getCoordsAt(i)->selectByTupleId2(cellPart[i].first,cellPart[i].second+1,1));
 353       ret->setCoordsAt(i,tmp);
 354     }
 355   return ret.retn();
 356 }
 357
 358 /*!
 359  * Return the space dimension of \a this. It only considers the arrays along X, Y and Z to deduce that.
 360  * This method throws exceptions if the not null arrays defining this are not contiguouly at the end. For example X!=0,Y==0,Z!=0 will throw.
 361  */
 362 int MEDCouplingCMesh::getSpaceDimension() const
 363 {
 364   static const char MSG[]="MEDCouplingCMesh::getSpaceDimension : mesh is invalid ! null vectors (X, Y or Z) must be put contiguously at the end !";
 365   int ret(0);
 366   bool isOK(true);
 367   if(_x_array)
 368     ret++;
 369   else
 370     isOK=false;
 371   if(_y_array)
 372     {
 373       if(!isOK)
 374         throw INTERP_KERNEL::Exception(MSG);
 375       ret++;
 376     }
 377   else
 378     isOK=false;
 379   if(_z_array)
 380     {
 381       if(!isOK)
 382         throw INTERP_KERNEL::Exception(MSG);
 383       ret++;
 384     }
 385   return ret;
 386 }
 387
 388 /*!
 389  * This method returns the mesh dimension of \a this. It can be different from space dimension in case of a not null dimension contains only one node.
 390  */
 391 int MEDCouplingCMesh::getMeshDimension() const
 392 {
 393   int ret(getSpaceDimension());
 394   if(_x_array)
 395     {
 396       if(_x_array->isAllocated())
 397         if(_x_array->getNumberOfTuples()==1)
 398           ret--;
 399     }
 400   if(_y_array)
 401     {
 402       if(_y_array->isAllocated())
 403         if(_y_array->getNumberOfTuples()==1)
 404           ret--;
 405     }
 406   if(_z_array)
 407     {
 408       if(_z_array->isAllocated())
 409         if(_z_array->getNumberOfTuples()==1)
 410           ret--;
 411     }
 412   return ret;
 413 }
 414
 415 void MEDCouplingCMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const
 416 {
 417   int tmp[3];
 418   int spaceDim=getSpaceDimension();
 419   getSplitNodeValues(tmp);
 420   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 421   int tmp2[3];
 422   GetPosFromId(nodeId,spaceDim,tmp,tmp2);
 423   for(int j=0;j<spaceDim;j++)
 424     if(tabs[j])
 425       coo.push_back(tabs[j]->getConstPointer()[tmp2[j]]);
 426 }
 427
 428 std::string MEDCouplingCMesh::simpleRepr() const
 429 {
 430   std::ostringstream ret;
 431   ret << "Cartesian mesh with name : \"" << getName() << "\"\n";
 432   ret << "Description of mesh : \"" << getDescription() << "\"\n";
 433   int tmpp1,tmpp2;
 434   double tt=getTime(tmpp1,tmpp2);
 435   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
 436   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
 437   ret << "Mesh and SpaceDimension dimension : " << getSpaceDimension() << "\n\nArrays :\n________\n\n";
 438   if(_x_array)
 439     {
 440       ret << "X Array :\n";
 441       _x_array->reprZipWithoutNameStream(ret);
 442     }
 443   if(_y_array)
 444     {
 445       ret << "Y Array :\n";
 446       _y_array->reprZipWithoutNameStream(ret);
 447     }
 448   if(_z_array)
 449     {
 450       ret << "Z Array :\n";
 451       _z_array->reprZipWithoutNameStream(ret);
 452     }
 453   return ret.str();
 454 }
 455
 456 std::string MEDCouplingCMesh::advancedRepr() const
 457 {
 458   return simpleRepr();
 459 }
 460
 461 /*!
 462  * Returns a DataArrayDouble holding positions of nodes along a given axis.
 463  * For more info on Cartesian meshes, see \ref MEDCouplingCMeshPage.
 464  *  \param [in] i - an index of axis, a value from [0,1,2].
 465  *  \return const DataArrayDouble * - a pointer to the data array of node coordinates
 466  *         referred by \a this mesh.
 467  *  \throw If \a i is not one of [0,1,2].
 468  *
 469  *  \if ENABLE_EXAMPLES
 470  *  \ref cpp_mccmesh_getCoordsAt "Here is a C++ example".<br>
 471  *  \ref  py_mccmesh_getCoordsAt "Here is a Python example".
 472  *  \endif
 473  */
 474 const DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i) const
 475 {
 476   switch(i)
 477   {
 478     case 0:
 479       return _x_array;
 480     case 1:
 481       return _y_array;
 482     case 2:
 483       return _z_array;
 484     default:
 485       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 486   }
 487 }
 488
 489 /*!
 490  * Returns a DataArrayDouble holding positions of nodes along a given axis.
 491  * For more info on Cartesian meshes, see \ref MEDCouplingCMeshPage.
 492  *  \param [in] i - an index of axis, a value from [0,1,2].
 493  *  \return const DataArrayDouble * - a pointer to the data array of node coordinates
 494  *         referred by \a this mesh.
 495  *  \throw If \a i is not one of [0,1,2].
 496  *
 497  *  \if ENABLE_EXAMPLES
 498  *  \ref cpp_mccmesh_getCoordsAt "Here is a C++ example".<br>
 499  *  \ref  py_mccmesh_getCoordsAt "Here is a Python example".
 500  *  \endif
 501  */
 502 DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i)
 503 {
 504   switch(i)
 505   {
 506     case 0:
 507       return _x_array;
 508     case 1:
 509       return _y_array;
 510     case 2:
 511       return _z_array;
 512     default:
 513       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 514   }
 515 }
 516
 517 /*!
 518  * Sets node coordinates along a given axis. For more info on Cartesian meshes, see
 519  * \ref MEDCouplingCMeshPage.
 520  *  \param [in] i - an index of axis, a value in range [0,1,2].
 521  *  \param [in] arr - DataArrayDouble holding positions of nodes along the i-th
 522  *         axis. It must be an array of one component.
 523  *  \throw If \a arr->getNumberOfComponents() != 1.
 524  *  \throw If \a i is not one of [0,1,2].
 525  *
 526  *  \if ENABLE_EXAMPLES
 527  *  \ref medcouplingcppexamplesCmeshStdBuild1 "Here is a C++ example".<br>
 528  *  \ref  medcouplingpyexamplesCmeshStdBuild1 "Here is a Python example".
 529  *  \endif
 530  */
 531 void MEDCouplingCMesh::setCoordsAt(int i, const DataArrayDouble *arr)
 532 {
 533   if(arr)
 534     arr->checkNbOfComps(1,"MEDCouplingCMesh::setCoordsAt");
 535   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
 536   if(i<0 || i>2)
 537     throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 538   if(arr!=*(thisArr[i]))
 539     {
 540       if(*(thisArr[i]))
 541         (*(thisArr[i]))->decrRef();
 542       (*(thisArr[i]))=const_cast<DataArrayDouble *>(arr);
 543       if(*(thisArr[i]))
 544         (*(thisArr[i]))->incrRef();
 545       declareAsNew();
 546     }
 547 }
 548
 549 /*!
 550  * Sets node coordinates along some of the tree axes. This method updates all the
 551  * three node coordinates arrays at once. For more info on Cartesian meshes, see
 552  * \ref MEDCouplingCMeshPage.
 553  *  \param [in] coordsX - DataArrayDouble holding positions of nodes along the X
 554  *         axis. It must be an array of one component or \c NULL.
 555  *  \param [in] coordsY - DataArrayDouble holding positions of nodes along the Y
 556  *         axis. It must be an array of one component or \c NULL.
 557  *  \param [in] coordsZ - DataArrayDouble holding positions of nodes along the Z
 558  *         axis. It must be an array of one component or \c NULL.
 559  *  \throw If \a coords*->getNumberOfComponents() != 1.
 560  *
 561  *  \if ENABLE_EXAMPLES
 562  *  \ref medcouplingcppexamplesCmeshStdBuild1 "Here is a C++ example".<br>
 563  *  \ref  medcouplingpyexamplesCmeshStdBuild1 "Here is a Python example".
 564  *  \endif
 565  */
 566 void MEDCouplingCMesh::setCoords(const DataArrayDouble *coordsX, const DataArrayDouble *coordsY, const DataArrayDouble *coordsZ)
 567 {
 568   if(coordsX)
 569     coordsX->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsX");
 570   if(coordsY)
 571     coordsY->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsY");
 572   if(coordsZ)
 573     coordsZ->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsZ");
 574   if(_x_array)
 575     _x_array->decrRef();
 576   _x_array=const_cast<DataArrayDouble *>(coordsX);
 577   if(_x_array)
 578     _x_array->incrRef();
 579   if(_y_array)
 580     _y_array->decrRef();
 581   _y_array=const_cast<DataArrayDouble *>(coordsY);
 582   if(_y_array)
 583     _y_array->incrRef();
 584   if(_z_array)
 585     _z_array->decrRef();
 586   _z_array=const_cast<DataArrayDouble *>(coordsZ);
 587   if(_z_array)
 588     _z_array->incrRef();
 589   declareAsNew();
 590 }
 591
 592 void MEDCouplingCMesh::getBoundingBox(double *bbox) const
 593 {
 594   int dim=getSpaceDimension();
 595   int j=0;
 596   for (int idim=0; idim<dim; idim++)
 597     {
 598       const DataArrayDouble *c=getCoordsAt(idim);
 599       if(c)
 600         {
 601           const double *coords=c->getConstPointer();
 602           int nb=c->getNbOfElems();
 603           bbox[2*j]=coords[0];
 604           bbox[2*j+1]=coords[nb-1];
 605           j++;
 606         }
 607     }
 608 }
 609
 610 /*!
 611  * Returns a new MEDCouplingFieldDouble containing volumes of cells constituting \a this
 612  * mesh.<br>
 613  * For 1D cells, the returned field contains lengths.<br>
 614  * For 2D cells, the returned field contains areas.<br>
 615  * For 3D cells, the returned field contains volumes.
 616  *  \param [in] isAbs - a not used parameter.
 617  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on cells
 618  *         and one time . The caller is to delete this field using decrRef() as it is no
 619  *         more needed.
 620  */
 621 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const
 622 {
 623   std::string name="MeasureOfMesh_";
 624   name+=getName();
 625   int nbelem=getNumberOfCells();
 626   MEDCouplingFieldDouble *field=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
 627   field->setName(name);
 628   DataArrayDouble* array=DataArrayDouble::New();
 629   array->alloc(nbelem,1);
 630   double *area_vol=array->getPointer();
 631   field->setArray(array) ;
 632   array->decrRef();
 633   field->setMesh(const_cast<MEDCouplingCMesh *>(this));
 634   field->synchronizeTimeWithMesh();
 635   int tmp[3];
 636   getSplitCellValues(tmp);
 637   int dim=getSpaceDimension();
 638   const double **thisArr=new const double *[dim];
 639   const DataArrayDouble *thisArr2[3]={_x_array,_y_array,_z_array};
 640   for(int i=0;i<dim;i++)
 641     thisArr[i]=thisArr2[i]->getConstPointer();
 642   for(int icell=0;icell<nbelem;icell++)
 643     {
 644       int tmp2[3];
 645       GetPosFromId(icell,dim,tmp,tmp2);
 646       area_vol[icell]=1.;
 647       for(int i=0;i<dim;i++)
 648         area_vol[icell]*=thisArr[i][tmp2[i]+1]-thisArr[i][tmp2[i]];
 649     }
 650   delete [] thisArr;
 651   return field;
 652 }
 653
 654 /*!
 655  * not implemented yet !
 656  */
 657 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureFieldOnNode(bool isAbs) const
 658 {
 659   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getMeasureFieldOnNode : not implemented yet !");
 660   //return 0;
 661 }
 662
 663 int MEDCouplingCMesh::getCellContainingPoint(const double *pos, double eps) const
 664 {
 665   int dim=getSpaceDimension();
 666   int ret=0;
 667   int coeff=1;
 668   for(int i=0;i<dim;i++)
 669     {
 670       const double *d=getCoordsAt(i)->getConstPointer();
 671       int nbOfNodes=getCoordsAt(i)->getNbOfElems();
 672       double ref=pos[i];
 673       const double *w=std::find_if(d,d+nbOfNodes,std::bind2nd(std::greater_equal<double>(),ref));
 674       int w2=(int)std::distance(d,w);
 675       if(w2<nbOfNodes)
 676         {
 677           if(w2==0)
 678             {
 679               if(ref>d[0]-eps)
 680                 w2=1;
 681               else
 682                 return -1;
 683             }
 684           ret+=coeff*(w2-1);
 685           coeff*=nbOfNodes-1;
 686         }
 687       else
 688         return -1;
 689     }
 690   return ret;
 691 }
 692
 693 void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
 694 {
 695   throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !");
 696 }
 697
 698 /*!
 699  * Translates all nodes of \a this mesh by a given vector. Actually, it adds each
 700  * component of the \a vector to all node coordinates of a corresponding axis.
 701  *  \param [in] vector - the translation vector whose size must be not less than \a
 702  *         this->getSpaceDimension().
 703  */
 704 void MEDCouplingCMesh::translate(const double *vector)
 705 {
 706   if(_x_array)
 707     std::transform(_x_array->getConstPointer(),_x_array->getConstPointer()+_x_array->getNbOfElems(),
 708                    _x_array->getPointer(),std::bind2nd(std::plus<double>(),vector[0]));
 709   if(_y_array)
 710     std::transform(_y_array->getConstPointer(),_y_array->getConstPointer()+_y_array->getNbOfElems(),
 711                    _y_array->getPointer(),std::bind2nd(std::plus<double>(),vector[1]));
 712   if(_z_array)
 713     std::transform(_z_array->getConstPointer(),_z_array->getConstPointer()+_z_array->getNbOfElems(),
 714                    _z_array->getPointer(),std::bind2nd(std::plus<double>(),vector[2]));
 715 }
 716
 717 /*!
 718  * Applies scaling transformation to all nodes of \a this mesh.
 719  *  \param [in] point - coordinates of a scaling center. This array is to be of
 720  *         size \a this->getSpaceDimension() at least.
 721  *  \param [in] factor - a scale factor.
 722  */
 723 void MEDCouplingCMesh::scale(const double *point, double factor)
 724 {
 725   for(int i=0;i<3;i++)
 726     {
 727       DataArrayDouble *c=getCoordsAt(i);
 728       if(c)
 729         {
 730           double *coords=c->getPointer();
 731           int lgth=c->getNbOfElems();
 732           std::transform(coords,coords+lgth,coords,std::bind2nd(std::minus<double>(),point[i]));
 733           std::transform(coords,coords+lgth,coords,std::bind2nd(std::multiplies<double>(),factor));
 734           std::transform(coords,coords+lgth,coords,std::bind2nd(std::plus<double>(),point[i]));
 735           c->declareAsNew();
 736         }
 737     }
 738   updateTime();
 739 }
 740
 741 MEDCouplingMesh *MEDCouplingCMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
 742 {
 743   //not implemented yet !
 744   return 0;
 745 }
 746
 747 /*!
 748  * Returns a new DataArrayDouble holding coordinates of all nodes of \a this mesh.
 749  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
 750  *          this->getNumberOfNodes() tuples per \a this->getSpaceDimension()
 751  *          components. The caller is to delete this array using decrRef() as it is
 752  *          no more needed.
 753  */
 754 DataArrayDouble *MEDCouplingCMesh::getCoordinatesAndOwner() const
 755 {
 756   DataArrayDouble *ret=DataArrayDouble::New();
 757   int spaceDim=getSpaceDimension();
 758   int nbNodes=getNumberOfNodes();
 759   ret->alloc(nbNodes,spaceDim);
 760   double *pt=ret->getPointer();
 761   int tmp[3];
 762   getSplitNodeValues(tmp);
 763   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 764   const double *tabsPtr[3];
 765   for(int j=0;j<spaceDim;j++)
 766     {
 767       tabsPtr[j]=tabs[j]->getConstPointer();
 768       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0));
 769     }
 770   int tmp2[3];
 771   for(int i=0;i<nbNodes;i++)
 772     {
 773       GetPosFromId(i,spaceDim,tmp,tmp2);
 774       for(int j=0;j<spaceDim;j++)
 775         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 776     }
 777   return ret;
 778 }
 779
 780 /*!
 781  * Returns a new DataArrayDouble holding barycenters of all cells. The barycenter is
 782  * computed by averaging coordinates of cell nodes.
 783  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
 784  *          this->getNumberOfCells() tuples per \a this->getSpaceDimension()
 785  *          components. The caller is to delete this array using decrRef() as it is
 786  *          no more needed.
 787  */
 788 DataArrayDouble *MEDCouplingCMesh::getBarycenterAndOwner() const
 789 {
 790   DataArrayDouble *ret=DataArrayDouble::New();
 791   int spaceDim=getSpaceDimension();
 792   int nbCells=getNumberOfCells();
 793   ret->alloc(nbCells,spaceDim);
 794   double *pt=ret->getPointer();
 795   int tmp[3];
 796   getSplitCellValues(tmp);
 797   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 798   std::vector<double> tabsPtr[3];
 799   for(int j=0;j<spaceDim;j++)
 800     {
 801       int sz=tabs[j]->getNbOfElems()-1;
 802       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0));
 803       const double *srcPtr=tabs[j]->getConstPointer();
 804       tabsPtr[j].insert(tabsPtr[j].end(),srcPtr,srcPtr+sz);
 805       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),srcPtr+1,tabsPtr[j].begin(),std::plus<double>());
 806       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),tabsPtr[j].begin(),std::bind2nd(std::multiplies<double>(),0.5));
 807     }
 808   int tmp2[3];
 809   for(int i=0;i<nbCells;i++)
 810     {
 811       GetPosFromId(i,spaceDim,tmp,tmp2);
 812       for(int j=0;j<spaceDim;j++)
 813         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 814     }
 815   return ret;
 816 }
 817
 818 DataArrayDouble *MEDCouplingCMesh::computeIsoBarycenterOfNodesPerCell() const
 819 {
 820   return MEDCouplingCMesh::getBarycenterAndOwner();
 821 }
 822
 823 void MEDCouplingCMesh::renumberCells(const int *old2NewBg, bool check)
 824 {
 825   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CMesh !");
 826 }
 827
 828 void MEDCouplingCMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 829 {
 830   int it,order;
 831   double time=getTime(it,order);
 832   tinyInfo.clear();
 833   tinyInfoD.clear();
 834   littleStrings.clear();
 835   littleStrings.push_back(getName());
 836   littleStrings.push_back(getDescription());
 837   littleStrings.push_back(getTimeUnit());
 838   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 839   for(int i=0;i<3;i++)
 840     {
 841       int val=-1;
 842       std::string st;
 843       if(thisArr[i])
 844         {
 845           val=thisArr[i]->getNumberOfTuples();
 846           st=thisArr[i]->getInfoOnComponent(0);
 847         }
 848       tinyInfo.push_back(val);
 849       littleStrings.push_back(st);
 850     }
 851   tinyInfo.push_back(it);
 852   tinyInfo.push_back(order);
 853   tinyInfoD.push_back(time);
 854 }
 855
 856 void MEDCouplingCMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 857 {
 858   a1->alloc(0,1);
 859   int sum=0;
 860   for(int i=0;i<3;i++)
 861     if(tinyInfo[i]!=-1)
 862       sum+=tinyInfo[i];
 863   a2->alloc(sum,1);
 864 }
 865
 866 void MEDCouplingCMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
 867 {
 868   a1=DataArrayInt::New();
 869   a1->alloc(0,1);
 870   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 871   int sz=0;
 872   for(int i=0;i<3;i++)
 873     {
 874       if(thisArr[i])
 875         sz+=thisArr[i]->getNumberOfTuples();
 876     }
 877   a2=DataArrayDouble::New();
 878   a2->alloc(sz,1);
 879   double *a2Ptr=a2->getPointer();
 880   for(int i=0;i<3;i++)
 881     if(thisArr[i])
 882       a2Ptr=std::copy(thisArr[i]->getConstPointer(),thisArr[i]->getConstPointer()+thisArr[i]->getNumberOfTuples(),a2Ptr);
 883 }
 884
 885 void MEDCouplingCMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
 886                                        const std::vector<std::string>& littleStrings)
 887 {
 888   setName(littleStrings[0]);
 889   setDescription(littleStrings[1]);
 890   setTimeUnit(littleStrings[2]);
 891   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
 892   const double *data=a2->getConstPointer();
 893   for(int i=0;i<3;i++)
 894     {
 895       if(tinyInfo[i]!=-1)
 896         {
 897           (*(thisArr[i]))=DataArrayDouble::New();
 898           (*(thisArr[i]))->alloc(tinyInfo[i],1);
 899           (*(thisArr[i]))->setInfoOnComponent(0,littleStrings[i+3]);
 900           std::copy(data,data+tinyInfo[i],(*(thisArr[i]))->getPointer());
 901           data+=tinyInfo[i];
 902         }
 903     }
 904   setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 905 }
 906
 907 void MEDCouplingCMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
 908 {
 909   std::ostringstream extent;
 910   DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 911   for(int i=0;i<3;i++)
 912     {
 913       if(thisArr[i])
 914         { extent << "0 " <<  thisArr[i]->getNumberOfTuples()-1 << " "; }
 915       else
 916         { extent << "0 0 "; }
 917     }
 918   ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\">\n";
 919   ofs << "    <Piece Extent=\"" << extent.str() << "\">\n";
 920   ofs << "      <PointData>\n" << pointData << std::endl;
 921   ofs << "      </PointData>\n";
 922   ofs << "      <CellData>\n" << cellData << std::endl;
 923   ofs << "      </CellData>\n";
 924   ofs << "      <Coordinates>\n";
 925   for(int i=0;i<3;i++)
 926     {
 927       if(thisArr[i])
 928         thisArr[i]->writeVTK(ofs,8,"Array",byteData);
 929       else
 930         {
 931           MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=DataArrayDouble::New(); coo->alloc(1,1);
 932           coo->setIJ(0,0,0.);
 933           coo->writeVTK(ofs,8,"Array",byteData);
 934         }
 935     }
 936   ofs << "      </Coordinates>\n";
 937   ofs << "    </Piece>\n";
 938   ofs << "  </" << getVTKDataSetType() << ">\n";
 939 }
 940
 941 void MEDCouplingCMesh::reprQuickOverview(std::ostream& stream) const
 942 {
 943   stream << "MEDCouplingCMesh C++ instance at " << this << ". Name : \"" << getName() << "\".";
 944   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 945   std::ostringstream stream2[3];
 946   bool isDef[3];
 947   int nbOfCells=1,nbOfNodes=1;
 948   for(int i=0;i<3;i++)
 949     {
 950       isDef[i]=thisArr[i]!=0;
 951       if(isDef[i])
 952         {
 953           char tmp='X'+i;
 954           stream2[i] << tmp << " positions array ";
 955           if(!thisArr[i]->isAllocated())
 956             stream2[i] << "set but not allocated.";
 957           else
 958             {
 959               int nbCompo=thisArr[i]->getNumberOfComponents();
 960               if(nbCompo==1)
 961                 {
 962                   int nbTuples=thisArr[i]->getNumberOfTuples();
 963                   if(nbTuples<1)
 964                     { stream2[i] << "set and allocated - WARNING number of elements < 1 !"; nbOfCells=-1; nbOfNodes=-1; }
 965                   else
 966                     {
 967                       stream2[i] << "(length=" << nbTuples << ")" << ": ";
 968                       thisArr[i]->reprQuickOverviewData(stream2[i],200);
 969                       if(nbOfCells!=-1)
 970                         { nbOfNodes*=nbTuples; nbOfCells*=nbTuples-1; }
 971                     }
 972                 }
 973               else
 974                 { stream2[i] << "set and allocated - WARNING number of components != 1 !"; nbOfCells=-1; nbOfNodes=-1; }
 975             }
 976         }
 977     }
 978   if(!isDef[0] && !isDef[1] && !isDef[2])
 979     { stream << " No arrays set !"; return; }
 980   if(nbOfCells>=0)
 981     { stream << std::endl << "Number of cells : " << nbOfCells << ". Number of nodes : " << nbOfNodes << "."; }
 982   for(int i=0;i<3;i++)
 983     {
 984       if(isDef[i])
 985         stream << std::endl << stream2[i].str();
 986     }
 987 }
 988
 989 std::string MEDCouplingCMesh::getVTKDataSetType() const
 990 {
 991   return std::string("RectilinearGrid");
 992 }