src/MEDCoupling/MEDCouplingCMesh.cxx

   1 // Copyright (C) 2007-2012  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
  20 #include "MEDCouplingCMesh.hxx"
  21 #include "MEDCouplingUMesh.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):MEDCouplingMesh(other)
  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 /*!
 110  * This method copyies all tiny strings from other (name and components name).
 111  * @throw if other and this have not same mesh type.
 112  */
 113 void MEDCouplingCMesh::copyTinyStringsFrom(const MEDCouplingMesh *other) throw(INTERP_KERNEL::Exception)
 114 {
 115   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
 116   if(!otherC)
 117     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::copyTinyStringsFrom : meshes have not same type !");
 118   MEDCouplingMesh::copyTinyStringsFrom(other);
 119   if(_x_array && otherC->_x_array)
 120     _x_array->copyStringInfoFrom(*otherC->_x_array);
 121   if(_y_array && otherC->_y_array)
 122     _y_array->copyStringInfoFrom(*otherC->_y_array);
 123   if(_z_array && otherC->_z_array)
 124     _z_array->copyStringInfoFrom(*otherC->_z_array);
 125 }
 126
 127 bool MEDCouplingCMesh::isEqualIfNotWhy(const MEDCouplingMesh *other, double prec, std::string& reason) const throw(INTERP_KERNEL::Exception)
 128 {
 129   if(!other)
 130     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::isEqualIfNotWhy : input other pointer is null !");
 131   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
 132   if(!otherC)
 133     {
 134       reason="mesh given in input is not castable in MEDCouplingCMesh !";
 135       return false;
 136     }
 137   if(!MEDCouplingMesh::isEqualIfNotWhy(other,prec,reason))
 138     return false;
 139   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 140   const DataArrayDouble *otherArr[3]={otherC->_x_array,otherC->_y_array,otherC->_z_array};
 141   std::ostringstream oss; oss.precision(15);
 142   for(int i=0;i<3;i++)
 143     {
 144       if((thisArr[i]!=0 && otherArr[i]==0) || (thisArr[i]==0 && otherArr[i]!=0))
 145         {
 146           oss << "Only one CMesh between the two this and other has its coordinates of rank" << i << " defined !";
 147           reason=oss.str();
 148           return false;
 149         }
 150       if(thisArr[i])
 151         if(!thisArr[i]->isEqualIfNotWhy(*otherArr[i],prec,reason))
 152           {
 153             oss << "Coordinates DataArrayDouble of rank #" << i << " differ :";
 154             reason.insert(0,oss.str());
 155             return false;
 156           }
 157     }
 158   return true;
 159 }
 160
 161 bool MEDCouplingCMesh::isEqualWithoutConsideringStr(const MEDCouplingMesh *other, double prec) const
 162 {
 163   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
 164   if(!otherC)
 165     return false;
 166   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 167   const DataArrayDouble *otherArr[3]={otherC->_x_array,otherC->_y_array,otherC->_z_array};
 168   for(int i=0;i<3;i++)
 169     {
 170       if((thisArr[i]!=0 && otherArr[i]==0) || (thisArr[i]==0 && otherArr[i]!=0))
 171         return false;
 172       if(thisArr[i])
 173         if(!thisArr[i]->isEqualWithoutConsideringStr(*otherArr[i],prec))
 174           return false;
 175     }
 176   return true;
 177 }
 178
 179 void MEDCouplingCMesh::checkDeepEquivalWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
 180                                             DataArrayInt *&cellCor, DataArrayInt *&nodeCor) const throw(INTERP_KERNEL::Exception)
 181 {
 182   if(!isEqualWithoutConsideringStr(other,prec))
 183     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkDeepEquivalWith : Meshes are not the same !");
 184 }
 185
 186 /*!
 187  * Nothing is done here (except to check that the other is a ParaMEDMEM::MEDCouplingCMesh instance too).
 188  * The user intend that the nodes are the same, so by construction of ParaMEDMEM::MEDCouplingCMesh, 'this' and 'other' are the same !
 189  */
 190 void MEDCouplingCMesh::checkDeepEquivalOnSameNodesWith(const MEDCouplingMesh *other, int cellCompPol, double prec,
 191                                                        DataArrayInt *&cellCor) const throw(INTERP_KERNEL::Exception)
 192 {
 193   const MEDCouplingCMesh *otherC=dynamic_cast<const MEDCouplingCMesh *>(other);
 194   if(!otherC)
 195     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkDeepEquivalOnSameNodesWith : other is NOT a cartesian mesh ! Impossible to check equivalence !");
 196 }
 197
 198 void MEDCouplingCMesh::checkCoherency() const throw(INTERP_KERNEL::Exception)
 199 {
 200   const char msg0[]="Invalid ";
 201   const char msg1[]=" array ! Must contain more than 1 element.";
 202   const char msg2[]=" array ! Must be with only one component.";
 203   if(_x_array)
 204     {
 205       if(_x_array->getNbOfElems()<2)
 206         {
 207           std::ostringstream os; os << msg0 << 'X' << msg1;
 208           throw INTERP_KERNEL::Exception(os.str().c_str());
 209         }
 210       if(_x_array->getNumberOfComponents()!=1)
 211         {
 212           std::ostringstream os; os << msg0 << 'X' << msg2;
 213           throw INTERP_KERNEL::Exception(os.str().c_str());
 214         }
 215     }
 216   if(_y_array)
 217     {
 218       if(_y_array->getNbOfElems()<2)
 219         {
 220           std::ostringstream os; os << msg0 << 'Y' << msg1;
 221           throw INTERP_KERNEL::Exception(os.str().c_str());
 222         }
 223       if(_y_array->getNumberOfComponents()!=1)
 224         {
 225           std::ostringstream os; os << msg0 << 'Y' << msg2;
 226           throw INTERP_KERNEL::Exception(os.str().c_str());
 227         }
 228
 229     }
 230   if(_z_array)
 231     {
 232       if(_z_array->getNbOfElems()<2)
 233         {
 234           std::ostringstream os; os << msg0 << 'Z' << msg1;
 235           throw INTERP_KERNEL::Exception(os.str().c_str());
 236         }
 237       if(_z_array->getNumberOfComponents()!=1)
 238         {
 239           std::ostringstream os; os << msg0 << 'Z' << msg2;
 240           throw INTERP_KERNEL::Exception(os.str().c_str());
 241         }
 242     }
 243 }
 244
 245 void MEDCouplingCMesh::checkCoherency1(double eps) const throw(INTERP_KERNEL::Exception)
 246 {
 247   checkCoherency();
 248   if(_x_array)
 249     _x_array->checkMonotonic(true, eps);
 250   if(_y_array)
 251     _y_array->checkMonotonic(true, eps);
 252   if(_z_array)
 253     _z_array->checkMonotonic(true, eps);
 254 }
 255
 256 void MEDCouplingCMesh::checkCoherency2(double eps) const throw(INTERP_KERNEL::Exception)
 257 {
 258   checkCoherency1(eps);
 259 }
 260
 261 int MEDCouplingCMesh::getNumberOfCells() const
 262 {
 263   int ret=1;
 264   if(_x_array)
 265     ret*=_x_array->getNbOfElems()-1;
 266   if(_y_array)
 267     ret*=_y_array->getNbOfElems()-1;
 268   if(_z_array)
 269     ret*=_z_array->getNbOfElems()-1;
 270   return ret;
 271 }
 272
 273 int MEDCouplingCMesh::getNumberOfNodes() const
 274 {
 275   int ret=1;
 276   if(_x_array)
 277     ret*=_x_array->getNbOfElems();
 278   if(_y_array)
 279     ret*=_y_array->getNbOfElems();
 280   if(_z_array)
 281     ret*=_z_array->getNbOfElems();
 282   return ret;
 283 }
 284
 285 void MEDCouplingCMesh::getSplitCellValues(int *res) const
 286 {
 287   int spaceDim=getSpaceDimension();
 288   for(int l=0;l<spaceDim;l++)
 289     {
 290       int val=1;
 291       for(int p=0;p<spaceDim-l-1;p++)
 292         val*=getCoordsAt(p)->getNbOfElems()-1;
 293       res[spaceDim-l-1]=val;
 294     }
 295 }
 296
 297 void MEDCouplingCMesh::getSplitNodeValues(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();
 305       res[spaceDim-l-1]=val;
 306     }
 307 }
 308
 309 int MEDCouplingCMesh::getCellIdFromPos(int i, int j, int k) const
 310 {
 311   int tmp[3]={i,j,k};
 312   int tmp2[3];
 313   int spaceDim=getSpaceDimension();
 314   getSplitCellValues(tmp2);
 315   std::transform(tmp,tmp+spaceDim,tmp2,tmp,std::multiplies<int>());
 316   return std::accumulate(tmp,tmp+spaceDim,0);
 317 }
 318
 319 int MEDCouplingCMesh::getNodeIdFromPos(int i, int j, int k) const
 320 {
 321   int tmp[3]={i,j,k};
 322   int tmp2[3];
 323   int spaceDim=getSpaceDimension();
 324   getSplitNodeValues(tmp2);
 325   std::transform(tmp,tmp+spaceDim,tmp2,tmp,std::multiplies<int>());
 326   return std::accumulate(tmp,tmp+spaceDim,0);
 327 }
 328
 329 void MEDCouplingCMesh::GetPosFromId(int nodeId, int spaceDim, const int *split, int *res)
 330 {
 331   int work=nodeId;
 332   for(int i=spaceDim-1;i>=0;i--)
 333     {
 334       int pos=work/split[i];
 335       work=work%split[i];
 336       res[i]=pos;
 337     }
 338 }
 339
 340 int MEDCouplingCMesh::getSpaceDimension() const
 341 {
 342   int ret=0;
 343   if(_x_array)
 344     ret++;
 345   if(_y_array)
 346     ret++;
 347   if(_z_array)
 348     ret++;
 349   return ret;
 350 }
 351
 352 int MEDCouplingCMesh::getMeshDimension() const
 353 {
 354   return getSpaceDimension();
 355 }
 356
 357 INTERP_KERNEL::NormalizedCellType MEDCouplingCMesh::getTypeOfCell(int cellId) const
 358 {
 359   switch(getMeshDimension())
 360     {
 361     case 3:
 362       return INTERP_KERNEL::NORM_HEXA8;
 363     case 2:
 364       return INTERP_KERNEL::NORM_QUAD4;
 365     case 1:
 366       return INTERP_KERNEL::NORM_SEG2;
 367     default:
 368       throw INTERP_KERNEL::Exception("Unexpected dimension for MEDCouplingCMesh::getTypeOfCell !");
 369     }
 370 }
 371
 372 std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingCMesh::getAllGeoTypes() const
 373 {
 374   INTERP_KERNEL::NormalizedCellType ret;
 375   switch(getMeshDimension())
 376     {
 377     case 3:
 378       ret=INTERP_KERNEL::NORM_HEXA8;
 379       break;
 380     case 2:
 381       ret=INTERP_KERNEL::NORM_QUAD4;
 382       break;
 383     case 1:
 384       ret=INTERP_KERNEL::NORM_SEG2;
 385       break;
 386     default:
 387       throw INTERP_KERNEL::Exception("Unexpected dimension for MEDCouplingCMesh::getAllGeoTypes !");
 388     }
 389   std::set<INTERP_KERNEL::NormalizedCellType> ret2;
 390   ret2.insert(ret);
 391   return ret2;
 392 }
 393
 394 int MEDCouplingCMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
 395 {
 396   int ret=getNumberOfCells();
 397   int dim=getMeshDimension();
 398   switch(type)
 399     {
 400     case INTERP_KERNEL::NORM_HEXA8:
 401       if(dim==3)
 402         return ret;
 403     case INTERP_KERNEL::NORM_QUAD4:
 404       if(dim==2)
 405         return ret;
 406     case INTERP_KERNEL::NORM_SEG2:
 407       if(dim==1)
 408         return ret;
 409     default:
 410       throw INTERP_KERNEL::Exception("Unexpected dimension for MEDCouplingCMesh::getTypeOfCell !");
 411     }
 412   return 0;
 413 }
 414
 415 void MEDCouplingCMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
 416 {
 417   int spaceDim=getSpaceDimension();
 418   int tmpCell[3],tmpNode[3];
 419   getSplitCellValues(tmpCell);
 420   getSplitNodeValues(tmpNode);
 421   int tmp2[3];
 422   GetPosFromId(cellId,spaceDim,tmpCell,tmp2);
 423   switch(spaceDim)
 424     {
 425     case 1:
 426       conn.push_back(tmp2[0]); conn.push_back(tmp2[0]+1);
 427       break;
 428     case 2:
 429       conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]); conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+1);
 430       conn.push_back((tmp2[1]+1)*(tmpCell[1]+1)+tmp2[0]+1); conn.push_back((tmp2[1]+1)*(tmpCell[1]+1)+tmp2[0]);
 431       break;
 432     case 3:
 433       conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+tmp2[2]*tmpNode[2]); conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+1+tmp2[2]*tmpNode[2]);
 434       conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+1+tmp2[2]*tmpNode[2]); conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+tmp2[2]*tmpNode[2]);
 435       conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+(tmp2[2]+1)*tmpNode[2]); conn.push_back(tmp2[1]*tmpCell[1]+tmp2[0]+1+(tmp2[2]+1)*tmpNode[2]);
 436       conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+1+(tmp2[2]+1)*tmpNode[2]); conn.push_back((tmp2[1]+1)*tmpNode[1]+tmp2[0]+(tmp2[2]+1)*tmpNode[2]);
 437       break;
 438     default:
 439       throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getNodeIdsOfCell : big problem spacedim must be in 1,2 or 3 !");
 440     };
 441 }
 442
 443 void MEDCouplingCMesh::getCoordinatesOfNode(int nodeId, std::vector<double>& coo) const throw(INTERP_KERNEL::Exception)
 444 {
 445   int tmp[3];
 446   int spaceDim=getSpaceDimension();
 447   getSplitNodeValues(tmp);
 448   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 449   int tmp2[3];
 450   GetPosFromId(nodeId,spaceDim,tmp,tmp2);
 451   for(int j=0;j<spaceDim;j++)
 452     if(tabs[j])
 453       coo.push_back(tabs[j]->getConstPointer()[tmp2[j]]);
 454 }
 455
 456 std::string MEDCouplingCMesh::simpleRepr() const
 457 {
 458   std::ostringstream ret;
 459   ret << "Cartesian mesh with name : \"" << getName() << "\"\n";
 460   ret << "Description of mesh : \"" << getDescription() << "\"\n";
 461   int tmpp1,tmpp2;
 462   double tt=getTime(tmpp1,tmpp2);
 463   ret << "Time attached to the mesh [unit] : " << tt << " [" << getTimeUnit() << "]\n";
 464   ret << "Iteration : " << tmpp1  << " Order : " << tmpp2 << "\n";
 465   ret << "Mesh and SpaceDimension dimension : " << getSpaceDimension() << "\n\nArrays :\n________\n\n";
 466   if(_x_array)
 467     {
 468       ret << "X Array :\n";
 469       _x_array->reprZipWithoutNameStream(ret);
 470     }
 471   if(_y_array)
 472     {
 473       ret << "Y Array :\n";
 474       _y_array->reprZipWithoutNameStream(ret);
 475     }
 476   if(_z_array)
 477     {
 478       ret << "Z Array :\n";
 479       _z_array->reprZipWithoutNameStream(ret);
 480     }
 481   return ret.str();
 482 }
 483
 484 std::string MEDCouplingCMesh::advancedRepr() const
 485 {
 486   return simpleRepr();
 487 }
 488
 489 const DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i) const throw(INTERP_KERNEL::Exception)
 490 {
 491   switch(i)
 492     {
 493     case 0:
 494       return _x_array;
 495     case 1:
 496       return _y_array;
 497     case 2:
 498       return _z_array;
 499     default:
 500       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 501     }
 502 }
 503
 504 DataArrayDouble *MEDCouplingCMesh::getCoordsAt(int i) throw(INTERP_KERNEL::Exception)
 505 {
 506   switch(i)
 507     {
 508     case 0:
 509       return _x_array;
 510     case 1:
 511       return _y_array;
 512     case 2:
 513       return _z_array;
 514     default:
 515       throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 516     }
 517 }
 518
 519 void MEDCouplingCMesh::setCoordsAt(int i, const DataArrayDouble *arr) throw(INTERP_KERNEL::Exception)
 520 {
 521   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
 522   if(i<0 || i>2)
 523     throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 524   if(arr!=*(thisArr[i]))
 525     {
 526       if(*(thisArr[i]))
 527         (*(thisArr[i]))->decrRef();
 528       (*(thisArr[i]))=const_cast<DataArrayDouble *>(arr);
 529       if(*(thisArr[i]))
 530         (*(thisArr[i]))->incrRef();
 531       declareAsNew();
 532     }
 533 }
 534
 535 void MEDCouplingCMesh::setCoords(const DataArrayDouble *coordsX, const DataArrayDouble *coordsY, const DataArrayDouble *coordsZ)
 536 {
 537   if(_x_array)
 538     _x_array->decrRef();
 539   _x_array=const_cast<DataArrayDouble *>(coordsX);
 540   if(_x_array)
 541     _x_array->incrRef();
 542   if(_y_array)
 543     _y_array->decrRef();
 544   _y_array=const_cast<DataArrayDouble *>(coordsY);
 545   if(_y_array)
 546     _y_array->incrRef();
 547   if(_z_array)
 548     _z_array->decrRef();
 549   _z_array=const_cast<DataArrayDouble *>(coordsZ);
 550   if(_z_array)
 551     _z_array->incrRef();
 552   declareAsNew();
 553 }
 554
 555 /*!
 556  * See MEDCouplingUMesh::getDistributionOfTypes for more information
 557  */
 558 std::vector<int> MEDCouplingCMesh::getDistributionOfTypes() const throw(INTERP_KERNEL::Exception)
 559 {
 560   //only one type of cell
 561   std::vector<int> ret(3);
 562   ret[0]=getTypeOfCell(0);
 563   ret[1]=getNumberOfCells();
 564   ret[2]=0; //ret[3*k+2]==0 because it has no sense here
 565   return ret;
 566 }
 567
 568 /*!
 569  * See MEDCouplingUMesh::checkTypeConsistencyAndContig for more information
 570  */
 571 DataArrayInt *MEDCouplingCMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 572 {
 573   if(code.empty())
 574     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : code is empty, should not !");
 575   std::size_t sz=code.size();
 576   if(sz!=3)
 577     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : code should be of size 3 exactly !");
 578
 579   int nbCells=getNumberOfCellsWithType((INTERP_KERNEL::NormalizedCellType)code[0]);
 580   if(code[2]==-1)
 581     {
 582       if(code[1]==nbCells)
 583         return 0;
 584       else
 585         throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : number of cells mismatch !");
 586     }
 587   else
 588     {
 589       if(code[2]<-1)
 590         throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : code[2]<-1 mismatch !");
 591       if(code[2]>=(int)idsPerType.size())
 592         throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : code[2]>size idsPerType !");
 593       return idsPerType[code[2]]->deepCpy();
 594     }
 595 }
 596
 597 /*!
 598  * See MEDCouplingUMesh::splitProfilePerType for more information
 599  */
 600 void MEDCouplingCMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 601 {
 602   int nbCells=getNumberOfCells();
 603   code.resize(3);
 604   code[0]=(int)getTypeOfCell(0);
 605   code[1]=nbCells;
 606   code[2]=0;
 607   idsInPflPerType.push_back(profile->deepCpy());
 608   idsPerType.push_back(profile->deepCpy());
 609 }
 610
 611 MEDCouplingUMesh *MEDCouplingCMesh::buildUnstructured() const throw(INTERP_KERNEL::Exception)
 612 {
 613   int spaceDim=getSpaceDimension();
 614   MEDCouplingUMesh *ret=MEDCouplingUMesh::New(getName(),spaceDim);
 615   DataArrayDouble *coords=getCoordinatesAndOwner();
 616   ret->setCoords(coords);
 617   coords->decrRef();
 618   switch(spaceDim)
 619     {
 620     case 1:
 621       fill1DUnstructuredMesh(ret);
 622       break;
 623     case 2:
 624       fill2DUnstructuredMesh(ret);
 625       break;
 626     case 3:
 627       fill3DUnstructuredMesh(ret);
 628       break;
 629     default:
 630       throw INTERP_KERNEL::Exception("MEDCouplingCMesh::buildUnstructured : big problem spacedim must be in 1,2 or 3 !");
 631     };
 632   return ret;
 633 }
 634
 635 MEDCouplingMesh *MEDCouplingCMesh::buildPart(const int *start, const int *end) const
 636 {
 637   MEDCouplingUMesh *um=buildUnstructured();
 638   MEDCouplingMesh *ret=um->buildPart(start,end);
 639   um->decrRef();
 640   return ret;
 641 }
 642
 643 MEDCouplingMesh *MEDCouplingCMesh::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
 644 {
 645   MEDCouplingUMesh *um=buildUnstructured();
 646   MEDCouplingMesh *ret=um->buildPartAndReduceNodes(start,end,arr);
 647   um->decrRef();
 648   return ret;
 649 }
 650
 651 DataArrayInt *MEDCouplingCMesh::simplexize(int policy) throw(INTERP_KERNEL::Exception)
 652 {
 653   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::simplexize : not available for Cartesian mesh !");
 654 }
 655
 656 void MEDCouplingCMesh::getBoundingBox(double *bbox) const
 657 {
 658   int dim=getSpaceDimension();
 659   int j=0;
 660   for (int idim=0; idim<dim; idim++)
 661     {
 662       const DataArrayDouble *c=getCoordsAt(idim);
 663       if(c)
 664         {
 665           const double *coords=c->getConstPointer();
 666           int nb=c->getNbOfElems();
 667           bbox[2*j]=coords[0];
 668           bbox[2*j+1]=coords[nb-1];
 669           j++;
 670         }
 671     }
 672 }
 673
 674 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const
 675 {
 676   std::string name="MeasureOfMesh_";
 677   name+=getName();
 678   int nbelem=getNumberOfCells();
 679   MEDCouplingFieldDouble *field=MEDCouplingFieldDouble::New(ON_CELLS);
 680   field->setName(name.c_str());
 681   DataArrayDouble* array=DataArrayDouble::New();
 682   array->alloc(nbelem,1);
 683   double *area_vol=array->getPointer();
 684   field->setArray(array) ;
 685   array->decrRef();
 686   field->setMesh(const_cast<MEDCouplingCMesh *>(this));
 687   int tmp[3];
 688   getSplitCellValues(tmp);
 689   int dim=getSpaceDimension();
 690   const double **thisArr=new const double *[dim];
 691   const DataArrayDouble *thisArr2[3]={_x_array,_y_array,_z_array};
 692   for(int i=0;i<dim;i++)
 693     thisArr[i]=thisArr2[i]->getConstPointer();
 694   for(int icell=0;icell<nbelem;icell++)
 695     {
 696       int tmp2[3];
 697       GetPosFromId(icell,dim,tmp,tmp2);
 698       area_vol[icell]=1.;
 699       for(int i=0;i<dim;i++)
 700         area_vol[icell]*=thisArr[i][tmp2[i]+1]-thisArr[i][tmp2[i]];
 701     }
 702   delete [] thisArr;
 703   return field;
 704 }
 705
 706 /*!
 707  * not implemented yet !
 708  */
 709 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureFieldOnNode(bool isAbs) const
 710 {
 711   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getMeasureFieldOnNode : not implemented yet !");
 712   //return 0;
 713 }
 714
 715 MEDCouplingFieldDouble *MEDCouplingCMesh::buildOrthogonalField() const
 716 {
 717   if(getMeshDimension()!=2)
 718     throw INTERP_KERNEL::Exception("Expected a cmesh with meshDim == 2 !");
 719   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 720   DataArrayDouble *array=DataArrayDouble::New();
 721   int nbOfCells=getNumberOfCells();
 722   array->alloc(nbOfCells,3);
 723   double *vals=array->getPointer();
 724   for(int i=0;i<nbOfCells;i++)
 725     { vals[3*i]=0.; vals[3*i+1]=0.; vals[3*i+2]=1.; }
 726   ret->setArray(array);
 727   array->decrRef();
 728   ret->setMesh(this);
 729   return ret;
 730 }
 731
 732 int MEDCouplingCMesh::getCellContainingPoint(const double *pos, double eps) const
 733 {
 734   int dim=getSpaceDimension();
 735   int ret=0;
 736   int coeff=1;
 737   for(int i=0;i<dim;i++)
 738     {
 739       const double *d=getCoordsAt(i)->getConstPointer();
 740       int nbOfNodes=getCoordsAt(i)->getNbOfElems();
 741       double ref=pos[i];
 742       const double *w=std::find_if(d,d+nbOfNodes,std::bind2nd(std::greater_equal<double>(),ref));
 743       int w2=(int)std::distance(d,w);
 744       if(w2<nbOfNodes)
 745         {
 746           if(w2==0)
 747             {
 748               if(ref>d[0]-eps)
 749                 w2=1;
 750               else
 751                 return -1;
 752             }
 753           ret+=coeff*(w2-1);
 754           coeff*=nbOfNodes-1;
 755         }
 756       else
 757         return -1;
 758     }
 759   return ret;
 760 }
 761
 762 void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
 763 {
 764   throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !");
 765 }
 766
 767 void MEDCouplingCMesh::translate(const double *vector)
 768 {
 769   if(_x_array)
 770     std::transform(_x_array->getConstPointer(),_x_array->getConstPointer()+_x_array->getNbOfElems(),
 771                    _x_array->getPointer(),std::bind2nd(std::plus<double>(),vector[0]));
 772   if(_y_array)
 773     std::transform(_y_array->getConstPointer(),_y_array->getConstPointer()+_y_array->getNbOfElems(),
 774                    _y_array->getPointer(),std::bind2nd(std::plus<double>(),vector[1]));
 775   if(_z_array)
 776     std::transform(_z_array->getConstPointer(),_z_array->getConstPointer()+_z_array->getNbOfElems(),
 777                    _z_array->getPointer(),std::bind2nd(std::plus<double>(),vector[2]));
 778 }
 779
 780 void MEDCouplingCMesh::scale(const double *point, double factor)
 781 {
 782   for(int i=0;i<3;i++)
 783     {
 784       DataArrayDouble *c=getCoordsAt(i);
 785       if(c)
 786         {
 787           double *coords=c->getPointer();
 788           int lgth=c->getNbOfElems();
 789           std::transform(coords,coords+lgth,coords,std::bind2nd(std::minus<double>(),point[i]));
 790           std::transform(coords,coords+lgth,coords,std::bind2nd(std::multiplies<double>(),factor));
 791           std::transform(coords,coords+lgth,coords,std::bind2nd(std::plus<double>(),point[i]));
 792           c->declareAsNew();
 793         }
 794     }
 795   updateTime();
 796 }
 797
 798 MEDCouplingMesh *MEDCouplingCMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
 799 {
 800   //not implemented yet !
 801   return 0;
 802 }
 803
 804 DataArrayDouble *MEDCouplingCMesh::getCoordinatesAndOwner() const
 805 {
 806   DataArrayDouble *ret=DataArrayDouble::New();
 807   int spaceDim=getSpaceDimension();
 808   int nbNodes=getNumberOfNodes();
 809   ret->alloc(nbNodes,spaceDim);
 810   double *pt=ret->getPointer();
 811   int tmp[3];
 812   getSplitNodeValues(tmp);
 813   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 814   const double *tabsPtr[3];
 815   for(int j=0;j<spaceDim;j++)
 816     {
 817       tabsPtr[j]=tabs[j]->getConstPointer();
 818       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0).c_str());
 819     }
 820   int tmp2[3];
 821   for(int i=0;i<nbNodes;i++)
 822     {
 823       GetPosFromId(i,spaceDim,tmp,tmp2);
 824       for(int j=0;j<spaceDim;j++)
 825         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 826     }
 827   return ret;
 828 }
 829
 830 DataArrayDouble *MEDCouplingCMesh::getBarycenterAndOwner() const
 831 {
 832   DataArrayDouble *ret=DataArrayDouble::New();
 833   int spaceDim=getSpaceDimension();
 834   int nbCells=getNumberOfCells();
 835   ret->alloc(nbCells,spaceDim);
 836   double *pt=ret->getPointer();
 837   int tmp[3];
 838   getSplitCellValues(tmp);
 839   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 840   std::vector<double> tabsPtr[3];
 841   for(int j=0;j<spaceDim;j++)
 842     {
 843       int sz=tabs[j]->getNbOfElems()-1;
 844       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0).c_str());
 845       const double *srcPtr=tabs[j]->getConstPointer();
 846       tabsPtr[j].insert(tabsPtr[j].end(),srcPtr,srcPtr+sz);
 847       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),srcPtr+1,tabsPtr[j].begin(),std::plus<double>());
 848       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),tabsPtr[j].begin(),std::bind2nd(std::multiplies<double>(),0.5));
 849     }
 850   int tmp2[3];
 851   for(int i=0;i<nbCells;i++)
 852     {
 853       GetPosFromId(i,spaceDim,tmp,tmp2);
 854       for(int j=0;j<spaceDim;j++)
 855         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 856     }
 857   return ret;
 858 }
 859
 860 void MEDCouplingCMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 861 {
 862   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CMesh !");
 863 }
 864
 865 void MEDCouplingCMesh::fill1DUnstructuredMesh(MEDCouplingUMesh *m) const
 866 {
 867   const DataArrayDouble *c=getCoordsAt(0);
 868   int nbOfCells=c->getNbOfElems()-1;
 869   DataArrayInt *connI=DataArrayInt::New();
 870   connI->alloc(nbOfCells+1,1);
 871   int *ci=connI->getPointer();
 872   DataArrayInt *conn=DataArrayInt::New();
 873   conn->alloc(3*nbOfCells,1);
 874   ci[0]=0;
 875   int *cp=conn->getPointer();
 876   for(int i=0;i<nbOfCells;i++)
 877     {
 878       cp[3*i]=(int)INTERP_KERNEL::NORM_SEG2;
 879       cp[3*i+1]=i;
 880       cp[3*i+2]=i+1;
 881       ci[i+1]=3*(i+1);
 882     }
 883   m->setConnectivity(conn,connI,true);
 884   conn->decrRef();
 885   connI->decrRef();
 886 }
 887
 888 void MEDCouplingCMesh::fill2DUnstructuredMesh(MEDCouplingUMesh *m) const
 889 {
 890   const DataArrayDouble *c1=getCoordsAt(0);
 891   const DataArrayDouble *c2=getCoordsAt(1);
 892   int n1=c1->getNbOfElems()-1;
 893   int n2=c2->getNbOfElems()-1;
 894   DataArrayInt *connI=DataArrayInt::New();
 895   connI->alloc(n1*n2+1,1);
 896   int *ci=connI->getPointer();
 897   DataArrayInt *conn=DataArrayInt::New();
 898   conn->alloc(5*n1*n2,1);
 899   ci[0]=0;
 900   int *cp=conn->getPointer();
 901   int pos=0;
 902   for(int j=0;j<n2;j++)
 903     for(int i=0;i<n1;i++,pos++)
 904       {
 905         cp[5*pos]=(int)INTERP_KERNEL::NORM_QUAD4;
 906         cp[5*pos+1]=i+1+j*(n1+1);
 907         cp[5*pos+2]=i+j*(n1+1);
 908         cp[5*pos+3]=i+(j+1)*(n1+1);
 909         cp[5*pos+4]=i+1+(j+1)*(n1+1);
 910         ci[pos+1]=5*(pos+1);
 911     }
 912   m->setConnectivity(conn,connI,true);
 913   conn->decrRef();
 914   connI->decrRef();
 915 }
 916
 917 void MEDCouplingCMesh::fill3DUnstructuredMesh(MEDCouplingUMesh *m) const
 918 {
 919   const DataArrayDouble *c1=getCoordsAt(0);
 920   const DataArrayDouble *c2=getCoordsAt(1);
 921   const DataArrayDouble *c3=getCoordsAt(2);
 922   int n1=c1->getNbOfElems()-1;
 923   int n2=c2->getNbOfElems()-1;
 924   int n3=c3->getNbOfElems()-1;
 925   DataArrayInt *connI=DataArrayInt::New();
 926   connI->alloc(n1*n2*n3+1,1);
 927   int *ci=connI->getPointer();
 928   DataArrayInt *conn=DataArrayInt::New();
 929   conn->alloc(9*n1*n2*n3,1);
 930   ci[0]=0;
 931   int *cp=conn->getPointer();
 932   int pos=0;
 933   for(int k=0;k<n3;k++)
 934     for(int j=0;j<n2;j++)
 935       for(int i=0;i<n1;i++,pos++)
 936         {
 937           cp[9*pos]=(int)INTERP_KERNEL::NORM_HEXA8;
 938           int tmp=(n1+1)*(n2+1);
 939           cp[9*pos+1]=i+1+j*(n1+1)+k*tmp;
 940           cp[9*pos+2]=i+j*(n1+1)+k*tmp;
 941           cp[9*pos+3]=i+(j+1)*(n1+1)+k*tmp;
 942           cp[9*pos+4]=i+1+(j+1)*(n1+1)+k*tmp;
 943           cp[9*pos+5]=i+1+j*(n1+1)+(k+1)*tmp;
 944           cp[9*pos+6]=i+j*(n1+1)+(k+1)*tmp;
 945           cp[9*pos+7]=i+(j+1)*(n1+1)+(k+1)*tmp;
 946           cp[9*pos+8]=i+1+(j+1)*(n1+1)+(k+1)*tmp;
 947           ci[pos+1]=9*(pos+1);
 948         }
 949   m->setConnectivity(conn,connI,true);
 950   conn->decrRef();
 951   connI->decrRef();
 952 }
 953
 954 void MEDCouplingCMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 955 {
 956   int it,order;
 957   double time=getTime(it,order);
 958   tinyInfo.clear();
 959   tinyInfoD.clear();
 960   littleStrings.clear();
 961   littleStrings.push_back(getName());
 962   littleStrings.push_back(getDescription());
 963   littleStrings.push_back(getTimeUnit());
 964   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 965   for(int i=0;i<3;i++)
 966     {
 967       int val=-1;
 968       std::string st;
 969       if(thisArr[i])
 970         {
 971           val=thisArr[i]->getNumberOfTuples();
 972           st=thisArr[i]->getInfoOnComponent(0);
 973         }
 974       tinyInfo.push_back(val);
 975       littleStrings.push_back(st);
 976     }
 977   tinyInfo.push_back(it);
 978   tinyInfo.push_back(order);
 979   tinyInfoD.push_back(time);
 980 }
 981
 982 void MEDCouplingCMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 983 {
 984   a1->alloc(0,1);
 985   int sum=0;
 986   for(int i=0;i<3;i++)
 987     if(tinyInfo[i]!=-1)
 988       sum+=tinyInfo[i];
 989   a2->alloc(sum,1);
 990 }
 991
 992 void MEDCouplingCMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
 993 {
 994   a1=DataArrayInt::New();
 995   a1->alloc(0,1);
 996   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 997   int sz=0;
 998   for(int i=0;i<3;i++)
 999     {
1000       if(thisArr[i])
1001         sz+=thisArr[i]->getNumberOfTuples();
1002     }
1003   a2=DataArrayDouble::New();
1004   a2->alloc(sz,1);
1005   double *a2Ptr=a2->getPointer();
1006   for(int i=0;i<3;i++)
1007     if(thisArr[i])
1008       a2Ptr=std::copy(thisArr[i]->getConstPointer(),thisArr[i]->getConstPointer()+thisArr[i]->getNumberOfTuples(),a2Ptr);
1009 }
1010
1011 void MEDCouplingCMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
1012                                        const std::vector<std::string>& littleStrings)
1013 {
1014   setName(littleStrings[0].c_str());
1015   setDescription(littleStrings[1].c_str());
1016   setTimeUnit(littleStrings[2].c_str());
1017   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
1018   const double *data=a2->getConstPointer();
1019   for(int i=0;i<3;i++)
1020     {
1021       if(tinyInfo[i]!=-1)
1022         {
1023           (*(thisArr[i]))=DataArrayDouble::New();
1024           (*(thisArr[i]))->alloc(tinyInfo[i],1);
1025           (*(thisArr[i]))->setInfoOnComponent(0,littleStrings[i+3].c_str());
1026           std::copy(data,data+tinyInfo[i],(*(thisArr[i]))->getPointer());
1027           data+=tinyInfo[i];
1028         }
1029     }
1030   setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
1031 }
1032
1033 void MEDCouplingCMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData) const throw(INTERP_KERNEL::Exception)
1034 {
1035   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::writeVTKLL : not implemented yet !");
1036 }
1037
1038 std::string MEDCouplingCMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exception)
1039 {
1040   return std::string("RectilinearGrid");
1041 }