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