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)
7 // Copyright : CEA 2003
9 //=============================================================================
12 #include "NETGENPlugin_NETGEN_3D.hxx"
13 #include "SMESH_Gen.hxx"
14 #include "SMESH_Mesh.hxx"
15 #include "SMESH_subMesh.hxx"
17 #include "SMDS_MeshElement.hxx"
18 #include "SMDS_MeshNode.hxx"
19 #include "SMDS_FacePosition.hxx"
22 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
23 #include <TopTools_ListOfShape.hxx>
24 #include <TopTools_ListIteratorOfListOfShape.hxx>
25 #include <TColStd_ListIteratorOfListOfInteger.hxx>
27 #include <BRep_Tool.hxx>
28 #include <Geom_Surface.hxx>
29 #include <Geom_Curve.hxx>
30 #include <Geom2d_Curve.hxx>
32 #include "utilities.h"
40 //=============================================================================
44 //=============================================================================
46 NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D(int hypId, int studyId,
48 : SMESH_3D_Algo(hypId, studyId, gen)
50 MESSAGE("NETGENPlugin_NETGEN_3D::NETGENPlugin_NETGEN_3D");
52 // _shapeType = TopAbs_SOLID;
53 _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);// 1 bit /shape type
54 // MESSAGE("_shapeType octal " << oct << _shapeType);
55 _compatibleHypothesis.push_back("MaxElementVolume");
57 _maxElementVolume = 0.;
59 _hypMaxElementVolume = NULL;
62 //=============================================================================
66 //=============================================================================
68 NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D()
70 MESSAGE("NETGENPlugin_NETGEN_3D::~NETGENPlugin_NETGEN_3D");
73 //=============================================================================
77 //=============================================================================
79 bool NETGENPlugin_NETGEN_3D::CheckHypothesis
81 const TopoDS_Shape& aShape,
82 SMESH_Hypothesis::Hypothesis_Status& aStatus)
84 MESSAGE("NETGENPlugin_NETGEN_3D::CheckHypothesis");
86 _hypMaxElementVolume = NULL;
88 list<const SMESHDS_Hypothesis*>::const_iterator itl;
89 const SMESHDS_Hypothesis* theHyp;
91 const list<const SMESHDS_Hypothesis*>& hyps = GetUsedHypothesis(aMesh, aShape);
92 int nbHyp = hyps.size();
95 aStatus = SMESH_Hypothesis::HYP_MISSING;
96 return false; // can't work with no hypothesis
100 theHyp = (*itl); // use only the first hypothesis
102 string hypName = theHyp->GetName();
103 // int hypId = theHyp->GetID();
108 if (hypName == "MaxElementVolume")
110 _hypMaxElementVolume = static_cast<const StdMeshers_MaxElementVolume*> (theHyp);
111 ASSERT(_hypMaxElementVolume);
112 _maxElementVolume = _hypMaxElementVolume->GetMaxVolume();
114 aStatus = SMESH_Hypothesis::HYP_OK;
117 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
122 //=============================================================================
124 *Here we are going to use the NETGEN mesher
126 //=============================================================================
128 bool NETGENPlugin_NETGEN_3D::Compute(SMESH_Mesh& aMesh,
129 const TopoDS_Shape& aShape)
131 MESSAGE("NETGENPlugin_NETGEN_3D::Compute with maxElmentsize = " << _maxElementVolume);
134 SMESHDS_Mesh* meshDS = aMesh.GetMeshDS();
135 SMESH_subMesh* theSubMesh = aMesh.GetSubMesh(aShape);
136 //const Handle(SMESHDS_SubMesh)& subMeshDS = theSubMesh->GetSubMeshDS();
138 map<int, const SMDS_MeshNode*> netgenToDS;
140 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Checking the mesh Faces");
142 // check if all faces were meshed by a triangle mesher (here MESFISTO_2D)
144 vector<SMESH_subMesh*> meshFaces;
145 vector<TopoDS_Shape> shapeFaces;
147 for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next())
149 TopoDS_Shape aShapeFace = exp.Current();
150 SMESH_subMesh* aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
152 int internal_size = meshFaces.size();
154 for (int i = 0;i<internal_size;i++)
156 if (aSubMesh == meshFaces[i]) index = 1;
158 if (index == 0) meshFaces.push_back(aSubMesh);
160 internal_size = shapeFaces.size();
162 for (int i = 0;i<internal_size;i++)
164 if (aShapeFace == shapeFaces[i]) index = 1;
166 if (index == 0) shapeFaces.push_back(aShapeFace);
169 int numberOfFaces = meshFaces.size();
170 int numberOfShapeFaces = shapeFaces.size();
172 SCRUTE(numberOfFaces);
173 SCRUTE(numberOfShapeFaces);
178 int NbTotOfNodesFaces = 0;
180 for (int i=0; i<numberOfFaces; i++)
182 TopoDS_Shape aShapeFace = meshFaces[i]->GetSubShape();
183 TopoDS_Shape aFace = shapeFaces[i];
184 SMESH_Algo* algoFace = _gen->GetAlgo(aMesh, aShapeFace);
185 string algoFaceName = algoFace->GetName();
186 SCRUTE(algoFaceName);
187 if (algoFaceName != "MEFISTO_2D")
189 SCRUTE(algoFaceName);
194 bool orientationMeshFace = (aFace.Orientation() == aShapeFace.Orientation());
196 const SMESHDS_SubMesh* aSubMeshDSFace = meshFaces[i]->GetSubMeshDS();
197 SCRUTE(aSubMeshDSFace);
199 int nbNodes = aSubMeshDSFace->NbNodes();
200 NbTotOfNodesFaces += nbNodes;
201 int nbTria = aSubMeshDSFace->NbElements();
202 NbTotOfTria += nbTria;
205 MESSAGE("NETGENPlugin_NETGEN_3D::Compute The mesh Face " << (i+1) << " has " << nbNodes << " face internal Nodes, " << nbTria << " triangles");
207 SCRUTE(orientationMeshFace);
209 if (orientationMeshFace)
211 MESSAGE("The mesh and face have the same orientation");
215 MESSAGE("The mesh and face have different orientations");
218 SMDS_NodeIteratorPtr iteratorNodes = aSubMeshDSFace->GetNodes();
221 while(iteratorNodes->more())
224 const SMDS_MeshNode * node = iteratorNodes->next();
225 // int nodeId = node->GetID();
226 // double nodeX = node->X();
227 // double nodeY = node->Y();
228 // double nodeZ = node->Z();
229 // MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
234 SMDS_ElemIteratorPtr iteratorTriangle = aSubMeshDSFace->GetElements();
238 int numberOfDegeneratedTriangle = 0;
239 while(iteratorTriangle->more())
242 const SMDS_MeshElement * triangle = iteratorTriangle->next();
243 int triangleId = triangle->GetID();
245 SMDS_ElemIteratorPtr triangleNodesIt = triangle->nodesIterator();
247 const SMDS_MeshNode * node1 = static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
248 double node1X = node1->X();
249 double node1Y = node1->Y();
250 double node1Z = node1->Z();
252 const SMDS_MeshNode * node2 = static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
253 double node2X = node2->X();
254 double node2Y = node2->Y();
255 double node2Z = node2->Z();
257 const SMDS_MeshNode * node3 = static_cast<const SMDS_MeshNode *>(triangleNodesIt->next());
258 double node3X = node3->X();
259 double node3Y = node3->Y();
260 double node3Z = node3->Z();
262 int triangleNode1 = node1->GetID();
263 int triangleNode2 = node2->GetID();
264 int triangleNode3 = node3->GetID();
266 // Compute the triangle surface
268 double vect1 = ((node2Y - node1Y)*(node3Z - node1Z) - (node2Z - node1Z)*(node3Y - node1Y));
269 double vect2 = - ((node2X - node1X)*(node3Z - node1Z) - (node2Z - node1Z)*(node3X - node1X));
270 double vect3 = ((node2X - node1X)*(node3Y - node1Y) - (node2Y - node1Y)*(node3X - node1X));
271 double epsilon = 1.0e-6;
273 bool triangleIsDegenerated = ((abs(vect1)<epsilon) && (abs(vect2)<epsilon) && (abs(vect3)<epsilon));
275 if (triangleIsDegenerated)
277 // MESSAGE("TRIANGLE -> ID = " << triangleId << " N1 = " << triangleNode1 << " N2 = " << triangleNode2 << " N3 = " << triangleNode3 << " is degenerated");
278 // MESSAGE("NODE -> ID = " << triangleNode1 << " X = " << node1X << " Y = " << node1Y << " Z = " << node1Z);
279 // MESSAGE("NODE -> ID = " << triangleNode2 << " X = " << node2X << " Y = " << node2Y << " Z = " << node2Z);
280 // MESSAGE("NODE -> ID = " << triangleNode3 << " X = " << node3X << " Y = " << node3Y << " Z = " << node3Z);
281 numberOfDegeneratedTriangle++;
285 // MESSAGE("TRIANGLE -> ID = " << triangleId << " N1 = " << triangleNode1 << " N2 = " << triangleNode2 << " N3 = " << triangleNode3 << " is normal");
289 if (numberOfDegeneratedTriangle > 0)
290 MESSAGE("WARNING THERE IS(ARE) " << numberOfDegeneratedTriangle << " degenerated triangle on this face");
298 SCRUTE(NbTotOfNodesFaces);
300 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Checking the mesh Edges");
302 // check if all edges were meshed by a edge mesher (here Regular_1D)
304 vector<SMESH_subMesh*> meshEdges;
305 for (TopExp_Explorer exp(aShape,TopAbs_EDGE);exp.More();exp.Next())
307 SMESH_subMesh* aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
309 int internal_size = meshEdges.size();
311 for (int i = 0;i<internal_size;i++)
313 if (aSubMesh == meshEdges[i]) index = 1;
315 if (index == 0) meshEdges.push_back(aSubMesh);
318 int numberOfEdges = meshEdges.size();
319 SCRUTE(numberOfEdges);
323 int NbTotOfNodesEdges = 0;
326 for (int i=0; i<numberOfEdges; i++)
328 TopoDS_Shape aShapeEdge = meshEdges[i]->GetSubShape();
329 SMESH_Algo* algoEdge = _gen->GetAlgo(aMesh, aShapeEdge);
330 string algoEdgeName = algoEdge->GetName();
331 SCRUTE(algoEdgeName);
332 if (algoEdgeName != "Regular_1D")
334 SCRUTE(algoEdgeName);
339 const SMESHDS_SubMesh* aSubMeshDSEdge = meshEdges[i]->GetSubMeshDS();
340 SCRUTE(aSubMeshDSEdge);
342 int nbNodes = aSubMeshDSEdge->NbNodes();
343 NbTotOfNodesEdges += nbNodes;
344 int nbSegs = aSubMeshDSEdge->NbElements();
345 NbTotOfSegs += nbSegs;
347 MESSAGE("NETGENPlugin_NETGEN_3D::Compute The mesh Edge " << (i+1) << " has " << nbNodes << " edge internal Nodes, " << nbSegs << " segments");
349 SMDS_NodeIteratorPtr iteratorNodes = aSubMeshDSEdge->GetNodes();
352 while(iteratorNodes->more())
355 const SMDS_MeshNode * node = iteratorNodes->next();
356 // int nodeId = node->GetID();
357 // double nodeX = node->X();
358 // double nodeY = node->Y();
359 // double nodeZ = node->Z();
360 // MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
366 SCRUTE(NbTotOfNodesEdges);
369 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Checking the mesh Vertices");
371 vector<SMESH_subMesh*> meshVertices;
372 for (TopExp_Explorer exp(aShape,TopAbs_VERTEX);exp.More();exp.Next())
374 SMESH_subMesh* aSubMesh = aMesh.GetSubMeshContaining(exp.Current());
376 int internal_size = meshVertices.size();
378 for (int i = 0;i<internal_size;i++)
380 if (aSubMesh == meshVertices[i]) index = 1;
382 if (index == 0) meshVertices.push_back(aSubMesh);
385 int numberOfVertices = meshVertices.size();
386 SCRUTE(numberOfVertices);
390 int NbTotOfNodesVertices = 0;
392 for (int i=0; i<numberOfVertices; i++)
394 TopoDS_Shape aShapeVertex = meshVertices[i]->GetSubShape();
396 const SMESHDS_SubMesh * aSubMeshDSVertex = meshVertices[i]->GetSubMeshDS();
397 SCRUTE(aSubMeshDSVertex);
399 int nbNodes = aSubMeshDSVertex->NbNodes();
400 NbTotOfNodesVertices += nbNodes;
402 MESSAGE("NETGENPlugin_NETGEN_3D::Compute The mesh Vertex " << (i+1) << " has " << nbNodes << " Nodes");
404 SMDS_NodeIteratorPtr iteratorNodes = aSubMeshDSVertex->GetNodes();
407 while(iteratorNodes->more())
410 const SMDS_MeshNode * node = iteratorNodes->next();
411 // int nodeId = node->GetID();
412 // double nodeX = node->X();
413 // double nodeY = node->Y();
414 // double nodeZ = node->Z();
415 // MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
421 SCRUTE(NbTotOfNodesVertices);
423 MESSAGE("NETGENPlugin_NETGEN_3D::Compute --> Analysis of all shell mesh");
425 vector<SMESH_subMesh*> meshShells;
428 for (TopExp_Explorer exp(aShape,TopAbs_SHELL);exp.More();exp.Next())
430 SMESH_subMesh* aSubMesh = aMesh.GetSubMesh(exp.Current());
433 aShell = TopoDS::Shell(exp.Current());
434 meshShells.push_back(aSubMesh);
437 int numberOfShells = meshShells.size();
438 SCRUTE(numberOfShells);
440 if (numberOfShells == 1)
442 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Only one shell --> generation of the mesh using directly Netgen");
445 Prepare the Netgen surface mesh from the SMESHDS
448 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Prepare the Netgen surface mesh from the SMESHDS");
450 int spaceDimension = 3;
451 int nbNodesByTri = 3;
452 int nbNodesByTetra = 4;
454 int Netgen_NbOfNodes = NbTotOfNodesFaces +
456 NbTotOfNodesVertices;
457 int Netgen_NbOfTria = NbTotOfTria;
458 int Netgen_param2ndOrder = 0;
459 double Netgen_paramFine = 1.;
460 double Netgen_paramSize = _maxElementVolume;
462 SCRUTE(Netgen_NbOfNodes);
463 SCRUTE(Netgen_NbOfTria);
465 double * Netgen_Coordinates = new double [spaceDimension*
467 int * listNodeCoresNetgenSmesh = new int [Netgen_NbOfNodes];
468 int * Netgen_Connectivity = new int [nbNodesByTri*Netgen_NbOfTria];
469 double * Netgen_point = new double [spaceDimension];
470 int * Netgen_triangle = new int [nbNodesByTri];
471 int * Netgen_tetrahedron = new int [nbNodesByTetra];
473 for (int i=0; i<Netgen_NbOfTria; i++)
475 for (int j=0; j<nbNodesByTri; j++)
476 Netgen_Connectivity[i*nbNodesByTri+j] = 0;
479 double bigNumber = 1.e20;
481 for (int i=0; i<Netgen_NbOfNodes; i++)
483 listNodeCoresNetgenSmesh[i] = 0;
484 for (int j=0; j<spaceDimension; j++)
485 Netgen_Coordinates[i*spaceDimension+j] = bigNumber;
489 for (int i=0; i<numberOfVertices; i++)
491 const SMESHDS_SubMesh * aSubMeshDSVertex =
492 meshVertices[i]->GetSubMeshDS();
494 SMDS_NodeIteratorPtr iteratorNodes = aSubMeshDSVertex->GetNodes();
496 while(iteratorNodes->more())
498 const SMDS_MeshNode * node = iteratorNodes->next();
499 int nodeId = node->GetID();
500 double nodeX = node->X();
501 double nodeY = node->Y();
502 double nodeZ = node->Z();
503 // MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
504 listNodeCoresNetgenSmesh[indexNodes] = nodeId;
505 int index = indexNodes*spaceDimension;
506 Netgen_Coordinates[index] = nodeX;
507 Netgen_Coordinates[index+1] = nodeY;
508 Netgen_Coordinates[index+2] = nodeZ;
509 netgenToDS[indexNodes] = node;
514 for (int i=0; i<numberOfEdges; i++)
516 const SMESHDS_SubMesh * aSubMeshDSEdge =
517 meshEdges[i]->GetSubMeshDS();
519 SMDS_NodeIteratorPtr iteratorNodes = aSubMeshDSEdge->GetNodes();
521 while(iteratorNodes->more())
523 const SMDS_MeshNode * node = iteratorNodes->next();
524 int nodeId = node->GetID();
525 double nodeX = node->X();
526 double nodeY = node->Y();
527 double nodeZ = node->Z();
528 // MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
529 listNodeCoresNetgenSmesh[indexNodes] = node->GetID();
530 int index = indexNodes*spaceDimension;
531 Netgen_Coordinates[index] = node->X();
532 Netgen_Coordinates[index+1] = node->Y();
533 Netgen_Coordinates[index+2] = node->Z();
534 netgenToDS[indexNodes] = node;
539 for (int i=0; i<numberOfFaces; i++)
541 const SMESHDS_SubMesh * aSubMeshDSFace =
542 meshFaces[i]->GetSubMeshDS();
544 SMDS_NodeIteratorPtr iteratorNodes = aSubMeshDSFace->GetNodes();
546 while(iteratorNodes->more())
548 const SMDS_MeshNode * node = iteratorNodes->next();
549 int nodeId = node->GetID();
550 double nodeX = node->X();
551 double nodeY = node->Y();
552 double nodeZ = node->Z();
553 // MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
554 listNodeCoresNetgenSmesh[indexNodes] = nodeId;
555 int index = indexNodes*spaceDimension;
556 Netgen_Coordinates[index] = nodeX;
557 Netgen_Coordinates[index+1] = nodeY;
558 Netgen_Coordinates[index+2] = nodeZ;
559 netgenToDS[indexNodes] = node;
566 for (int i=0; i<Netgen_NbOfNodes; i++)
568 ASSERT(listNodeCoresNetgenSmesh[i] != 0);
570 for (int j=0; j<Netgen_NbOfNodes && j!=i; j++)
571 ASSERT(listNodeCoresNetgenSmesh[i] != listNodeCoresNetgenSmesh[j]);
573 for (int j=0; j<spaceDimension; j++)
574 ASSERT(Netgen_Coordinates[i*spaceDimension+j] != bigNumber);
578 for (int i=0; i<numberOfFaces; i++)
580 const SMESHDS_SubMesh * aSubMeshDSFace =
581 meshFaces[i]->GetSubMeshDS();
583 TopoDS_Shape aFace = shapeFaces[i];
585 SMDS_ElemIteratorPtr iteratorTriangle = aSubMeshDSFace->GetElements();
587 TopoDS_Shape aShapeFace = meshFaces[i]->GetSubShape();
589 bool orientationMeshFace = (aFace.Orientation() == aShapeFace.Orientation());
591 SCRUTE(orientationMeshFace);
593 if (orientationMeshFace)
595 MESSAGE("The mesh and face have the same orientation");
597 while(iteratorTriangle->more())
599 const SMDS_MeshElement * triangle = iteratorTriangle->next();
600 int triangleId = triangle->GetID();
602 SMDS_ElemIteratorPtr triangleNodesIt = triangle->nodesIterator();
604 int triangleNode1 = (triangleNodesIt->next())->GetID();
605 int triangleNode2 = (triangleNodesIt->next())->GetID();
606 int triangleNode3 = (triangleNodesIt->next())->GetID();
608 // MESSAGE("TRIANGLE -> ID = " << triangleId << " N1 = " << triangleNode1 << " N2 = " << triangleNode2 << " N3 = " << triangleNode3);
613 int index = indexTrias*nbNodesByTri;
615 for (int j=0; j<Netgen_NbOfNodes; j++)
619 if (triangleNode1 == listNodeCoresNetgenSmesh[j])
621 else if (triangleNode2 == listNodeCoresNetgenSmesh[j])
623 else if (triangleNode3 == listNodeCoresNetgenSmesh[j])
627 triangleNode1 = N1New;
628 triangleNode2 = N2New;
629 triangleNode3 = N3New;
631 Netgen_Connectivity[index] = triangleNode1;
632 Netgen_Connectivity[index+1] = triangleNode2;
633 Netgen_Connectivity[index+2] = triangleNode3;
640 MESSAGE("The mesh and face have different orientations");
642 while(iteratorTriangle->more())
644 const SMDS_MeshElement * triangle = iteratorTriangle->next();
645 int triangleId = triangle->GetID();
647 SMDS_ElemIteratorPtr triangleNodesIt = triangle->nodesIterator();
649 int triangleNode1 = (triangleNodesIt->next())->GetID();
650 int triangleNode3 = (triangleNodesIt->next())->GetID();
651 int triangleNode2 = (triangleNodesIt->next())->GetID();
653 // MESSAGE("TRIANGLE -> ID = " << triangleId << " N1 = " << triangleNode1 << " N2 = " << triangleNode2 << " N3 = " << triangleNode3);
658 int index = indexTrias*nbNodesByTri;
660 for (int j=0; j<Netgen_NbOfNodes; j++)
664 if (triangleNode1 == listNodeCoresNetgenSmesh[j])
666 else if (triangleNode2 == listNodeCoresNetgenSmesh[j])
668 else if (triangleNode3 == listNodeCoresNetgenSmesh[j])
672 triangleNode1 = N1New;
673 triangleNode2 = N2New;
674 triangleNode3 = N3New;
676 Netgen_Connectivity[index] = triangleNode1;
677 Netgen_Connectivity[index+1] = triangleNode2;
678 Netgen_Connectivity[index+2] = triangleNode3;
687 int * nodesUsed = new int[Netgen_NbOfNodes];
689 for (int i=0; i<Netgen_NbOfNodes; i++) nodesUsed[i] = 0;
691 for (int i=0; i<Netgen_NbOfTria; i++)
692 for (int j=0; j<nbNodesByTri; j++)
694 int Nij = Netgen_Connectivity[i*nbNodesByTri+j];
696 ASSERT((Nij>=1) && (Nij<=Netgen_NbOfNodes));
698 nodesUsed[Nij-1] = 1;
699 Netgen_Connectivity[i*nbNodesByTri+j] = Nij;
702 for (int i=0; i<Netgen_NbOfNodes; i++)
704 ASSERT(nodesUsed[i] != 0);
710 Feed the Netgen surface mesh
713 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Feed the Netgen surface mesh");
715 Ng_Mesh * Netgen_mesh;
719 Netgen_mesh = Ng_NewMesh();
721 Ng_Meshing_Parameters Netgen_param;
723 for (int i=0; i<Netgen_NbOfNodes; i++)
725 for (int j=0; j<spaceDimension; j++)
726 Netgen_point[j] = Netgen_Coordinates[i*spaceDimension+j];
728 Ng_AddPoint(Netgen_mesh, Netgen_point);
731 for (int i=0; i<Netgen_NbOfTria; i++)
733 for (int j=0; j<nbNodesByTri; j++)
734 Netgen_triangle[j] = Netgen_Connectivity[i*nbNodesByTri+j];
736 Ng_AddSurfaceElement(Netgen_mesh, NG_TRIG, Netgen_triangle);
739 SCRUTE(Netgen_paramSize);
741 Netgen_param.secondorder = Netgen_param2ndOrder;
742 Netgen_param.fineness = Netgen_paramFine;
743 Netgen_param.maxh = Netgen_paramSize;
746 Generate the volume mesh
749 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Generate the volume mesh");
751 SCRUTE(Netgen_NbOfNodes);
752 SCRUTE(Netgen_NbOfTria);
754 SCRUTE(Ng_GetNP(Netgen_mesh));
755 SCRUTE(Ng_GetNE(Netgen_mesh));
756 SCRUTE(Ng_GetNSE(Netgen_mesh));
758 ASSERT(Netgen_NbOfNodes == Ng_GetNP(Netgen_mesh));
759 ASSERT(Ng_GetNE(Netgen_mesh) == 0);
760 ASSERT(Netgen_NbOfTria == Ng_GetNSE(Netgen_mesh));
765 status = Ng_GenerateVolumeMesh(Netgen_mesh, &Netgen_param);
767 MESSAGE("An exception has been caught during the Volume Mesh Generation ...");
768 status = NG_VOLUME_FAILURE;
773 int Netgen_NbOfNodesNew = Ng_GetNP(Netgen_mesh);
775 int Netgen_NbOfTetra = Ng_GetNE(Netgen_mesh);
777 SCRUTE(Netgen_NbOfNodesNew);
779 SCRUTE(Netgen_NbOfTetra);
781 if ((status != NG_OK) ||
782 (Netgen_NbOfNodesNew <= Netgen_NbOfNodes) ||
783 (Netgen_NbOfTetra <= 0))
785 MESSAGE("NETGENPlugin_NETGEN_3D::Compute The Volume Mesh Generation has failed ...");
789 Free the memory needed by to generate the Netgen Mesh
792 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Free the memory needed by to generate the Netgen Mesh");
794 delete [] Netgen_Coordinates;
795 delete [] Netgen_Connectivity;
796 delete [] Netgen_point;
797 delete [] Netgen_triangle;
798 delete [] Netgen_tetrahedron;
800 delete [] listNodeCoresNetgenSmesh;
802 Ng_DeleteMesh(Netgen_mesh);
808 MESSAGE("NETGENPlugin_NETGEN_3D::Compute End of Volume Mesh Generation");
811 double * Netgen_CoordinatesNew = new double [spaceDimension*Netgen_NbOfNodesNew];
812 int * Netgen_ConnectivityNew = new int [nbNodesByTetra*Netgen_NbOfTetra];
814 for (int i=0; i<Netgen_NbOfNodesNew; i++)
816 Ng_GetPoint(Netgen_mesh, (i+1), Netgen_point);
818 for (int j=0; j<spaceDimension; j++)
819 Netgen_CoordinatesNew[i*spaceDimension+j] = Netgen_point[j];
822 for (int i=0; i<Netgen_NbOfNodes; i++)
823 for (int j=0; j<spaceDimension; j++)
824 ASSERT(Netgen_CoordinatesNew[i*spaceDimension+j] == Netgen_Coordinates[i*spaceDimension+j])
826 for (int i=0; i<Netgen_NbOfTetra; i++)
828 Ng_GetVolumeElement(Netgen_mesh, (i+1), Netgen_tetrahedron);
830 for (int j=0; j<nbNodesByTetra; j++)
831 Netgen_ConnectivityNew[i*nbNodesByTetra+j] = Netgen_tetrahedron[j];
835 Feed back the SMESHDS with the generated Nodes and Volume Elements
838 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Feed back the SMESHDS with the generated Nodes and Volume Elements");
840 int NbTotOfNodesShell = Netgen_NbOfNodesNew - Netgen_NbOfNodes;
842 SCRUTE(NbTotOfNodesShell);
844 int * listNodeShellCoresNetgenSmesh = new int [NbTotOfNodesShell];
846 for (int i=0; i<NbTotOfNodesShell; i++)
847 listNodeShellCoresNetgenSmesh[i] = 0;
849 MESSAGE("NETGENPlugin_NETGEN_3D::Compute --> Adding the New Nodes to SMESHDS");
851 for (int i=0; i<NbTotOfNodesShell; i++)
853 int index = (i+Netgen_NbOfNodes)*spaceDimension;
855 SMDS_MeshNode * node =
856 meshDS->AddNode(Netgen_CoordinatesNew[index],
857 Netgen_CoordinatesNew[index+1],
858 Netgen_CoordinatesNew[index+2]);
860 meshDS->SetNodeInVolume(node, aShell);
862 index = i+Netgen_NbOfNodes;
863 netgenToDS[index] = node;
865 listNodeShellCoresNetgenSmesh[i] = node->GetID();
868 SCRUTE(Netgen_NbOfNodesNew);
870 SCRUTE(netgenToDS.size());
872 for (int i=0; i<NbTotOfNodesShell; i++)
874 ASSERT(listNodeShellCoresNetgenSmesh[i] != 0);
876 for (int j=0; j<NbTotOfNodesShell && j!=i; j++)
877 ASSERT(listNodeShellCoresNetgenSmesh[i] != listNodeShellCoresNetgenSmesh[j]);
880 MESSAGE("NETGENPlugin_NETGEN_3D::Compute --> Adding the New elements (Tetrahedrons) to the SMESHDS");
882 for (int i=0; i<Netgen_NbOfTetra; i++)
884 int index = i*nbNodesByTetra;
885 int tetraNode1 = Netgen_ConnectivityNew[index];
886 int tetraNode2 = Netgen_ConnectivityNew[index+1];
887 int tetraNode3 = Netgen_ConnectivityNew[index+2];
888 int tetraNode4 = Netgen_ConnectivityNew[index+3];
890 const SMDS_MeshNode * node1 = netgenToDS[tetraNode1-1];
891 const SMDS_MeshNode * node2 = netgenToDS[tetraNode2-1];
892 const SMDS_MeshNode * node3 = netgenToDS[tetraNode3-1];
893 const SMDS_MeshNode * node4 = netgenToDS[tetraNode4-1];
896 if (index <= Netgen_NbOfNodes)
897 tetraNode1 = listNodeCoresNetgenSmesh[index-1];
899 tetraNode1 = listNodeShellCoresNetgenSmesh[index-Netgen_NbOfNodes-1];
902 if (index <= Netgen_NbOfNodes)
903 tetraNode2 = listNodeCoresNetgenSmesh[index-1];
905 tetraNode2 = listNodeShellCoresNetgenSmesh[index-Netgen_NbOfNodes-1];
908 if (index <= Netgen_NbOfNodes)
909 tetraNode3 = listNodeCoresNetgenSmesh[index-1];
911 tetraNode3 = listNodeShellCoresNetgenSmesh[index-Netgen_NbOfNodes-1];
914 if (index <= Netgen_NbOfNodes)
915 tetraNode4 = listNodeCoresNetgenSmesh[index-1];
917 tetraNode4 = listNodeShellCoresNetgenSmesh[index-Netgen_NbOfNodes-1];
919 SMDS_MeshVolume * elt =
920 meshDS->AddVolume(node1,node2,node3,node4);
922 meshDS->SetMeshElementOnShape(elt, aShell);
926 Free the memory needed by to generate the Netgen Mesh
929 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Free the memory needed by to generate the Netgen Mesh");
931 delete [] Netgen_Coordinates;
932 delete [] Netgen_Connectivity;
933 delete [] Netgen_CoordinatesNew;
934 delete [] Netgen_ConnectivityNew;
935 delete [] Netgen_point;
936 delete [] Netgen_triangle;
937 delete [] Netgen_tetrahedron;
939 delete [] listNodeCoresNetgenSmesh;
940 delete [] listNodeShellCoresNetgenSmesh;
942 Ng_DeleteMesh(Netgen_mesh);
950 MESSAGE("NETGENPlugin_NETGEN_3D::Compute Verification of the Shell mesh");
952 TopoDS_Shape aShapeShell = meshShells[0]->GetSubShape();
953 SMESH_Algo* algoShell = _gen->GetAlgo(aMesh, aShapeShell);
954 string algoShellName = algoShell->GetName();
955 SCRUTE(algoShellName);
956 if (algoShellName != "NETGEN_3D")
958 SCRUTE(algoShellName);
963 const SMESHDS_SubMesh * aSubMeshDSShell = meshShells[0]->GetSubMeshDS();
964 SCRUTE(&aSubMeshDSShell);
966 int nbNodes = aSubMeshDSShell->NbNodes();
967 int nbTetra = aSubMeshDSShell->NbElements();
969 MESSAGE("NETGENPlugin_NETGEN_3D::Compute The mesh Shell has " << nbNodes << " shell internal Nodes, " << nbTetra << " tetrahedrons");
971 SMDS_NodeIteratorPtr iteratorNodes = aSubMeshDSShell->GetNodes();
979 while(iteratorNodes->more())
982 const SMDS_MeshNode * node = iteratorNodes->next();
983 int nodeId = node->GetID();
984 double nodeX = node->X();
985 double nodeY = node->Y();
986 double nodeZ = node->Z();
987 // MESSAGE("NODE -> ID = " << nodeId << " X = " << nodeX << " Y = " << nodeY << " Z = " << nodeZ);
992 SMDS_ElemIteratorPtr iteratorTetra = aSubMeshDSShell->GetElements();
997 while(iteratorTetra->more())
1000 const SMDS_MeshElement * tetra = iteratorTetra->next();
1001 int tetraId = tetra->GetID();
1003 SMDS_ElemIteratorPtr tetraNodesIt = tetra->nodesIterator();
1005 int tetraNode1 = (tetraNodesIt->next())->GetID();
1006 int tetraNode2 = (tetraNodesIt->next())->GetID();
1007 int tetraNode3 = (tetraNodesIt->next())->GetID();
1008 int tetraNode4 = (tetraNodesIt->next())->GetID();
1010 // MESSAGE("TETRAHEDRON -> ID = " << tetraId << " N1 = " << tetraNode1 << " N2 = " << tetraNode2 << " N3 = " << tetraNode3 << " N4 = " << tetraNode4);
1019 SCRUTE(numberOfShells);
1020 MESSAGE("NETGENPlugin_NETGEN_3D::Compute ERROR More than one shell ????? ");
1028 //=============================================================================
1032 //=============================================================================
1034 ostream & NETGENPlugin_NETGEN_3D::SaveTo(ostream & save)
1039 //=============================================================================
1043 //=============================================================================
1045 istream & NETGENPlugin_NETGEN_3D::LoadFrom(istream & load)
1050 //=============================================================================
1054 //=============================================================================
1056 ostream & operator << (ostream & save, NETGENPlugin_NETGEN_3D & hyp)
1058 return hyp.SaveTo( save );
1061 //=============================================================================
1065 //=============================================================================
1067 istream & operator >> (istream & load, NETGENPlugin_NETGEN_3D & hyp)
1069 return hyp.LoadFrom( load );