src/NETGENPlugin/NETGENPlugin_NETGEN_3D.cxx

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