src/MEDCoupling/MEDCouplingCMesh.cxx

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