Salome HOME
Update version number: 3.1.0a2
[modules/smesh.git] / 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 }