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