CDMATH/mesh/src/Mesh.cxx

   1 /*
   2  * mesh.cxx
   3  *
   4  *  Created on: 22 janv. 2012
   5  *      Authors: CDMATH
   6  */
   7
   8 #include "Mesh.hxx"
   9 #include "Node.hxx"
  10 #include "Cell.hxx"
  11 #include "Face.hxx"
  12
  13 #include "CdmathException.hxx"
  14 #include "MEDFileMesh.hxx"
  15 #include "MEDLoader.hxx"
  16
  17 #include <iostream>
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20 #include <iterator>
  21 #include <algorithm>
  22 #include <string.h>
  23
  24 using namespace MEDCoupling;
  25 using namespace std;
  26
  27 //----------------------------------------------------------------------
  28 Mesh::Mesh( void )
  29 //----------------------------------------------------------------------
  30 {
  31         _mesh=NULL;
  32         _cells=NULL;
  33         _nodes=NULL;
  34         _faces=NULL;
  35         _spaceDim = 0 ;
  36         _meshDim  = 0 ;
  37         _numberOfNodes = 0;
  38         _numberOfFaces = 0;
  39         _numberOfCells = 0;
  40         _numberOfEdges = 0;
  41     _isStructured=false;
  42         _xMin=0.;
  43         _xMax=0.;
  44         _yMin=0.;
  45         _yMax=0.;
  46         _zMin=0.;
  47         _zMax=0.;
  48     _nxyz.resize(0);
  49     _dxyz.resize(0.);
  50         _faceGroupNames.resize(0);
  51         _faceGroups.resize(0);
  52         _nodeGroupNames.resize(0);
  53         _nodeGroups.resize(0);
  54     _indexFacePeriodicSet=false;
  55     _name="";
  56 }
  57
  58 //----------------------------------------------------------------------
  59 Mesh::~Mesh( void )
  60 //----------------------------------------------------------------------
  61 {
  62         delete [] _cells;
  63         delete [] _nodes;
  64         delete [] _faces;
  65         //for(int i=0; i< _faceGroups.size(); i++)
  66         //      _faceGroups[i]->decrRef();
  67         //      _nodeGroups[i]->decrRef();
  68 }
  69
  70 std::string
  71 Mesh::getName( void ) const
  72 {
  73     return _name;
  74 }
  75
  76 Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh )
  77 {
  78         _spaceDim=mesh->getSpaceDimension();
  79         _meshDim=mesh->getMeshDimension();
  80     _isStructured=true;
  81         _dxyz=mesh->getDXYZ();
  82         _nxyz=mesh->getCellGridStructure();
  83         double* Box0=new double[2*_spaceDim];
  84         mesh->getBoundingBox(Box0);
  85     _name=mesh->getName();
  86     _indexFacePeriodicSet=false;
  87
  88         _xMin=Box0[0];
  89         _xMax=Box0[1];
  90         if (_spaceDim>=2)
  91         {
  92                 _yMin=Box0[2];
  93                 _yMax=Box0[3];
  94         }
  95         if (_spaceDim>=3)
  96         {
  97                 _zMin=Box0[4];
  98                 _zMax=Box0[5];
  99         }
 100
 101         double *originPtr = new double[_spaceDim];
 102         double *dxyzPtr = new double[_spaceDim];
 103         mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
 104
 105         for(int i=0;i<_spaceDim;i++)
 106         {
 107                 originPtr[i]=Box0[2*i];
 108                 nodeStrctPtr[i]=_nxyz[i]+1;
 109                 dxyzPtr[i]=_dxyz[i];
 110         }
 111         _mesh=MEDCouplingIMesh::New(_name,
 112                         _spaceDim,
 113                         nodeStrctPtr,
 114                         nodeStrctPtr+_spaceDim,
 115                         originPtr,
 116                         originPtr+_spaceDim,
 117                         dxyzPtr,
 118                         dxyzPtr+_spaceDim);
 119         delete [] originPtr;
 120         delete [] dxyzPtr;
 121         delete [] nodeStrctPtr;
 122         delete [] Box0 ;
 123         setMesh();
 124 }
 125
 126 Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh )
 127 {
 128         _spaceDim=mesh->getSpaceDimension();
 129         _meshDim=mesh->getMeshDimension();
 130     _isStructured=false;
 131         double* Box0=new double[2*_spaceDim];
 132         mesh->getBoundingBox(Box0);
 133     _name=mesh->getName();
 134     _indexFacePeriodicSet=false;
 135
 136         _xMin=Box0[0];
 137         _xMax=Box0[1];
 138         if (_spaceDim>=2)
 139         {
 140                 _yMin=Box0[2];
 141                 _yMax=Box0[3];
 142         }
 143         if (_spaceDim>=3)
 144         {
 145                 _zMin=Box0[4];
 146                 _zMax=Box0[5];
 147         }
 148
 149         _mesh=mesh->deepCopy();
 150         delete [] Box0 ;
 151         setMesh();
 152 }
 153
 154 //----------------------------------------------------------------------
 155 Mesh::Mesh( const Mesh& mesh )
 156 //----------------------------------------------------------------------
 157 {
 158         _spaceDim = mesh.getSpaceDimension() ;
 159         _meshDim = mesh.getMeshDimension() ;
 160     _name=mesh.getName();
 161     _xMax=mesh.getXMax();
 162     _yMin=mesh.getYMin();
 163     _yMax=mesh.getYMax();
 164     _zMin=mesh.getZMin();
 165     _zMax=mesh.getZMax();
 166     _isStructured=mesh.isStructured();
 167     if(_isStructured)
 168     {
 169         _nxyz = mesh.getCellGridStructure() ;
 170         _dxyz=mesh.getDXYZ();
 171     }
 172         _numberOfNodes = mesh.getNumberOfNodes();
 173         _numberOfFaces = mesh.getNumberOfFaces();
 174         _numberOfCells = mesh.getNumberOfCells();
 175         _numberOfEdges = mesh.getNumberOfEdges();
 176
 177         _faceGroupNames = mesh.getNameOfFaceGroups() ;
 178         _faceGroups = mesh.getFaceGroups() ;
 179         _nodeGroupNames = mesh.getNameOfNodeGroups() ;
 180         _nodeGroups = mesh.getNodeGroups() ;
 181
 182         _nodes   = new Node[_numberOfNodes] ;
 183         _faces   = new Face[_numberOfFaces] ;
 184         _cells   = new Cell[_numberOfCells] ;
 185
 186     //memcpy(_nodes,mesh.getNodes(),sizeof(*mesh.getNodes())) ;
 187     //memcpy(_cells,mesh.getCells(),sizeof(*mesh.getCells())) ;
 188     //memcpy(_faces,mesh.getFaces(),sizeof(*mesh.getFaces())) ;
 189
 190         for (int i=0;i<_numberOfNodes;i++)
 191                 _nodes[i]=mesh.getNode(i);
 192
 193         for (int i=0;i<_numberOfFaces;i++)
 194                 _faces[i]=mesh.getFace(i);
 195
 196         for (int i=0;i<_numberOfCells;i++)
 197                 _cells[i]=mesh.getCell(i);
 198
 199     _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
 200     if(_indexFacePeriodicSet)
 201         _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
 202
 203     _boundaryFaceIds=mesh.getBoundaryFaceIds();
 204     _boundaryNodeIds=mesh.getBoundaryNodeIds();
 205
 206     _eltsTypes=mesh.getElementTypes();
 207     _eltsTypesNames=mesh.getElementTypesNames();
 208
 209         MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
 210         _mesh=m1;
 211 }
 212
 213 //----------------------------------------------------------------------
 214 Mesh::Mesh( const std::string filename, int meshLevel )
 215 //----------------------------------------------------------------------
 216 {
 217         readMeshMed(filename, meshLevel);
 218 }
 219
 220 //----------------------------------------------------------------------
 221 void
 222 Mesh::readMeshMed( const std::string filename, const int meshLevel)
 223 //----------------------------------------------------------------------
 224 {
 225         MEDFileUMesh *m=MEDFileUMesh::New(filename.c_str());//reads the first mesh encountered in the file, otherwise call New (const char *fileName, const char *mName, int dt=-1, int it=-1)
 226         _mesh=m->getMeshAtLevel(meshLevel);
 227     _mesh->checkConsistencyLight();
 228         _mesh->setName(_mesh->getName());
 229         _meshDim=_mesh->getMeshDimension();
 230         _spaceDim=_mesh->getSpaceDimension();
 231     _name=_mesh->getName();
 232     _indexFacePeriodicSet=false;
 233     MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingIMesh*> (_mesh.retn());
 234     if(structuredMesh)
 235     {
 236         _isStructured=true;
 237         _dxyz=structuredMesh->getDXYZ();
 238         _nxyz=structuredMesh->getCellGridStructure();
 239         double* Box0=new double[2*_spaceDim];
 240         structuredMesh->getBoundingBox(Box0);
 241
 242         _xMin=Box0[0];
 243         _xMax=Box0[1];
 244         if (_spaceDim>=2)
 245         {
 246             _yMin=Box0[2];
 247             _yMax=Box0[3];
 248         }
 249         if (_spaceDim>=3)
 250         {
 251             _zMin=Box0[4];
 252             _zMax=Box0[5];
 253         }
 254     }
 255     else
 256         _isStructured=false;
 257
 258         MEDCouplingUMesh*  mu = setMesh();
 259         setGroups(m, mu);
 260
 261         cout<<endl<< "Loaded file "<< filename<<endl;
 262     cout<<"Mesh name= "<<m->getName()<<", mesh dim="<< _meshDim<< ", space dim="<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
 263
 264         m->decrRef();
 265         mu->decrRef();
 266 }
 267
 268 void
 269 Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName)
 270 {
 271         int nbFace=getNumberOfFaces();
 272         bool flag=false;
 273         for (int iface=0;iface<nbFace;iface++)
 274         {
 275                 double FX=_faces[iface].x();
 276                 double FY=_faces[iface].y();
 277                 double FZ=_faces[iface].z();
 278                 if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
 279                 {
 280                         _faces[iface].setGroupName(groupName);
 281                         std::vector< int > nodesID= _faces[iface].getNodesId();
 282                         int nbNodes = _faces[iface].getNumberOfNodes();
 283                         for(int inode=0 ; inode<nbNodes ; inode++)
 284                                 _nodes[nodesID[inode]].setGroupName(groupName);
 285
 286                         flag=true;
 287                 }
 288         }
 289         if (flag)
 290     {
 291                 _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
 292                 _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
 293         //To do : update _faceGroups and _nodeGroups
 294     }
 295 }
 296
 297 void
 298 Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName)
 299 {
 300         int nbFace=getNumberOfFaces();
 301         bool flag=false;
 302         for (int iface=0;iface<nbFace;iface++)
 303         {
 304                 double cord=_faces[iface].getBarryCenter()[direction];
 305                 if (abs(cord-value)<eps)
 306                 {
 307                         _faces[iface].setGroupName(groupName);
 308                         std::vector< int > nodesID= _faces[iface].getNodesId();
 309                         int nbNodes = _faces[iface].getNumberOfNodes();
 310                         for(int inode=0 ; inode<nbNodes ; inode++)
 311                 {
 312                                 _nodes[nodesID[inode]].setGroupName(groupName);
 313                 }
 314
 315                         flag=true;
 316                 }
 317         }
 318         if (flag)
 319     {
 320                 _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
 321                 _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
 322         //To do : update _faceGroups, _nodeGroups
 323     }
 324 }
 325
 326 void
 327 Mesh::setBoundaryNodesFromFaces()
 328 {
 329     for (int iface=0;iface<_boundaryFaceIds.size();iface++)
 330     {
 331         std::vector< int > nodesID= _faces[_boundaryFaceIds[iface]].getNodesId();
 332         int nbNodes = _faces[_boundaryFaceIds[iface]].getNumberOfNodes();
 333         for(int inode=0 ; inode<nbNodes ; inode++)
 334         {
 335             std::vector<int>::const_iterator  it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]);
 336             if( it != _boundaryNodeIds.end() )
 337                 _boundaryNodeIds.push_back(nodesID[inode]);
 338         }
 339     }
 340 }
 341
 342 std::map<int,int>
 343 Mesh::getIndexFacePeriodic( void ) const
 344 {
 345     return _indexFacePeriodicMap;
 346 }
 347
 348 void
 349 Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
 350 {
 351     if(_indexFacePeriodicSet)
 352         return;
 353
 354     double eps=1.E-10;
 355
 356     for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
 357     {
 358         Face my_face=_faces[_boundaryFaceIds[indexFace]];
 359         int iface_perio=-1;
 360         if(_meshDim==1)
 361         {
 362             for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
 363                 if(iface!=indexFace)
 364                 {
 365                     iface_perio=_boundaryFaceIds[iface];
 366                     break;
 367                 }
 368         }
 369         else if(_meshDim==2)
 370         {
 371             double x=my_face.x();
 372             double y=my_face.y();
 373
 374             for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
 375             {
 376                 Face face_i=_faces[_boundaryFaceIds[iface]];
 377                 double xi=face_i.x();
 378                 double yi=face_i.y();
 379                 if (   (abs(y-yi)<eps || abs(x-xi)<eps )// Case of a square geometry
 380                     && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
 381                     && ( !use_central_inversion || abs(y+yi) + abs(x+xi)<eps ) // Case of a central inversion
 382                     && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
 383                     && fabs(my_face.getXN()      + face_i.getXN())<eps
 384                     && fabs(my_face.getYN()      + face_i.getYN())<eps
 385                     && fabs(my_face.getZN()      + face_i.getZN())<eps )
 386                 {
 387                     iface_perio=_boundaryFaceIds[iface];
 388                     break;
 389                 }
 390             }
 391         }
 392         else if(_meshDim==3)
 393         {
 394             double x=my_face.x();
 395             double y=my_face.y();
 396             double z=my_face.z();
 397
 398             for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
 399             {
 400                 Face face_i=_faces[_boundaryFaceIds[iface]];
 401                 double xi=face_i.x();
 402                 double yi=face_i.y();
 403                 double zi=face_i.z();
 404                 if ( ((abs(y-yi)<eps && abs(x-xi)<eps) || (abs(x-xi)<eps && abs(z-zi)<eps) || (abs(y-yi)<eps && abs(z-zi)<eps))// Case of a cube geometry
 405                     && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
 406                     && ( !use_central_inversion || abs(y+yi) + abs(x+xi) + abs(z+zi)<eps )// Case of a central inversion
 407                     && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
 408                     && fabs(my_face.getXN()      + face_i.getXN())<eps
 409                     && fabs(my_face.getYN()      + face_i.getYN())<eps
 410                     && fabs(my_face.getZN()      + face_i.getZN())<eps )
 411                 {
 412                     iface_perio=_boundaryFaceIds[iface];
 413                     break;
 414                 }
 415             }
 416         }
 417         else
 418             throw CdmathException("Mesh::setPeriodicFaces: Mesh dimension should be 1, 2 or 3");
 419
 420         if (iface_perio==-1)
 421             throw CdmathException("Mesh::setPeriodicFaces: periodic face not found, iface_perio==-1 " );
 422         else
 423             _indexFacePeriodicMap[_boundaryFaceIds[indexFace]]=iface_perio;
 424     }
 425     _indexFacePeriodicSet=true;
 426 }
 427
 428 int
 429 Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion)
 430 {
 431         if (!_faces[indexFace].isBorder())
 432         {
 433             cout<<"Pb with indexFace= "<<indexFace<<endl;
 434             throw CdmathException("Mesh::getIndexFacePeriodic: not a border face" );
 435         }
 436
 437     if(!_indexFacePeriodicSet)
 438         setPeriodicFaces(check_groups, use_central_inversion);
 439
 440     std::map<int,int>::const_iterator  it = _indexFacePeriodicMap.find(indexFace);
 441     if( it != _indexFacePeriodicMap.end() )
 442         return it->second;
 443     else
 444     {
 445         cout<<"Pb with indexFace= "<<indexFace<<endl;
 446         throw CdmathException("Mesh::getIndexFacePeriodic: not a periodic face" );
 447     }
 448 }
 449
 450 bool
 451 Mesh::isBorderNode(int nodeid) const
 452 {
 453         return _nodes[nodeid].isBorder();
 454 }
 455
 456 bool
 457 Mesh::isBorderFace(int faceid) const
 458 {
 459         return _faces[faceid].isBorder();
 460 }
 461
 462 std::vector< int >
 463 Mesh::getBoundaryFaceIds() const
 464 {
 465     return _boundaryFaceIds;
 466 }
 467
 468 std::vector< int >
 469 Mesh::getBoundaryNodeIds() const
 470 {
 471     return _boundaryNodeIds;
 472 }
 473
 474 bool
 475 Mesh::isTriangular() const
 476 {
 477         return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TRI3;
 478 }
 479 bool
 480 Mesh::isTetrahedral() const
 481 {
 482         return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TETRA4;
 483 }
 484 bool
 485 Mesh::isQuadrangular() const
 486 {
 487         return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_QUAD4;
 488 }
 489 bool
 490 Mesh::isHexahedral() const
 491 {
 492         return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_HEXA8;
 493 }
 494 bool
 495 Mesh::isStructured() const
 496 {
 497         return _isStructured;
 498 }
 499
 500 std::vector< INTERP_KERNEL::NormalizedCellType >
 501 Mesh::getElementTypes() const
 502 {
 503   return _eltsTypes;
 504 }
 505
 506 std::vector< string >
 507 Mesh::getElementTypesNames() const
 508 {
 509     std::vector< string > result(0);
 510     for(int i=0; i< _eltsTypes.size(); i++)
 511     {
 512         if( _eltsTypes[i]==INTERP_KERNEL::NORM_POINT1)
 513             result.push_back("Points");
 514         else if( _eltsTypes[i]==INTERP_KERNEL::NORM_SEG2)
 515             result.push_back("Segments");
 516         else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TRI3)
 517             result.push_back("Triangles");
 518         else if( _eltsTypes[i]==INTERP_KERNEL::NORM_QUAD4)
 519             result.push_back("Quadrangles");
 520         else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYGON)
 521             result.push_back("Polygons");
 522         else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TETRA4)
 523             result.push_back("Tetrahedra");
 524         else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PYRA5)
 525             result.push_back("Pyramids");
 526         else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PENTA6)
 527             result.push_back("Pentahedra");
 528         else if( _eltsTypes[i]==INTERP_KERNEL::NORM_HEXA8)
 529             result.push_back("Hexahedra");
 530         else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYHED)
 531             result.push_back("Polyhedrons");
 532         else
 533                 {
 534                         cout<< "Mesh " + _name + " contains an element of type " <<endl;
 535                         cout<< _eltsTypes[i]<<endl;
 536                         throw CdmathException("Mesh::getElementTypes : recognised cell med types are NORM_POINT1, NORM_SEG2, NORM_TRI3, NORM_QUAD4, NORM_TETRA4, NORM_PYRA5, NORM_PENTA6, NORM_HEXA8, NORM_POLYGON, NORM_POLYHED");
 537         }
 538     }
 539     return result;
 540 }
 541
 542 void
 543 Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh*  mu)
 544 {
 545         //Searching for face groups
 546         vector<string> faceGroups=medmesh->getGroupsNames() ;
 547
 548         for (unsigned int i=0;i<faceGroups.size();i++ )
 549         {
 550                 string groupName=faceGroups[i];
 551                 vector<mcIdType> nonEmptyGrp(medmesh->getGrpNonEmptyLevels(groupName));
 552                 //We check if the group has a relative dimension equal to -1
 553                 //before call to the function getGroup(-1,groupName.c_str())
 554                 vector<mcIdType>::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1);
 555                 if (it != nonEmptyGrp.end())
 556                 {
 557                         cout<<"Boundary face group named "<< groupName << " found"<<endl;
 558                         MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str());
 559             mu->unPolyze();
 560                         _faceGroups.insert(_faceGroups.begin(),m);
 561                         _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
 562                         DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
 563                         const double *coorBary=baryCell->getConstPointer();
 564
 565                         int nbCellsSubMesh=m->getNumberOfCells();
 566                         for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
 567                         {
 568                                 vector<double> coorBaryXyz(3,0);
 569                                 for (int d=0; d<_spaceDim; d++)
 570                                         coorBaryXyz[d] = coorBary[k+d];
 571                                 Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
 572
 573                                 int flag=0;
 574                 /* double min=INFINITY;
 575                 Point p3;
 576                 int closestFaceNb;*/
 577                                 for (int iface=0;iface<_numberOfFaces;iface++ )
 578                                 {
 579                                         Point p2=_faces[iface].getBarryCenter();
 580                     /*if(groupName=="Wall" and p1.distance(p2)<min)
 581                     {
 582                         min=p1.distance(p2);
 583                         p3=p2;
 584                         closestFaceNb=iface;
 585                     }*/
 586                                         if(p1.distance(p2)<1.E-10)
 587                                         {
 588                                                 _faces[iface].setGroupName(groupName);
 589                                                 flag=1;
 590                                                 break;
 591                                         }
 592                                 }
 593                 /* if(groupName=="Wall" and min>1.E-10)
 594                     {
 595                         cout.precision(15);
 596                         cout<<"Subcell number "<< ic <<" Min distance to Wall = "<<min <<" p= "<<p1[0]<<" , "<<p1[1]<<" , "<<p1[2]<<endl;
 597                         cout<<" attained at face "<< closestFaceNb << " p_min= "<<p3[0]<<" , "<<p3[1]<<" , "<<p3[2]<<endl;
 598                     }*/
 599                                 if (flag==0)
 600                                         throw CdmathException("No face belonging to group " + groupName + " found");
 601                         }
 602                         baryCell->decrRef();
 603                         //m->decrRef();
 604                 }
 605         }
 606
 607         //Searching for node groups
 608         vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
 609
 610         for (unsigned int i=0;i<nodeGroups.size();i++ )
 611         {
 612                 string groupName=nodeGroups[i];
 613                 DataArrayIdType * nodeGroup=medmesh->getNodeGroupArr( groupName );
 614                 const mcIdType *nodeids=nodeGroup->getConstPointer();
 615
 616                 if(nodeids!=NULL)
 617                 {
 618                         cout<<"Boundary node group named "<< groupName << " found"<<endl;
 619
 620                         _nodeGroups.insert(_nodeGroups.begin(),nodeGroup);
 621                         _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
 622
 623                         int nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
 624
 625                         DataArrayDouble *coo = mu->getCoords() ;
 626                         const double *cood=coo->getConstPointer();
 627
 628                         for (int ic(0); ic<nbNodesSubMesh; ic++)
 629                         {
 630                                 vector<double> coorP(3,0);
 631                                 for (int d=0; d<_spaceDim; d++)
 632                                         coorP[d] = cood[nodeids[ic]*_spaceDim+d];
 633                                 Point p1(coorP[0],coorP[1],coorP[2]) ;
 634
 635                                 int flag=0;
 636                                 for (int inode=0;inode<_numberOfNodes;inode++ )
 637                                 {
 638                                         Point p2=_nodes[inode].getPoint();
 639                                         if(p1.distance(p2)<1.E-10)
 640                                         {
 641                                                 _nodes[inode].setGroupName(groupName);
 642                                                 flag=1;
 643                                                 break;
 644                                         }
 645                                 }
 646                                 if (flag==0)
 647                                         throw CdmathException("No node belonging to group " + groupName + " found");
 648                         }
 649                 }
 650         }
 651 }
 652
 653 //----------------------------------------------------------------------
 654 MEDCouplingUMesh*
 655 Mesh::setMesh( void )
 656 //----------------------------------------------------------------------
 657 {
 658         DataArrayIdType *desc  = DataArrayIdType::New();
 659         DataArrayIdType *descI = DataArrayIdType::New();
 660         DataArrayIdType *revDesc  = DataArrayIdType::New();
 661         DataArrayIdType *revDescI = DataArrayIdType::New();
 662         MEDCouplingUMesh* mu = _mesh->buildUnstructured();
 663
 664         mu->unPolyze();
 665         /* Save the boundary mesh for later use*/
 666         _boundaryMesh = mu->computeSkin();
 667
 668         MEDCouplingUMesh* mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);//mesh of dimension N-1 containing the cell interfaces
 669
 670     const mcIdType *tmp = desc->getConstPointer();
 671     const mcIdType *tmpI=descI->getConstPointer();
 672
 673         const mcIdType *tmpA =revDesc->getConstPointer();
 674         const mcIdType *tmpAI=revDescI->getConstPointer();
 675
 676     //const int *work=tmp+tmpI[id];//corresponds to buildDescendingConnectivity
 677
 678         //Test du type d'éléments contenus dans le maillage afin d'éviter les éléments contenant des points de gauss
 679         _eltsTypes=mu->getAllGeoTypesSorted();
 680         for(int i=0; i<_eltsTypes.size();i++)
 681         {
 682                 if(
 683                                    _eltsTypes[i]!= INTERP_KERNEL::NORM_POINT1 && _eltsTypes[i]!= INTERP_KERNEL::NORM_SEG2
 684                                 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TRI3   && _eltsTypes[i]!= INTERP_KERNEL::NORM_QUAD4
 685                                 && _eltsTypes[i]!= INTERP_KERNEL::NORM_TETRA4 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PYRA5
 686                                 && _eltsTypes[i]!= INTERP_KERNEL::NORM_PENTA6 && _eltsTypes[i]!= INTERP_KERNEL::NORM_HEXA8
 687                                 && _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYGON&& _eltsTypes[i]!= INTERP_KERNEL::NORM_POLYHED
 688                 )
 689                 {
 690                         cout<< "Mesh " + mu->getName() + " contains an element of type " <<endl;
 691                         cout<< _eltsTypes[i]<<endl;
 692                         throw CdmathException("Mesh::setMesh : in order to avoid gauss points, mesh should contain elements of type NORM_POINT1, NORM_SEG2, NORM_TRI3, NORM_QUAD4, NORM_TETRA4, NORM_PYRA5, NORM_PENTA6, NORM_HEXA8, NORM_POLYGON, NORM_POLYHED");
 693                 }
 694         }
 695
 696         DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
 697         const double *coorBary=baryCell->getConstPointer();
 698
 699         MEDCouplingFieldDouble* fields=mu->getMeasureField(true);
 700         DataArrayDouble *surface = fields->getArray();
 701         const double *surf=surface->getConstPointer();
 702
 703         DataArrayDouble *coo = mu->getCoords() ;
 704         const double    *cood=coo->getConstPointer();
 705
 706         DataArrayIdType *revNode =DataArrayIdType::New();
 707         DataArrayIdType *revNodeI=DataArrayIdType::New();
 708         mu->getReverseNodalConnectivity(revNode,revNodeI) ;
 709         const mcIdType *tmpN =revNode->getConstPointer();
 710         const mcIdType *tmpNI=revNodeI->getConstPointer();
 711
 712         DataArrayIdType *revCell =DataArrayIdType::New();
 713         DataArrayIdType *revCellI=DataArrayIdType::New();
 714         mu2->getReverseNodalConnectivity(revCell,revCellI) ;
 715         const mcIdType *tmpC =revCell->getConstPointer();
 716         const mcIdType *tmpCI=revCellI->getConstPointer();
 717
 718         const DataArrayIdType *nodal  = mu2->getNodalConnectivity() ;
 719         const DataArrayIdType *nodalI = mu2->getNodalConnectivityIndex() ;
 720         const mcIdType *tmpNE =nodal->getConstPointer();
 721         const mcIdType *tmpNEI=nodalI->getConstPointer();
 722
 723         _numberOfCells = mu->getNumberOfCells() ;
 724         _cells      = new Cell[_numberOfCells] ;
 725
 726         _numberOfNodes = mu->getNumberOfNodes() ;
 727         _nodes      = new Node[_numberOfNodes] ;
 728
 729         _numberOfFaces = mu2->getNumberOfCells();
 730         _faces       = new Face[_numberOfFaces] ;
 731
 732     _indexFacePeriodicSet=false;
 733
 734     //Definition used if _meshDim =3 to determine the edges
 735     DataArrayIdType *desc2 =DataArrayIdType::New();
 736     DataArrayIdType *descI2=DataArrayIdType::New();
 737     DataArrayIdType *revDesc2 =DataArrayIdType::New();
 738     DataArrayIdType *revDescI2=DataArrayIdType::New();
 739     DataArrayIdType *revNode2 =DataArrayIdType::New();
 740     DataArrayIdType *revNodeI2=DataArrayIdType::New();
 741     const mcIdType *tmpN2 ;
 742     const mcIdType *tmpNI2;
 743     MEDCouplingUMesh* mu3;
 744
 745         if (_meshDim == 1)
 746         _numberOfEdges = mu->getNumberOfCells();
 747     else if (_meshDim == 2)
 748         _numberOfEdges = mu2->getNumberOfCells();
 749     else
 750     {
 751         mu3=mu2->buildDescendingConnectivity(desc2,descI2,revDesc2,revDescI2);//1D mesh of segments
 752         _numberOfEdges = mu3->getNumberOfCells();
 753         mu3->getReverseNodalConnectivity(revNode2,revNodeI2) ;
 754         tmpN2 =revNode2->getConstPointer();
 755         tmpNI2=revNodeI2->getConstPointer();
 756     }
 757
 758         // _cells, _nodes and _faces initialization:
 759         if (_spaceDim == 1)
 760         {
 761                 for( int id=0;id<_numberOfCells;id++ )
 762                 {
 763                         const mcIdType *work=tmp+tmpI[id];
 764             int nbFaces=tmpI[id+1]-tmpI[id];
 765
 766                         int nbVertices=mu->getNumberOfNodesInCell(id) ;
 767
 768                         Cell ci( nbVertices, nbFaces, surf[id], Point(coorBary[id], 0.0, 0.0) ) ;
 769
 770                         std::vector<mcIdType> nodeIdsOfCell ;
 771                         mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
 772                         for( int el=0;el<nbVertices;el++ )
 773             {
 774                                 ci.addNodeId(el,nodeIdsOfCell[el]) ;
 775                 ci.addFaceId(el,nodeIdsOfCell[el]) ;
 776             }
 777
 778             //This assumes that in a 1D mesh the cell numbers form a consecutive sequence
 779             //going either from left to right (xn=-1) or right to left (xn=1)
 780                         double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
 781
 782                         for( int el=0;el<nbFaces;el++ )
 783                         {
 784                                 ci.addNormalVector(el,xn,0.0,0.0) ;
 785                                 int indexFace=abs(work[el])-1;
 786                 ci.addFaceId(el,indexFace) ;
 787                                 xn = - xn;
 788                         }
 789                         _cells[id] = ci ;
 790                 }
 791
 792                 for( int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
 793                 {
 794                         Point p(cood[k], 0.0, 0.0) ;
 795                         const mcIdType *workc=tmpN+tmpNI[id];
 796                         mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
 797
 798                         const mcIdType *workf=tmpC+tmpCI[id];
 799                         mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
 800             const mcIdType *workn=tmpA+tmpAI[id];
 801             mcIdType nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
 802                         Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
 803
 804                         for( int el=0;el<nbCells;el++ )
 805                                 vi.addCellId(el,workc[el]) ;
 806                         for( int el=0;el<nbFaces;el++ )
 807                                 vi.addFaceId(el,workf[el]) ;
 808                         for( int el=0;el<nbNeighbourNodes;el++ )
 809                                 vi.addNeighbourNodeId(el,workn[el]) ;
 810                         _nodes[id] = vi ;
 811                 }
 812         _boundaryNodeIds.push_back(0);
 813         _boundaryNodeIds.push_back(_numberOfNodes-1);
 814
 815                 for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
 816                 {
 817                         Point p(cood[k], 0.0, 0.0) ;
 818                         const mcIdType *workc=tmpA+tmpAI[id];
 819                         mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
 820
 821                         const mcIdType *workv=tmpNE+tmpNEI[id]+1;
 822                         Face fi( 1, nbCells, 1.0, p, 1.0, 0.0, 0.0) ;
 823                         fi.addNodeId(0,workv[0]) ;
 824
 825                         for(int idCell=0; idCell<nbCells; idCell++)
 826                                 fi.addCellId(idCell,workc[idCell]) ;
 827
 828                         _faces[id] = fi ;
 829                 }
 830         _boundaryFaceIds.push_back(0);
 831         _boundaryFaceIds.push_back(_numberOfFaces-1);
 832         }
 833         else if(_spaceDim==2  || _spaceDim==3)
 834         {
 835                 DataArrayDouble *barySeg = mu2->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
 836                 const double *coorBarySeg=barySeg->getConstPointer();
 837
 838                 MEDCouplingFieldDouble* fieldn;
 839                 DataArrayDouble *normal;
 840                 const double *tmpNormal;
 841
 842                 if(_spaceDim==_meshDim)
 843                         fieldn = mu2->buildOrthogonalField();//Compute the normal to each cell interface
 844                 else
 845                         fieldn = mu->buildOrthogonalField();//compute the 3D normal vector to the 2D cell
 846
 847                 normal = fieldn->getArray();
 848                 tmpNormal = normal->getConstPointer();
 849
 850                 /*Building mesh cells */
 851                 for(int id(0), k(0); id<_numberOfCells; id++, k+=_spaceDim)
 852                 {
 853             const mcIdType *work=tmp+tmpI[id];
 854                         mcIdType nbFaces=tmpI[id+1]-tmpI[id];
 855
 856                         mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;
 857
 858                         vector<double> coorBaryXyz(3,0);
 859                         for (int d=0; d<_spaceDim; d++)
 860                                 coorBaryXyz[d] = coorBary[k+d];
 861
 862                         Point p(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
 863                         Cell ci( nbVertices, nbFaces, surf[id], p ) ;
 864
 865                         /* Filling cell nodes */
 866                         std::vector<mcIdType> nodeIdsOfCell ;
 867                         mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
 868                         for( int el=0;el<nbVertices;el++ )
 869                                 ci.addNodeId(el,nodeIdsOfCell[el]) ;
 870
 871                         /* Filling cell faces */
 872                         if(_spaceDim==_meshDim)//use the normal field generated by buildOrthogonalField()
 873                                 for( int el=0;el<nbFaces;el++ )
 874                                 {
 875                     mcIdType faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
 876                                         vector<double> xyzn(3,0);//Outer normal to the cell
 877                                         if (work[el]<0)
 878                                                 for (int d=0; d<_spaceDim; d++)
 879                                                         xyzn[d] = -tmpNormal[_spaceDim*faceIndex+d];
 880                                         else
 881                                                 for (int d=0; d<_spaceDim; d++)
 882                                                         xyzn[d] = +tmpNormal[_spaceDim*faceIndex+d];
 883                                         ci.addNormalVector(el,xyzn[0],xyzn[1],xyzn[2]) ;
 884                                         ci.addFaceId(el,faceIndex) ;
 885                                 }
 886                         else//build normals associated to the couple (cell id, face el)
 887                         {
 888                                 if(_meshDim==1)//we know in this case there are only two faces around the cell id, each face is composed of a single node
 889                                 {//work[0]= first face global number, work[1]= second face global number
 890                     mcIdType indexFace0=abs(work[0])-1;//=work[0] since Fortran type numbering was used, and negative sign means anticlockwise numbering
 891                     mcIdType indexFace1=abs(work[1])-1;//=work[1] since Fortran type numbering was used, and negative sign means anticlockwise numbering
 892                                         mcIdType idNodeA=(tmpNE+tmpNEI[indexFace0]+1)[0];//global number of the first  face node work[0]=(abs(work[0])-1)
 893                                         mcIdType idNodeB=(tmpNE+tmpNEI[indexFace1]+1)[0];//global number of the second face node work[1]=(abs(work[1])-1)
 894                                         Vector vecAB(3);
 895                                         for(int i=0;i<_spaceDim;i++)
 896                                                 vecAB[i]=coo->getIJ(idNodeB,i) - coo->getIJ(idNodeA,i);
 897                                         vecAB/=vecAB.norm();
 898                                         ci.addNormalVector(0,-vecAB[0],-vecAB[1],-vecAB[2]) ;
 899                                         ci.addNormalVector(1,vecAB[0],vecAB[1],vecAB[2]) ;
 900                                         ci.addFaceId(0,indexFace0) ;
 901                                         ci.addFaceId(1,indexFace1) ;
 902                                 }
 903                                 else//_meshDim==2, number of faces around the cell id is variable, each face is composed of two nodes
 904                                 {
 905                                         Vector xyzn(3);
 906                                         for (int d=0; d<_spaceDim; d++)
 907                                                 xyzn[d] = tmpNormal[_spaceDim*id+d];
 908                                         for( int el=0;el<nbFaces;el++ )
 909                                         {
 910                         int faceIndex=(abs(work[el])-1);//=work[el] since Fortran type numbering was used, and negative sign means anticlockwise numbering
 911                                                 const mcIdType *workv=tmpNE+tmpNEI[faceIndex]+1;
 912                                                 mcIdType nbNodes= tmpNEI[faceIndex+1]-tmpNEI[faceIndex]-1;
 913                                                 if(nbNodes!=2)//We want to compute the normal to a straight line, not a curved interface composed of more thant 2 points
 914                                                 {
 915                                                         cout<<"Mesh name "<< mu->getName()<< " space dim= "<< _spaceDim <<" mesh dim= "<< _meshDim <<endl;
 916                                                         cout<<"For cell id "<<id<<" and local face number "<<el<<", the number of nodes is "<< nbNodes<< ", total number of faces is "<< nbFaces <<endl;
 917                                                         throw CdmathException("Mesh::setMesh number of nodes around a face should be 2");
 918                                                 }
 919
 920                                                 mcIdType idNodeA=workv[0];
 921                                                 mcIdType idNodeB=workv[1];
 922                                                 vector<double> nodeA(_spaceDim), nodeB(_spaceDim), nodeP(_spaceDim);
 923                                                 for(int i=0;i<_spaceDim;i++)
 924                                                 {
 925                                                         nodeA[i]=coo->getIJ(idNodeA,i);
 926                                                         nodeB[i]=coo->getIJ(idNodeB,i);
 927                                                         nodeP[i]=coorBary[_spaceDim*id+i];
 928                                                 }
 929                                                 //Let P be the barycenter of the cell id
 930                                                 Vector vecAB(3), vecPA(3);
 931                                                 for(int i=0;i<_spaceDim;i++)
 932                                                 {
 933                                                         vecAB[i]=coo->getIJ(idNodeB,i)       - coo->getIJ(idNodeA,i);
 934                                                         vecPA[i]=coo->getIJ(idNodeA,i) - coorBary[_spaceDim*id+i];
 935                                                 }
 936
 937                                                 Vector normale = xyzn % vecAB;//Normal to the edge
 938                                                 normale/=normale.norm();
 939
 940                                                 if(normale*vecPA<0)
 941                                                         ci.addNormalVector(el,normale[0],normale[1],normale[2]) ;
 942                                                 else
 943                                                         ci.addNormalVector(el,-normale[0],-normale[1],-normale[2]) ;
 944                                                 ci.addFaceId(el,faceIndex) ;
 945                                         }
 946                                 }
 947                         }
 948                         _cells[id] = ci ;
 949                 }
 950
 951                 MEDCouplingFieldDouble* fieldl=mu2->getMeasureField(true);
 952                 DataArrayDouble *longueur = fieldl->getArray();
 953                 const double *lon=longueur->getConstPointer();
 954
 955                 if(_spaceDim!=_meshDim)
 956                 {
 957                         /* Since spaceDim!=meshDim, don't build normal to faces */
 958                         fieldn->decrRef();
 959             normal=NULL;
 960             tmpNormal=NULL;
 961                 }
 962
 963                 /*Building mesh faces */
 964                 for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
 965                 {
 966                         vector<double> coorBarySegXyz(3,0);
 967                         for (int d=0; d<_spaceDim; d++)
 968                                 coorBarySegXyz[d] = coorBarySeg[k+d];
 969                         Point p(coorBarySegXyz[0],coorBarySegXyz[1],coorBarySegXyz[2]) ;
 970                         const mcIdType *workc=tmpA+tmpAI[id];
 971                         mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
 972
 973             if (nbCells>2 && _spaceDim==_meshDim)
 974             {
 975                 cout<<"Warning : nbCells>2, numberOfFaces="<<_numberOfFaces<<endl;
 976                 cout<<"nbCells= "<<nbCells<<", _spaceDim="<<_spaceDim<<", _meshDim="<<_meshDim<<endl;
 977                 for(int icell=0; icell<nbCells; icell++)
 978                     cout<<workc[icell]<<", ";
 979                 cout<<endl;
 980                 throw CdmathException("Wrong mesh : nbCells>2 and spaceDim==meshDim");
 981             }
 982             if (nbCells==1)
 983                 _boundaryFaceIds.push_back(id);
 984
 985                         const mcIdType *workv=tmpNE+tmpNEI[id]+1;
 986                         mcIdType nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;
 987
 988                         Face fi;
 989                         if(_spaceDim==_meshDim)//Euclidean flat mesh geometry
 990                 if(_spaceDim==2)
 991                     fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], 0.0) ;
 992                 else
 993                     fi=Face( nbNodes, nbCells, lon[id], p, tmpNormal[k], tmpNormal[k+1], tmpNormal[k+2]) ;
 994                         else//Curved mesh geometry
 995                                 fi=Face( nbNodes, nbCells, lon[id], p, 0.0, 0.0, 0.0) ;//Since spaceDim!=meshDim, normal to face is not defined
 996
 997                         for(int node_id=0; node_id<nbNodes;node_id++)
 998                                 fi.addNodeId(node_id,workv[node_id]) ;
 999
1000                         fi.addCellId(0,workc[0]) ;
1001                         for(int cell_id=1; cell_id<nbCells;cell_id++)
1002             {
1003                 int cell_idx=0;
1004                 if (workc[cell_id]!=workc[cell_id-1])//For some meshes (bad ones) the same cell can appear several times
1005                     {
1006                     fi.addCellId(cell_idx+1,workc[cell_id]) ;
1007                     cell_idx++;
1008                     }
1009             }
1010
1011                         _faces[id] = fi ;
1012                 }
1013
1014                 /*Building mesh nodes, should be done after building mesh faces in order to detect boundary nodes*/
1015                 for(int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
1016                 {
1017                         vector<double> coorP(3,0);
1018                         for (int d=0; d<_spaceDim; d++)
1019                                 coorP[d] = cood[k+d];
1020                         Point p(coorP[0],coorP[1],coorP[2]) ;
1021
1022                         const mcIdType *workc=tmpN+tmpNI[id];
1023                         mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
1024                         const mcIdType *workf=tmpC+tmpCI[id];
1025                         mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
1026                         const mcIdType *workn;
1027                         mcIdType nbNeighbourNodes;
1028             if (_meshDim == 1)
1029             {
1030                 workn=tmpA+tmpAI[id];
1031                 nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
1032             }
1033             else if (_meshDim == 2)
1034             {
1035                 workn=tmpC+tmpCI[id];
1036                 nbNeighbourNodes=tmpCI[id+1]-tmpCI[id];
1037             }
1038             else//_meshDim == 3
1039             {
1040                 workn=tmpN2+tmpNI2[id];
1041                 nbNeighbourNodes=tmpNI2[id+1]-tmpNI2[id];
1042             }
1043                         Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1044
1045                         for( int el=0;el<nbCells;el++ )
1046                                 vi.addCellId(el,workc[el]) ;
1047                         for( int el=0;el<nbNeighbourNodes;el++ )
1048                                 vi.addNeighbourNodeId(el,workn[el]) ;
1049             //Detection of border nodes
1050             bool isBorder=false;
1051                         for( int el=0;el<nbFaces;el++ )
1052             {
1053                                 vi.addFaceId(el,workf[el],_faces[workf[el]].isBorder()) ;
1054                 isBorder= isBorder || _faces[workf[el]].isBorder();
1055             }
1056             if(isBorder)
1057                 _boundaryNodeIds.push_back(id);
1058                         _nodes[id] = vi ;
1059                 }
1060
1061                 if(_spaceDim==_meshDim)
1062                         fieldn->decrRef();
1063                 fieldl->decrRef();
1064                 barySeg->decrRef();
1065         }
1066         else
1067                 throw CdmathException("Mesh::setMesh space dimension should be 1, 2 or 3");
1068
1069     //definition of the bounding box for unstructured meshes
1070     if(!_isStructured)//Structured case is treated in function readMeshMed
1071     {
1072         double Box0[2*_spaceDim];
1073         coo->getMinMaxPerComponent(Box0);
1074
1075         _xMin=Box0[0];
1076         _xMax=Box0[1];
1077         if (_spaceDim>=2)
1078         {
1079             _yMin=Box0[2];
1080             _yMax=Box0[3];
1081         }
1082         if (_spaceDim>=3)
1083         {
1084             _zMin=Box0[4];
1085             _zMax=Box0[5];
1086         }
1087     }
1088
1089     //Set boundary groups
1090         _faceGroupNames.push_back("Boundary");
1091     _nodeGroupNames.push_back("Boundary");
1092
1093     desc->decrRef();
1094         descI->decrRef();
1095         revDesc->decrRef();
1096         revDescI->decrRef();
1097         mu2->decrRef();
1098         baryCell->decrRef();
1099         fields->decrRef();
1100         revNode->decrRef();
1101         revNodeI->decrRef();
1102         revCell->decrRef();
1103         revCellI->decrRef();
1104
1105     if (_meshDim == 3)
1106     {
1107         revNode2->decrRef();
1108         revNodeI2->decrRef();
1109         desc2->decrRef();
1110         descI2->decrRef();
1111         revDesc2->decrRef();
1112         revDescI2->decrRef();
1113         mu3->decrRef();
1114     }
1115
1116     return mu;
1117 }
1118
1119 //----------------------------------------------------------------------
1120 Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
1121 //----------------------------------------------------------------------
1122 {
1123         if(nx<=0)
1124                 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0");
1125         if(xmin>=xmax)
1126                 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax");
1127
1128         double dx = (xmax - xmin)/nx ;
1129
1130         _spaceDim = 1 ;
1131         _meshDim  = 1 ;
1132     _name=meshName;
1133     _isStructured = true;
1134     _indexFacePeriodicSet=false;
1135
1136         _xMin=xmin;
1137         _xMax=xmax;
1138         _yMin=0.;
1139         _yMax=0.;
1140         _zMin=0.;
1141         _zMax=0.;
1142
1143         _dxyz.resize(_spaceDim);
1144         _dxyz[0]=dx;
1145         _nxyz.resize(_spaceDim);
1146         _nxyz[0]=nx;
1147
1148         double *originPtr = new double[_spaceDim];
1149         double *dxyzPtr = new double[_spaceDim];
1150         mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
1151
1152         originPtr[0]=xmin;
1153         nodeStrctPtr[0]=nx+1;
1154         dxyzPtr[0]=dx;
1155
1156         _mesh=MEDCouplingIMesh::New(meshName,
1157                         _spaceDim,
1158                         nodeStrctPtr,
1159                         nodeStrctPtr+_spaceDim,
1160                         originPtr,
1161                         originPtr+_spaceDim,
1162                         dxyzPtr,
1163                         dxyzPtr+_spaceDim);
1164         delete [] originPtr;
1165         delete [] dxyzPtr;
1166         delete [] nodeStrctPtr;
1167
1168         DataArrayIdType *desc=DataArrayIdType::New();
1169         DataArrayIdType *descI=DataArrayIdType::New();
1170         DataArrayIdType *revDesc=DataArrayIdType::New();
1171         DataArrayIdType *revDescI=DataArrayIdType::New();
1172         MEDCouplingUMesh* mu=_mesh->buildUnstructured();
1173         MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
1174
1175         const mcIdType *tmp=desc->getConstPointer();
1176         const mcIdType *tmpI=descI->getConstPointer();
1177
1178         const mcIdType *tmpA =revDesc->getConstPointer();
1179         const mcIdType *tmpAI=revDescI->getConstPointer();
1180
1181         _eltsTypes=mu->getAllGeoTypesSorted();
1182
1183         DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
1184         const double *coorBary=baryCell->getConstPointer();
1185
1186         _numberOfCells = _mesh->getNumberOfCells() ;
1187         _cells    = new Cell[_numberOfCells] ;
1188
1189         _numberOfNodes = mu->getNumberOfNodes() ;
1190         _nodes    = new Node[_numberOfNodes] ;
1191
1192         _numberOfFaces = _numberOfNodes;
1193         _faces    = new Face[_numberOfFaces] ;
1194
1195     _numberOfEdges = _numberOfCells;
1196
1197         MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
1198         DataArrayDouble *longueur = fieldl->getArray();
1199         const double *lon=longueur->getConstPointer();
1200
1201         DataArrayDouble *coo = mu->getCoords() ;
1202         const double *cood=coo->getConstPointer();
1203
1204         int comp=0;
1205         for( int id=0;id<_numberOfCells;id++ )
1206         {
1207                 int nbVertices=mu->getNumberOfNodesInCell(id) ;
1208                 Point p(coorBary[id],0.0,0.0) ;
1209                 Cell ci( nbVertices, nbVertices, lon[id], p ) ;
1210
1211                 std::vector<mcIdType> nodeIdsOfCell ;
1212                 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
1213                 for( int el=0;el<nbVertices;el++ )
1214                 {
1215                         ci.addNodeId(el,nodeIdsOfCell[el]) ;
1216                         ci.addFaceId(el,nodeIdsOfCell[el]) ;
1217                 }
1218
1219         double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
1220
1221         mcIdType nbFaces=tmpI[id+1]-tmpI[id];
1222         const mcIdType *work=tmp+tmpI[id];
1223
1224         for( int el=0;el<nbFaces;el++ )
1225                 {
1226                         ci.addNormalVector(el,xn,0.0,0.0) ;
1227                         ci.addFaceId(el,work[el]) ;
1228                         xn = - xn;
1229                 }
1230
1231                 _cells[id] = ci ;
1232
1233                 comp=comp+2;
1234         }
1235
1236     //Suppress the following since tmpN=tmpA
1237         DataArrayIdType *revNode=DataArrayIdType::New();
1238         DataArrayIdType *revNodeI=DataArrayIdType::New();
1239         mu->getReverseNodalConnectivity(revNode,revNodeI) ;
1240         const mcIdType *tmpN=revNode->getConstPointer();
1241         const mcIdType *tmpNI=revNodeI->getConstPointer();
1242
1243         for( int id=0;id<_numberOfNodes;id++ )
1244         {
1245                 std::vector<double> coo ;
1246                 mu->getCoordinatesOfNode(id,coo);
1247                 Point p(coo[0],0.0,0.0) ;
1248                 const mcIdType *workc=tmpN+tmpNI[id];
1249                 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
1250                 mcIdType nbFaces=1;
1251         const mcIdType *workn=tmpA+tmpAI[id];
1252         mcIdType nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
1253
1254                 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1255         for( int el=0;el<nbCells;el++ )
1256                         vi.addCellId(el,workc[el]) ;
1257                 for( int el=0;el<nbFaces;el++ )
1258                         vi.addFaceId(el,id) ;
1259         for( int el=0;el<nbNeighbourNodes;el++ )
1260                         vi.addNeighbourNodeId(el,workn[el]) ;
1261                 _nodes[id] = vi ;
1262
1263
1264                 int nbVertices=1;
1265                 Face fi( nbVertices, nbCells, 1.0, p, 1., 0., 0. ) ;
1266         for( int el=0;el<nbVertices;el++ )
1267                         fi.addNodeId(el,id) ;
1268
1269                 for( int el=0;el<nbCells;el++ )
1270                         fi.addCellId(el,workc[el]) ;
1271                 _faces[id] = fi ;
1272         }
1273
1274     //Set boundary groups
1275         _faceGroupNames.push_back("Boundary");
1276     _nodeGroupNames.push_back("Boundary");
1277     _boundaryFaceIds.push_back(0);
1278     _boundaryFaceIds.push_back(_numberOfFaces-1);
1279     _boundaryNodeIds.push_back(0);
1280     _boundaryNodeIds.push_back(_numberOfNodes-1);
1281
1282         fieldl->decrRef();
1283         baryCell->decrRef();
1284         desc->decrRef();
1285         descI->decrRef();
1286         revDesc->decrRef();
1287         revDescI->decrRef();
1288         revNode->decrRef();
1289         revNodeI->decrRef();
1290         mu2->decrRef();
1291         mu->decrRef();
1292 }
1293
1294 //----------------------------------------------------------------------
1295 Mesh::Mesh( std::vector<double> points, std::string meshName )
1296 //----------------------------------------------------------------------
1297 {
1298     int nx=points.size();
1299
1300         if(nx<2)
1301                 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : nx < 2, vector should contain at least two values");
1302     int i=0;
1303     while( i<nx-1 && points[i+1]>points[i] )
1304         i++;
1305         if( i!=nx-1 )
1306     {
1307         //cout<< points << endl;
1308                 throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : vector values should be sorted");
1309     }
1310
1311         _spaceDim = 1 ;
1312         _meshDim  = 1 ;
1313     _name=meshName;
1314         _xMin=points[0];
1315         _xMax=points[nx-1];
1316         _yMin=0.;
1317         _yMax=0.;
1318         _zMin=0.;
1319         _zMax=0.;
1320
1321     _isStructured = false;
1322     _indexFacePeriodicSet=false;
1323
1324     MEDCouplingUMesh * mesh1d = MEDCouplingUMesh::New(meshName, 1);
1325     mesh1d->allocateCells(nx - 1);
1326     double * coords = new double[nx];
1327     mcIdType * nodal_con = new mcIdType[2];
1328     coords[0]=points[0];
1329     for(int i=0; i<nx- 1 ; i++)
1330     {
1331         nodal_con[0]=i;
1332         nodal_con[1]=i+1;
1333         mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con);
1334         coords[i+1]=points[i + 1];
1335     }
1336     mesh1d->finishInsertingCells();
1337
1338     DataArrayDouble * coords_arr = DataArrayDouble::New();
1339     coords_arr->alloc(nx,1);
1340     std::copy(coords,coords+nx,coords_arr->getPointer());
1341     mesh1d->setCoords(coords_arr);
1342
1343     delete [] coords, nodal_con;
1344     coords_arr->decrRef();
1345
1346     _mesh=mesh1d;
1347
1348         DataArrayIdType *desc=DataArrayIdType::New();
1349         DataArrayIdType *descI=DataArrayIdType::New();
1350         DataArrayIdType *revDesc=DataArrayIdType::New();
1351         DataArrayIdType *revDescI=DataArrayIdType::New();
1352     MEDCouplingUMesh* mu=_mesh->buildUnstructured();
1353         MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
1354
1355         DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
1356         const double *coorBary=baryCell->getConstPointer();
1357
1358         _numberOfCells = _mesh->getNumberOfCells() ;
1359         _cells    = new Cell[_numberOfCells] ;
1360
1361     _numberOfEdges = _numberOfCells;
1362
1363         _eltsTypes=mu->getAllGeoTypesSorted();
1364
1365         MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
1366         DataArrayDouble *longueur = fieldl->getArray();
1367         const double *lon=longueur->getConstPointer();
1368
1369         int comp=0;
1370         for( int id=0;id<_numberOfCells;id++ )
1371         {
1372                 int nbVertices=mu->getNumberOfNodesInCell(id) ;
1373                 Point p(coorBary[id],0.0,0.0) ;
1374                 Cell ci( nbVertices, nbVertices, lon[id], p ) ;
1375
1376                 std::vector<mcIdType> nodeIdsOfCell ;
1377                 mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
1378                 for( int el=0;el<nbVertices;el++ )
1379                 {
1380                         ci.addNodeId(el,nodeIdsOfCell[el]) ;
1381                         ci.addFaceId(el,nodeIdsOfCell[el]) ;
1382                 }
1383                 _cells[id] = ci ;
1384                 comp=comp+2;
1385         }
1386
1387
1388     //Suppress the following since tmpN=tmpA
1389         DataArrayIdType *revNode=DataArrayIdType::New();
1390         DataArrayIdType *revNodeI=DataArrayIdType::New();
1391         mu->getReverseNodalConnectivity(revNode,revNodeI) ;
1392         const mcIdType *tmpN=revNode->getConstPointer();
1393         const mcIdType *tmpNI=revNodeI->getConstPointer();
1394
1395         _numberOfNodes = mu->getNumberOfNodes() ;
1396         _nodes    = new Node[_numberOfNodes] ;
1397         _numberOfFaces = _numberOfNodes;
1398         _faces    = new Face[_numberOfFaces] ;
1399
1400         for( int id=0;id<_numberOfNodes;id++ )
1401         {
1402                 std::vector<double> coo ;
1403                 mu->getCoordinatesOfNode(id,coo);
1404                 Point p(coo[0],0.0,0.0) ;
1405                 const mcIdType *workc=tmpN+tmpNI[id];
1406                 mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
1407                 mcIdType nbFaces=1;
1408         const mcIdType *workn=tmpN+tmpNI[id];
1409         mcIdType nbNeighbourNodes=tmpNI[id+1]-tmpNI[id];
1410                 Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
1411                 int nbVertices=1;
1412                 /* provisoire !!!!!!!!!!!!*/
1413                 //        Point pf(0.0,0.0,0.0) ;
1414                 Face fi( nbVertices, nbCells, 0.0, p, 0., 0., 0. ) ;
1415
1416                 for( int el=0;el<nbCells;el++ )
1417                         vi.addCellId(el,workc[el]) ;
1418                 for( int el=0;el<nbFaces;el++ )
1419                         vi.addFaceId(el,id) ;
1420         for( int el=0;el<nbNeighbourNodes;el++ )
1421                         vi.addNeighbourNodeId(el,workn[el]) ;
1422                 _nodes[id] = vi ;
1423
1424                 for( int el=0;el<nbVertices;el++ )
1425                         fi.addNodeId(el,id) ;
1426
1427                 for( int el=0;el<nbCells;el++ )
1428                         fi.addCellId(el,workc[el]) ;
1429                 _faces[id] = fi ;
1430         }
1431     //Set boundary groups
1432         _faceGroupNames.push_back("Boundary");
1433     _nodeGroupNames.push_back("Boundary");
1434     _boundaryFaceIds.push_back(0);
1435     _boundaryFaceIds.push_back(_numberOfFaces-1);
1436     _boundaryNodeIds.push_back(0);
1437     _boundaryNodeIds.push_back(_numberOfNodes-1);
1438
1439         fieldl->decrRef();
1440         baryCell->decrRef();
1441         desc->decrRef();
1442         descI->decrRef();
1443         revDesc->decrRef();
1444         revDescI->decrRef();
1445         revNode->decrRef();
1446         revNodeI->decrRef();
1447         mu2->decrRef();
1448         mu->decrRef();
1449 }
1450 //----------------------------------------------------------------------
1451 Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
1452 //----------------------------------------------------------------------
1453 {
1454         if(nx<=0 || ny<=0)
1455                 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : nx <= 0 or ny <= 0");
1456         if(xmin>=xmax)
1457                 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : xmin >= xmax");
1458         if(ymin>=ymax)
1459                 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax");
1460
1461         _xMin=xmin;
1462         _xMax=xmax;
1463         _yMin=ymin;
1464         _yMax=ymax;
1465         _zMin=0.;
1466         _zMax=0.;
1467
1468
1469         double dx = (xmax - xmin)/nx ;
1470         double dy = (ymax - ymin)/ny ;
1471
1472         _spaceDim = 2 ;
1473         _meshDim  = 2 ;
1474     _name=meshName;
1475     _indexFacePeriodicSet=false;
1476     _isStructured = true;
1477         _nxyz.resize(_spaceDim);
1478         _nxyz[0]=nx;
1479         _nxyz[1]=ny;
1480
1481         _dxyz.resize(_spaceDim);
1482         _dxyz[0]=dx;
1483         _dxyz[1]=dy;
1484
1485         double *originPtr = new double[_spaceDim];
1486         double *dxyzPtr   = new double[_spaceDim];
1487         mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
1488
1489         originPtr[0]=xmin;
1490         originPtr[1]=ymin;
1491         nodeStrctPtr[0]=nx+1;
1492         nodeStrctPtr[1]=ny+1;
1493         dxyzPtr[0]=dx;
1494         dxyzPtr[1]=dy;
1495
1496         _mesh=MEDCouplingIMesh::New(meshName,
1497                         _spaceDim,
1498                         nodeStrctPtr,
1499                         nodeStrctPtr+_spaceDim,
1500                         originPtr,
1501                         originPtr+_spaceDim,
1502                         dxyzPtr,
1503                         dxyzPtr+_spaceDim);
1504
1505     if(split_to_triangles_policy==0 || split_to_triangles_policy==1)
1506         {
1507             _mesh=_mesh->buildUnstructured();
1508             _mesh->simplexize(split_to_triangles_policy);
1509         }
1510     else if (split_to_triangles_policy != -1)
1511         {
1512             cout<< "split_to_triangles_policy = "<< split_to_triangles_policy << endl;
1513             throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : Unknown splitting policy");
1514         }
1515
1516         delete [] originPtr;
1517         delete [] dxyzPtr;
1518         delete [] nodeStrctPtr;
1519
1520         setMesh();
1521 }
1522
1523 //----------------------------------------------------------------------
1524 Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz, int split_to_tetrahedra_policy, std::string meshName)
1525 //----------------------------------------------------------------------
1526 {
1527         if(nx<=0 || ny<=0 || nz<=0)
1528                 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : nx <= 0 or ny <= 0 or nz <= 0");
1529         if(xmin>=xmax)
1530                 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : xmin >= xmax");
1531         if(ymin>=ymax)
1532                 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : ymin >= ymax");
1533         if(zmin>=zmax)
1534                 throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : zmin >= zmax");
1535
1536         _spaceDim = 3;
1537         _meshDim  = 3;
1538     _name=meshName;
1539     _indexFacePeriodicSet=false;
1540         _xMin=xmin;
1541         _xMax=xmax;
1542         _yMin=ymin;
1543         _yMax=ymax;
1544         _zMin=zmin;
1545         _zMax=zmax;
1546
1547         double dx = (xmax - xmin)/nx ;
1548         double dy = (ymax - ymin)/ny ;
1549         double dz = (zmax - zmin)/nz ;
1550
1551     _isStructured = true;
1552         _dxyz.resize(_spaceDim);
1553         _dxyz[0]=dx;
1554         _dxyz[1]=dy;
1555         _dxyz[2]=dz;
1556
1557         _nxyz.resize(_spaceDim);
1558         _nxyz[0]=nx;
1559         _nxyz[1]=ny;
1560         _nxyz[2]=nz;
1561
1562         double *originPtr = new double[_spaceDim];
1563         double *dxyzPtr = new double[_spaceDim];
1564         mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
1565
1566         originPtr[0]=xmin;
1567         originPtr[1]=ymin;
1568         originPtr[2]=zmin;
1569         nodeStrctPtr[0]=nx+1;
1570         nodeStrctPtr[1]=ny+1;
1571         nodeStrctPtr[2]=nz+1;
1572         dxyzPtr[0]=dx;
1573         dxyzPtr[1]=dy;
1574         dxyzPtr[2]=dz;
1575
1576         _mesh=MEDCouplingIMesh::New(meshName,
1577                         _spaceDim,
1578                         nodeStrctPtr,
1579                         nodeStrctPtr+_spaceDim,
1580                         originPtr,
1581                         originPtr+_spaceDim,
1582                         dxyzPtr,
1583                         dxyzPtr+_spaceDim);
1584
1585     if( split_to_tetrahedra_policy == 0 )
1586         {
1587             _mesh=_mesh->buildUnstructured();
1588             _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5);
1589         }
1590     else if( split_to_tetrahedra_policy == 1 )
1591         {
1592             _mesh=_mesh->buildUnstructured();
1593             _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6);
1594         }
1595     else if ( split_to_tetrahedra_policy != -1 )
1596         {
1597             cout<< "split_to_tetrahedra_policy = "<< split_to_tetrahedra_policy << endl;
1598             throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : splitting policy value should be 0 or 1");
1599         }
1600
1601         delete [] originPtr;
1602         delete [] dxyzPtr;
1603         delete [] nodeStrctPtr;
1604
1605         setMesh();
1606 }
1607
1608 //----------------------------------------------------------------------
1609 int
1610 Mesh::getSpaceDimension( void )  const
1611 //----------------------------------------------------------------------
1612 {
1613         return _spaceDim ;
1614 }
1615
1616 //----------------------------------------------------------------------
1617 int
1618 Mesh::getMeshDimension( void )  const
1619 //----------------------------------------------------------------------
1620 {
1621         return _meshDim ;
1622 }
1623
1624 std::vector<double>
1625 Mesh::getDXYZ() const
1626 {
1627     if(!_isStructured)
1628                 throw CdmathException("std::vector<double> Mesh::getDXYZ() : dx,dy and dz are defined only for structured meshes !");
1629
1630         return _dxyz;
1631 }
1632
1633 std::vector<mcIdType>
1634 Mesh::getCellGridStructure() const
1635 {
1636     if(!_isStructured)
1637                 throw CdmathException("std::vector<int> Mesh::getCellGridStructure() : nx, ny and nz are defined only for structured meshes !");
1638
1639         return _nxyz;
1640 }
1641
1642 //----------------------------------------------------------------------
1643 int
1644 Mesh::getNx( void )  const
1645 //----------------------------------------------------------------------
1646 {
1647     if(!_isStructured)
1648                 throw CdmathException("int Mesh::getNx( void ) : Nx is defined only for structured meshes !");
1649
1650         return _nxyz[0];
1651 }
1652
1653 //----------------------------------------------------------------------
1654 int
1655 Mesh::getNy( void )  const
1656 //----------------------------------------------------------------------
1657 {
1658     if(!_isStructured)
1659                 throw CdmathException("int Mesh::getNy( void ) : Ny is defined only for structured meshes !");
1660         if(_spaceDim < 2)
1661                 throw CdmathException("int double& Field::operator[ielem] : Ny is not defined in dimension < 2!");
1662
1663         return _nxyz[1];
1664 }
1665
1666 //----------------------------------------------------------------------
1667 int
1668 Mesh::getNz( void )  const
1669 //----------------------------------------------------------------------
1670 {
1671     if(!_isStructured)
1672                 throw CdmathException("int Mesh::getNz( void ) : Nz is defined only for structured meshes !");
1673         if(_spaceDim < 3)
1674                 throw CdmathException("int Mesh::getNz( void ) : Nz is not defined in dimension < 3!");
1675
1676         return _nxyz[2];
1677 }
1678
1679 //----------------------------------------------------------------------
1680 double
1681 Mesh::getXMin( void )  const
1682 //----------------------------------------------------------------------
1683 {
1684         return _xMin ;
1685 }
1686
1687 //----------------------------------------------------------------------
1688 double
1689 Mesh::getXMax( void )  const
1690 //----------------------------------------------------------------------
1691 {
1692         return _xMax ;
1693 }
1694
1695 //----------------------------------------------------------------------
1696 double
1697 Mesh::getYMin( void )  const
1698 //----------------------------------------------------------------------
1699 {
1700         return _yMin ;
1701 }
1702
1703 //----------------------------------------------------------------------
1704 double
1705 Mesh::getYMax( void )  const
1706 //----------------------------------------------------------------------
1707 {
1708         return _yMax ;
1709 }
1710
1711 //----------------------------------------------------------------------
1712 double
1713 Mesh::getZMin( void )  const
1714 //----------------------------------------------------------------------
1715 {
1716         return _zMin ;
1717 }
1718
1719 //----------------------------------------------------------------------
1720 double
1721 Mesh::getZMax( void )  const
1722 //----------------------------------------------------------------------
1723 {
1724         return _zMax ;
1725 }
1726
1727 //----------------------------------------------------------------------
1728 MCAuto<MEDCouplingMesh>
1729 Mesh::getMEDCouplingMesh( void )  const
1730 //----------------------------------------------------------------------
1731 {
1732         return _mesh ;
1733 }
1734
1735 //----------------------------------------------------------------------
1736 int
1737 Mesh::getNumberOfNodes ( void ) const
1738 //----------------------------------------------------------------------
1739 {
1740         return _numberOfNodes ;
1741 }
1742
1743 //----------------------------------------------------------------------
1744 int
1745 Mesh::getNumberOfCells ( void ) const
1746 //----------------------------------------------------------------------
1747 {
1748         return _numberOfCells ;
1749 }
1750
1751 //----------------------------------------------------------------------
1752 int
1753 Mesh::getNumberOfFaces ( void ) const
1754 //----------------------------------------------------------------------
1755 {
1756         return _numberOfFaces ;
1757 }
1758
1759 //----------------------------------------------------------------------
1760 int
1761 Mesh::getNumberOfEdges ( void ) const
1762 //----------------------------------------------------------------------
1763 {
1764         return _numberOfEdges ;
1765 }
1766
1767 //----------------------------------------------------------------------
1768 Face*
1769 Mesh::getFaces ( void )  const
1770 //----------------------------------------------------------------------
1771 {
1772         return _faces ;
1773 }
1774
1775 //----------------------------------------------------------------------
1776 Cell*
1777 Mesh::getCells ( void ) const
1778 //----------------------------------------------------------------------
1779 {
1780         return _cells ;
1781 }
1782
1783 //----------------------------------------------------------------------
1784 Cell&
1785 Mesh::getCell ( int i ) const
1786 //----------------------------------------------------------------------
1787 {
1788         return _cells[i] ;
1789 }
1790
1791 //----------------------------------------------------------------------
1792 Face&
1793 Mesh::getFace ( int i ) const
1794 //----------------------------------------------------------------------
1795 {
1796         return _faces[i] ;
1797 }
1798
1799 //----------------------------------------------------------------------
1800 Node&
1801 Mesh::getNode ( int i ) const
1802 //----------------------------------------------------------------------
1803 {
1804         return _nodes[i] ;
1805 }
1806
1807 //----------------------------------------------------------------------
1808 Node*
1809 Mesh::getNodes ( void )  const
1810 //----------------------------------------------------------------------
1811 {
1812         return _nodes ;
1813 }
1814
1815 vector<string>
1816 Mesh::getNameOfFaceGroups( void )  const
1817 {
1818         return _faceGroupNames;
1819 }
1820
1821 vector<MEDCoupling::MEDCouplingUMesh *>
1822 Mesh::getFaceGroups( void )  const
1823 {
1824         return _faceGroups;
1825 }
1826
1827 vector<string>
1828 Mesh::getNameOfNodeGroups( void )  const
1829 {
1830         return _nodeGroupNames;
1831 }
1832
1833 vector<MEDCoupling::DataArrayIdType *>
1834 Mesh::getNodeGroups( void )  const
1835 {
1836         return _nodeGroups;
1837 }
1838
1839 //----------------------------------------------------------------------
1840 const Mesh&
1841 Mesh::operator= ( const Mesh& mesh )
1842 //----------------------------------------------------------------------
1843 {
1844         _spaceDim = mesh.getSpaceDimension() ;
1845         _meshDim  = mesh.getMeshDimension() ;
1846     _name = mesh.getName();
1847         _numberOfNodes = mesh.getNumberOfNodes();
1848         _numberOfFaces = mesh.getNumberOfFaces();
1849         _numberOfCells = mesh.getNumberOfCells();
1850         _numberOfEdges = mesh.getNumberOfEdges();
1851
1852     _isStructured = mesh.isStructured();
1853     if(_isStructured)
1854     {
1855         _nxyz = mesh.getCellGridStructure() ;
1856         _dxyz=mesh.getDXYZ();
1857         _xMin=mesh.getXMin();
1858         _xMax=mesh.getXMax();
1859         _yMin=mesh.getYMin();
1860         _yMax=mesh.getYMax();
1861         _zMin=mesh.getZMin();
1862         _zMax=mesh.getZMax();
1863     }
1864         _faceGroupNames = mesh.getNameOfFaceGroups() ;
1865         _faceGroups = mesh.getFaceGroups() ;
1866         _nodeGroupNames = mesh.getNameOfNodeGroups() ;
1867         _nodeGroups = mesh.getNodeGroups() ;
1868
1869         if (_nodes)
1870         {
1871                 delete [] _nodes ;
1872                 _nodes=NULL;
1873         }
1874         if (_faces)
1875         {
1876                 delete [] _faces ;
1877                 _faces=NULL;
1878         }
1879         if (_cells)
1880         {
1881                 delete [] _cells ;
1882                 _cells=NULL;
1883         }
1884
1885         _nodes   = new Node[_numberOfNodes] ;
1886         _faces   = new Face[_numberOfFaces] ;
1887         _cells   = new Cell[_numberOfCells] ;
1888
1889         for (int i=0;i<_numberOfNodes;i++)
1890                 _nodes[i]=mesh.getNode(i);
1891
1892         for (int i=0;i<_numberOfFaces;i++)
1893                 _faces[i]=mesh.getFace(i);
1894
1895         for (int i=0;i<_numberOfCells;i++)
1896                 _cells[i]=mesh.getCell(i);
1897
1898     _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
1899     if(_indexFacePeriodicSet)
1900         _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
1901
1902     _boundaryFaceIds=mesh.getBoundaryFaceIds();
1903     _boundaryNodeIds=mesh.getBoundaryNodeIds();
1904
1905     _eltsTypes=mesh.getElementTypes();
1906     _eltsTypesNames=mesh.getElementTypesNames();
1907
1908         MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
1909         _mesh=m1;
1910         return *this;
1911 }
1912
1913 bool Mesh::isIndexFacePeriodicSet() const
1914 {
1915  return    _indexFacePeriodicSet;
1916 }
1917 //----------------------------------------------------------------------
1918 double
1919 Mesh::minRatioVolSurf() const
1920 {
1921     double dx_min  = 1e30;
1922     for(int i=0; i<_numberOfCells; i++)
1923     {
1924         Cell Ci = getCell(i);
1925         if (_meshDim > 1)
1926         {
1927             double perimeter=0;
1928             for(int k=0; k< Ci.getNumberOfFaces(); k++)
1929             {
1930                 int indexFace=Ci.getFacesId()[k];
1931                 Face Fk = getFace(indexFace);
1932                 perimeter+=Fk.getMeasure();
1933             }
1934             dx_min = min(dx_min,Ci.getMeasure()/perimeter);
1935         }
1936         else
1937             dx_min = min(dx_min,Ci.getMeasure());
1938     }
1939
1940     return dx_min;
1941 }
1942
1943 Mesh
1944 Mesh::getBoundaryMesh ( void )  const
1945 {
1946         return Mesh(_boundaryMesh);
1947 }
1948
1949 int
1950 Mesh::getMaxNbNeighbours(EntityType type) const
1951 {
1952     double result=0;
1953
1954     if (type==CELLS)
1955         {
1956         for(int i=0; i<_numberOfCells; i++)
1957             if(result < _cells[i].getNumberOfFaces())
1958                 result=_cells[i].getNumberOfFaces();
1959         }
1960     else if(type==NODES)
1961         {
1962         for(int i=0; i<_numberOfNodes; i++)
1963             if(result < _nodes[i].getNumberOfEdges())
1964                 result=_nodes[i].getNumberOfEdges();
1965         }
1966     else
1967                 throw CdmathException("Mesh::getMaxNbNeighbours : entity type is not accepted. Should be CELLS or NODES");
1968
1969     return result;
1970 }
1971 //----------------------------------------------------------------------
1972 void
1973 Mesh::writeVTK ( const std::string fileName ) const
1974 //----------------------------------------------------------------------
1975 {
1976         string fname=fileName+".vtu";
1977         _mesh->writeVTK(fname.c_str()) ;
1978 }
1979
1980 //----------------------------------------------------------------------
1981 void
1982 Mesh::writeMED ( const std::string fileName ) const
1983 //----------------------------------------------------------------------
1984 {
1985         string fname=fileName+".med";
1986         MEDCouplingUMesh* mu=_mesh->buildUnstructured();
1987         MEDCoupling::WriteUMesh(fname.c_str(),mu,true);
1988
1989         //MEDFileUMesh meshMEDFile;
1990         //meshMEDFile.setMeshAtLevel(0,mu);
1991         //for(int i=0; i< _faceGroups.size(); i++)
1992         //meshMEDFile.setMeshAtLevel(-1,_faceGroups[i]);
1993         //if (fromScratch)
1994         //MEDCoupling::meshMEDFile.write(fname.c_str(),2)       ;
1995         //else
1996         //MEDCoupling::meshMEDFile.write(fname.c_str(),1)       ;
1997
1998
1999         mu->decrRef();
2000 }
2001
2002 std::vector< int >
2003 Mesh::getGroupFaceIds(std::string groupName) const
2004 {
2005     if( std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName) == _faceGroupNames.end() )
2006     {
2007         cout<<"Mesh::getGroupFaceIds Error : face group " << groupName << " does not exist"<<endl;
2008         throw CdmathException("Required face group does not exist");
2009     }
2010     else
2011     {
2012         std::vector< int > result(0);
2013         for(int i=0; i<_boundaryFaceIds.size(); i++)
2014         {
2015             vector< std::string > groupNames = getFace(_boundaryFaceIds[i]).getGroupNames();
2016             if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() )
2017                 result.push_back(_boundaryFaceIds[i]);
2018         }
2019         return result;
2020     }
2021 }
2022
2023 std::vector< int >
2024 Mesh::getGroupNodeIds(std::string groupName) const
2025 {
2026     if( std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName) == _nodeGroupNames.end() )
2027     {
2028         cout<<"Mesh::getGroupNodeIds Error : node group " << groupName << " does not exist"<<endl;
2029         throw CdmathException("Required node group does not exist");
2030     }
2031     else
2032     {
2033         std::vector< int > result(0);
2034         for(int i=0; i<_boundaryFaceIds.size(); i++)
2035         {
2036             vector< std::string > groupNames = getNode(_boundaryFaceIds[i]).getGroupNames();
2037             if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() )
2038                 result.push_back(_boundaryFaceIds[i]);
2039         }
2040         return result;
2041     }
2042 }
2043
2044 void
2045 Mesh::setFaceGroupByIds(std::vector< int > faceIds, std::string groupName)
2046 {
2047     for(int i=0; i< faceIds.size(); i++)
2048         getFace(faceIds[i]).setGroupName(groupName);
2049 }
2050
2051 void
2052 Mesh::setNodeGroupByIds(std::vector< int > nodeIds, std::string groupName)
2053 {
2054     for(int i=0; i< nodeIds.size(); i++)
2055         getNode(nodeIds[i]).setGroupName(groupName);
2056 }
2057