src/MEDCoupling/MEDCouplingCMesh.cxx

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