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  *  \if ENABLE_EXAMPLES
 471  *  \ref cpp_mccmesh_getCoordsAt "Here is a C++ example".<br>
 472  *  \ref  py_mccmesh_getCoordsAt "Here is a Python example".
 473  *  \endif
 474  */
 475 const DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i) const
 476 {
 477   switch(i)
 478     {
 479     case 0:
 480       return _x_array;
 481     case 1:
 482       return _y_array;
 483     case 2:
 484       return _z_array;
 485     default:
 486       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 487     }
 488 }
 489
 490 /*!
 491  * Returns a DataArrayDouble holding positions of nodes along a given axis.
 492  * For more info on Cartesian meshes, see \ref MEDCouplingCMeshPage.
 493  *  \param [in] i - an index of axis, a value from [0,1,2].
 494  *  \return const DataArrayDouble * - a pointer to the data array of node coordinates
 495  *         referred by \a this mesh.
 496  *  \throw If \a i is not one of [0,1,2].
 497  *
 498  *  \if ENABLE_EXAMPLES
 499  *  \ref cpp_mccmesh_getCoordsAt "Here is a C++ example".<br>
 500  *  \ref  py_mccmesh_getCoordsAt "Here is a Python example".
 501  *  \endif
 502  */
 503 DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i)
 504 {
 505   switch(i)
 506     {
 507     case 0:
 508       return _x_array;
 509     case 1:
 510       return _y_array;
 511     case 2:
 512       return _z_array;
 513     default:
 514       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 515     }
 516 }
 517
 518 /*!
 519  * Sets node coordinates along a given axis. For more info on Cartesian meshes, see
 520  * \ref MEDCouplingCMeshPage.
 521  *  \param [in] i - an index of axis, a value in range [0,1,2].
 522  *  \param [in] arr - DataArrayDouble holding positions of nodes along the i-th
 523  *         axis. It must be an array of one component.
 524  *  \throw If \a arr->getNumberOfComponents() != 1.
 525  *  \throw If \a i is not one of [0,1,2].
 526  *
 527  *  \if ENABLE_EXAMPLES
 528  *  \ref medcouplingcppexamplesCmeshStdBuild1 "Here is a C++ example".<br>
 529  *  \ref  medcouplingpyexamplesCmeshStdBuild1 "Here is a Python example".
 530  *  \endif
 531  */
 532 void MEDCouplingCMesh::setCoordsAt(int i, const DataArrayDouble *arr)
 533 {
 534   if(arr)
 535     arr->checkNbOfComps(1,"MEDCouplingCMesh::setCoordsAt");
 536   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
 537   if(i<0 || i>2)
 538     throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 539   if(arr!=*(thisArr[i]))
 540     {
 541       if(*(thisArr[i]))
 542         (*(thisArr[i]))->decrRef();
 543       (*(thisArr[i]))=const_cast<DataArrayDouble *>(arr);
 544       if(*(thisArr[i]))
 545         (*(thisArr[i]))->incrRef();
 546       declareAsNew();
 547     }
 548 }
 549
 550 /*!
 551  * Sets node coordinates along some of the tree axes. This method updates all the
 552  * three node coordinates arrays at once. For more info on Cartesian meshes, see
 553  * \ref MEDCouplingCMeshPage.
 554  *  \param [in] coordsX - DataArrayDouble holding positions of nodes along the X
 555  *         axis. It must be an array of one component or \c NULL.
 556  *  \param [in] coordsY - DataArrayDouble holding positions of nodes along the Y
 557  *         axis. It must be an array of one component or \c NULL.
 558  *  \param [in] coordsZ - DataArrayDouble holding positions of nodes along the Z
 559  *         axis. It must be an array of one component or \c NULL.
 560  *  \throw If \a coords*->getNumberOfComponents() != 1.
 561  *
 562  *  \if ENABLE_EXAMPLES
 563  *  \ref medcouplingcppexamplesCmeshStdBuild1 "Here is a C++ example".<br>
 564  *  \ref  medcouplingpyexamplesCmeshStdBuild1 "Here is a Python example".
 565  *  \endif
 566  */
 567 void MEDCouplingCMesh::setCoords(const DataArrayDouble *coordsX, const DataArrayDouble *coordsY, const DataArrayDouble *coordsZ)
 568 {
 569   if(coordsX)
 570     coordsX->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsX");
 571   if(coordsY)
 572     coordsY->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsY");
 573   if(coordsZ)
 574     coordsZ->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsZ");
 575   if(_x_array)
 576     _x_array->decrRef();
 577   _x_array=const_cast<DataArrayDouble *>(coordsX);
 578   if(_x_array)
 579     _x_array->incrRef();
 580   if(_y_array)
 581     _y_array->decrRef();
 582   _y_array=const_cast<DataArrayDouble *>(coordsY);
 583   if(_y_array)
 584     _y_array->incrRef();
 585   if(_z_array)
 586     _z_array->decrRef();
 587   _z_array=const_cast<DataArrayDouble *>(coordsZ);
 588   if(_z_array)
 589     _z_array->incrRef();
 590   declareAsNew();
 591 }
 592
 593 void MEDCouplingCMesh::getBoundingBox(double *bbox) const
 594 {
 595   int dim=getSpaceDimension();
 596   int j=0;
 597   for (int idim=0; idim<dim; idim++)
 598     {
 599       const DataArrayDouble *c=getCoordsAt(idim);
 600       if(c)
 601         {
 602           const double *coords=c->getConstPointer();
 603           int nb=c->getNbOfElems();
 604           bbox[2*j]=coords[0];
 605           bbox[2*j+1]=coords[nb-1];
 606           j++;
 607         }
 608     }
 609 }
 610
 611 /*!
 612  * Returns a new MEDCouplingFieldDouble containing volumes of cells constituting \a this
 613  * mesh.<br>
 614  * For 1D cells, the returned field contains lengths.<br>
 615  * For 2D cells, the returned field contains areas.<br>
 616  * For 3D cells, the returned field contains volumes.
 617  *  \param [in] isAbs - a not used parameter.
 618  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on cells
 619  *         and one time . The caller is to delete this field using decrRef() as it is no
 620  *         more needed.
 621  */
 622 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const
 623 {
 624   std::string name="MeasureOfMesh_";
 625   name+=getName();
 626   int nbelem=getNumberOfCells();
 627   MEDCouplingFieldDouble *field=MEDCouplingFieldDouble::New(ON_CELLS,ONE_TIME);
 628   field->setName(name);
 629   DataArrayDouble* array=DataArrayDouble::New();
 630   array->alloc(nbelem,1);
 631   double *area_vol=array->getPointer();
 632   field->setArray(array) ;
 633   array->decrRef();
 634   field->setMesh(const_cast<MEDCouplingCMesh *>(this));
 635   field->synchronizeTimeWithMesh();
 636   int tmp[3];
 637   getSplitCellValues(tmp);
 638   int dim=getSpaceDimension();
 639   const double **thisArr=new const double *[dim];
 640   const DataArrayDouble *thisArr2[3]={_x_array,_y_array,_z_array};
 641   for(int i=0;i<dim;i++)
 642     thisArr[i]=thisArr2[i]->getConstPointer();
 643   for(int icell=0;icell<nbelem;icell++)
 644     {
 645       int tmp2[3];
 646       GetPosFromId(icell,dim,tmp,tmp2);
 647       area_vol[icell]=1.;
 648       for(int i=0;i<dim;i++)
 649         area_vol[icell]*=thisArr[i][tmp2[i]+1]-thisArr[i][tmp2[i]];
 650     }
 651   delete [] thisArr;
 652   return field;
 653 }
 654
 655 /*!
 656  * not implemented yet !
 657  */
 658 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureFieldOnNode(bool isAbs) const
 659 {
 660   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getMeasureFieldOnNode : not implemented yet !");
 661   //return 0;
 662 }
 663
 664 int MEDCouplingCMesh::getCellContainingPoint(const double *pos, double eps) const
 665 {
 666   int dim=getSpaceDimension();
 667   int ret=0;
 668   int coeff=1;
 669   for(int i=0;i<dim;i++)
 670     {
 671       const double *d=getCoordsAt(i)->getConstPointer();
 672       int nbOfNodes=getCoordsAt(i)->getNbOfElems();
 673       double ref=pos[i];
 674       const double *w=std::find_if(d,d+nbOfNodes,std::bind2nd(std::greater_equal<double>(),ref));
 675       int w2=(int)std::distance(d,w);
 676       if(w2<nbOfNodes)
 677         {
 678           if(w2==0)
 679             {
 680               if(ref>d[0]-eps)
 681                 w2=1;
 682               else
 683                 return -1;
 684             }
 685           ret+=coeff*(w2-1);
 686           coeff*=nbOfNodes-1;
 687         }
 688       else
 689         return -1;
 690     }
 691   return ret;
 692 }
 693
 694 void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
 695 {
 696   throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !");
 697 }
 698
 699 /*!
 700  * Translates all nodes of \a this mesh by a given vector. Actually, it adds each
 701  * component of the \a vector to all node coordinates of a corresponding axis.
 702  *  \param [in] vector - the translation vector whose size must be not less than \a
 703  *         this->getSpaceDimension().
 704  */
 705 void MEDCouplingCMesh::translate(const double *vector)
 706 {
 707   if(_x_array)
 708     std::transform(_x_array->getConstPointer(),_x_array->getConstPointer()+_x_array->getNbOfElems(),
 709                    _x_array->getPointer(),std::bind2nd(std::plus<double>(),vector[0]));
 710   if(_y_array)
 711     std::transform(_y_array->getConstPointer(),_y_array->getConstPointer()+_y_array->getNbOfElems(),
 712                    _y_array->getPointer(),std::bind2nd(std::plus<double>(),vector[1]));
 713   if(_z_array)
 714     std::transform(_z_array->getConstPointer(),_z_array->getConstPointer()+_z_array->getNbOfElems(),
 715                    _z_array->getPointer(),std::bind2nd(std::plus<double>(),vector[2]));
 716 }
 717
 718 /*!
 719  * Applies scaling transformation to all nodes of \a this mesh.
 720  *  \param [in] point - coordinates of a scaling center. This array is to be of
 721  *         size \a this->getSpaceDimension() at least.
 722  *  \param [in] factor - a scale factor.
 723  */
 724 void MEDCouplingCMesh::scale(const double *point, double factor)
 725 {
 726   for(int i=0;i<3;i++)
 727     {
 728       DataArrayDouble *c=getCoordsAt(i);
 729       if(c)
 730         {
 731           double *coords=c->getPointer();
 732           int lgth=c->getNbOfElems();
 733           std::transform(coords,coords+lgth,coords,std::bind2nd(std::minus<double>(),point[i]));
 734           std::transform(coords,coords+lgth,coords,std::bind2nd(std::multiplies<double>(),factor));
 735           std::transform(coords,coords+lgth,coords,std::bind2nd(std::plus<double>(),point[i]));
 736           c->declareAsNew();
 737         }
 738     }
 739   updateTime();
 740 }
 741
 742 MEDCouplingMesh *MEDCouplingCMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
 743 {
 744   //not implemented yet !
 745   return 0;
 746 }
 747
 748 /*!
 749  * Returns a new DataArrayDouble holding coordinates of all nodes of \a this mesh.
 750  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
 751  *          this->getNumberOfNodes() tuples per \a this->getSpaceDimension()
 752  *          components. The caller is to delete this array using decrRef() as it is
 753  *          no more needed.
 754  */
 755 DataArrayDouble *MEDCouplingCMesh::getCoordinatesAndOwner() const
 756 {
 757   DataArrayDouble *ret=DataArrayDouble::New();
 758   int spaceDim=getSpaceDimension();
 759   int nbNodes=getNumberOfNodes();
 760   ret->alloc(nbNodes,spaceDim);
 761   double *pt=ret->getPointer();
 762   int tmp[3];
 763   getSplitNodeValues(tmp);
 764   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 765   const double *tabsPtr[3];
 766   for(int j=0;j<spaceDim;j++)
 767     {
 768       tabsPtr[j]=tabs[j]->getConstPointer();
 769       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0));
 770     }
 771   int tmp2[3];
 772   for(int i=0;i<nbNodes;i++)
 773     {
 774       GetPosFromId(i,spaceDim,tmp,tmp2);
 775       for(int j=0;j<spaceDim;j++)
 776         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 777     }
 778   return ret;
 779 }
 780
 781 /*!
 782  * Returns a new DataArrayDouble holding barycenters of all cells. The barycenter is
 783  * computed by averaging coordinates of cell nodes.
 784  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
 785  *          this->getNumberOfCells() tuples per \a this->getSpaceDimension()
 786  *          components. The caller is to delete this array using decrRef() as it is
 787  *          no more needed.
 788  */
 789 DataArrayDouble *MEDCouplingCMesh::getBarycenterAndOwner() const
 790 {
 791   DataArrayDouble *ret=DataArrayDouble::New();
 792   int spaceDim=getSpaceDimension();
 793   int nbCells=getNumberOfCells();
 794   ret->alloc(nbCells,spaceDim);
 795   double *pt=ret->getPointer();
 796   int tmp[3];
 797   getSplitCellValues(tmp);
 798   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 799   std::vector<double> tabsPtr[3];
 800   for(int j=0;j<spaceDim;j++)
 801     {
 802       int sz=tabs[j]->getNbOfElems()-1;
 803       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0));
 804       const double *srcPtr=tabs[j]->getConstPointer();
 805       tabsPtr[j].insert(tabsPtr[j].end(),srcPtr,srcPtr+sz);
 806       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),srcPtr+1,tabsPtr[j].begin(),std::plus<double>());
 807       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),tabsPtr[j].begin(),std::bind2nd(std::multiplies<double>(),0.5));
 808     }
 809   int tmp2[3];
 810   for(int i=0;i<nbCells;i++)
 811     {
 812       GetPosFromId(i,spaceDim,tmp,tmp2);
 813       for(int j=0;j<spaceDim;j++)
 814         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 815     }
 816   return ret;
 817 }
 818
 819 DataArrayDouble *MEDCouplingCMesh::computeIsoBarycenterOfNodesPerCell() const
 820 {
 821   return MEDCouplingCMesh::getBarycenterAndOwner();
 822 }
 823
 824 void MEDCouplingCMesh::renumberCells(const int *old2NewBg, bool check)
 825 {
 826   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CMesh !");
 827 }
 828
 829 void MEDCouplingCMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 830 {
 831   int it,order;
 832   double time=getTime(it,order);
 833   tinyInfo.clear();
 834   tinyInfoD.clear();
 835   littleStrings.clear();
 836   littleStrings.push_back(getName());
 837   littleStrings.push_back(getDescription());
 838   littleStrings.push_back(getTimeUnit());
 839   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 840   for(int i=0;i<3;i++)
 841     {
 842       int val=-1;
 843       std::string st;
 844       if(thisArr[i])
 845         {
 846           val=thisArr[i]->getNumberOfTuples();
 847           st=thisArr[i]->getInfoOnComponent(0);
 848         }
 849       tinyInfo.push_back(val);
 850       littleStrings.push_back(st);
 851     }
 852   tinyInfo.push_back(it);
 853   tinyInfo.push_back(order);
 854   tinyInfoD.push_back(time);
 855 }
 856
 857 void MEDCouplingCMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 858 {
 859   a1->alloc(0,1);
 860   int sum=0;
 861   for(int i=0;i<3;i++)
 862     if(tinyInfo[i]!=-1)
 863       sum+=tinyInfo[i];
 864   a2->alloc(sum,1);
 865 }
 866
 867 void MEDCouplingCMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
 868 {
 869   a1=DataArrayInt::New();
 870   a1->alloc(0,1);
 871   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 872   int sz=0;
 873   for(int i=0;i<3;i++)
 874     {
 875       if(thisArr[i])
 876         sz+=thisArr[i]->getNumberOfTuples();
 877     }
 878   a2=DataArrayDouble::New();
 879   a2->alloc(sz,1);
 880   double *a2Ptr=a2->getPointer();
 881   for(int i=0;i<3;i++)
 882     if(thisArr[i])
 883       a2Ptr=std::copy(thisArr[i]->getConstPointer(),thisArr[i]->getConstPointer()+thisArr[i]->getNumberOfTuples(),a2Ptr);
 884 }
 885
 886 void MEDCouplingCMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
 887                                        const std::vector<std::string>& littleStrings)
 888 {
 889   setName(littleStrings[0]);
 890   setDescription(littleStrings[1]);
 891   setTimeUnit(littleStrings[2]);
 892   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
 893   const double *data=a2->getConstPointer();
 894   for(int i=0;i<3;i++)
 895     {
 896       if(tinyInfo[i]!=-1)
 897         {
 898           (*(thisArr[i]))=DataArrayDouble::New();
 899           (*(thisArr[i]))->alloc(tinyInfo[i],1);
 900           (*(thisArr[i]))->setInfoOnComponent(0,littleStrings[i+3]);
 901           std::copy(data,data+tinyInfo[i],(*(thisArr[i]))->getPointer());
 902           data+=tinyInfo[i];
 903         }
 904     }
 905   setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
 906 }
 907
 908 void MEDCouplingCMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData, DataArrayByte *byteData) const
 909 {
 910   std::ostringstream extent;
 911   DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 912   for(int i=0;i<3;i++)
 913     {
 914       if(thisArr[i])
 915         { extent << "0 " <<  thisArr[i]->getNumberOfTuples()-1 << " "; }
 916       else
 917         { extent << "0 0 "; }
 918     }
 919   ofs << "  <" << getVTKDataSetType() << " WholeExtent=\"" << extent.str() << "\">\n";
 920   ofs << "    <Piece Extent=\"" << extent.str() << "\">\n";
 921   ofs << "      <PointData>\n" << pointData << std::endl;
 922   ofs << "      </PointData>\n";
 923   ofs << "      <CellData>\n" << cellData << std::endl;
 924   ofs << "      </CellData>\n";
 925   ofs << "      <Coordinates>\n";
 926   for(int i=0;i<3;i++)
 927     {
 928       if(thisArr[i])
 929         thisArr[i]->writeVTK(ofs,8,"Array",byteData);
 930       else
 931         {
 932           MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> coo=DataArrayDouble::New(); coo->alloc(1,1);
 933           coo->setIJ(0,0,0.);
 934           coo->writeVTK(ofs,8,"Array",byteData);
 935         }
 936     }
 937   ofs << "      </Coordinates>\n";
 938   ofs << "    </Piece>\n";
 939   ofs << "  </" << getVTKDataSetType() << ">\n";
 940 }
 941
 942 void MEDCouplingCMesh::reprQuickOverview(std::ostream& stream) const
 943 {
 944   stream << "MEDCouplingCMesh C++ instance at " << this << ". Name : \"" << getName() << "\".";
 945   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 946   std::ostringstream stream2[3];
 947   bool isDef[3];
 948   int nbOfCells=1,nbOfNodes=1;
 949   for(int i=0;i<3;i++)
 950     {
 951       isDef[i]=thisArr[i]!=0;
 952       if(isDef[i])
 953         {
 954           char tmp='X'+i;
 955           stream2[i] << tmp << " positions array ";
 956           if(!thisArr[i]->isAllocated())
 957             stream2[i] << "set but not allocated.";
 958           else
 959             {
 960               int nbCompo=thisArr[i]->getNumberOfComponents();
 961               if(nbCompo==1)
 962                 {
 963                   int nbTuples=thisArr[i]->getNumberOfTuples();
 964                   if(nbTuples<1)
 965                     { stream2[i] << "set and allocated - WARNING number of elements < 1 !"; nbOfCells=-1; nbOfNodes=-1; }
 966                   else
 967                     {
 968                       stream2[i] << "(length=" << nbTuples << ")" << ": ";
 969                       thisArr[i]->reprQuickOverviewData(stream2[i],200);
 970                       if(nbOfCells!=-1)
 971                         { nbOfNodes*=nbTuples; nbOfCells*=nbTuples-1; }
 972                     }
 973                 }
 974               else
 975                 { stream2[i] << "set and allocated - WARNING number of components != 1 !"; nbOfCells=-1; nbOfNodes=-1; }
 976             }
 977         }
 978     }
 979   if(!isDef[0] && !isDef[1] && !isDef[2])
 980     { stream << " No arrays set !"; return; }
 981   if(nbOfCells>=0)
 982     { stream << std::endl << "Number of cells : " << nbOfCells << ". Number of nodes : " << nbOfNodes << "."; }
 983   for(int i=0;i<3;i++)
 984     {
 985       if(isDef[i])
 986         stream << std::endl << stream2[i].str();
 987     }
 988
 989 }
 990
 991 std::string MEDCouplingCMesh::getVTKDataSetType() const
 992 {
 993   return std::string("RectilinearGrid");
 994 }