0020982: EDF 1547 SMESH: Creation of non-conformal quadratic pyramids

move check of NO_FixQuadraticElements environment var to SMESH_MesherHelper

src/SMESH/SMESH_MeshEditor.cxx

@@ -9366,7 +9366,7 @@ void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d)
     }
   }
 
-  if ( !theForce3d  && !getenv("NO_FixQuadraticElements"))
+  if ( !theForce3d )
   { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
     aHelper.SetSubShape(0); // apply FixQuadraticElements() to the whole mesh
     aHelper.FixQuadraticElements();

src/SMESH/SMESH_MesherHelper.cxx

@@ -1830,6 +1830,57 @@ double SMESH_MesherHelper::GetOtherParam(const double param) const
   return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i];
 }
 
+namespace {
+
+  //=======================================================================
+  /*!
+   * \brief Iterator on ancestors of the given type
+   */
+  //=======================================================================
+
+  struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
+  {
+    TopTools_ListIteratorOfListOfShape _ancIter;
+    TopAbs_ShapeEnum                   _type;
+    TopTools_MapOfShape                _encountered;
+    TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type)
+      : _ancIter( ancestors ), _type( type )
+    {
+      if ( _ancIter.More() ) {
+        if ( _ancIter.Value().ShapeType() != _type ) next();
+        else _encountered.Add( _ancIter.Value() );
+      }
+    }
+    virtual bool more()
+    {
+      return _ancIter.More();
+    }
+    virtual const TopoDS_Shape* next()
+    {
+      const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0;
+      if ( _ancIter.More() )
+        for ( _ancIter.Next();  _ancIter.More(); _ancIter.Next())
+          if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() ))
+            break;
+      return s;
+    }
+  };
+
+} // namespace
+
+//=======================================================================
+/*!
+ * \brief Return iterator on ancestors of the given type
+ */
+//=======================================================================
+
+PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape,
+                                                   const SMESH_Mesh&   mesh,
+                                                   TopAbs_ShapeEnum    ancestorType)
+{
+  return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType));
+}
+
 //#include <Perf_Meter.hxx>
 
 //=======================================================================
@@ -1882,7 +1933,7 @@ namespace { // Structures used by FixQuadraticElements()
     void Move(const gp_Vec& move, bool sum=false) const
     { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; }
     gp_XYZ Move() const { return _nodeMove.XYZ() / _nbMoves; }
-    bool IsMoved() const { return (_nbMoves > 0 && !IsStraight()); }
+    bool IsMoved() const { return (_nbMoves > 0 /*&& !IsStraight()*/); }
     bool IsStraight() const
     { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(),
                              _nodeMove.SquareMagnitude());
@@ -1928,6 +1979,8 @@ namespace { // Structures used by FixQuadraticElements()
     const QLink* operator->() const { return _qlink; }
 
     gp_Vec Normal() const;
+
+    bool IsStraight() const;
   };
   // --------------------------------------------------------------------
   typedef list< TChainLink > TChain;
@@ -2105,9 +2158,10 @@ namespace { // Structures used by FixQuadraticElements()
 //               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 )
+            if ( face->_sides[i]->MediumPos() == pos )
               if ( const QFace* contFace = face->_sides[i]->GetContinuesFace( face ))
-                faces.push_back( contFace );
+                if ( contFace->_sides.size() == 3 )
+                  faces.push_back( contFace );
           }
         }
         faces.pop_front();
@@ -2130,10 +2184,11 @@ namespace { // Structures used by FixQuadraticElements()
     // propagate from quadrangle to neighbour faces
     if ( link->MediumPos() >= pos ) {
       int nbLinkFaces = link->_faces.size();
-      if ( nbLinkFaces == 4 || (nbLinkFaces < 4 && link->OnBoundary())) {
+      if ( nbLinkFaces == 4 || (/*nbLinkFaces < 4 && */link->OnBoundary())) {
         // hexahedral mesh or boundary quadrangles - goto a continous face
         if ( const QFace* f = link->GetContinuesFace( this ))
-          return f->GetLinkChain( *chLink, chain, pos, error );
+          if ( f->_sides.size() == 4 )
+            return f->GetLinkChain( *chLink, chain, pos, error );
       }
       else {
         TChainLink chLink(link); // side face of prismatic mesh - visit all faces of iSide
@@ -2327,7 +2382,7 @@ namespace { // Structures used by FixQuadraticElements()
     gp_Vec linkDir2(0,0,0);
     try {
       OCC_CATCH_SIGNALS;
-      if ( f1 )
+      if ( f1 && theLink->MediumPos() <= (*link1)->MediumPos() )
         len1 = f1->MoveByBoundary
           ( *link1, theRefVec, theLinks, theFaceHelper, len1, theStep-1, &linkDir1, theSign);
       else
@@ -2338,7 +2393,7 @@ namespace { // Structures used by FixQuadraticElements()
     }
     try {
       OCC_CATCH_SIGNALS;
-      if ( f2 )
+      if ( f2 && theLink->MediumPos() <= (*link2)->MediumPos() )
         len2 = f2->MoveByBoundary
           ( *link2, theRefVec, theLinks, theFaceHelper, len2, theStep-1, &linkDir2, theSign);
       else
@@ -2417,7 +2472,9 @@ namespace { // Structures used by FixQuadraticElements()
 
     if ( _faces.empty() )
       return;
-    int iFaceCont = -1;
+    int iFaceCont = -1, nbBoundary = 0, iBoundary[2]={-1,-1};
+    if ( _faces[0]->IsBoundary() )
+      iBoundary[ nbBoundary++ ] = 0;
     for ( int iF = 1; iFaceCont < 0 && iF < _faces.size(); ++iF )
     {
       // look for a face bounding none of volumes bound by _faces[0]
@@ -2428,8 +2485,21 @@ namespace { // Structures used by FixQuadraticElements()
                     _faces[iF]->_volumes[iV] == _faces[0]->_volumes[1]);
       if ( !sameVol )
         iFaceCont = iF;
+      if ( _faces[iF]->IsBoundary() )
+        iBoundary[ nbBoundary++ ] = iF;
     }
-    if ( iFaceCont > 0 ) // continues faces found, set one by the other
+    // Set continues faces: arrange _faces to have
+    // _faces[0] continues to _faces[1]
+    // _faces[2] continues to _faces[3]
+    if ( nbBoundary == 2 ) // bnd faces are continues
+    {
+      if (( iBoundary[0] < 2 ) != ( iBoundary[1] < 2 ))
+      {
+        int iNear0 = iBoundary[0] < 2 ? 1-iBoundary[0] : 5-iBoundary[0];
+        std::swap( _faces[ iBoundary[1] ], _faces[iNear0] );
+      }
+    }
+    else if ( iFaceCont > 0 ) // continues faces found
     {
       if ( iFaceCont != 1 )
         std::swap( _faces[1], _faces[iFaceCont] );
@@ -2479,6 +2549,27 @@ namespace { // Structures used by FixQuadraticElements()
     if (_qfaces[1]) norm += _qfaces[1]->_normal;
     return norm;
   }
+  //================================================================================
+  /*!
+   * \brief Test link curvature taking into account size of faces
+   */
+  //================================================================================
+
+  bool TChainLink::IsStraight() const
+  {
+    bool isStraight = _qlink->IsStraight();
+    if ( isStraight && _qfaces[0] && !_qfaces[1] )
+    {
+      int i = _qfaces[0]->LinkIndex( _qlink );
+      int iOpp = ( i + 2 ) % _qfaces[0]->_sides.size();
+      gp_XYZ mid1 = _qlink->MiddlePnt();
+      gp_XYZ mid2 = _qfaces[0]->_sides[ iOpp ]->MiddlePnt();
+      double faceSize2 = (mid1-mid2).SquareModulus();
+      isStraight = _qlink->_nodeMove.SquareMagnitude() < 1/3./3. * faceSize2;
+    }
+    return isStraight;
+  }
+  
   //================================================================================
   /*!
    * \brief Move medium nodes of vertical links of pentahedrons adjacent by side faces
@@ -2497,7 +2588,7 @@ namespace { // Structures used by FixQuadraticElements()
         bndLinks1.insert( lnk->_qlink );
       else
         interLinks.insert( lnk->_qlink );
-      isCurved = isCurved || !(*lnk)->IsStraight();
+      isCurved = isCurved || !lnk->IsStraight();
     }
     if ( !isCurved )
       return; // no need to move
@@ -2546,7 +2637,7 @@ namespace { // Structures used by FixQuadraticElements()
 
     for ( linkIt = linkSet.begin(); linkIt != linksEnd; ++linkIt)
     {
-      if ( linkIt->IsBoundary() && !(*linkIt)->IsStraight() && linkIt->_qfaces[0])
+      if ( linkIt->IsBoundary() && !linkIt->IsStraight() && linkIt->_qfaces[0])
       {
         // move iff a boundary link is bent towards inside of a face (issue 0021084)
         const QFace* face = linkIt->_qfaces[0];
@@ -2757,6 +2848,10 @@ namespace { // Structures used by FixQuadraticElements()
 
 void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
 {
+  // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
+  if ( getenv("NO_FixQuadraticElements") )
+    return;
+
   // 0. Apply algorithm to solids or geom faces
   // ----------------------------------------------
   if ( myShape.IsNull() ) {
@@ -2789,7 +2884,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
     }
     // fix nodes on geom faces
 #ifdef _DEBUG_
-    //int nbfaces = faces.Extent();
+    int nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; 
 #endif
     for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) {
       MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
@@ -2831,15 +2926,18 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
   bool isCurved = false;
   //bool hasRectFaces = false;
   //set<int> nbElemNodeSet;
+  SMDS_VolumeTool volTool;
+
+  TIDSortedNodeSet apexOfPyramid;
+  const int apexIndex = 4;
 
   if ( elemType == SMDSAbs_Volume )
   {
-    SMDS_VolumeTool volTool;
     while ( elemIt->more() ) // loop on volumes
     {
       const SMDS_MeshElement* vol = elemIt->next();
       if ( !vol->IsQuadratic() || !volTool.Set( vol ))
-        return; //continue;
+        return;
       for ( int iF = 0; iF < volTool.NbFaces(); ++iF ) // loop on faces of volume
       {
         int nbN = volTool.NbFaceNodes( iF );
@@ -2873,6 +2971,9 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
                                                faceNodes[4],faceNodes[6] );
 #endif
       }
+      // collect pyramid apexes for further correction
+      if ( vol->NbCornerNodes() == 5 )
+        apexOfPyramid.insert( vol->GetNode( apexIndex ));
     }
     set< QLink >::iterator pLink = links.begin();
     for ( ; pLink != links.end(); ++pLink )
@@ -2907,9 +3008,9 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
     return; // no curved edges of faces
 
   // 3. Compute displacement of medium nodes
-  // -------------------------------------
+  // ---------------------------------------
 
-  // two loops on faces: the first is to treat boundary links, the second is for internal ones
+  // two loops on QFaces: the first is to treat boundary links, the second is for internal ones
   TopLoc_Location loc;
   // not treat boundary of volumic submesh
   int isInside = ( elemType == SMDSAbs_Volume && volumeOnly ) ? 1 : 0;
@@ -2921,7 +3022,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
     for ( pFace = faces.begin(); pFace != faces.end(); ++pFace ) {
       if ( bool(isInside) == pFace->IsBoundary() )
         continue;
-      for ( int dir = 0; dir < 2; ++dir ) // 2 directions of propagation from quadrangle
+      for ( int dir = 0; dir < 2; ++dir ) // 2 directions of propagation from the quadrangle
       {
         MSG( "CHAIN");
         // make chain of links connected via continues faces
@@ -2954,12 +3055,12 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
         {
           TChain& chain = chains[iC];
           if ( chain.empty() ) continue;
-          if ( chain.front()->IsStraight() && chain.back()->IsStraight() ) {
+          if ( chain.front().IsStraight() && chain.back().IsStraight() ) {
             MSG("3D straight - ignore");
             continue;
           }
           if ( chain.front()->MediumPos() > bndPos ||
-               chain.back()->MediumPos() > bndPos ) {
+               chain.back() ->MediumPos() > bndPos ) {
             MSG("Internal chain - ignore");
             continue;
           }
@@ -2989,9 +3090,11 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
 
           TopoDS_Face face;
           bool checkUV = true;
-          if ( !isInside ) {
+          if ( !isInside )
-            // compute node displacement of end links in parametric space of face
+          {
-            const SMDS_MeshNode* nodeOnFace = (*(++chain.begin()))->_mediumNode;
+            // compute node displacement of end links of chain in parametric space of face
+            TChainLink& linkOnFace = *(++chain.begin());
+            const SMDS_MeshNode* nodeOnFace = linkOnFace->_mediumNode;
             TopoDS_Shape f = GetSubShapeByNode( nodeOnFace, GetMeshDS() );
             if ( !f.IsNull() && f.ShapeType() == TopAbs_FACE )
             {
@@ -3010,14 +3113,24 @@ 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());
+                isStraight[is1] = isStraightLink( (uv2-uv1).SquareModulus(),
+                                                  10 * uvMove.SquareModulus());
               }
-//               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
               }
+
+              // check if a chain is already fixed
+              gp_XY uvm = GetNodeUV( face, linkOnFace->_mediumNode, 0, &checkUV);
+              gp_XY uv1 = GetNodeUV( face, linkOnFace->node1(), nodeOnFace, &checkUV);
+              gp_XY uv2 = GetNodeUV( face, linkOnFace->node2(), nodeOnFace, &checkUV);
+              gp_XY uv12 = GetMiddleUV( surf, uv1, uv2);
+              if (( uvm - uv12 ).SquareModulus() > 1e-10 )
+              {
+                MSG("Already fixed - ignore");
+                continue;
+              }
             }
           }
           gp_Trsf trsf;
@@ -3087,60 +3200,105 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
   }
 
   // 4. Move nodes
-  // -----------
+  // -------------
 
+//   vector<const SMDS_MeshElement*> vols( 100 );
+//   vector<double>                  volSize( 100 );
+//   int nbVols;
+//   bool ok;
   for ( pLink = links.begin(); pLink != links.end(); ++pLink ) {
     if ( pLink->IsMoved() ) {
-      //gp_Pnt p = pLink->MediumPnt() + pLink->Move();
       gp_Pnt p = pLink->MiddlePnt() + pLink->Move();
       GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z());
+      //
+//       gp_Pnt pNew = pLink->MiddlePnt() + pLink->Move();
+//       if ( pLink->MediumPos() != SMDS_TOP_3DSPACE )
+//       {
+//         // avoid making distorted volumes near boundary
+//         SMDS_ElemIteratorPtr volIt =
+//           (*pLink)._mediumNode->GetInverseElementIterator( SMDSAbs_Volume );
+//         for ( nbVols = 0; volIt->more() && volTool.Set( volIt->next() ); ++nbVols )
+//         {
+//           vols   [ nbVols ] = volTool.Element();
+//           volSize[ nbVols ] = volTool.GetSize();
+//         }
+//         gp_Pnt pOld = pLink->MediumPnt();
+//         const_cast<SMDS_MeshNode*>( pLink->_mediumNode )->setXYZ( pNew.X(), pNew.Y(), pNew.Z() );
+//         ok = true;
+//         while ( nbVols-- && ok )
+//         {
+//           volTool.Set( vols[ nbVols ]);
+//           ok = ( volSize[ nbVols ] * volTool.GetSize() > 1e-20 ); 
+//         }
+//         if ( !ok )
+//         {
+//           const_cast<SMDS_MeshNode*>( pLink->_mediumNode )->setXYZ( pOld.X(), pOld.Y(), pOld.Z() );
+//           MSG( "Do NOT move \t" << pLink->_mediumNode->GetID()
+//                << " because of distortion of volume " << vols[ nbVols+1 ]->GetID());
+//           continue;
+//         }
+//       }
+//       GetMeshDS()->MoveNode( pLink->_mediumNode, pNew.X(), pNew.Y(), pNew.Z() );
+    }
+  }
+
+  //return;
+
+  // issue 0020982
+  // Move the apex of pyramid together with the most curved link
+
+  TIDSortedNodeSet::iterator apexIt = apexOfPyramid.begin();
+  for ( ; apexIt != apexOfPyramid.end(); ++apexIt )
+  {
+    SMESH_TNodeXYZ apex = *apexIt;
+
+    gp_Vec maxMove( 0,0,0 );
+    double maxMoveSize2 = 0;
+
+    // shift of node index to get medium nodes between the base nodes
+    const int base2MediumShift = 5;
+
+    // find maximal movement of medium node
+    SMDS_ElemIteratorPtr volIt = apex._node->GetInverseElementIterator( SMDSAbs_Volume );
+    vector< const SMDS_MeshElement* > pyramids;
+    while ( volIt->more() )
+    {
+      const SMDS_MeshElement* pyram = volIt->next();
+      if ( pyram->GetEntityType() != SMDSEntity_Quad_Pyramid ) continue;
+      pyramids.push_back( pyram );
+
+      for ( int iBase = 0; iBase < apexIndex; ++iBase )
+      {
+        SMESH_TNodeXYZ medium = pyram->GetNode( iBase + base2MediumShift );
+        if ( medium._node->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE )
+        {
+          SMESH_TNodeXYZ n1 = pyram->GetNode( iBase );
+          SMESH_TNodeXYZ n2 = pyram->GetNode( ( iBase+1 ) % 4 );
+          gp_Pnt middle = 0.5 * ( n1 + n2 );
+          gp_Vec move( middle, medium );
+          double moveSize2 = move.SquareMagnitude();
+          if ( moveSize2 > maxMoveSize2 )
+            maxMove = move, maxMoveSize2 = moveSize2;
+        }
+      }
+    }
+
+    // move the apex
+    if ( maxMoveSize2 > 1e-20 )
+    {
+      apex += maxMove.XYZ();
+      GetMeshDS()->MoveNode( apex._node, apex.X(), apex.Y(), apex.Z());
+
+      // move medium nodes neighboring the apex to the middle
+      const int base2MediumShift_2 = 9;
+      for ( unsigned i = 0; i < pyramids.size(); ++i )
+        for ( int iBase = 0; iBase < apexIndex; ++iBase )
+        {
+          SMESH_TNodeXYZ         base = pyramids[i]->GetNode( iBase );
+          const SMDS_MeshNode* medium = pyramids[i]->GetNode( iBase + base2MediumShift_2 );
+          gp_XYZ middle = 0.5 * ( apex + base );
+          GetMeshDS()->MoveNode( medium, middle.X(), middle.Y(), middle.Z());
+        }
     }
   }
 }
-
-//=======================================================================
-/*!
- * \brief Iterator on ancestors of the given type
- */
-//=======================================================================
-
-struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
-{
-  TopTools_ListIteratorOfListOfShape _ancIter;
-  TopAbs_ShapeEnum                   _type;
-  TopTools_MapOfShape                _encountered;
-  TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type)
-    : _ancIter( ancestors ), _type( type )
-  {
-    if ( _ancIter.More() ) {
-      if ( _ancIter.Value().ShapeType() != _type ) next();
-      else _encountered.Add( _ancIter.Value() );
-    }
-  }
-  virtual bool more()
-  {
-    return _ancIter.More();
-  }
-  virtual const TopoDS_Shape* next()
-  {
-    const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0;
-    if ( _ancIter.More() )
-      for ( _ancIter.Next();  _ancIter.More(); _ancIter.Next())
-        if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() ))
-          break;
-    return s;
-  }
-};
-
-//=======================================================================
-/*!
- * \brief Return iterator on ancestors of the given type
- */
-//=======================================================================
-
-PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape,
-                                                   const SMESH_Mesh&   mesh,
-                                                   TopAbs_ShapeEnum    ancestorType)
-{
-  return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType));
-}