src/SMESH/SMESH_NETGEN_3D.cxx

   1 //=============================================================================
   2 // File      : SMESH_NETGEN_3D.cxx
   3 // Created   : lundi 27 Janvier 2003
   4 // Author    : Nadir BOUHAMOU (CEA)
   5 // Project   : SALOME
   6 // Copyright : CEA 2003
   7 // $Header$
   8 //=============================================================================
   9 using namespace std;
  10
  11 #include "SMESH_NETGEN_3D.hxx"
  12 #include "SMESH_MEFISTO_2D.hxx"
  13 #include "SMESH_Gen.hxx"
  14 #include "SMESH_Mesh.hxx"
  15
  16 #include "SMDS_MeshElement.hxx"
  17 #include "SMDS_MeshNode.hxx"
  18 #include "SMDS_FacePosition.hxx"
  19
  20 #include <TopExp.hxx>
  21 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  22 #include <TopTools_ListOfShape.hxx>
  23 #include <TopTools_ListIteratorOfListOfShape.hxx>
  24 #include <TColStd_ListIteratorOfListOfInteger.hxx>
  25
  26 #include <BRep_Tool.hxx>
  27 #include <Geom_Surface.hxx>
  28 #include <Geom_Curve.hxx>
  29 #include <Geom2d_Curve.hxx>
  30
  31 #include "utilities.h"
  32
  33 /*
  34   Netgen include files
  35 */
  36
  37 #include "nglib.h"
  38
  39 //=============================================================================
  40 /*!
  41  *
  42  */
  43 //=============================================================================
  44
  45 SMESH_NETGEN_3D::SMESH_NETGEN_3D(int hypId, int studyId,
  46                              SMESH_Gen* gen)
  47   : SMESH_3D_Algo(hypId, studyId, gen)
  48 {
  49   MESSAGE("SMESH_NETGEN_3D::SMESH_NETGEN_3D");
  50   _name = "NETGEN_3D";
  51 //   _shapeType = TopAbs_SOLID;
  52   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
  53 //   MESSAGE("_shapeType octal " << oct << _shapeType);
  54   _compatibleHypothesis.push_back("MaxElementVolume");
  55
  56   _maxElementVolume = 0.;
  57
  58   _hypMaxElementVolume = NULL;
  59 }
  60
  61 //=============================================================================
  62 /*!
  63  *
  64  */
  65 //=============================================================================
  66
  67 SMESH_NETGEN_3D::~SMESH_NETGEN_3D()
  68 {
  69   MESSAGE("SMESH_NETGEN_3D::~SMESH_NETGEN_3D");
  70 }
  71
  72 //=============================================================================
  73 /*!
  74  *
  75  */
  76 //=============================================================================
  77
  78 bool SMESH_NETGEN_3D::CheckHypothesis(SMESH_Mesh& aMesh,
  79                                      const TopoDS_Shape& aShape)
  80 {
  81   MESSAGE("SMESH_NETGEN_3D::CheckHypothesis");
  82
  83   _hypMaxElementVolume = NULL;
  84
  85   list<const SMESHDS_Hypothesis*>::const_iterator itl;
  86   const SMESHDS_Hypothesis* theHyp;
  87
  88   const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
  89   int nbHyp = hyps.size();
  90   if (nbHyp != 1) return false;  // only one compatible hypothesis allowed
  91
  92   itl = hyps.begin();
  93   theHyp = (*itl);
  94
  95   string hypName = theHyp->GetName();
  96   int hypId = theHyp->GetID();
  97   SCRUTE(hypName);
  98
  99   bool isOk = false;
 100
 101   if (hypName == "MaxElementVolume")
 102     {
 103       _hypMaxElementVolume = static_cast<const SMESH_MaxElementVolume*> (theHyp);
 104       ASSERT(_hypMaxElementVolume);
 105       _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
 106       isOk =true;
 107     }
 108
 109   return isOk;
 110 }
 111
 112 //=============================================================================
 113 /*!
 114  *Here we are going to use the NETGEN mesher
 115  */
 116 //=============================================================================
 117
 118 bool SMESH_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
 119                              const TopoDS_Shape& aShape)
 120   throw (SALOME_Exception)
 121 {
 122   MESSAGE("SMESH_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
 123
 124   bool isOk = false;
 125   SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
 126   SMESH_subMesh* theSubMesh = aMesh.GetSubMesh(aShape);
 127   //const Handle(SMESHDS_SubMesh)& subMeshDS = theSubMesh->GetSubMeshDS();
 128
 129   map<int, const SMDS_MeshNode*> netgenToDS;
 130
 131   MESSAGE("SMESH_NETGEN_3D::Compute Checking the mesh Faces");
 132
 133   // check if all faces were meshed by a triangle mesher (here MESFISTO_2D)
 134
 135   vector<SMESH_subMesh*> meshFaces;
 136   vector<TopoDS_Shape> shapeFaces;
 137
 138   for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
 139     {
 140       TopoDS_Shape aShapeFace = exp.Current();
 141       SMESH_subMesh* aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
 142       ASSERT (aSubMesh);
 143       int internal_size = meshFaces.size();
 144       int index = 0;
 145       for (int i = 0;i<internal_size;i++)
 146         {
 147           if (aSubMesh == meshFaces[i]) index = 1;
 148         }
 149       if (index == 0) meshFaces.push_back(aSubMesh);
 150
 151       internal_size = shapeFaces.size();
 152       index = 0;
 153       for (int i = 0;i<internal_size;i++)
 154         {
 155           if (aShapeFace == shapeFaces[i]) index = 1;
 156         }
 157       if (index == 0) shapeFaces.push_back(aShapeFace);
 158     }
 159
 160   int numberOfFaces = meshFaces.size();
 161   int numberOfShapeFaces = shapeFaces.size();
 162
 163   SCRUTE(numberOfFaces);
 164   SCRUTE(numberOfShapeFaces);
 165
 166   MESSAGE("---");
 167
 168   int NbTotOfTria = 0;
 169   int NbTotOfNodesFaces = 0;
 170
 171   for (int i=0; i<numberOfFaces; i++)
 172     {
 173       TopoDS_Shape aShapeFace = meshFaces[i]->GetSubShape();
 174       TopoDS_Shape aFace = shapeFaces[i];
 175       SMESH_Algo* algoFace = _gen->GetAlgo(aMesh, aShapeFace);
 176       string algoFaceName = algoFace->GetName();
 177       SCRUTE(algoFaceName);
 178       if (algoFaceName != "MEFISTO_2D")
 179         {
 180           SCRUTE(algoFaceName);
 181           ASSERT(0);
 182           return false;
 183         }
 184
 185       bool orientationMeshFace = (aFace.Orientation() == aShapeFace.Orientation());
 186
 187       const SMESHDS_SubMesh* aSubMeshDSFace = meshFaces[i]->GetSubMeshDS();
 188       SCRUTE(aSubMeshDSFace);
 189
 190       int nbNodes = aSubMeshDSFace->NbNodes();
 191       NbTotOfNodesFaces += nbNodes;
 192       int nbTria = aSubMeshDSFace->NbElements();
 193       NbTotOfTria += nbTria;
 194       int index = 0;
 195
 196       MESSAGE("SMESH_NETGEN_3D::Compute The mesh Face " << (i+1) << " has " << nbNodes << " face internal Nodes, " << nbTria << " triangles");
 197
 198       SCRUTE(orientationMeshFace);
 199
 200       if (orientationMeshFace)
 201         {
 202           MESSAGE("The mesh and face have the same orientation");
 203         }
 204       else
 205         {
 206           MESSAGE("The mesh and face have different orientations");
 207         }
 208
 209       SMDS_Iterator<const SMDS_MeshNode *> * iteratorNodes = aSubMeshDSFace->GetNodes();
 210       SCRUTE(nbNodes);
 211       index = 0;
 212       while(iteratorNodes->more())
 213         {
 214           index++;
 215           const SMDS_MeshNode * node = iteratorNodes->next();
 216 //        int nodeId = node->GetID();
 217 //        double nodeX = node->X();
 218 //        double nodeY = node->Y();
 219 //        double nodeZ = node->Z();
 220 //        MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
 221         }
 222       delete iteratorNodes;
 223
 224       SCRUTE(index);
 225
 226       SMDS_Iterator<const SMDS_MeshElement *> * iteratorTriangle = aSubMeshDSFace->GetElements();
 227
 228       SCRUTE(nbTria);
 229       index = 0;
 230       int numberOfDegeneratedTriangle = 0;
 231       while(iteratorTriangle->more())
 232         {
 233           index++;
 234           const SMDS_MeshElement * triangle = iteratorTriangle->next();
 235           int triangleId = triangle->GetID();
 236
 237           SMDS_Iterator<const SMDS_MeshElement *> * triangleNodesIt = triangle->nodesIterator();
 238
 239           const SMDS_MeshNode * node1 = static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 240           double node1X = node1->X();
 241           double node1Y = node1->Y();
 242           double node1Z = node1->Z();
 243
 244           const SMDS_MeshNode * node2 = static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 245           double node2X = node2->X();
 246           double node2Y = node2->Y();
 247           double node2Z = node2->Z();
 248
 249           const SMDS_MeshNode * node3 = static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
 250           double node3X = node3->X();
 251           double node3Y = node3->Y();
 252           double node3Z = node3->Z();
 253
 254           int triangleNode1 = node1->GetID();
 255           int triangleNode2 = node2->GetID();
 256           int triangleNode3 = node3->GetID();
 257
 258           // Compute the triangle surface
 259
 260           double vect1 = ((node2Y - node1Y)*(node3Z - node1Z) - (node2Z - node1Z)*(node3Y - node1Y));
 261           double vect2 = - ((node2X - node1X)*(node3Z - node1Z) - (node2Z - node1Z)*(node3X - node1X));
 262           double vect3 = ((node2X - node1X)*(node3Y - node1Y) - (node2Y - node1Y)*(node3X - node1X));
 263           double epsilon = 1.0e-6;
 264
 265           bool triangleIsDegenerated = ((abs(vect1)<epsilon) && (abs(vect2)<epsilon) && (abs(vect3)<epsilon));
 266
 267           if (triangleIsDegenerated)
 268             {
 269 //            MESSAGE("TRIANGLE -> ID = " << triangleId << " N1 = " << triangleNode1 << " N2 = " << triangleNode2 << " N3 = " << triangleNode3 << " is degenerated");
 270 //            MESSAGE("NODE -> ID = " << triangleNode1 << " X = " << node1X << " Y = " << node1Y << " Z = " << node1Z);
 271 //            MESSAGE("NODE -> ID = " << triangleNode2 << " X = " << node2X << " Y = " << node2Y << " Z = " << node2Z);
 272 //            MESSAGE("NODE -> ID = " << triangleNode3 << " X = " << node3X << " Y = " << node3Y << " Z = " << node3Z);
 273               numberOfDegeneratedTriangle++;
 274             }
 275           else
 276             {
 277 //            MESSAGE("TRIANGLE -> ID = " << triangleId << " N1 = " << triangleNode1 << " N2 = " << triangleNode2 << " N3 = " << triangleNode3 << " is normal");
 278             }
 279         }
 280
 281       delete iteratorTriangle;
 282
 283       if (numberOfDegeneratedTriangle > 0)
 284         MESSAGE("WARNING THERE IS(ARE) " << numberOfDegeneratedTriangle << " degenerated triangle on this face");
 285
 286       SCRUTE(index);
 287     }
 288
 289
 290
 291   SCRUTE(NbTotOfTria);
 292   SCRUTE(NbTotOfNodesFaces);
 293
 294   MESSAGE("SMESH_NETGEN_3D::Compute Checking the mesh Edges");
 295
 296   // check if all edges were meshed by a edge mesher (here Regular_1D)
 297
 298   vector<SMESH_subMesh*> meshEdges;
 299   for (TopExp_Explorer exp(aShape,TopAbs_EDGE);exp.More();exp.Next())
 300     {
 301       SMESH_subMesh* aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
 302       ASSERT (aSubMesh);
 303       int internal_size = meshEdges.size();
 304       int index = 0;
 305       for (int i = 0;i<internal_size;i++)
 306         {
 307           if (aSubMesh == meshEdges[i]) index = 1;
 308         }
 309       if (index == 0) meshEdges.push_back(aSubMesh);
 310     }
 311
 312   int numberOfEdges = meshEdges.size();
 313   SCRUTE(numberOfEdges);
 314
 315   MESSAGE("---");
 316
 317   int NbTotOfNodesEdges = 0;
 318   int NbTotOfSegs = 0;
 319
 320   for (int i=0; i<numberOfEdges; i++)
 321     {
 322       TopoDS_Shape aShapeEdge = meshEdges[i]->GetSubShape();
 323       SMESH_Algo* algoEdge = _gen->GetAlgo(aMesh, aShapeEdge);
 324       string algoEdgeName = algoEdge->GetName();
 325       SCRUTE(algoEdgeName);
 326       if (algoEdgeName != "Regular_1D")
 327         {
 328           SCRUTE(algoEdgeName);
 329           ASSERT(0);
 330           return false;
 331         }
 332
 333       const SMESHDS_SubMesh* aSubMeshDSEdge = meshEdges[i]->GetSubMeshDS();
 334       SCRUTE(aSubMeshDSEdge);
 335
 336       int nbNodes = aSubMeshDSEdge->NbNodes();
 337       NbTotOfNodesEdges += nbNodes;
 338       int nbSegs = aSubMeshDSEdge->NbElements();
 339       NbTotOfSegs += nbSegs;
 340
 341       MESSAGE("SMESH_NETGEN_3D::Compute The mesh Edge " << (i+1) << " has " << nbNodes << " edge internal Nodes, " << nbSegs << " segments");
 342
 343       SMDS_Iterator<const SMDS_MeshNode *> * iteratorNodes = aSubMeshDSEdge->GetNodes();
 344       SCRUTE(nbNodes);
 345       int index = 0;
 346       while(iteratorNodes->more())
 347         {
 348           index++;
 349           const SMDS_MeshNode * node = iteratorNodes->next();
 350 //        int nodeId = node->GetID();
 351 //        double nodeX = node->X();
 352 //        double nodeY = node->Y();
 353 //        double nodeZ = node->Z();
 354 //        MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
 355         }
 356       delete iteratorNodes;
 357
 358       SCRUTE(index);
 359     }
 360
 361   SCRUTE(NbTotOfNodesEdges);
 362   SCRUTE(NbTotOfSegs);
 363
 364   MESSAGE("SMESH_NETGEN_3D::Compute Checking the mesh Vertices");
 365
 366   vector<SMESH_subMesh*> meshVertices;
 367   for (TopExp_Explorer exp(aShape,TopAbs_VERTEX);exp.More();exp.Next())
 368     {
 369       SMESH_subMesh* aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
 370       ASSERT (aSubMesh);
 371       int internal_size = meshVertices.size();
 372       int index = 0;
 373       for (int i = 0;i<internal_size;i++)
 374         {
 375           if (aSubMesh == meshVertices[i]) index = 1;
 376         }
 377       if (index == 0) meshVertices.push_back(aSubMesh);
 378     }
 379
 380   int numberOfVertices = meshVertices.size();
 381   SCRUTE(numberOfVertices);
 382
 383   MESSAGE("---");
 384
 385   int NbTotOfNodesVertices = 0;
 386
 387   for (int i=0; i<numberOfVertices; i++)
 388     {
 389       TopoDS_Shape aShapeVertex = meshVertices[i]->GetSubShape();
 390
 391       const SMESHDS_SubMesh * aSubMeshDSVertex = meshVertices[i]->GetSubMeshDS();
 392       SCRUTE(aSubMeshDSVertex);
 393
 394       int nbNodes = aSubMeshDSVertex->NbNodes();
 395       NbTotOfNodesVertices += nbNodes;
 396
 397       MESSAGE("SMESH_NETGEN_3D::Compute The mesh Vertex " << (i+1) << " has " << nbNodes << " Nodes");
 398
 399       SMDS_Iterator<const SMDS_MeshNode *> * iteratorNodes = aSubMeshDSVertex->GetNodes();
 400       SCRUTE(nbNodes);
 401       int index = 0;
 402       while(iteratorNodes->more())
 403         {
 404           index++;
 405           const SMDS_MeshNode * node = iteratorNodes->next();
 406 //        int nodeId = node->GetID();
 407 //        double nodeX = node->X();
 408 //        double nodeY = node->Y();
 409 //        double nodeZ = node->Z();
 410 //        MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
 411         }
 412       delete iteratorNodes;
 413
 414       SCRUTE(index);
 415     }
 416
 417   SCRUTE(NbTotOfNodesVertices);
 418
 419   MESSAGE("SMESH_NETGEN_3D::Compute --> Analysis of all shell mesh");
 420
 421   vector<SMESH_subMesh*> meshShells;
 422   TopoDS_Shell aShell;
 423
 424   for (TopExp_Explorer exp(aShape,TopAbs_SHELL);exp.More();exp.Next())
 425     {
 426       SMESH_subMesh* aSubMesh = aMesh.GetSubMesh(exp.Current());
 427       ASSERT(aSubMesh);
 428       SCRUTE(aSubMesh);
 429       aShell = TopoDS::Shell(exp.Current());
 430       meshShells.push_back(aSubMesh);
 431     }
 432
 433   int numberOfShells = meshShells.size();
 434   SCRUTE(numberOfShells);
 435
 436   if (numberOfShells == 1)
 437     {
 438       MESSAGE("SMESH_NETGEN_3D::Compute Only one shell --> generation of the mesh using directly Netgen");
 439
 440       /*
 441         Prepare the Netgen surface mesh from the SMESHDS
 442       */
 443
 444       MESSAGE("SMESH_NETGEN_3D::Compute Prepare the Netgen surface mesh from the SMESHDS");
 445
 446       int spaceDimension = 3;
 447       int nbNodesByTri = 3;
 448       int nbNodesByTetra = 4;
 449
 450       int Netgen_NbOfNodes = NbTotOfNodesFaces +
 451                              NbTotOfNodesEdges +
 452                              NbTotOfNodesVertices;
 453       int Netgen_NbOfTria = NbTotOfTria;
 454       int Netgen_param2ndOrder = 0;
 455       double Netgen_paramFine = 1.;
 456       double Netgen_paramSize = _maxElementVolume;
 457
 458       SCRUTE(Netgen_NbOfNodes);
 459       SCRUTE(Netgen_NbOfTria);
 460
 461       double * Netgen_Coordinates = new double [spaceDimension*
 462                                                 Netgen_NbOfNodes];
 463       int * listNodeCoresNetgenSmesh = new int [Netgen_NbOfNodes];
 464       int * Netgen_Connectivity = new int [nbNodesByTri*Netgen_NbOfTria];
 465       double * Netgen_point = new double [spaceDimension];
 466       int * Netgen_triangle = new int [nbNodesByTri];
 467       int * Netgen_tetrahedron = new int [nbNodesByTetra];
 468
 469       for (int i=0; i<Netgen_NbOfTria; i++)
 470         {
 471           for (int j=0; j<nbNodesByTri; j++)
 472             Netgen_Connectivity[i*nbNodesByTri+j] = 0;
 473         }
 474
 475       double bigNumber = 1.e20;
 476
 477       for (int i=0; i<Netgen_NbOfNodes; i++)
 478         {
 479           listNodeCoresNetgenSmesh[i] = 0;
 480           for (int j=0; j<spaceDimension; j++)
 481             Netgen_Coordinates[i*spaceDimension+j] = bigNumber;
 482         }
 483
 484       int indexNodes = 0;
 485       for (int i=0; i<numberOfVertices; i++)
 486         {
 487           const SMESHDS_SubMesh * aSubMeshDSVertex =
 488             meshVertices[i]->GetSubMeshDS();
 489
 490           SMDS_Iterator<const SMDS_MeshNode *> * iteratorNodes = aSubMeshDSVertex->GetNodes();
 491
 492           while(iteratorNodes->more())
 493             {
 494               const SMDS_MeshNode * node = iteratorNodes->next();
 495               int nodeId = node->GetID();
 496               double nodeX = node->X();
 497               double nodeY = node->Y();
 498               double nodeZ = node->Z();
 499 //            MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
 500               listNodeCoresNetgenSmesh[indexNodes] = nodeId;
 501               int index = indexNodes*spaceDimension;
 502               Netgen_Coordinates[index] = nodeX;
 503               Netgen_Coordinates[index+1] = nodeY;
 504               Netgen_Coordinates[index+2] = nodeZ;
 505               netgenToDS[indexNodes] = node;
 506               indexNodes++;
 507             }
 508           delete iteratorNodes;
 509         }
 510
 511       for (int i=0; i<numberOfEdges; i++)
 512         {
 513           const SMESHDS_SubMesh *  aSubMeshDSEdge =
 514             meshEdges[i]->GetSubMeshDS();
 515
 516           SMDS_Iterator<const SMDS_MeshNode *> * iteratorNodes = aSubMeshDSEdge->GetNodes();
 517
 518           while(iteratorNodes->more())
 519             {
 520               const SMDS_MeshNode * node = iteratorNodes->next();
 521               int nodeId = node->GetID();
 522               double nodeX = node->X();
 523               double nodeY = node->Y();
 524               double nodeZ = node->Z();
 525 //            MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
 526               listNodeCoresNetgenSmesh[indexNodes] = node->GetID();
 527               int index = indexNodes*spaceDimension;
 528               Netgen_Coordinates[index] = node->X();
 529               Netgen_Coordinates[index+1] = node->Y();
 530               Netgen_Coordinates[index+2] = node->Z();
 531               netgenToDS[indexNodes] = node;
 532               indexNodes++;
 533             }
 534           delete iteratorNodes;
 535         }
 536
 537       for (int i=0; i<numberOfFaces; i++)
 538         {
 539           const SMESHDS_SubMesh * aSubMeshDSFace =
 540             meshFaces[i]->GetSubMeshDS();
 541
 542           SMDS_Iterator<const SMDS_MeshNode *> * iteratorNodes = aSubMeshDSFace->GetNodes();
 543
 544           while(iteratorNodes->more())
 545             {
 546               const SMDS_MeshNode * node = iteratorNodes->next();
 547               int nodeId = node->GetID();
 548               double nodeX = node->X();
 549               double nodeY = node->Y();
 550               double nodeZ = node->Z();
 551 //            MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
 552               listNodeCoresNetgenSmesh[indexNodes] = nodeId;
 553               int index = indexNodes*spaceDimension;
 554               Netgen_Coordinates[index] = nodeX;
 555               Netgen_Coordinates[index+1] = nodeY;
 556               Netgen_Coordinates[index+2] = nodeZ;
 557               netgenToDS[indexNodes] = node;
 558               indexNodes++;
 559             }
 560           delete iteratorNodes;
 561         }
 562
 563       SCRUTE(indexNodes);
 564
 565       for (int i=0; i<Netgen_NbOfNodes; i++)
 566         {
 567           ASSERT(listNodeCoresNetgenSmesh[i] != 0);
 568
 569           for (int j=0; j<Netgen_NbOfNodes && j!=i; j++)
 570             ASSERT(listNodeCoresNetgenSmesh[i] != listNodeCoresNetgenSmesh[j]);
 571
 572           for (int j=0; j<spaceDimension; j++)
 573             ASSERT(Netgen_Coordinates[i*spaceDimension+j] != bigNumber);
 574         }
 575
 576       int indexTrias = 0;
 577       for (int i=0; i<numberOfFaces; i++)
 578         {
 579           const SMESHDS_SubMesh * aSubMeshDSFace =
 580             meshFaces[i]->GetSubMeshDS();
 581
 582           TopoDS_Shape aFace = shapeFaces[i];
 583
 584           SMDS_Iterator<const SMDS_MeshElement *> * iteratorTriangle = aSubMeshDSFace->GetElements();
 585
 586           TopoDS_Shape aShapeFace = meshFaces[i]->GetSubShape();
 587
 588           bool orientationMeshFace = (aFace.Orientation() == aShapeFace.Orientation());
 589
 590           SCRUTE(orientationMeshFace);
 591
 592           if (orientationMeshFace)
 593             {
 594               MESSAGE("The mesh and face have the same orientation");
 595
 596               while(iteratorTriangle->more())
 597                 {
 598                   const SMDS_MeshElement * triangle = iteratorTriangle->next();
 599                   int triangleId = triangle->GetID();
 600
 601                   SMDS_Iterator<const SMDS_MeshElement *> * triangleNodesIt = triangle->nodesIterator();
 602
 603                   int triangleNode1 = (triangleNodesIt->next())->GetID();
 604                   int triangleNode2 = (triangleNodesIt->next())->GetID();
 605                   int triangleNode3 = (triangleNodesIt->next())->GetID();
 606
 607 //                MESSAGE("TRIANGLE -> ID = " << triangleId << " N1 = " << triangleNode1 << " N2 = " << triangleNode2 << " N3 = " << triangleNode3);
 608
 609                   int N1New = 0;
 610                   int N2New = 0;
 611                   int N3New = 0;
 612                   int index = indexTrias*nbNodesByTri;
 613
 614                   for (int j=0; j<Netgen_NbOfNodes; j++)
 615                     {
 616                       int jp1 = j+1;
 617
 618                       if (triangleNode1 == listNodeCoresNetgenSmesh[j])
 619                         N1New = jp1;
 620                       else if (triangleNode2 == listNodeCoresNetgenSmesh[j])
 621                         N2New = jp1;
 622                       else if (triangleNode3 == listNodeCoresNetgenSmesh[j])
 623                         N3New = jp1;
 624                     }
 625
 626                   triangleNode1 = N1New;
 627                   triangleNode2 = N2New;
 628                   triangleNode3 = N3New;
 629
 630                   Netgen_Connectivity[index] = triangleNode1;
 631                   Netgen_Connectivity[index+1] = triangleNode2;
 632                   Netgen_Connectivity[index+2] = triangleNode3;
 633
 634                   indexTrias++;
 635                 }
 636               delete iteratorTriangle;
 637             }
 638           else
 639             {
 640               MESSAGE("The mesh and face have different orientations");
 641
 642               while(iteratorTriangle->more())
 643                 {
 644                   const SMDS_MeshElement * triangle = iteratorTriangle->next();
 645                   int triangleId = triangle->GetID();
 646
 647                   SMDS_Iterator<const SMDS_MeshElement *> * triangleNodesIt = triangle->nodesIterator();
 648
 649                   int triangleNode1 = (triangleNodesIt->next())->GetID();
 650                   int triangleNode3 = (triangleNodesIt->next())->GetID();
 651                   int triangleNode2 = (triangleNodesIt->next())->GetID();
 652
 653 //                MESSAGE("TRIANGLE -> ID = " << triangleId << " N1 = " << triangleNode1 << " N2 = " << triangleNode2 << " N3 = " << triangleNode3);
 654
 655                   int N1New = 0;
 656                   int N2New = 0;
 657                   int N3New = 0;
 658                   int index = indexTrias*nbNodesByTri;
 659
 660                   for (int j=0; j<Netgen_NbOfNodes; j++)
 661                     {
 662                       int jp1 = j+1;
 663
 664                       if (triangleNode1 == listNodeCoresNetgenSmesh[j])
 665                         N1New = jp1;
 666                       else if (triangleNode2 == listNodeCoresNetgenSmesh[j])
 667                         N2New = jp1;
 668                       else if (triangleNode3 == listNodeCoresNetgenSmesh[j])
 669                         N3New = jp1;
 670                     }
 671
 672                   triangleNode1 = N1New;
 673                   triangleNode2 = N2New;
 674                   triangleNode3 = N3New;
 675
 676                   Netgen_Connectivity[index] = triangleNode1;
 677                   Netgen_Connectivity[index+1] = triangleNode2;
 678                   Netgen_Connectivity[index+2] = triangleNode3;
 679
 680                   indexTrias++;
 681                 }
 682               delete iteratorTriangle;
 683             }
 684         }
 685
 686       SCRUTE(indexTrias);
 687
 688       int * nodesUsed = new int[Netgen_NbOfNodes];
 689
 690       for (int i=0; i<Netgen_NbOfNodes; i++) nodesUsed[i] = 0;
 691
 692       for (int i=0; i<Netgen_NbOfTria; i++)
 693         for (int j=0; j<nbNodesByTri; j++)
 694           {
 695             int Nij = Netgen_Connectivity[i*nbNodesByTri+j];
 696
 697             ASSERT((Nij>=1) && (Nij<=Netgen_NbOfNodes));
 698
 699             nodesUsed[Nij-1] = 1;
 700             Netgen_Connectivity[i*nbNodesByTri+j] = Nij;
 701           }
 702
 703       for (int i=0; i<Netgen_NbOfNodes; i++)
 704         {
 705           ASSERT(nodesUsed[i] != 0);
 706         }
 707
 708       delete [] nodesUsed;
 709
 710       /*
 711         Feed the Netgen surface mesh
 712       */
 713
 714       MESSAGE("SMESH_NETGEN_3D::Compute Feed the Netgen surface mesh");
 715
 716       Ng_Mesh * Netgen_mesh;
 717
 718       Ng_Init();
 719
 720       Netgen_mesh = Ng_NewMesh();
 721
 722       Ng_Meshing_Parameters Netgen_param;
 723
 724       for (int i=0; i<Netgen_NbOfNodes; i++)
 725         {
 726           for (int j=0; j<spaceDimension; j++)
 727             Netgen_point[j] = Netgen_Coordinates[i*spaceDimension+j];
 728
 729           Ng_AddPoint(Netgen_mesh, Netgen_point);
 730         }
 731
 732       for (int i=0; i<Netgen_NbOfTria; i++)
 733         {
 734           for (int j=0; j<nbNodesByTri; j++)
 735             Netgen_triangle[j] = Netgen_Connectivity[i*nbNodesByTri+j];
 736
 737           Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
 738         }
 739
 740       SCRUTE(Netgen_paramSize);
 741
 742       Netgen_param.secondorder = Netgen_param2ndOrder;
 743       Netgen_param.fineness = Netgen_paramFine;
 744       Netgen_param.maxh = Netgen_paramSize;
 745
 746       /*
 747         Generate the volume mesh
 748       */
 749
 750       MESSAGE("SMESH_NETGEN_3D::Compute Generate the volume mesh");
 751
 752       SCRUTE(Netgen_NbOfNodes);
 753       SCRUTE(Netgen_NbOfTria);
 754
 755       SCRUTE(Ng_GetNP(Netgen_mesh));
 756       SCRUTE(Ng_GetNE(Netgen_mesh));
 757       SCRUTE(Ng_GetNSE(Netgen_mesh));
 758
 759       ASSERT(Netgen_NbOfNodes == Ng_GetNP(Netgen_mesh));
 760       ASSERT(Ng_GetNE(Netgen_mesh) == 0);
 761       ASSERT(Netgen_NbOfTria == Ng_GetNSE(Netgen_mesh));
 762
 763       Ng_Result status;
 764
 765       status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
 766
 767       SCRUTE(status);
 768
 769       int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
 770
 771       int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
 772
 773       SCRUTE(Netgen_NbOfNodesNew);
 774
 775       SCRUTE(Netgen_NbOfTetra);
 776
 777       if ((status != NG_OK) ||
 778           (Netgen_NbOfNodesNew <= Netgen_NbOfNodes) ||
 779           (Netgen_NbOfTetra <= 0))
 780         {
 781           MESSAGE("SMESH_NETGEN_3D::Compute The Volume Mesh Generation has failed ...");
 782           SCRUTE(status);
 783
 784           /*
 785             Free the memory needed by to generate the Netgen Mesh
 786           */
 787
 788           MESSAGE("SMESH_NETGEN_3D::Compute Free the memory needed by to generate the Netgen Mesh");
 789
 790           delete [] Netgen_Coordinates;
 791           delete [] Netgen_Connectivity;
 792           delete [] Netgen_point;
 793           delete [] Netgen_triangle;
 794           delete [] Netgen_tetrahedron;
 795
 796           delete [] listNodeCoresNetgenSmesh;
 797
 798           Ng_DeleteMesh(Netgen_mesh);
 799           Ng_Exit();
 800
 801           return false;
 802         }
 803
 804       MESSAGE("SMESH_NETGEN_3D::Compute End of Volume Mesh Generation");
 805       SCRUTE(status);
 806
 807       double * Netgen_CoordinatesNew = new double [spaceDimension*Netgen_NbOfNodesNew];
 808       int * Netgen_ConnectivityNew = new int [nbNodesByTetra*Netgen_NbOfTetra];
 809
 810       for (int i=0; i<Netgen_NbOfNodesNew; i++)
 811         {
 812           Ng_GetPoint(Netgen_mesh, (i+1), Netgen_point);
 813
 814           for (int j=0; j<spaceDimension; j++)
 815             Netgen_CoordinatesNew[i*spaceDimension+j] = Netgen_point[j];
 816         }
 817
 818       for (int i=0; i<Netgen_NbOfNodes; i++)
 819         for (int j=0; j<spaceDimension; j++)
 820           ASSERT(Netgen_CoordinatesNew[i*spaceDimension+j] == Netgen_Coordinates[i*spaceDimension+j])
 821
 822       for (int i=0; i<Netgen_NbOfTetra; i++)
 823         {
 824           Ng_GetVolumeElement(Netgen_mesh, (i+1), Netgen_tetrahedron);
 825
 826           for (int j=0; j<nbNodesByTetra; j++)
 827             Netgen_ConnectivityNew[i*nbNodesByTetra+j] = Netgen_tetrahedron[j];
 828         }
 829
 830       /*
 831         Feed back the SMESHDS with the generated Nodes and Volume Elements
 832       */
 833
 834       MESSAGE("SMESH_NETGEN_3D::Compute Feed back the SMESHDS with the generated Nodes and Volume Elements");
 835
 836       int NbTotOfNodesShell = Netgen_NbOfNodesNew - Netgen_NbOfNodes;
 837
 838       SCRUTE(NbTotOfNodesShell);
 839
 840       int * listNodeShellCoresNetgenSmesh = new int [NbTotOfNodesShell];
 841
 842       for (int i=0; i<NbTotOfNodesShell; i++)
 843         listNodeShellCoresNetgenSmesh[i] = 0;
 844
 845       MESSAGE("SMESH_NETGEN_3D::Compute --> Adding the New Nodes to SMESHDS");
 846
 847       for (int i=0; i<NbTotOfNodesShell; i++)
 848         {
 849           int index = (i+Netgen_NbOfNodes)*spaceDimension;
 850
 851           SMDS_MeshNode * node =
 852             meshDS->AddNode(Netgen_CoordinatesNew[index],
 853                             Netgen_CoordinatesNew[index+1],
 854                             Netgen_CoordinatesNew[index+2]);
 855
 856           meshDS->SetNodeInVolume(node, aShell);
 857
 858           index = i+Netgen_NbOfNodes;
 859           netgenToDS[index] = node;
 860
 861           listNodeShellCoresNetgenSmesh[i] = node->GetID();
 862         }
 863
 864       SCRUTE(Netgen_NbOfNodesNew);
 865
 866       SCRUTE(netgenToDS.size());
 867
 868       for (int i=0; i<NbTotOfNodesShell; i++)
 869         {
 870           ASSERT(listNodeShellCoresNetgenSmesh[i] != 0);
 871
 872           for (int j=0; j<NbTotOfNodesShell && j!=i; j++)
 873             ASSERT(listNodeShellCoresNetgenSmesh[i] != listNodeShellCoresNetgenSmesh[j]);
 874         }
 875
 876       MESSAGE("SMESH_NETGEN_3D::Compute --> Adding the New elements (Tetrahedrons) to the SMESHDS");
 877
 878       for (int i=0; i<Netgen_NbOfTetra; i++)
 879         {
 880           int index = i*nbNodesByTetra;
 881           int tetraNode1 = Netgen_ConnectivityNew[index];
 882           int tetraNode2 = Netgen_ConnectivityNew[index+1];
 883           int tetraNode3 = Netgen_ConnectivityNew[index+2];
 884           int tetraNode4 = Netgen_ConnectivityNew[index+3];
 885
 886           const SMDS_MeshNode * node1 = netgenToDS[tetraNode1-1];
 887           const SMDS_MeshNode * node2 = netgenToDS[tetraNode2-1];
 888           const SMDS_MeshNode * node3 = netgenToDS[tetraNode3-1];
 889           const SMDS_MeshNode * node4 = netgenToDS[tetraNode4-1];
 890
 891           index = tetraNode1;
 892           if (index <= Netgen_NbOfNodes)
 893             tetraNode1 = listNodeCoresNetgenSmesh[index-1];
 894           else
 895             tetraNode1 = listNodeShellCoresNetgenSmesh[index-Netgen_NbOfNodes-1];
 896
 897           index = tetraNode2;
 898           if (index <= Netgen_NbOfNodes)
 899             tetraNode2 = listNodeCoresNetgenSmesh[index-1];
 900           else
 901             tetraNode2 = listNodeShellCoresNetgenSmesh[index-Netgen_NbOfNodes-1];
 902
 903           index = tetraNode3;
 904           if (index <= Netgen_NbOfNodes)
 905             tetraNode3 = listNodeCoresNetgenSmesh[index-1];
 906           else
 907             tetraNode3 = listNodeShellCoresNetgenSmesh[index-Netgen_NbOfNodes-1];
 908
 909           index = tetraNode4;
 910           if (index <= Netgen_NbOfNodes)
 911             tetraNode4 = listNodeCoresNetgenSmesh[index-1];
 912           else
 913             tetraNode4 = listNodeShellCoresNetgenSmesh[index-Netgen_NbOfNodes-1];
 914
 915           SMDS_MeshVolume * elt =
 916             meshDS->AddVolume(node1,node2,node3,node4);
 917
 918           meshDS->SetMeshElementOnShape(elt, aShell);
 919         }
 920
 921       /*
 922         Free the memory needed by to generate the Netgen Mesh
 923       */
 924
 925       MESSAGE("SMESH_NETGEN_3D::Compute Free the memory needed by to generate the Netgen Mesh");
 926
 927       delete [] Netgen_Coordinates;
 928       delete [] Netgen_Connectivity;
 929       delete [] Netgen_CoordinatesNew;
 930       delete [] Netgen_ConnectivityNew;
 931       delete [] Netgen_point;
 932       delete [] Netgen_triangle;
 933       delete [] Netgen_tetrahedron;
 934
 935       delete [] listNodeCoresNetgenSmesh;
 936       delete [] listNodeShellCoresNetgenSmesh;
 937
 938       Ng_DeleteMesh(Netgen_mesh);
 939       Ng_Exit();
 940
 941       /*
 942         Verification
 943       */
 944
 945       {
 946         MESSAGE("SMESH_NETGEN_3D::Compute Verification of the Shell mesh");
 947
 948         TopoDS_Shape aShapeShell = meshShells[0]->GetSubShape();
 949         SMESH_Algo* algoShell = _gen->GetAlgo(aMesh, aShapeShell);
 950         string algoShellName = algoShell->GetName();
 951         SCRUTE(algoShellName);
 952         if (algoShellName != "NETGEN_3D")
 953           {
 954             SCRUTE(algoShellName);
 955             ASSERT(0);
 956             return false;
 957           }
 958
 959         const SMESHDS_SubMesh * aSubMeshDSShell = meshShells[0]->GetSubMeshDS();
 960         SCRUTE(&aSubMeshDSShell);
 961
 962         int nbNodes = aSubMeshDSShell->NbNodes();
 963         int nbTetra = aSubMeshDSShell->NbElements();
 964
 965         MESSAGE("SMESH_NETGEN_3D::Compute The mesh Shell has " << nbNodes << " shell internal Nodes, " << nbTetra << " tetrahedrons");
 966
 967         SMDS_Iterator<const SMDS_MeshNode *> * iteratorNodes = aSubMeshDSShell->GetNodes();
 968
 969         SCRUTE(nbNodes);
 970
 971         int index;
 972
 973         index = 0;
 974
 975         while(iteratorNodes->more())
 976           {
 977             index++;
 978             const SMDS_MeshNode * node = iteratorNodes->next();
 979             int nodeId = node->GetID();
 980             double nodeX = node->X();
 981             double nodeY = node->Y();
 982             double nodeZ = node->Z();
 983 //          MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
 984           }
 985         delete iteratorNodes;
 986
 987         SCRUTE(index);
 988
 989         SMDS_Iterator<const SMDS_MeshElement *> * iteratorTetra = aSubMeshDSShell->GetElements();
 990
 991         SCRUTE(nbTetra);
 992
 993         index = 0;
 994         while(iteratorTetra->more())
 995           {
 996             index++;
 997             const SMDS_MeshElement * tetra = iteratorTetra->next();
 998             int tetraId = tetra->GetID();
 999
1000             SMDS_Iterator<const SMDS_MeshElement *> * tetraNodesIt = tetra->nodesIterator();
1001
1002             int tetraNode1 = (tetraNodesIt->next())->GetID();
1003             int tetraNode2 = (tetraNodesIt->next())->GetID();
1004             int tetraNode3 = (tetraNodesIt->next())->GetID();
1005             int tetraNode4 = (tetraNodesIt->next())->GetID();
1006
1007 //          MESSAGE("TETRAHEDRON -> ID = " << tetraId << " N1 = " << tetraNode1 << " N2 = " << tetraNode2 << " N3 = " << tetraNode3 << " N4 = " << tetraNode4);
1008
1009           }
1010         delete iteratorTetra;
1011
1012         SCRUTE(index);
1013       }
1014     }
1015   else
1016     {
1017       SCRUTE(numberOfShells);
1018       MESSAGE("SMESH_NETGEN_3D::Compute ERROR More than one shell ????? ");
1019       return false;
1020     }
1021
1022   return true;
1023 }
1024
1025
1026 //=============================================================================
1027 /*!
1028  *
1029  */
1030 //=============================================================================
1031
1032 ostream & SMESH_NETGEN_3D::SaveTo(ostream & save)
1033 {
1034   return save << this;
1035 }
1036
1037 //=============================================================================
1038 /*!
1039  *
1040  */
1041 //=============================================================================
1042
1043 istream & SMESH_NETGEN_3D::LoadFrom(istream & load)
1044 {
1045   return load >> (*this);
1046 }
1047
1048 //=============================================================================
1049 /*!
1050  *
1051  */
1052 //=============================================================================
1053
1054 ostream & operator << (ostream & save, SMESH_NETGEN_3D & hyp)
1055 {
1056   return save;
1057 }
1058
1059 //=============================================================================
1060 /*!
1061  *
1062  */
1063 //=============================================================================
1064
1065 istream & operator >> (istream & load, SMESH_NETGEN_3D & hyp)
1066 {
1067   return load;
1068 }