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