CDMATH/mesh/src/Mesh.cxx

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