Salome HOME
0020986: EDF 1557 SMESH: Convert to quadratic with medium node on geometry fails...
authoreap <eap@opencascade.com>
Wed, 22 Sep 2010 11:24:21 +0000 (11:24 +0000)
committereap <eap@opencascade.com>
Wed, 22 Sep 2010 11:24:21 +0000 (11:24 +0000)
   Optimize FixQuadraticElements()

src/SMESH/SMESH_MesherHelper.cxx

index bdf3690b60cd9172525f12916cd60a101bcf4933..05a27b16ac3077b3f1951e99f87ef5bebb47656a 100644 (file)
@@ -479,7 +479,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
     if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
     if ( Precision::IsInfinite( uv.X() ) ||
          Precision::IsInfinite( uv.Y() ) ||
-         nodePnt.Distance( surface->Value( uv.X(), uv.Y() )) > tol )
+         nodePnt.Distance( surface->Value( uv.X(), uv.Y() )) > 2*tol )
     {
       // uv incorrect, project the node to surface
       GeomAPI_ProjectPointOnSurf& projector = GetProjector( F, loc, tol );
@@ -492,12 +492,11 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
       Quantity_Parameter U,V;
       projector.LowerDistanceParameters(U,V);
       uv.SetCoord( U,V );
-      if ( nodePnt.Distance( surface->Value( U, V )) > tol )
+      if ( nodePnt.Distance( surface->Value( U, V )) > 2*tol )
       {
         MESSAGE( "SMESH_MesherHelper::CheckNodeUV(), invalid projection" );
         return false;
       }
-      //uv.SetCoord( U,V );
     }
     else if ( uv.Modulus() > numeric_limits<double>::min() )
     {
@@ -1390,6 +1389,8 @@ double SMESH_MesherHelper::GetOtherParam(const double param) const
   return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i];
 }
 
+//#include <Perf_Meter.hxx>
+
 //=======================================================================
 namespace { // Structures used by FixQuadraticElements()
 //=======================================================================
@@ -1399,7 +1400,12 @@ namespace { // Structures used by FixQuadraticElements()
 #define MSG(txt) __DMP__(txt<<endl)
 #define MSGBEG(txt) __DMP__(txt)
 
-  const double straightTol2 = 1e-33; // to detect straing links
+  //const double straightTol2 = 1e-33; // to detect straing links
+  bool isStraightLink(double linkLen2, double middleNodeMove2)
+  {
+    // straight if <node move> < 1/15 * <link length>
+    return middleNodeMove2 < 1/15./15. * linkLen2;
+  }
 
   struct QFace;
   // ---------------------------------------
@@ -1436,8 +1442,10 @@ namespace { // Structures used by FixQuadraticElements()
     { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; }
     gp_XYZ Move() const { return _nodeMove.XYZ() / _nbMoves; }
     bool IsMoved() const { return (_nbMoves > 0 && !IsStraight()); }
-    bool IsStraight() const { return _nodeMove.SquareMagnitude() <= straightTol2; }
-
+    bool IsStraight() const
+    { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(),
+                             _nodeMove.SquareMagnitude());
+    }
     bool operator<(const QLink& other) const {
       return (node1()->GetID() == other.node1()->GetID() ?
               node2()->GetID() < other.node2()->GetID() :
@@ -1459,7 +1467,7 @@ namespace { // Structures used by FixQuadraticElements()
     TChainLink(const QLink* qlink=0):_qlink(qlink) {
       _qfaces[0] = _qfaces[1] = 0;
     }
-    void SetFace(const QFace* face) { int iF = _qfaces[0] ? 1 : 0; _qfaces[iF]=face; }
+    void SetFace(const QFace* face) const { int iF = _qfaces[0] ? 1 : 0; _qfaces[iF]=face; }
 
     bool IsBoundary() const { return !_qfaces[1]; }
 
@@ -1637,6 +1645,7 @@ namespace { // Structures used by FixQuadraticElements()
 
     if ( _sides.size() != 4 ) { // triangle - visit all my continous faces
       MSGBEG( *this );
+      TLinkSet links;
       list< const QFace* > faces( 1, this );
       while ( !faces.empty() ) {
         const QFace* face = faces.front();
@@ -1644,12 +1653,13 @@ namespace { // Structures used by FixQuadraticElements()
           if ( !face->_sideIsAdded[i] && face->_sides[i] ) {
             face->_sideIsAdded[i] = true;
             // find a face side in the chain
-            TChain::iterator chLink = chain.begin();
-            for ( ; chLink != chain.end(); ++chLink )
-              if ( chLink->_qlink == face->_sides[i] )
-                break;
-            if ( chLink == chain.end() )
-              chLink = chain.insert( chain.begin(), TChainLink(face->_sides[i]));
+            TLinkInSet chLink = links.insert( TChainLink(face->_sides[i])).first;
+//             TChain::iterator chLink = chain.begin();
+//             for ( ; chLink != chain.end(); ++chLink )
+//               if ( chLink->_qlink == face->_sides[i] )
+//                 break;
+//             if ( chLink == chain.end() )
+//               chLink = chain.insert( chain.begin(), TChainLink(face->_sides[i]));
             // add a face to a chained link and put a continues face in the queue
             chLink->SetFace( face );
             if ( face->_sides[i]->MediumPos() >= pos )
@@ -1661,6 +1671,7 @@ namespace { // Structures used by FixQuadraticElements()
       }
       if ( error < ERR_TRI )
         error = ERR_TRI;
+      chain.insert( chain.end(), links.begin(),links.end() );
       return false;
     }
     _sideIsAdded[iSide] = true; // not to add this link to chain again
@@ -2293,37 +2304,49 @@ namespace { // Structures used by FixQuadraticElements()
 
 void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
 {
-  // apply algorithm to solids or geom faces
+  // 0. Apply algorithm to solids or geom faces
   // ----------------------------------------------
   if ( myShape.IsNull() ) {
     if ( !myMesh->HasShapeToMesh() ) return;
     SetSubShape( myMesh->GetShapeToMesh() );
 
+#ifdef _DEBUG_
+    int nbSolids = 0;
+    TopTools_IndexedMapOfShape solids;
+    TopExp::MapShapes(myShape,TopAbs_SOLID,solids);
+    nbSolids = solids.Extent();
+#endif
     TopTools_MapOfShape faces; // faces not in solid or in not meshed solid
     for ( TopExp_Explorer f(myShape,TopAbs_FACE,TopAbs_SOLID); f.More(); f.Next() ) {
-      faces.Add( f.Current() );
+      faces.Add( f.Current() ); // not in solid
     }
     for ( TopExp_Explorer s(myShape,TopAbs_SOLID); s.More(); s.Next() ) {
       if ( myMesh->GetSubMesh( s.Current() )->IsEmpty() ) { // get faces of solid
         for ( TopExp_Explorer f( s.Current(), TopAbs_FACE); f.More(); f.Next() )
-          faces.Add( f.Current() );
+          faces.Add( f.Current() ); // in not meshed solid
       }
       else { // fix nodes in the solid and its faces
+        MSG("FIX SOLID " << nbSolids-- << " #" << GetMeshDS()->ShapeToIndex(s.Current()));
         SMESH_MesherHelper h(*myMesh);
         h.SetSubShape( s.Current() );
         h.FixQuadraticElements(false);
       }
     }
     // fix nodes on geom faces
+#ifdef _DEBUG_
+    int nbfaces = faces.Extent();
+#endif
     for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) {
+      MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
       SMESH_MesherHelper h(*myMesh);
       h.SetSubShape( fIt.Key() );
       h.FixQuadraticElements(true);
     }
+    //perf_print_all_meters(1);
     return;
   }
 
-  // Find out type of elements and get iterator on them
+  // 1. Find out type of elements and get iterator on them
   // ---------------------------------------------------
 
   SMDS_ElemIteratorPtr elemIt;
@@ -2342,7 +2365,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
   if ( !elemIt || !elemIt->more() || elemType < SMDSAbs_Face )
     return;
 
-  // Fill in auxiliary data structures
+  // 2. Fill in auxiliary data structures
   // ----------------------------------
 
   set< QLink > links;
@@ -2351,7 +2374,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
   set< QFace >::iterator pFace;
 
   bool isCurved = false;
-  bool hasRectFaces = false;
+  //bool hasRectFaces = false;
   set<int> nbElemNodeSet;
 
   if ( elemType == SMDSAbs_Volume )
@@ -2384,9 +2407,9 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
         if ( pFace->NbVolumes() == 0 )
           pFace->AddSelfToLinks();
         pFace->SetVolume( vol );
-        hasRectFaces = hasRectFaces ||
-          ( volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_HEXA ||
-            volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_PENTA );
+//         hasRectFaces = hasRectFaces ||
+//           ( volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_HEXA ||
+//             volTool.GetVolumeType() == SMDS_VolumeTool::QUAD_PENTA );
 #ifdef _DEBUG_
         if ( nbN == 6 )
           pFace->_face = GetMeshDS()->FindFace(faceNodes[0],faceNodes[2],faceNodes[4]);
@@ -2422,13 +2445,13 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
       // store QFace
       pFace = faces.insert( QFace( faceLinks )).first;
       pFace->AddSelfToLinks();
-      hasRectFaces = ( hasRectFaces || nbN == 4 );
+      //hasRectFaces = ( hasRectFaces || nbN == 4 );
     }
   }
   if ( !isCurved )
     return; // no curved edges of faces
 
-  // Compute displacement of medium nodes
+  // 3. Compute displacement of medium nodes
   // -------------------------------------
 
   // two loops on faces: the first is to treat boundary links, the second is for internal ones
@@ -2519,6 +2542,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
             {
               face = TopoDS::Face( f );
               Handle(Geom_Surface) surf = BRep_Tool::Surface(face,loc);
+              bool isStraight[2];
               for ( int is1 = 0; is1 < 2; ++is1 ) // move0 or move1
               {
                 TChainLink& link = is1 ? chain.back() : chain.front();
@@ -2531,9 +2555,11 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
                 ( is1 ? move1 : move0 ).SetCoord( uvMove.X(), uvMove.Y(), 0 );
                 if ( !is1 ) // correct nodeOnFace for move1 (issue 0020919)
                   nodeOnFace = (*(++chain.rbegin()))->_mediumNode;
+                isStraight[is1] = isStraightLink( (uv2-uv1).SquareModulus(),uvMove.SquareModulus());
               }
-              if ( move0.SquareMagnitude() < straightTol2 &&
-                   move1.SquareMagnitude() < straightTol2 ) {
+//               if ( move0.SquareMagnitude() < straightTol2 &&
+//                    move1.SquareMagnitude() < straightTol2 ) {
+              if ( isStraight[0] && isStraight[1] ) {
                 MSG("2D straight - ignore");
                 continue; // straight - no need to move nodes of internal links
               }
@@ -2605,7 +2631,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
     } // loop on faces
   }
 
-  // Move nodes
+  // 4. Move nodes
   // -----------
 
   for ( pLink = links.begin(); pLink != links.end(); ++pLink ) {