CDMATH/mesh/src/Mesh.cxx

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