src/MEDCoupling/MEDCouplingCMesh.cxx

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