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