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   if(arr)
 523     arr->checkNbOfComps(1,"MEDCouplingCMesh::setCoordsAt");
 524   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
 525   if(i<0 || i>2)
 526     throw INTERP_KERNEL::Exception("Invalid rank specified must be 0 or 1 or 2.");
 527   if(arr!=*(thisArr[i]))
 528     {
 529       if(*(thisArr[i]))
 530         (*(thisArr[i]))->decrRef();
 531       (*(thisArr[i]))=const_cast<DataArrayDouble *>(arr);
 532       if(*(thisArr[i]))
 533         (*(thisArr[i]))->incrRef();
 534       declareAsNew();
 535     }
 536 }
 537
 538 void MEDCouplingCMesh::setCoords(const DataArrayDouble *coordsX, const DataArrayDouble *coordsY, const DataArrayDouble *coordsZ)
 539 {
 540   if(coordsX)
 541     coordsX->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsX");
 542   if(coordsY)
 543     coordsY->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsY");
 544   if(coordsZ)
 545     coordsZ->checkNbOfComps(1,"MEDCouplingCMesh::setCoords : coordsZ");
 546   if(_x_array)
 547     _x_array->decrRef();
 548   _x_array=const_cast<DataArrayDouble *>(coordsX);
 549   if(_x_array)
 550     _x_array->incrRef();
 551   if(_y_array)
 552     _y_array->decrRef();
 553   _y_array=const_cast<DataArrayDouble *>(coordsY);
 554   if(_y_array)
 555     _y_array->incrRef();
 556   if(_z_array)
 557     _z_array->decrRef();
 558   _z_array=const_cast<DataArrayDouble *>(coordsZ);
 559   if(_z_array)
 560     _z_array->incrRef();
 561   declareAsNew();
 562 }
 563
 564 /*!
 565  * See MEDCouplingUMesh::getDistributionOfTypes for more information
 566  */
 567 std::vector<int> MEDCouplingCMesh::getDistributionOfTypes() const throw(INTERP_KERNEL::Exception)
 568 {
 569   //only one type of cell
 570   std::vector<int> ret(3);
 571   ret[0]=getTypeOfCell(0);
 572   ret[1]=getNumberOfCells();
 573   ret[2]=0; //ret[3*k+2]==0 because it has no sense here
 574   return ret;
 575 }
 576
 577 /*!
 578  * See MEDCouplingUMesh::checkTypeConsistencyAndContig for more information
 579  */
 580 DataArrayInt *MEDCouplingCMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 581 {
 582   if(code.empty())
 583     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : code is empty, should not !");
 584   std::size_t sz=code.size();
 585   if(sz!=3)
 586     throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : code should be of size 3 exactly !");
 587
 588   int nbCells=getNumberOfCellsWithType((INTERP_KERNEL::NormalizedCellType)code[0]);
 589   if(code[2]==-1)
 590     {
 591       if(code[1]==nbCells)
 592         return 0;
 593       else
 594         throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : number of cells mismatch !");
 595     }
 596   else
 597     {
 598       if(code[2]<-1)
 599         throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : code[2]<-1 mismatch !");
 600       if(code[2]>=(int)idsPerType.size())
 601         throw INTERP_KERNEL::Exception("MEDCouplingCMesh::checkTypeConsistencyAndContig : code[2]>size idsPerType !");
 602       return idsPerType[code[2]]->deepCpy();
 603     }
 604 }
 605
 606 /*!
 607  * See MEDCouplingUMesh::splitProfilePerType for more information
 608  */
 609 void MEDCouplingCMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const throw(INTERP_KERNEL::Exception)
 610 {
 611   int nbCells=getNumberOfCells();
 612   code.resize(3);
 613   code[0]=(int)getTypeOfCell(0);
 614   code[1]=nbCells;
 615   code[2]=0;
 616   idsInPflPerType.push_back(profile->deepCpy());
 617   idsPerType.push_back(profile->deepCpy());
 618 }
 619
 620 MEDCouplingUMesh *MEDCouplingCMesh::buildUnstructured() const throw(INTERP_KERNEL::Exception)
 621 {
 622   int spaceDim=getSpaceDimension();
 623   MEDCouplingUMesh *ret=MEDCouplingUMesh::New(getName(),spaceDim);
 624   DataArrayDouble *coords=getCoordinatesAndOwner();
 625   ret->setCoords(coords);
 626   coords->decrRef();
 627   switch(spaceDim)
 628     {
 629     case 1:
 630       fill1DUnstructuredMesh(ret);
 631       break;
 632     case 2:
 633       fill2DUnstructuredMesh(ret);
 634       break;
 635     case 3:
 636       fill3DUnstructuredMesh(ret);
 637       break;
 638     default:
 639       throw INTERP_KERNEL::Exception("MEDCouplingCMesh::buildUnstructured : big problem spacedim must be in 1,2 or 3 !");
 640     };
 641   return ret;
 642 }
 643
 644 MEDCouplingMesh *MEDCouplingCMesh::buildPart(const int *start, const int *end) const
 645 {
 646   MEDCouplingUMesh *um=buildUnstructured();
 647   MEDCouplingMesh *ret=um->buildPart(start,end);
 648   um->decrRef();
 649   return ret;
 650 }
 651
 652 MEDCouplingMesh *MEDCouplingCMesh::buildPartAndReduceNodes(const int *start, const int *end, DataArrayInt*& arr) const
 653 {
 654   MEDCouplingUMesh *um=buildUnstructured();
 655   MEDCouplingMesh *ret=um->buildPartAndReduceNodes(start,end,arr);
 656   um->decrRef();
 657   return ret;
 658 }
 659
 660 DataArrayInt *MEDCouplingCMesh::simplexize(int policy) throw(INTERP_KERNEL::Exception)
 661 {
 662   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::simplexize : not available for Cartesian mesh !");
 663 }
 664
 665 void MEDCouplingCMesh::getBoundingBox(double *bbox) const
 666 {
 667   int dim=getSpaceDimension();
 668   int j=0;
 669   for (int idim=0; idim<dim; idim++)
 670     {
 671       const DataArrayDouble *c=getCoordsAt(idim);
 672       if(c)
 673         {
 674           const double *coords=c->getConstPointer();
 675           int nb=c->getNbOfElems();
 676           bbox[2*j]=coords[0];
 677           bbox[2*j+1]=coords[nb-1];
 678           j++;
 679         }
 680     }
 681 }
 682
 683 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureField(bool isAbs) const
 684 {
 685   std::string name="MeasureOfMesh_";
 686   name+=getName();
 687   int nbelem=getNumberOfCells();
 688   MEDCouplingFieldDouble *field=MEDCouplingFieldDouble::New(ON_CELLS);
 689   field->setName(name.c_str());
 690   DataArrayDouble* array=DataArrayDouble::New();
 691   array->alloc(nbelem,1);
 692   double *area_vol=array->getPointer();
 693   field->setArray(array) ;
 694   array->decrRef();
 695   field->setMesh(const_cast<MEDCouplingCMesh *>(this));
 696   int tmp[3];
 697   getSplitCellValues(tmp);
 698   int dim=getSpaceDimension();
 699   const double **thisArr=new const double *[dim];
 700   const DataArrayDouble *thisArr2[3]={_x_array,_y_array,_z_array};
 701   for(int i=0;i<dim;i++)
 702     thisArr[i]=thisArr2[i]->getConstPointer();
 703   for(int icell=0;icell<nbelem;icell++)
 704     {
 705       int tmp2[3];
 706       GetPosFromId(icell,dim,tmp,tmp2);
 707       area_vol[icell]=1.;
 708       for(int i=0;i<dim;i++)
 709         area_vol[icell]*=thisArr[i][tmp2[i]+1]-thisArr[i][tmp2[i]];
 710     }
 711   delete [] thisArr;
 712   return field;
 713 }
 714
 715 /*!
 716  * not implemented yet !
 717  */
 718 MEDCouplingFieldDouble *MEDCouplingCMesh::getMeasureFieldOnNode(bool isAbs) const
 719 {
 720   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::getMeasureFieldOnNode : not implemented yet !");
 721   //return 0;
 722 }
 723
 724 MEDCouplingFieldDouble *MEDCouplingCMesh::buildOrthogonalField() const
 725 {
 726   if(getMeshDimension()!=2)
 727     throw INTERP_KERNEL::Exception("Expected a cmesh with meshDim == 2 !");
 728   MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(ON_CELLS,NO_TIME);
 729   DataArrayDouble *array=DataArrayDouble::New();
 730   int nbOfCells=getNumberOfCells();
 731   array->alloc(nbOfCells,3);
 732   double *vals=array->getPointer();
 733   for(int i=0;i<nbOfCells;i++)
 734     { vals[3*i]=0.; vals[3*i+1]=0.; vals[3*i+2]=1.; }
 735   ret->setArray(array);
 736   array->decrRef();
 737   ret->setMesh(this);
 738   return ret;
 739 }
 740
 741 int MEDCouplingCMesh::getCellContainingPoint(const double *pos, double eps) const
 742 {
 743   int dim=getSpaceDimension();
 744   int ret=0;
 745   int coeff=1;
 746   for(int i=0;i<dim;i++)
 747     {
 748       const double *d=getCoordsAt(i)->getConstPointer();
 749       int nbOfNodes=getCoordsAt(i)->getNbOfElems();
 750       double ref=pos[i];
 751       const double *w=std::find_if(d,d+nbOfNodes,std::bind2nd(std::greater_equal<double>(),ref));
 752       int w2=(int)std::distance(d,w);
 753       if(w2<nbOfNodes)
 754         {
 755           if(w2==0)
 756             {
 757               if(ref>d[0]-eps)
 758                 w2=1;
 759               else
 760                 return -1;
 761             }
 762           ret+=coeff*(w2-1);
 763           coeff*=nbOfNodes-1;
 764         }
 765       else
 766         return -1;
 767     }
 768   return ret;
 769 }
 770
 771 void MEDCouplingCMesh::rotate(const double *center, const double *vector, double angle)
 772 {
 773   throw INTERP_KERNEL::Exception("No rotation available on CMesh : Traduce it to StructuredMesh to apply it !");
 774 }
 775
 776 void MEDCouplingCMesh::translate(const double *vector)
 777 {
 778   if(_x_array)
 779     std::transform(_x_array->getConstPointer(),_x_array->getConstPointer()+_x_array->getNbOfElems(),
 780                    _x_array->getPointer(),std::bind2nd(std::plus<double>(),vector[0]));
 781   if(_y_array)
 782     std::transform(_y_array->getConstPointer(),_y_array->getConstPointer()+_y_array->getNbOfElems(),
 783                    _y_array->getPointer(),std::bind2nd(std::plus<double>(),vector[1]));
 784   if(_z_array)
 785     std::transform(_z_array->getConstPointer(),_z_array->getConstPointer()+_z_array->getNbOfElems(),
 786                    _z_array->getPointer(),std::bind2nd(std::plus<double>(),vector[2]));
 787 }
 788
 789 void MEDCouplingCMesh::scale(const double *point, double factor)
 790 {
 791   for(int i=0;i<3;i++)
 792     {
 793       DataArrayDouble *c=getCoordsAt(i);
 794       if(c)
 795         {
 796           double *coords=c->getPointer();
 797           int lgth=c->getNbOfElems();
 798           std::transform(coords,coords+lgth,coords,std::bind2nd(std::minus<double>(),point[i]));
 799           std::transform(coords,coords+lgth,coords,std::bind2nd(std::multiplies<double>(),factor));
 800           std::transform(coords,coords+lgth,coords,std::bind2nd(std::plus<double>(),point[i]));
 801           c->declareAsNew();
 802         }
 803     }
 804   updateTime();
 805 }
 806
 807 MEDCouplingMesh *MEDCouplingCMesh::mergeMyselfWith(const MEDCouplingMesh *other) const
 808 {
 809   //not implemented yet !
 810   return 0;
 811 }
 812
 813 DataArrayDouble *MEDCouplingCMesh::getCoordinatesAndOwner() const
 814 {
 815   DataArrayDouble *ret=DataArrayDouble::New();
 816   int spaceDim=getSpaceDimension();
 817   int nbNodes=getNumberOfNodes();
 818   ret->alloc(nbNodes,spaceDim);
 819   double *pt=ret->getPointer();
 820   int tmp[3];
 821   getSplitNodeValues(tmp);
 822   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 823   const double *tabsPtr[3];
 824   for(int j=0;j<spaceDim;j++)
 825     {
 826       tabsPtr[j]=tabs[j]->getConstPointer();
 827       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0).c_str());
 828     }
 829   int tmp2[3];
 830   for(int i=0;i<nbNodes;i++)
 831     {
 832       GetPosFromId(i,spaceDim,tmp,tmp2);
 833       for(int j=0;j<spaceDim;j++)
 834         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 835     }
 836   return ret;
 837 }
 838
 839 DataArrayDouble *MEDCouplingCMesh::getBarycenterAndOwner() const
 840 {
 841   DataArrayDouble *ret=DataArrayDouble::New();
 842   int spaceDim=getSpaceDimension();
 843   int nbCells=getNumberOfCells();
 844   ret->alloc(nbCells,spaceDim);
 845   double *pt=ret->getPointer();
 846   int tmp[3];
 847   getSplitCellValues(tmp);
 848   const DataArrayDouble *tabs[3]={getCoordsAt(0),getCoordsAt(1),getCoordsAt(2)};
 849   std::vector<double> tabsPtr[3];
 850   for(int j=0;j<spaceDim;j++)
 851     {
 852       int sz=tabs[j]->getNbOfElems()-1;
 853       ret->setInfoOnComponent(j,tabs[j]->getInfoOnComponent(0).c_str());
 854       const double *srcPtr=tabs[j]->getConstPointer();
 855       tabsPtr[j].insert(tabsPtr[j].end(),srcPtr,srcPtr+sz);
 856       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),srcPtr+1,tabsPtr[j].begin(),std::plus<double>());
 857       std::transform(tabsPtr[j].begin(),tabsPtr[j].end(),tabsPtr[j].begin(),std::bind2nd(std::multiplies<double>(),0.5));
 858     }
 859   int tmp2[3];
 860   for(int i=0;i<nbCells;i++)
 861     {
 862       GetPosFromId(i,spaceDim,tmp,tmp2);
 863       for(int j=0;j<spaceDim;j++)
 864         pt[i*spaceDim+j]=tabsPtr[j][tmp2[j]];
 865     }
 866   return ret;
 867 }
 868
 869 void MEDCouplingCMesh::renumberCells(const int *old2NewBg, bool check) throw(INTERP_KERNEL::Exception)
 870 {
 871   throw INTERP_KERNEL::Exception("Functionnality of renumbering cell not available for CMesh !");
 872 }
 873
 874 void MEDCouplingCMesh::fill1DUnstructuredMesh(MEDCouplingUMesh *m) const
 875 {
 876   const DataArrayDouble *c=getCoordsAt(0);
 877   int nbOfCells=c->getNbOfElems()-1;
 878   DataArrayInt *connI=DataArrayInt::New();
 879   connI->alloc(nbOfCells+1,1);
 880   int *ci=connI->getPointer();
 881   DataArrayInt *conn=DataArrayInt::New();
 882   conn->alloc(3*nbOfCells,1);
 883   ci[0]=0;
 884   int *cp=conn->getPointer();
 885   for(int i=0;i<nbOfCells;i++)
 886     {
 887       cp[3*i]=(int)INTERP_KERNEL::NORM_SEG2;
 888       cp[3*i+1]=i;
 889       cp[3*i+2]=i+1;
 890       ci[i+1]=3*(i+1);
 891     }
 892   m->setConnectivity(conn,connI,true);
 893   conn->decrRef();
 894   connI->decrRef();
 895 }
 896
 897 void MEDCouplingCMesh::fill2DUnstructuredMesh(MEDCouplingUMesh *m) const
 898 {
 899   const DataArrayDouble *c1=getCoordsAt(0);
 900   const DataArrayDouble *c2=getCoordsAt(1);
 901   int n1=c1->getNbOfElems()-1;
 902   int n2=c2->getNbOfElems()-1;
 903   DataArrayInt *connI=DataArrayInt::New();
 904   connI->alloc(n1*n2+1,1);
 905   int *ci=connI->getPointer();
 906   DataArrayInt *conn=DataArrayInt::New();
 907   conn->alloc(5*n1*n2,1);
 908   ci[0]=0;
 909   int *cp=conn->getPointer();
 910   int pos=0;
 911   for(int j=0;j<n2;j++)
 912     for(int i=0;i<n1;i++,pos++)
 913       {
 914         cp[5*pos]=(int)INTERP_KERNEL::NORM_QUAD4;
 915         cp[5*pos+1]=i+1+j*(n1+1);
 916         cp[5*pos+2]=i+j*(n1+1);
 917         cp[5*pos+3]=i+(j+1)*(n1+1);
 918         cp[5*pos+4]=i+1+(j+1)*(n1+1);
 919         ci[pos+1]=5*(pos+1);
 920     }
 921   m->setConnectivity(conn,connI,true);
 922   conn->decrRef();
 923   connI->decrRef();
 924 }
 925
 926 void MEDCouplingCMesh::fill3DUnstructuredMesh(MEDCouplingUMesh *m) const
 927 {
 928   const DataArrayDouble *c1=getCoordsAt(0);
 929   const DataArrayDouble *c2=getCoordsAt(1);
 930   const DataArrayDouble *c3=getCoordsAt(2);
 931   int n1=c1->getNbOfElems()-1;
 932   int n2=c2->getNbOfElems()-1;
 933   int n3=c3->getNbOfElems()-1;
 934   DataArrayInt *connI=DataArrayInt::New();
 935   connI->alloc(n1*n2*n3+1,1);
 936   int *ci=connI->getPointer();
 937   DataArrayInt *conn=DataArrayInt::New();
 938   conn->alloc(9*n1*n2*n3,1);
 939   ci[0]=0;
 940   int *cp=conn->getPointer();
 941   int pos=0;
 942   for(int k=0;k<n3;k++)
 943     for(int j=0;j<n2;j++)
 944       for(int i=0;i<n1;i++,pos++)
 945         {
 946           cp[9*pos]=(int)INTERP_KERNEL::NORM_HEXA8;
 947           int tmp=(n1+1)*(n2+1);
 948           cp[9*pos+1]=i+1+j*(n1+1)+k*tmp;
 949           cp[9*pos+2]=i+j*(n1+1)+k*tmp;
 950           cp[9*pos+3]=i+(j+1)*(n1+1)+k*tmp;
 951           cp[9*pos+4]=i+1+(j+1)*(n1+1)+k*tmp;
 952           cp[9*pos+5]=i+1+j*(n1+1)+(k+1)*tmp;
 953           cp[9*pos+6]=i+j*(n1+1)+(k+1)*tmp;
 954           cp[9*pos+7]=i+(j+1)*(n1+1)+(k+1)*tmp;
 955           cp[9*pos+8]=i+1+(j+1)*(n1+1)+(k+1)*tmp;
 956           ci[pos+1]=9*(pos+1);
 957         }
 958   m->setConnectivity(conn,connI,true);
 959   conn->decrRef();
 960   connI->decrRef();
 961 }
 962
 963 void MEDCouplingCMesh::getTinySerializationInformation(std::vector<double>& tinyInfoD, std::vector<int>& tinyInfo, std::vector<std::string>& littleStrings) const
 964 {
 965   int it,order;
 966   double time=getTime(it,order);
 967   tinyInfo.clear();
 968   tinyInfoD.clear();
 969   littleStrings.clear();
 970   littleStrings.push_back(getName());
 971   littleStrings.push_back(getDescription());
 972   littleStrings.push_back(getTimeUnit());
 973   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
 974   for(int i=0;i<3;i++)
 975     {
 976       int val=-1;
 977       std::string st;
 978       if(thisArr[i])
 979         {
 980           val=thisArr[i]->getNumberOfTuples();
 981           st=thisArr[i]->getInfoOnComponent(0);
 982         }
 983       tinyInfo.push_back(val);
 984       littleStrings.push_back(st);
 985     }
 986   tinyInfo.push_back(it);
 987   tinyInfo.push_back(order);
 988   tinyInfoD.push_back(time);
 989 }
 990
 991 void MEDCouplingCMesh::resizeForUnserialization(const std::vector<int>& tinyInfo, DataArrayInt *a1, DataArrayDouble *a2, std::vector<std::string>& littleStrings) const
 992 {
 993   a1->alloc(0,1);
 994   int sum=0;
 995   for(int i=0;i<3;i++)
 996     if(tinyInfo[i]!=-1)
 997       sum+=tinyInfo[i];
 998   a2->alloc(sum,1);
 999 }
1000
1001 void MEDCouplingCMesh::serialize(DataArrayInt *&a1, DataArrayDouble *&a2) const
1002 {
1003   a1=DataArrayInt::New();
1004   a1->alloc(0,1);
1005   const DataArrayDouble *thisArr[3]={_x_array,_y_array,_z_array};
1006   int sz=0;
1007   for(int i=0;i<3;i++)
1008     {
1009       if(thisArr[i])
1010         sz+=thisArr[i]->getNumberOfTuples();
1011     }
1012   a2=DataArrayDouble::New();
1013   a2->alloc(sz,1);
1014   double *a2Ptr=a2->getPointer();
1015   for(int i=0;i<3;i++)
1016     if(thisArr[i])
1017       a2Ptr=std::copy(thisArr[i]->getConstPointer(),thisArr[i]->getConstPointer()+thisArr[i]->getNumberOfTuples(),a2Ptr);
1018 }
1019
1020 void MEDCouplingCMesh::unserialization(const std::vector<double>& tinyInfoD, const std::vector<int>& tinyInfo, const DataArrayInt *a1, DataArrayDouble *a2,
1021                                        const std::vector<std::string>& littleStrings)
1022 {
1023   setName(littleStrings[0].c_str());
1024   setDescription(littleStrings[1].c_str());
1025   setTimeUnit(littleStrings[2].c_str());
1026   DataArrayDouble **thisArr[3]={&_x_array,&_y_array,&_z_array};
1027   const double *data=a2->getConstPointer();
1028   for(int i=0;i<3;i++)
1029     {
1030       if(tinyInfo[i]!=-1)
1031         {
1032           (*(thisArr[i]))=DataArrayDouble::New();
1033           (*(thisArr[i]))->alloc(tinyInfo[i],1);
1034           (*(thisArr[i]))->setInfoOnComponent(0,littleStrings[i+3].c_str());
1035           std::copy(data,data+tinyInfo[i],(*(thisArr[i]))->getPointer());
1036           data+=tinyInfo[i];
1037         }
1038     }
1039   setTime(tinyInfoD[0],tinyInfo[3],tinyInfo[4]);
1040 }
1041
1042 void MEDCouplingCMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData) const throw(INTERP_KERNEL::Exception)
1043 {
1044   throw INTERP_KERNEL::Exception("MEDCouplingCMesh::writeVTKLL : not implemented yet !");
1045 }
1046
1047 std::string MEDCouplingCMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exception)
1048 {
1049   return std::string("RectilinearGrid");
1050 }