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 throw(INTERP_KERNEL::Exception)
 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 throw(INTERP_KERNEL::Exception)
 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     }
 248   if(_z_array)
 249     {
 250       if(_z_array->getNbOfElems()<2)
 251         {
 252           std::ostringstream os; os << msg0 << 'Z' << msg1;
 253           throw INTERP_KERNEL::Exception(os.str().c_str());
 254         }
 255       if(_z_array->getNumberOfComponents()!=1)
 256         {
 257           std::ostringstream os; os << msg0 << 'Z' << msg2;
 258           throw INTERP_KERNEL::Exception(os.str().c_str());
 259         }
 260     }
 261 }
 262
 263 void MEDCouplingCMesh::checkCoherency1(double eps) const
 264 {
 265   checkCoherency();
 266   if(_x_array)
 267     _x_array->checkMonotonic(true, eps);
 268   if(_y_array)
 269     _y_array->checkMonotonic(true, eps);
 270   if(_z_array)
 271     _z_array->checkMonotonic(true, eps);
 272 }
 273
 274 void MEDCouplingCMesh::checkCoherency2(double eps) const
 275 {
 276   checkCoherency1(eps);
 277 }
 278
 279 int MEDCouplingCMesh::getNumberOfCells() const
 280 {
 281   int ret=1;
 282   if(_x_array)
 283     ret*=_x_array->getNbOfElems()-1;
 284   if(_y_array)
 285     ret*=_y_array->getNbOfElems()-1;
 286   if(_z_array)
 287     ret*=_z_array->getNbOfElems()-1;
 288   return ret;
 289 }
 290
 291 int MEDCouplingCMesh::getNumberOfNodes() const
 292 {
 293   int ret=1;
 294   if(_x_array)
 295     ret*=_x_array->getNbOfElems();
 296   if(_y_array)
 297     ret*=_y_array->getNbOfElems();
 298   if(_z_array)
 299     ret*=_z_array->getNbOfElems();
 300   return ret;
 301 }
 302
 303 void MEDCouplingCMesh::getSplitCellValues(int *res) const
 304 {
 305   int meshDim(getMeshDimension());
 306   for(int l=0;l<meshDim;l++)
 307     {
 308       int val=1;
 309       for(int p=0;p<meshDim-l-1;p++)
 310         val*=getCoordsAt(p)->getNbOfElems()-1;
 311       res[meshDim-l-1]=val;
 312     }
 313 }
 314
 315 void MEDCouplingCMesh::getSplitNodeValues(int *res) const
 316 {
 317   int spaceDim(getSpaceDimension());
 318   for(int l=0;l<spaceDim;l++)
 319     {
 320       int val=1;
 321       for(int p=0;p<spaceDim-l-1;p++)
 322         val*=getCoordsAt(p)->getNbOfElems();
 323       res[spaceDim-l-1]=val;
 324     }
 325 }
 326
 327 void MEDCouplingCMesh::getNodeGridStructure(int *res) const
 328 {
 329   int spaceDim(getSpaceDimension());
 330   for(int i=0;i<spaceDim;i++)
 331     res[i]=getCoordsAt(i)->getNbOfElems();
 332 }
 333
 334 std::vector<int> MEDCouplingCMesh::getNodeGridStructure() const
 335 {
 336   std::vector<int> ret(getSpaceDimension());
 337   getNodeGridStructure(&ret[0]);
 338   return ret;
 339 }
 340
 341 MEDCouplingStructuredMesh *MEDCouplingCMesh::buildStructuredSubPart(const std::vector< std::pair<int,int> >& cellPart) const
 342 {
 343   checkCoherency();
 344   int dim(getMeshDimension());
 345   if(dim!=(int)cellPart.size())
 346     {
 347       std::ostringstream oss; oss << "MEDCouplingCMesh::buildStructuredSubPart : the mesh dimension is " << dim << " and cell part size is " << cellPart.size() << " !";
 348       throw INTERP_KERNEL::Exception(oss.str().c_str());
 349     }
 350   MEDCouplingAutoRefCountObjectPtr<MEDCouplingCMesh> ret(dynamic_cast<MEDCouplingCMesh *>(deepCpy()));
 351   for(int i=0;i<dim;i++)
 352     {
 353       MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> tmp(ret->getCoordsAt(i)->selectByTupleId2(cellPart[i].first,cellPart[i].second+1,1));
 354       ret->setCoordsAt(i,tmp);
 355     }
 356   return ret.retn();
 357 }
 358
 359 /*!
 360  * Return the space dimension of \a this. It only considers the arrays along X, Y and Z to deduce that.
 361  * 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.
 362  */
 363 int MEDCouplingCMesh::getSpaceDimension() const
 364 {
 365   static const char MSG[]="MEDCouplingCMesh::getSpaceDimension : mesh is invalid ! null vectors (X, Y or Z) must be put contiguously at the end !";
 366   int ret(0);
 367   bool isOK(true);
 368   if(_x_array)
 369     ret++;
 370   else
 371     isOK=false;
 372   if(_y_array)
 373     {
 374       if(!isOK)
 375         throw INTERP_KERNEL::Exception(MSG);
 376       ret++;
 377     }
 378   else
 379     isOK=false;
 380   if(_z_array)
 381     {
 382       if(!isOK)
 383         throw INTERP_KERNEL::Exception(MSG);
 384       ret++;
 385     }
 386   return ret;
 387 }
 388
 389 /*!
 390  * 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.
 391  */
 392 int MEDCouplingCMesh::getMeshDimension() const
 393 {
 394   int ret(getSpaceDimension());
 395   if(_x_array)
 396     {
 397       if(_x_array->isAllocated())
 398         if(_x_array->getNumberOfTuples()==1)
 399           ret--;
 400     }
 401   if(_y_array)
 402     {
 403       if(_y_array->isAllocated())
 404         if(_y_array->getNumberOfTuples()==1)
 405           ret--;
 406     }
 407   if(_z_array)
 408     {
 409       if(_z_array->isAllocated())
 410         if(_z_array->getNumberOfTuples()==1)
 411           ret--;
 412     }
 413   return ret;
 414 }
 415
 416 void MEDCouplingCMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const
 417 {
 418   int tmp[3];
 419   int spaceDim=getSpaceDimension();
 420   getSplitNodeValues(tmp);
 421   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 422   int tmp2[3];
 423   GetPosFromId(nodeId,spaceDim,tmp,tmp2);
 424   for(int j=0;j<spaceDim;j++)
 425     if(tabs[j])
 426       coo.push_back(tabs[j]->getConstPointer()[tmp2[j]]);
 427 }
 428
 429 std::string MEDCouplingCMesh::simpleRepr() const
 430 {
 431   std::ostringstream ret;
 432   ret << "Cartesian mesh with name : \"" << getName() << "\"\n";
 433   ret << "Description of mesh : \"" << getDescription() << "\"\n";
 434   int tmpp1,tmpp2;
 435   double tt=getTime(tmpp1,tmpp2);
 436   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
 437   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
 438   ret << "Mesh and SpaceDimension dimension : " << getSpaceDimension() << "\n\nArrays :\n________\n\n";
 439   if(_x_array)
 440     {
 441       ret << "X Array :\n";
 442       _x_array->reprZipWithoutNameStream(ret);
 443     }
 444   if(_y_array)
 445     {
 446       ret << "Y Array :\n";
 447       _y_array->reprZipWithoutNameStream(ret);
 448     }
 449   if(_z_array)
 450     {
 451       ret << "Z Array :\n";
 452       _z_array->reprZipWithoutNameStream(ret);
 453     }
 454   return ret.str();
 455 }
 456
 457 std::string MEDCouplingCMesh::advancedRepr() const
 458 {
 459   return simpleRepr();
 460 }
 461
 462 /*!
 463  * Returns a DataArrayDouble holding positions of nodes along a given axis.
 464  * For more info on Cartesian meshes, see \ref MEDCouplingCMeshPage.
 465  *  \param [in] i - an index of axis, a value from [0,1,2].
 466  *  \return const DataArrayDouble * - a pointer to the data array of node coordinates
 467  *         referred by \a this mesh.
 468  *  \throw If \a i is not one of [0,1,2].
 469  *
 470  *  \ref cpp_mccmesh_getCoordsAt "Here is a C++ example".<br>
 471  *  \ref  py_mccmesh_getCoordsAt "Here is a Python example".
 472  */
 473 const DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i) const
 474 {
 475   switch(i)
 476     {
 477     case 0:
 478       return _x_array;
 479     case 1:
 480       return _y_array;
 481     case 2:
 482       return _z_array;
 483     default:
 484       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 485     }
 486 }
 487
 488 /*!
 489  * Returns a DataArrayDouble holding positions of nodes along a given axis.
 490  * For more info on Cartesian meshes, see \ref MEDCouplingCMeshPage.
 491  *  \param [in] i - an index of axis, a value from [0,1,2].
 492  *  \return const DataArrayDouble * - a pointer to the data array of node coordinates
 493  *         referred by \a this mesh.
 494  *  \throw If \a i is not one of [0,1,2].
 495  *
 496  *  \ref cpp_mccmesh_getCoordsAt "Here is a C++ example".<br>
 497  *  \ref  py_mccmesh_getCoordsAt "Here is a Python example".
 498  */
 499 DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i)
 500 {
 501   switch(i)
 502     {
 503     case 0:
 504       return _x_array;
 505     case 1:
 506       return _y_array;
 507     case 2:
 508       return _z_array;
 509     default:
 510       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 511     }
 512 }
 513
 514 /*!
 515  * Sets node coordinates along a given axis. For more info on Cartesian meshes, see
 516  * \ref MEDCouplingCMeshPage.
 517  *  \param [in] i - an index of axis, a value in range [0,1,2].
 518  *  \param [in] arr - DataArrayDouble holding positions of nodes along the i-th
 519  *         axis. It must be an array of one component.
 520  *  \throw If \a arr->getNumberOfComponents() != 1.
 521  *  \throw If \a i is not one of [0,1,2].
 522  *
 523  *  \ref medcouplingcppexamplesCmeshStdBuild1 "Here is a C++ example".<br>
 524  *  \ref  medcouplingpyexamplesCmeshStdBuild1 "Here is a Python example".
 525  */
 526 void MEDCouplingCMesh::setCoordsAt(int i, const DataArrayDouble *arr)
 527 {
 528   if(arr)
 529     arr->checkNbOfComps(1,"MEDCouplingCMesh::setCoordsAt");
 530   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
 531   if(i<0 || i>2)
 532     throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 533   if(arr!=*(thisArr[i]))
 534     {
 535       if(*(thisArr[i]))
 536         (*(thisArr[i]))->decrRef();
 537       (*(thisArr[i]))=const_cast<DataArrayDouble *>(arr);
 538       if(*(thisArr[i]))
 539         (*(thisArr[i]))->incrRef();
 540       declareAsNew();
 541     }
 542 }
 543
 544 /*!
 545  * Sets node coordinates along some of the tree axes. This method updates all the
 546  * three node coordinates arrays at once. For more info on Cartesian meshes, see
 547  * \ref MEDCouplingCMeshPage.
 548  *  \param [in] coordsX - DataArrayDouble holding positions of nodes along the X
 549  *         axis. It must be an array of one component or \c NULL.
 550  *  \param [in] coordsY - DataArrayDouble holding positions of nodes along the Y
 551  *         axis. It must be an array of one component or \c NULL.
 552  *  \param [in] coordsZ - DataArrayDouble holding positions of nodes along the Z
 553  *         axis. It must be an array of one component or \c NULL.
 554  *  \throw If \a coords*->getNumberOfComponents() != 1.
 555  *
 556  *  \ref medcouplingcppexamplesCmeshStdBuild1 "Here is a C++ example".<br>
 557  *  \ref  medcouplingpyexamplesCmeshStdBuild1 "Here is a Python example".
 558  */
 559 void MEDCouplingCMesh::setCoords(const DataArrayDouble *coordsX, const DataArrayDouble *coordsY, const DataArrayDouble *coordsZ)
 560 {
 561   if(coordsX)
 562     coordsX->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsX");
 563   if(coordsY)
 564     coordsY->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsY");
 565   if(coordsZ)
 566     coordsZ->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsZ");
 567   if(_x_array)
 568     _x_array->decrRef();
 569   _x_array=const_cast<DataArrayDouble *>(coordsX);
 570   if(_x_array)
 571     _x_array->incrRef();
 572   if(_y_array)
 573     _y_array->decrRef();
 574   _y_array=const_cast<DataArrayDouble *>(coordsY);
 575   if(_y_array)
 576     _y_array->incrRef();
 577   if(_z_array)
 578     _z_array->decrRef();
 579   _z_array=const_cast<DataArrayDouble *>(coordsZ);
 580   if(_z_array)
 581     _z_array->incrRef();
 582   declareAsNew();
 583 }
 584
 585 void MEDCouplingCMesh::getBoundingBox(double *bbox) const
 586 {
 587   int dim=getSpaceDimension();
 588   int j=0;
 589   for (int idim=0; idim<dim; idim++)
 590     {
 591       const DataArrayDouble *c=getCoordsAt(idim);
 592       if(c)
 593         {
 594           const double *coords=c->getConstPointer();
 595           int nb=c->getNbOfElems();
 596           bbox[2*j]=coords[0];
 597           bbox[2*j+1]=coords[nb-1];
 598           j++;
 599         }
 600     }
 601 }
 602
 603 /*!
 604  * Returns a new MEDCouplingFieldDouble containing volumes of cells constituting \a this
 605  * mesh.<br>
 606  * For 1D cells, the returned field contains lengths.<br>
 607  * For 2D cells, the returned field contains areas.<br>
 608  * For 3D cells, the returned field contains volumes.
 609  *  \param [in] isAbs - a not used parameter.
 610  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on cells
 611  *         and one time . The caller is to delete this field using decrRef() as it is no
 612  *         more needed.
 613  */
 614 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const
 615 {
 616   std::string name="MeasureOfMesh_";
 617   name+=getName();
 618   int nbelem=getNumberOfCells();
 619   MEDCouplingFieldDouble *field=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
 620   field->setName(name);
 621   DataArrayDouble* array=DataArrayDouble::New();
 622   array->alloc(nbelem,1);
 623   double *area_vol=array->getPointer();
 624   field->setArray(array) ;
 625   array->decrRef();
 626   field->setMesh(const_cast<MEDCouplingCMesh *>(this));
 627   field->synchronizeTimeWithMesh();
 628   int tmp[3];
 629   getSplitCellValues(tmp);
 630   int dim=getSpaceDimension();
 631   const double **thisArr=new const double *[dim];
 632   const DataArrayDouble *thisArr2[3]={_x_array,_y_array,_z_array};
 633   for(int i=0;i<dim;i++)
 634     thisArr[i]=thisArr2[i]->getConstPointer();
 635   for(int icell=0;icell<nbelem;icell++)
 636     {
 637       int tmp2[3];
 638       GetPosFromId(icell,dim,tmp,tmp2);
 639       area_vol[icell]=1.;
 640       for(int i=0;i<dim;i++)
 641         area_vol[icell]*=thisArr[i][tmp2[i]+1]-thisArr[i][tmp2[i]];
 642     }
 643   delete [] thisArr;
 644   return field;
 645 }
 646
 647 /*!
 648  * not implemented yet !
 649  */
 650 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureFieldOnNode(bool isAbs) const
 651 {
 652   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getMeasureFieldOnNode : not implemented yet !");
 653   //return 0;
 654 }
 655
 656 int MEDCouplingCMesh::getCellContainingPoint(const double *pos, double eps) const
 657 {
 658   int dim=getSpaceDimension();
 659   int ret=0;
 660   int coeff=1;
 661   for(int i=0;i<dim;i++)
 662     {
 663       const double *d=getCoordsAt(i)->getConstPointer();
 664       int nbOfNodes=getCoordsAt(i)->getNbOfElems();
 665       double ref=pos[i];
 666       const double *w=std::find_if(d,d+nbOfNodes,std::bind2nd(std::greater_equal<double>(),ref));
 667       int w2=(int)std::distance(d,w);
 668       if(w2<nbOfNodes)
 669         {
 670           if(w2==0)
 671             {
 672               if(ref>d[0]-eps)
 673                 w2=1;
 674               else
 675                 return -1;
 676             }
 677           ret+=coeff*(w2-1);
 678           coeff*=nbOfNodes-1;
 679         }
 680       else
 681         return -1;
 682     }
 683   return ret;
 684 }
 685
 686 void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
 687 {
 688   throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !");
 689 }
 690
 691 /*!
 692  * Translates all nodes of \a this mesh by a given vector. Actually, it adds each
 693  * component of the \a vector to all node coordinates of a corresponding axis.
 694  *  \param [in] vector - the translation vector whose size must be not less than \a
 695  *         this->getSpaceDimension().
 696  */
 697 void MEDCouplingCMesh::translate(const double *vector)
 698 {
 699   if(_x_array)
 700     std::transform(_x_array->getConstPointer(),_x_array->getConstPointer()+_x_array->getNbOfElems(),
 701                    _x_array->getPointer(),std::bind2nd(std::plus<double>(),vector[0]));
 702   if(_y_array)
 703     std::transform(_y_array->getConstPointer(),_y_array->getConstPointer()+_y_array->getNbOfElems(),
 704                    _y_array->getPointer(),std::bind2nd(std::plus<double>(),vector[1]));
 705   if(_z_array)
 706     std::transform(_z_array->getConstPointer(),_z_array->getConstPointer()+_z_array->getNbOfElems(),
 707                    _z_array->getPointer(),std::bind2nd(std::plus<double>(),vector[2]));
 708 }
 709
 710 /*!
 711  * Applies scaling transformation to all nodes of \a this mesh.
 712  *  \param [in] point - coordinates of a scaling center. This array is to be of
 713  *         size \a this->getSpaceDimension() at least.
 714  *  \param [in] factor - a scale factor.
 715  */
 716 void MEDCouplingCMesh::scale(const double *point, double factor)
 717 {
 718   for(int i=0;i<3;i++)
 719     {
 720       DataArrayDouble *c=getCoordsAt(i);
 721       if(c)
 722         {
 723           double *coords=c->getPointer();
 724           int lgth=c->getNbOfElems();
 725           std::transform(coords,coords+lgth,coords,std::bind2nd(std::minus<double>(),point[i]));
 726           std::transform(coords,coords+lgth,coords,std::bind2nd(std::multiplies<double>(),factor));
 727           std::transform(coords,coords+lgth,coords,std::bind2nd(std::plus<double>(),point[i]));
 728           c->declareAsNew();
 729         }
 730     }
 731   updateTime();
 732 }
 733
 734 MEDCouplingMesh *MEDCouplingCMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
 735 {
 736   //not implemented yet !
 737   return 0;
 738 }
 739
 740 /*!
 741  * Returns a new DataArrayDouble holding coordinates of all nodes of \a this mesh.
 742  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
 743  *          this->getNumberOfNodes() tuples per \a this->getSpaceDimension()
 744  *          components. The caller is to delete this array using decrRef() as it is
 745  *          no more needed.
 746  */
 747 DataArrayDouble *MEDCouplingCMesh::getCoordinatesAndOwner() const
 748 {
 749   DataArrayDouble *ret=DataArrayDouble::New();
 750   int spaceDim=getSpaceDimension();
 751   int nbNodes=getNumberOfNodes();
 752   ret->alloc(nbNodes,spaceDim);
 753   double *pt=ret->getPointer();
 754   int tmp[3];
 755   getSplitNodeValues(tmp);
 756   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 757   const double *tabsPtr[3];
 758   for(int j=0;j<spaceDim;j++)
 759     {
 760       tabsPtr[j]=tabs[j]->getConstPointer();
 761       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0));
 762     }
 763   int tmp2[3];
 764   for(int i=0;i<nbNodes;i++)
 765     {
 766       GetPosFromId(i,spaceDim,tmp,tmp2);
 767       for(int j=0;j<spaceDim;j++)
 768         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 769     }
 770   return ret;
 771 }
 772
 773 /*!
 774  * Returns a new DataArrayDouble holding barycenters of all cells. The barycenter is
 775  * computed by averaging coordinates of cell nodes.
 776  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
 777  *          this->getNumberOfCells() tuples per \a this->getSpaceDimension()
 778  *          components. The caller is to delete this array using decrRef() as it is
 779  *          no more needed.
 780  */
 781 DataArrayDouble *MEDCouplingCMesh::getBarycenterAndOwner() const
 782 {
 783   DataArrayDouble *ret=DataArrayDouble::New();
 784   int spaceDim=getSpaceDimension();
 785   int nbCells=getNumberOfCells();
 786   ret->alloc(nbCells,spaceDim);
 787   double *pt=ret->getPointer();
 788   int tmp[3];
 789   getSplitCellValues(tmp);
 790   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 791   std::vector<double> tabsPtr[3];
 792   for(int j=0;j<spaceDim;j++)
 793     {
 794       int sz=tabs[j]->getNbOfElems()-1;
 795       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0));
 796       const double *srcPtr=tabs[j]->getConstPointer();
 797       tabsPtr[j].insert(tabsPtr[j].end(),srcPtr,srcPtr+sz);
 798       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),srcPtr+1,tabsPtr[j].begin(),std::plus<double>());
 799       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),tabsPtr[j].begin(),std::bind2nd(std::multiplies<double>(),0.5));
 800     }
 801   int tmp2[3];
 802   for(int i=0;i<nbCells;i++)
 803     {
 804       GetPosFromId(i,spaceDim,tmp,tmp2);
 805       for(int j=0;j<spaceDim;j++)
 806         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 807     }
 808   return ret;
 809 }
 810
 811 DataArrayDouble *MEDCouplingCMesh::computeIsoBarycenterOfNodesPerCell() const
 812 {
 813   return MEDCouplingCMesh::getBarycenterAndOwner();
 814 }
 815
 816 void MEDCouplingCMesh::renumberCells(const int *old2NewBg, bool check)
 817 {
 818   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CMesh !");
 819 }
 820
 821 void MEDCouplingCMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 822 {
 823   int it,order;
 824   double time=getTime(it,order);
 825   tinyInfo.clear();
 826   tinyInfoD.clear();
 827   littleStrings.clear();
 828   littleStrings.push_back(getName());
 829   littleStrings.push_back(getDescription());
 830   littleStrings.push_back(getTimeUnit());
 831   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 832   for(int i=0;i<3;i++)
 833     {
 834       int val=-1;
 835       std::string st;
 836       if(thisArr[i])
 837         {
 838           val=thisArr[i]->getNumberOfTuples();
 839           st=thisArr[i]->getInfoOnComponent(0);
 840         }
 841       tinyInfo.push_back(val);
 842       littleStrings.push_back(st);
 843     }
 844   tinyInfo.push_back(it);
 845   tinyInfo.push_back(order);
 846   tinyInfoD.push_back(time);
 847 }
 848
 849 void MEDCouplingCMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 850 {
 851   a1->alloc(0,1);
 852   int sum=0;
 853   for(int i=0;i<3;i++)
 854     if(tinyInfo[i]!=-1)
 855       sum+=tinyInfo[i];
 856   a2->alloc(sum,1);
 857 }
 858
 859 void MEDCouplingCMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
 860 {
 861   a1=DataArrayInt::New();
 862   a1->alloc(0,1);
 863   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 864   int sz=0;
 865   for(int i=0;i<3;i++)
 866     {
 867       if(thisArr[i])
 868         sz+=thisArr[i]->getNumberOfTuples();
 869     }
 870   a2=DataArrayDouble::New();
 871   a2->alloc(sz,1);
 872   double *a2Ptr=a2->getPointer();
 873   for(int i=0;i<3;i++)
 874     if(thisArr[i])
 875       a2Ptr=std::copy(thisArr[i]->getConstPointer(),thisArr[i]->getConstPointer()+thisArr[i]->getNumberOfTuples(),a2Ptr);
 876 }
 877
 878 void MEDCouplingCMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
 879                                        const std::vector<std::string>& littleStrings)
 880 {
 881   setName(littleStrings[0]);
 882   setDescription(littleStrings[1]);
 883   setTimeUnit(littleStrings[2]);
 884   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
 885   const double *data=a2->getConstPointer();
 886   for(int i=0;i<3;i++)
 887     {
 888       if(tinyInfo[i]!=-1)
 889         {
 890           (*(thisArr[i]))=DataArrayDouble::New();
 891           (*(thisArr[i]))->alloc(tinyInfo[i],1);
 892           (*(thisArr[i]))->setInfoOnComponent(0,littleStrings[i+3]);
 893           std::copy(data,data+tinyInfo[i],(*(thisArr[i]))->getPointer());
 894           data+=tinyInfo[i];
 895         }
 896     }
 897   setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 898 }
 899
 900 void MEDCouplingCMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
 901 {
 902   std::ostringstream extent;
 903   DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 904   for(int i=0;i<3;i++)
 905     {
 906       if(thisArr[i])
 907         { extent << "0 " <<  thisArr[i]->getNumberOfTuples()-1 << " "; }
 908       else
 909         { extent << "0 0 "; }
 910     }
 911   ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\">\n";
 912   ofs << "    <Piece Extent=\"" << extent.str() << "\">\n";
 913   ofs << "      <PointData>\n" << pointData << std::endl;
 914   ofs << "      </PointData>\n";
 915   ofs << "      <CellData>\n" << cellData << std::endl;
 916   ofs << "      </CellData>\n";
 917   ofs << "      <Coordinates>\n";
 918   for(int i=0;i<3;i++)
 919     {
 920       if(thisArr[i])
 921         thisArr[i]->writeVTK(ofs,8,"Array",byteData);
 922       else
 923         {
 924           MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=DataArrayDouble::New(); coo->alloc(1,1);
 925           coo->setIJ(0,0,0.);
 926           coo->writeVTK(ofs,8,"Array",byteData);
 927         }
 928     }
 929   ofs << "      </Coordinates>\n";
 930   ofs << "    </Piece>\n";
 931   ofs << "  </" << getVTKDataSetType() << ">\n";
 932 }
 933
 934 void MEDCouplingCMesh::reprQuickOverview(std::ostream& stream) const
 935 {
 936   stream << "MEDCouplingCMesh C++ instance at " << this << ". Name : \"" << getName() << "\".";
 937   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 938   std::ostringstream stream2[3];
 939   bool isDef[3];
 940   int nbOfCells=1,nbOfNodes=1;
 941   for(int i=0;i<3;i++)
 942     {
 943       isDef[i]=thisArr[i]!=0;
 944       if(isDef[i])
 945         {
 946           char tmp='X'+i;
 947           stream2[i] << tmp << " positions array ";
 948           if(!thisArr[i]->isAllocated())
 949             stream2[i] << "set but not allocated.";
 950           else
 951             {
 952               int nbCompo=thisArr[i]->getNumberOfComponents();
 953               if(nbCompo==1)
 954                 {
 955                   int nbTuples=thisArr[i]->getNumberOfTuples();
 956                   if(nbTuples<1)
 957                     { stream2[i] << "set and allocated - WARNING number of elements < 1 !"; nbOfCells=-1; nbOfNodes=-1; }
 958                   else
 959                     {
 960                       stream2[i] << "(length=" << nbTuples << ")" << ": ";
 961                       thisArr[i]->reprQuickOverviewData(stream2[i],200);
 962                       if(nbOfCells!=-1)
 963                         { nbOfNodes*=nbTuples; nbOfCells*=nbTuples-1; }
 964                     }
 965                 }
 966               else
 967                 { stream2[i] << "set and allocated - WARNING number of components != 1 !"; nbOfCells=-1; nbOfNodes=-1; }
 968             }
 969         }
 970     }
 971   if(!isDef[0] && !isDef[1] && !isDef[2])
 972     { stream << " No arrays set !"; return; }
 973   if(nbOfCells>=0)
 974     { stream << std::endl << "Number of cells : " << nbOfCells << ". Number of nodes : " << nbOfNodes << "."; }
 975   for(int i=0;i<3;i++)
 976     {
 977       if(isDef[i])
 978         stream << std::endl << stream2[i].str();
 979     }
 980
 981 }
 982
 983 std::string MEDCouplingCMesh::getVTKDataSetType() const
 984 {
 985   return std::string("RectilinearGrid");
 986 }