src/MEDCoupling/MEDCouplingCMesh.cxx

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