0021238: EDF 1817 SMESH: GHS3D on quadratic meshes

   Enable generation of quadratic tetrahedrons on input quadratic faces

diff --git a/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx b/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx
index abf3db6ca4efe8ec1522e31c87949af570fd6567..48caf143285bee60c14820de211edd2c79abc88e 100644 (file)
--- a/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx
+++ b/src/GHS3DPlugin/GHS3DPlugin_GHS3D.cxx
@@ -302,12 +302,12 @@ static bool writeFaces (ofstream &             theFile,
     while ( itOnSubMesh->more() )
     {
       aFace = itOnSubMesh->next();
-      nbNodes = aFace->NbNodes();
+      nbNodes = aFace->NbCornerNodes();
 
       theFile << space << nbNodes;
 
       itOnSubFace = aFace->nodesIterator();
-      while ( itOnSubFace->more() ) {
+      while ( itOnSubFace->more() && nbNodes-- ) {
         // find GHS3D ID
         aSmdsID = itOnSubFace->next()->GetID();
         itOnMap = theSmdsToGhs3dIdMap.find( aSmdsID );
@@ -321,6 +321,7 @@ static bool writeFaces (ofstream &             theFile,
       }
 
       // (NB_NODES + 1) times: DUMMY_INT
+      nbNodes = aFace->NbCornerNodes();
       for ( int j=0; j<=nbNodes; j++)
         theFile << space << dummyint;
 
@@ -375,12 +376,12 @@ static bool writeFaces (ofstream &                      theFile,
   {
     elem = eIt->next();
     // NB_NODES PER FACE
-    nbNodes = elem->NbNodes();
+    nbNodes = elem->NbCornerNodes();
     theFile << space << nbNodes;
 
     // NODE_NB_1 NODE_NB_2 ...
     nodeIt = elem->nodesIterator();
-    while ( nodeIt->more() )
+    while ( nodeIt->more() && nbNodes-- )
     {
       // find GHS3D ID
       const SMDS_MeshNode* node = castToNode( nodeIt->next() );
@@ -390,6 +391,7 @@ static bool writeFaces (ofstream &                      theFile,
     }
 
     // (NB_NODES + 1) times: DUMMY_INT
+    nbNodes = elem->NbCornerNodes();
     for ( int i=0; i<=nbNodes; i++)
       theFile << space << dummyint;
     theFile << std::endl;
@@ -427,8 +429,13 @@ static bool writePoints (ofstream &                       theFile,
   int nbNodes = theMesh->NbNodes();
   if ( nbNodes == 0 )
     return false;
+
   int nbEnforcedVertices = theEnforcedVertices.size();
 
+  int aGhs3dID = 1;
+  SMDS_NodeIteratorPtr it = theMesh->nodesIterator();
+  const SMDS_MeshNode* node;
+
   // Issue 020674: EDF 870 SMESH: Mesh generated by Netgen not usable by GHS3D
   // The problem is in nodes on degenerated edges, we need to skip them
   if ( theHelper.HasDegeneratedEdges() )
@@ -444,13 +451,26 @@ static bool writePoints (ofstream &                       theFile,
       }
     }
   }
+
+  const bool isQuadMesh = 
+    theHelper.GetMesh()->NbEdges( ORDER_QUADRATIC ) ||
+    theHelper.GetMesh()->NbFaces( ORDER_QUADRATIC ) ||
+    theHelper.GetMesh()->NbVolumes( ORDER_QUADRATIC );
+  if ( isQuadMesh )
+  {
+    // descrease nbNodes by nb of medium nodes
+    while ( it->more() )
+    {
+      node = it->next();
+      if ( !theHelper.IsDegenShape( node->getshapeId() ))
+        nbNodes -= int( theHelper.IsMedium( node ));
+    }
+    it = theMesh->nodesIterator();
+  }
+
   const char* space    = "  ";
   const int   dummyint = 0;
 
-  int aGhs3dID = 1;
-  SMDS_NodeIteratorPtr it = theMesh->nodesIterator();
-  const SMDS_MeshNode* node;
-
   // NB_NODES
   std::cout << std::endl;
   std::cout << "The initial 2D mesh contains :" << std::endl;
@@ -467,8 +487,8 @@ static bool writePoints (ofstream &                       theFile,
   while ( it->more() )
   {
     node = it->next();
-    if ( node->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE &&
-         theHelper.IsDegenShape( node->getshapeId() )) // Issue 020674
+    if (( isQuadMesh && theHelper.IsMedium( node )) || // Issue 0021238
+        theHelper.IsDegenShape( node->getshapeId() )) // Issue 0020674
       continue;
 
     theSmdsToGhs3dIdMap.insert( make_pair( node->GetID(), aGhs3dID ));
@@ -809,7 +829,7 @@ static bool readResultFile(const int                       fileOpen,
 #ifdef WITH_SMESH_CANCEL_COMPUTE
                            GHS3DPlugin_GHS3D*              theAlgo,
 #endif
-                           SMESH_Mesh&                     theMesh,
+                           SMESH_MesherHelper&             theHelper,
                            TopoDS_Shape                    tabShape[],
                            double**                        tabBox,
                            const int                       nbShape,
@@ -826,11 +846,11 @@ static bool readResultFile(const int                       fileOpen,
   char *tetraPtr;
   char *shapePtr;
 
-  SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
+  SMESHDS_Mesh* theMeshDS = theHelper.GetMeshDS();
 
   int fileStat;
   int nbElems, nbNodes, nbInputNodes;
-  int nodeId/*, triangleId*/;
+  int nodeId;
   int nbTriangle;
   int ID, shapeID, ghs3dShapeID;
   int IdShapeRef = 1;
@@ -842,7 +862,6 @@ static bool readResultFile(const int                       fileOpen,
   const SMDS_MeshNode **node;
 
   tab    = new int[3];
-  //tabID  = new int[nbShape];
   nodeID = new int[4];
   coord  = new double[3];
   node   = new const SMDS_MeshNode*[4];
@@ -931,7 +950,7 @@ static bool readResultFile(const int                       fileOpen,
       const SMDS_MeshNode* n3 = theGhs3dIdToNodeMap[ nodeId3 ];
       try {
         OCC_CATCH_SIGNALS;
-        tabID[i] = findShapeID( theMesh, n1, n2, n3, toMeshHoles );
+        tabID[i] = findShapeID( *theHelper.GetMesh(), n1, n2, n3, toMeshHoles );
         // -- 0020330: Pb with ghs3d as a submesh
         // check that found shape is to be meshed
         if ( tabID[i] > 0 ) {
@@ -1024,7 +1043,8 @@ static bool readResultFile(const int                       fileOpen,
       }
     }
     if ( toMeshHoles || shapeID != HOLE_ID ) {
-      aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
+      aTet = theHelper.AddVolume( node[1], node[0], node[2], node[3],
+                                  /*id=*/0, /*force3d=*/false);
       theMeshDS->SetMeshElementOnShape( aTet, shapeID );
     }
 #ifdef _DEBUG_
@@ -1086,12 +1106,12 @@ static bool readResultFile(const int                      fileOpen,
 #ifdef WITH_SMESH_CANCEL_COMPUTE
                            GHS3DPlugin_GHS3D*             theAlgo,
 #endif
-                           SMESH_Mesh&                    theMesh,
+                           SMESH_MesherHelper&            theHelper,
                            TopoDS_Shape                   aSolid,
                            vector <const SMDS_MeshNode*>& theNodeByGhs3dId,
                            int                            nbEnforcedVertices)
 {
-  SMESHDS_Mesh* theMeshDS = theMesh.GetMeshDS();
+  SMESHDS_Mesh* theMeshDS = theHelper.GetMeshDS();
 
   Kernel_Utils::Localizer loc;
   struct stat  status;
@@ -1152,8 +1172,8 @@ static bool readResultFile(const int                      fileOpen,
   // volumic elements already exist
   SMESH_ElementSearcher* elemSearcher = 0;
   vector< const SMDS_MeshElement* > foundVolumes;
-  if ( theMesh.NbVolumes() > 0 )
-    elemSearcher = SMESH_MeshEditor( &theMesh ).GetElementSearcher();
+  if ( theHelper.GetMesh()->NbVolumes() > 0 )
+    elemSearcher = SMESH_MeshEditor( theHelper.GetMesh() ).GetElementSearcher();
 
   // Reading the nodeCoord and update the nodeMap
   shapeID = theMeshDS->ShapeToIndex( aSolid );
@@ -1190,7 +1210,11 @@ static bool readResultFile(const int                      fileOpen,
 
   shapePtr = ptr;
 
+  if ( theHelper.IsQuadraticMesh() != SMESH_MesherHelper::LINEAR )
+    theHelper.IsQuadraticSubMesh( aSolid );
+
   // Associating the tetrahedrons to the shapes
+
   for (int iElem = 0; iElem < nbElems; iElem++) {
 #ifdef WITH_SMESH_CANCEL_COMPUTE
     if(theAlgo->computeCanceled())
@@ -1213,8 +1237,8 @@ static bool readResultFile(const int                      fileOpen,
                                              SMDSAbs_Volume, foundVolumes ))
         continue;
     }
-    aTet = theMeshDS->AddVolume( node[1], node[0], node[2], node[3] );
-    shapeID = theMeshDS->ShapeToIndex( aSolid );
+    aTet = theHelper.AddVolume( node[1], node[0], node[2], node[3],
+                                /*id=*/0, /*force3d=*/true);
     theMeshDS->SetMeshElementOnShape( aTet, shapeID );
   }
   if ( nbElems )
@@ -1419,6 +1443,10 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
   else {
     bool toMeshHoles =
       _hyp ? _hyp->GetToMeshHoles(true) : GHS3DPlugin_Hypothesis::DefaultMeshHoles();
+
+    helper.IsQuadraticSubMesh( theShape );
+    helper.SetElementsOnShape( false );
+
     Ok = readResultFile( fileOpen,
 #ifdef WNT
                          aResultFileName.ToCString(),
@@ -1426,7 +1454,7 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
 #ifdef WITH_SMESH_CANCEL_COMPUTE
                          this,
 #endif
-                         theMesh, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap,
+                         helper, tabShape, tabBox, _nbShape, aGhs3dIdToNodeMap,
                          toMeshHoles, nbEnforcedVertices );
   }
 
@@ -1485,12 +1513,12 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
  */
 //=============================================================================
 bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
-                                SMESH_MesherHelper* aHelper)
+                                SMESH_MesherHelper* theHelper)
 {
   MESSAGE("GHS3DPlugin_GHS3D::Compute()");
 
   //SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
-  TopoDS_Shape theShape = aHelper->GetSubShape();
+  TopoDS_Shape theShape = theHelper->GetSubShape();
 
   // a unique working file name
   // to avoid access to the same files by eg different users
@@ -1587,7 +1615,7 @@ bool GHS3DPlugin_GHS3D::Compute(SMESH_Mesh&         theMesh,
 #ifdef WITH_SMESH_CANCEL_COMPUTE
                          this,
 #endif
-                         theMesh, theShape ,aNodeByGhs3dId, nbEnforcedVertices );
+                         *theHelper, theShape ,aNodeByGhs3dId, nbEnforcedVertices );
   }
   
   // ---------------------