Salome HOME
IPAL0054631: NETGEN-1D2D3D fails on two adjacent boxes with Viscous Layers
[modules/smesh.git] / src / StdMeshers / StdMeshers_QuadToTriaAdaptor.cxx
index a06fde4952fab068e6da774cedde2a2c3cc27fa9..9afac27333671429736827228f0ba52d197f49b6 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2019  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -258,6 +258,7 @@ namespace
           meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false);
       }
     }
+    return;
   }
 
   //================================================================================
@@ -306,6 +307,37 @@ namespace
 
     return false; // == "algo fails"
   }
+
+  //================================================================================
+  /*!
+   * \brief Check if a face is in a SOLID
+   */
+  //================================================================================
+
+  bool isInSolid( vector<const SMDS_MeshNode*> & faceNodes,
+                  const int                      nbNodes,
+                  const int                      solidID )
+  {
+    if ( !faceNodes[0] )
+      return true; // NOT_QUAD
+    for ( int i = 0; i < nbNodes; ++i )
+    {
+      int shapeID = faceNodes[i]->GetShapeID();
+      if ( shapeID == solidID )
+        return true;
+    }
+    faceNodes.resize( nbNodes );
+    std::vector<const SMDS_MeshElement*> vols;
+    SMDS_Mesh::GetElementsByNodes( faceNodes, vols, SMDSAbs_Volume );
+    bool inSolid = false;
+    for ( size_t i = 0; i < vols.size() && !inSolid; ++i )
+    {
+      int shapeID = vols[i]->GetShapeID();
+      inSolid = ( shapeID == solidID );
+    }
+    faceNodes.push_back( faceNodes[0] );
+    return inSolid;
+  }
 }
 
 //================================================================================
@@ -319,7 +351,7 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
                                                   set<const SMDS_MeshNode*> & nodesToMove)
 {
   // cout << endl << "Merge " << PrmI->GetID() << " " << PrmJ->GetID() << " "
-  //      << PrmI->GetNode(4) << PrmJ->GetNode(4) << endl;
+  //      << PrmI->GetNode(4)->GetID() << " " << PrmJ->GetNode(4)->GetID() << endl;
   const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
   //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
   SMESH_TNodeXYZ Pj( Nrem );
@@ -338,6 +370,9 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
   typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
   TStdElemIterator itEnd;
 
+  typedef std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > TNNMap;
+  TNNMap mediumReplaceMap;
+
   // find and remove coincided faces of merged pyramids
   vector< const SMDS_MeshElement* > inverseElems
     // copy inverse elements to avoid iteration on changing container
@@ -355,6 +390,11 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
     }
     if ( FJEqual )
     {
+      if ( FJEqual->NbNodes() == 6 ) // find medium nodes to replace
+      {
+        mediumReplaceMap.insert( std::make_pair( FJEqual->GetNode(3), FI->GetNode(5) ));
+        mediumReplaceMap.insert( std::make_pair( FJEqual->GetNode(5), FI->GetNode(3) ));
+      }
       removeTmpElement( FI );
       removeTmpElement( FJEqual );
       myRemovedTrias.insert( FI );
@@ -370,12 +410,27 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
     const SMDS_MeshElement* elem = inverseElems[i];
     nodes.assign( elem->begin_nodes(), elem->end_nodes() );
     nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode;
+    if ( !mediumReplaceMap.empty() )
+      for ( size_t iN = elem->NbCornerNodes(); iN < nodes.size(); ++iN )
+      {
+        TNNMap::iterator n2n = mediumReplaceMap.find( nodes[iN] );
+        if ( n2n != mediumReplaceMap.end() )
+          nodes[iN] = n2n->second;
+      }
     GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size());
   }
   ASSERT( Nrem->NbInverseElements() == 0 );
   GetMeshDS()->RemoveFreeNode( Nrem,
                                GetMeshDS()->MeshElements( Nrem->getshapeId()),
                                /*fromGroups=*/false);
+  if ( !mediumReplaceMap.empty() )
+    for ( TNNMap::iterator n2n = mediumReplaceMap.begin(); n2n != mediumReplaceMap.end(); ++n2n )
+    {
+      const SMDS_MeshNode* remNode = n2n->first;
+      if ( !remNode->IsNull() && remNode->NbInverseElements() == 0 )
+        GetMeshDS()->RemoveFreeNode( remNode, 0, /*fromGroups=*/false);
+    }
+  return;
 }
 
 //================================================================================
@@ -390,7 +445,7 @@ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement*    PrmI
 {
   TIDSortedElemSet adjacentPyrams;
   bool mergedPyrams = false;
-  for ( int k=0; k<4; k++ ) // loop on 4 base nodes of PrmI
+  for ( int k = 0; k < 4; k++ ) // loop on 4 base nodes of PrmI
   {
     const SMDS_MeshNode*   n = PrmI->GetNode(k);
     SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
@@ -414,6 +469,7 @@ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement*    PrmI
     for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
       MergeAdjacent( *prm, nodesToMove, true );
   }
+  return;
 }
 
 //================================================================================
@@ -474,7 +530,7 @@ static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
 //           and a segment [PC,P]
 //=======================================================================
 
-static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
+static bool HasIntersection3(const gp_Pnt& P,  const gp_Pnt& PC, gp_Pnt&       Pint,
                              const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3)
 {
   const double EPSILON = 1e-6;
@@ -522,7 +578,27 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 
   Pint = orig + dir * t;
 
-  return ( t > 0.  &&  t < segLen );
+  bool hasInt = ( t > 0.  &&  t < segLen );
+
+  if ( hasInt && det < EPSILON ) // t is inaccurate, additionally check
+  {
+    gp_XYZ triNorm = edge1 ^ edge2;
+    gp_XYZ int0vec = Pint.XYZ() - vert0;
+    gp_XYZ in      = triNorm ^ edge1; // dir inside triangle from edge1
+    double dot     = int0vec * in;
+    if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON )
+      return false;
+    in  = edge2 ^ triNorm;
+    dot = int0vec * in;
+    if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON )
+      return false;
+    gp_XYZ int1vec = Pint.XYZ() - vert1;
+    in  = triNorm ^ ( vert2 - vert1 );
+    dot = int1vec * in;
+    if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON )
+      return false;
+  }
+  return hasInt;
 }
 
 //=======================================================================
@@ -805,8 +881,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
 {
   SMESH_ProxyMesh::setMesh( aMesh );
 
-  if ( aShape.ShapeType() != TopAbs_SOLID &&
-       aShape.ShapeType() != TopAbs_SHELL )
+  if ( aShape.ShapeType() != TopAbs_SOLID )
     return false;
 
   myShape = aShape;
@@ -815,9 +890,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
 
   const SMESHDS_SubMesh * aSubMeshDSFace;
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
-  SMESH_MesherHelper helper(aMesh);
-  helper.IsQuadraticSubMesh(aShape);
-  helper.SetElementsOnShape( true );
+  SMESH_MesherHelper helper1(aMesh);
+  helper1.IsQuadraticSubMesh(aShape);
 
   if ( myElemSearcher ) delete myElemSearcher;
   vector< SMDS_ElemIteratorPtr > itVec;
@@ -841,6 +915,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
   vector<const SMDS_MeshNode*> FNodes(5);
   gp_Pnt PC;
   gp_Vec VNorm;
+  const int solidID = meshDS->ShapeToIndex( aShape );
 
   for ( TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next() )
   {
@@ -849,22 +924,46 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
       aSubMeshDSFace = aProxyMesh->GetSubMesh( aShapeFace );
     else
       aSubMeshDSFace = meshDS->MeshElements( aShapeFace );
+    if ( !aSubMeshDSFace )
+      continue;
 
     vector<const SMDS_MeshElement*> trias, quads;
     bool hasNewTrias = false;
 
-    if ( aSubMeshDSFace )
+    const bool toCheckFaceInSolid =
+      aProxyMesh ? aProxyMesh->HasPrismsOnTwoSides( meshDS->MeshElements( aShapeFace )) : false;
+    if ( toCheckFaceInSolid && !dynamic_cast< const SMESH_ProxyMesh::SubMesh* >( aSubMeshDSFace ))
+      continue; // no room for pyramids as prisms are on both sides 
+
     {
-      bool isRev = false;
-      if ( helper.NbAncestors( aShapeFace, aMesh, aShape.ShapeType() ) > 1 )
-        isRev = helper.IsReversedSubMesh( TopoDS::Face(aShapeFace) );
+      bool isRevGlob = false;
+      SMESH_MesherHelper helper2( aMesh );
+      PShapeIteratorPtr sIt = helper2.GetAncestors( aShapeFace, aMesh, aShape.ShapeType() );
+      while ( const TopoDS_Shape * solid = sIt->next() )
+        if ( !solid->IsSame( aShape ))
+        {
+          isRevGlob = helper2.IsReversedSubMesh( TopoDS::Face( aShapeFace ));
+          if ( toCheckFaceInSolid )
+            helper2.IsQuadraticSubMesh( *solid );
+          break;
+        }
 
       SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
       while ( iteratorElem->more() ) // loop on elements on a geometrical face
       {
         const SMDS_MeshElement* face = iteratorElem->next();
+
         // preparation step to get face info
-        int stat = Preparation(face, PN, VN, FNodes, PC, VNorm);
+        int stat = Preparation( face, PN, VN, FNodes, PC, VNorm );
+
+        bool isRev = isRevGlob;
+        SMESH_MesherHelper* helper = &helper1;
+        if ( toCheckFaceInSolid && !isInSolid( FNodes, face->NbCornerNodes(), solidID ))
+        {
+          isRev  = !isRevGlob;
+          helper = &helper2;
+        }
+
         switch ( stat )
         {
         case NOT_QUAD:
@@ -877,10 +976,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
             // degenerate face
             // add triangles to result map
             SMDS_MeshFace* NewFace;
+            helper->SetElementsOnShape( false );
             if(!isRev)
-              NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
+              NewFace = helper->AddFace( FNodes[0], FNodes[1], FNodes[2] );
             else
-              NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
+              NewFace = helper->AddFace( FNodes[0], FNodes[2], FNodes[1] );
             storeTmpElement( NewFace );
             trias.push_back ( NewFace );
             quads.push_back( face );
@@ -916,20 +1016,25 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
               if ( !LimitHeight( PCbest, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/true, aShape ))
                 return false;
             }
-            // create node for PCbest
-            SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
+            // create node at PCbest
+            helper->SetElementsOnShape( true );
+            SMDS_MeshNode* NewNode = helper->AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
+
+            // create a pyramid
+            SMDS_MeshVolume* aPyram;
+            if ( isRev )
+              aPyram = helper->AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
+            else
+              aPyram = helper->AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
+            myPyramids.push_back(aPyram);
 
             // add triangles to result map
-            for(i=0; i<4; i++)
+            helper->SetElementsOnShape( false );
+            for ( i = 0; i < 4; i++ )
             {
-              trias.push_back ( meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ));
+              trias.push_back ( helper->AddFace( NewNode, FNodes[i], FNodes[i+1] ));
               storeTmpElement( trias.back() );
             }
-            // create a pyramid
-            if ( isRev ) swap( FNodes[1], FNodes[3]);
-            SMDS_MeshVolume* aPyram =
-              helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
-            myPyramids.push_back(aPyram);
 
             quads.push_back( face );
             hasNewTrias = true;
@@ -1018,7 +1123,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
   vector<const SMDS_MeshElement*> myPyramids;
   SMESH_MesherHelper helper(aMesh);
   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
-  helper.SetElementsOnShape( true );
 
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
   SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh();
@@ -1116,10 +1220,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
           IsRev = true;
         }
       }
+      helper.SetElementsOnShape( false );
       if(!IsRev)
-        NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
+        NewFace = helper.AddFace( FNodes[0], FNodes[1], FNodes[2] );
       else
-        NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
+        NewFace = helper.AddFace( FNodes[0], FNodes[2], FNodes[1] );
       storeTmpElement( NewFace );
       prxSubMesh->AddElement( NewFace );
       continue;
@@ -1238,25 +1343,28 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
         return false;
 
       // create node for Papex
+      helper.SetElementsOnShape( true );
       SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() );
 
+      // create a pyramid
+      SMDS_MeshVolume* aPyram;
+      if(isRev)
+        aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
+      else
+        aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
+      myPyramids.push_back(aPyram);
+
       // add triangles to result map
+      helper.SetElementsOnShape( false );
       for ( i = 0; i < 4; i++) {
         SMDS_MeshFace* NewFace;
         if(isRev)
-          NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] );
+          NewFace = helper.AddFace( NewNode, FNodes[i], FNodes[i+1] );
         else
-          NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] );
+          NewFace = helper.AddFace( NewNode, FNodes[i+1], FNodes[i] );
         storeTmpElement( NewFace );
         prxSubMesh->AddElement( NewFace );
       }
-      // create a pyramid
-      SMDS_MeshVolume* aPyram;
-      if(isRev)
-        aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
-      else
-        aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
-      myPyramids.push_back(aPyram);
     }
   } // end loop on all faces
 
@@ -1316,7 +1424,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
       while ( vIt->more() )
       {
         const SMDS_MeshElement* PrmJ = vIt->next();
-        if ( SMESH_MeshAlgos::GetCommonNodes( PrmI, PrmJ ).size() > 1 )
+        if ( SMESH_MeshAlgos::NbCommonNodes( PrmI, PrmJ ) > 1 )
           checkedPyrams.insert( PrmJ );
       }
     }