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