fix centroidal smoothing of quadratic elements

src/SMESH/SMESH_MeshEditor.cxx

@@ -207,6 +207,25 @@ int SMESH_MeshEditor::FindShape (const SMDS_MeshElement * theElem)
   return 0;
 }
 
+//=======================================================================
+//function : IsMedium
+//purpose  : 
+//=======================================================================
+
+bool SMESH_MeshEditor::IsMedium(const SMDS_MeshNode*      node,
+                                const SMDSAbs_ElementType typeToCheck)
+{
+  bool isMedium = false;
+  SMDS_ElemIteratorPtr it = node->GetInverseElementIterator();
+  while (it->more()) {
+    const SMDS_MeshElement* elem = it->next();
+    isMedium = elem->IsMediumNode(node);
+    if ( typeToCheck == SMDSAbs_All || elem->GetType() == typeToCheck )
+      break;
+  }
+  return isMedium;
+}
+
 //=======================================================================
 //function : ShiftNodesQuadTria
 //purpose  : auxilary
@@ -1739,34 +1758,23 @@ void laplacianSmooth(const SMDS_MeshNode*                 theNode,
     if ( elem->GetType() != SMDSAbs_Face )
       continue;
 
-    // put all nodes in array
+    for ( int i = 0; i < elem->NbNodes(); ++i ) {
-    int nbNodes = 0, iNode = 0;
+      if ( elem->GetNode( i ) == theNode ) {
-    vector< const SMDS_MeshNode*> aNodes( elem->NbNodes() );
+        // add linked nodes
-    if(!elem->IsQuadratic()) {
+        int iBefore = i - 1;
-      SMDS_ElemIteratorPtr itN = elem->nodesIterator();
+        int iAfter = i + 1;
-      while ( itN->more() ) {
+        if ( elem->IsQuadratic() ) {
-        aNodes[ nbNodes ] = static_cast<const SMDS_MeshNode*>( itN->next() );
+          int nbCorners = elem->NbNodes() / 2;
-        if ( aNodes[ nbNodes ] == theNode )
+          if ( iAfter >= nbCorners )
-          iNode = nbNodes; // index of theNode within aNodes
+            iAfter = 0; // elem->GetNode() wraps index
-        nbNodes++;
+          if ( iBefore == -1 )
+            iBefore = nbCorners - 1;
+        }
+        nodeSet.insert( elem->GetNode( iAfter ));
+        nodeSet.insert( elem->GetNode( iBefore ));
+        break;
       }
     }
-    else {
-      int nbn = elem->NbNodes()/2;
-      aNodes.resize(nbn);
-      SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-      while ( nbNodes<nbn ) {
-        aNodes[ nbNodes ] = static_cast<const SMDS_MeshNode*>( itN->next() );
-        if ( aNodes[ nbNodes ] == theNode )
-          iNode = nbNodes; // index of theNode within aNodes
-        nbNodes++;
-      }
-    }
-    // add linked nodes
-    int iAfter = ( iNode + 1 == nbNodes ) ? 0 : iNode + 1;
-    nodeSet.insert( aNodes[ iAfter ]);
-    int iBefore = ( iNode == 0 ) ? nbNodes - 1 : iNode - 1;
-    nodeSet.insert( aNodes[ iBefore ]);
   }
 
   // compute new coodrs
@@ -1855,12 +1863,11 @@ void centroidalSmooth(const SMDS_MeshNode*                 theNode,
     }
     double elemArea = anAreaFunc.GetValue( aNodePoints );
     totalArea += elemArea;
-    elemCenter /= elem->NbNodes();
+    elemCenter /= nn;
     aNewXYZ += elemCenter * elemArea;
   }
   aNewXYZ /= totalArea;
   if ( !theSurface.IsNull() ) {
-    ASSERT( theUVMap.find( theNode ) != theUVMap.end() );
     theUVMap[ theNode ]->SetCoord( aNewXYZ.X(), aNewXYZ.Y() );
     aNewXYZ = theSurface->Value( aNewXYZ.X(), aNewXYZ.Y() ).XYZ();
   }
@@ -1915,7 +1922,7 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
   if ( theTgtAspectRatio < 1.0 )
     theTgtAspectRatio = 1.0;
 
-  double disttol = 1.e-16;
+  const double disttol = 1.e-16;
 
   SMESH::Controls::AspectRatio aQualityFunc;
 
@@ -1948,8 +1955,6 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
   // smooth elements on each TopoDS_Face separately
   // ===============================================
 
-  set<const SMDS_MeshElement*> QuadElems;
-
   set< int >::reverse_iterator fId = faceIdSet.rbegin(); // treate 0 fId at the end
   for ( ; fId != faceIdSet.rend(); ++fId ) {
     // get face surface and submesh
@@ -1978,6 +1983,7 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
     // compute UV for them
     // ---------------------------------------------------------
     bool checkBoundaryNodes = false;
+    bool isQuadratic = false;
     set<const SMDS_MeshNode*> setMovableNodes;
     map< const SMDS_MeshNode*, gp_XY* > uvMap, uvMap2;
     list< gp_XY > listUV; // uvs the 2 uvMaps refer to
@@ -2003,23 +2009,20 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
         continue;
       }
       elemsOnFace.push_back( elem );
-      if(elem->IsQuadratic()) {
-        QuadElems.insert(elem);
-      }
       theElems.erase( itElem++ );
       nbElemOnFace++;
 
+      if ( !isQuadratic )
+        isQuadratic = elem->IsQuadratic();
+
       // get movable nodes of elem
       const SMDS_MeshNode* node;
       SMDS_TypeOfPosition posType;
       SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-      int nbn =  elem->NbNodes();
+      int nn = 0, nbn =  elem->NbNodes();
       if(elem->IsQuadratic())
         nbn = nbn/2;
-      int nn = 0;
+      while ( nn++ < nbn ) {
-      //while ( itN->more() ) {
-      while ( nn<nbn ) {
-        nn++;
         node = static_cast<const SMDS_MeshNode*>( itN->next() );
         const SMDS_PositionPtr& pos = node->GetPosition();
         posType = pos.get() ? pos->GetTypeOfPosition() : SMDS_TOP_3DSPACE;
@@ -2053,7 +2056,10 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
       // get nodes to check UV
       list< const SMDS_MeshNode* > uvCheckNodes;
       itN = elem->nodesIterator();
-      while ( itN->more() ) {
+      nn = 0; nbn =  elem->NbNodes();
+      if(elem->IsQuadratic())
+        nbn = nbn/2;
+      while ( nn++ < nbn ) {
         node = static_cast<const SMDS_MeshNode*>( itN->next() );
         if ( uvMap.find( node ) == uvMap.end() )
           uvCheckNodes.push_back( node );
@@ -2149,31 +2155,17 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
       list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
       for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
         const SMDS_MeshElement* elem = (*elemIt);
-        // put elem nodes in array
+        int nbn =  elem->NbNodes();
-        vector< const SMDS_MeshNode* > nodes;
+        if(elem->IsQuadratic())
-        if(!elem->IsQuadratic()) {
+          nbn = nbn/2;
-          nodes.reserve( elem->NbNodes() + 1 );
-          SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-          while ( itN->more() )
-            nodes.push_back( static_cast<const SMDS_MeshNode*>( itN->next() ));
-        }
-        else {
-          nodes.reserve( elem->NbNodes()/2 + 1 );
-          SMDS_ElemIteratorPtr itN = elem->nodesIterator();
-          int nbn = 0;
-          while ( nbn<elem->NbNodes()/2 ) {
-            nbn++;
-            nodes.push_back( static_cast<const SMDS_MeshNode*>( itN->next() ));
-          }
-        }
-        nodes.push_back( nodes.front() );
         // loop on elem links: insert them in linkNbMap
-        for ( int iN = 1; iN < nodes.size(); ++iN ) {
+        const SMDS_MeshNode* curNode, *prevNode = elem->GetNode( nbn );
+        for ( int iN = 0; iN < nbn; ++iN ) {
+          curNode = elem->GetNode( iN );
           TLink link;
-          if ( nodes[ iN-1 ]->GetID() < nodes[ iN ]->GetID() )
+          if ( curNode < prevNode ) link = make_pair( curNode , prevNode );
-            link = make_pair( nodes[ iN-1 ], nodes[ iN ] );
+          else                      link = make_pair( prevNode , curNode );
-          else
+          prevNode = curNode;
-            link = make_pair( nodes[ iN ], nodes[ iN-1 ] );
           link_nb = linkNbMap.find( link );
           if ( link_nb == linkNbMap.end() )
             linkNbMap.insert( make_pair ( link, 1 ));
@@ -2223,8 +2215,11 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
         // get nodes on seam and its vertices
         list< const SMDS_MeshNode* > seamNodes;
         SMDS_NodeIteratorPtr nSeamIt = sm->GetNodes();
-        while ( nSeamIt->more() )
+        while ( nSeamIt->more() ) {
-          seamNodes.push_back( nSeamIt->next() );
+          const SMDS_MeshNode* node = nSeamIt->next();
+          if ( !isQuadratic || !IsMedium( node ))
+            seamNodes.push_back( node );
+        }
         TopExp_Explorer vExp( edge, TopAbs_VERTEX );
         for ( ; vExp.More(); vExp.Next() ) {
           sm = aMesh->MeshElements( vExp.Current() );
@@ -2258,7 +2253,9 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
               continue;
             int nbUseMap1 = 0, nbUseMap2 = 0;
             SMDS_ElemIteratorPtr nIt = e->nodesIterator();
-            while ( nIt->more() )
+            int nn = 0, nbn =  e->NbNodes();
+            if(e->IsQuadratic()) nbn = nbn/2;
+            while ( nn++ < nbn )
             {
               const SMDS_MeshNode* n =
                 static_cast<const SMDS_MeshNode*>( nIt->next() );
@@ -2282,7 +2279,8 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
             // be on one side of a seam
             if ( theSmoothMethod == CENTROIDAL && nbUseMap1 && nbUseMap2 ) {
               SMDS_ElemIteratorPtr nIt = e->nodesIterator();
-              while ( nIt->more() ) {
+              nn = 0;
+              while ( nn++ < nbn ) {
                 const SMDS_MeshNode* n =
                   static_cast<const SMDS_MeshNode*>( nIt->next() );
                 setMovableNodes.erase( n );
@@ -2376,38 +2374,36 @@ void SMESH_MeshEditor::Smooth (set<const SMDS_MeshElement*> & theElems,
       }
     }
 
-  } // loop on face ids
+    // move medium nodes of quadratic elements
-
+    if ( isQuadratic )
-  set< const SMDS_MeshElement* >::iterator itq;
+    {
-  for ( itq = QuadElems.begin(); itq != QuadElems.end(); itq++ ) {
+      list< const SMDS_MeshElement* >::iterator elemIt = elemsOnFace.begin();
-    const SMDS_QuadraticFaceOfNodes* QF =
+      for ( ; elemIt != elemsOnFace.end(); ++elemIt ) {
-      dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (*itq);
+        const SMDS_QuadraticFaceOfNodes* QF =
-    if(QF) {
+          dynamic_cast<const SMDS_QuadraticFaceOfNodes*> (*elemIt);
-      vector<const SMDS_MeshNode*> Ns(QF->NbNodes()+1);
+        if(QF) {
-      SMDS_NodeIteratorPtr anIter = QF->interlacedNodesIterator();
+          vector<const SMDS_MeshNode*> Ns;
-      int nn = 0;
+          Ns.reserve(QF->NbNodes()+1);
-      while( anIter->more() ) {
+          SMDS_NodeIteratorPtr anIter = QF->interlacedNodesIterator();
-        Ns[nn++] = anIter->next();
+          while ( anIter->more() )
-        //Ns.push_back(anIter->next());
+            Ns.push_back( anIter->next() );
-      }
+          Ns.push_back( Ns[0] );
-      Ns[nn] = Ns[0];
+          for(int i=0; i<QF->NbNodes(); i=i+2) {
-      //Ns.push_back(Ns[0]);
+            double x = (Ns[i]->X() + Ns[i+2]->X())/2;
-      int i=0;
+            double y = (Ns[i]->Y() + Ns[i+2]->Y())/2;
-      for(; i<QF->NbNodes(); i=i+2) {
+            double z = (Ns[i]->Z() + Ns[i+2]->Z())/2;
-        double x = (Ns[i]->X() + Ns[i+2]->X())/2;
+            if( fabs( Ns[i+1]->X() - x ) > disttol ||
-        double y = (Ns[i]->Y() + Ns[i+2]->Y())/2;
+                fabs( Ns[i+1]->Y() - y ) > disttol ||
-        double z = (Ns[i]->Z() + Ns[i+2]->Z())/2;
+                fabs( Ns[i+1]->Z() - z ) > disttol ) {
-        if( fabs( Ns[i+1]->X() - x ) > disttol ||
+              // we have to move i+1 node
-           fabs( Ns[i+1]->Y() - y ) > disttol ||
+              aMesh->MoveNode( Ns[i+1], x, y, z );
-           fabs( Ns[i+1]->Z() - z ) > disttol ) {
+            }
-          // we have to move i+1 node
+          }
-          const_cast<SMDS_MeshNode*>(Ns[i+1])->setXYZ(x,y,z);
-          aMesh->MoveNode( Ns[i+1], x, y, z );
         }
       }
     }
-  }
+    
-  QuadElems.clear();
+  } // loop on face ids
 
 }
 

src/SMESH/SMESH_MeshEditor.hxx

@@ -364,6 +364,15 @@ class SMESH_MeshEditor {
   // - not in avoidSet,
   // - in elemSet provided that !elemSet.empty()
 
+  /*!
+   * \brief Returns true if given node is medium
+    * \param n - node to check
+    * \param typeToCheck - type of elements containing the node to ask about node status
+    * \retval bool - check result
+   */
+  static bool IsMedium(const SMDS_MeshNode*      node,
+                       const SMDSAbs_ElementType typeToCheck = SMDSAbs_All);
+
   int FindShape (const SMDS_MeshElement * theElem);
   // Return an index of the shape theElem is on
   // or zero if a shape not found