src/MEDCoupling/MEDCouplingCMesh.cxx

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