medtool/src/MEDCoupling/MEDCouplingCMesh.cxx

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