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