Salome HOME
Merge from V6_4_BR 05/12/2011
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index dd79438606bbf89049efc3740d35e5a8456aa865..34cef148deca4712e7e86afcb159b758061a334c 100644 (file)
@@ -229,7 +229,8 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
     TopLoc_Location loc;
     Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
 
-    if ( surface->IsUPeriodic() || surface->IsVPeriodic() )
+    if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
+         surface->IsUClosed()   || surface->IsVClosed() )
     {
       //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
       //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface();
@@ -892,6 +893,85 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
   return true;
 }
 
+//=======================================================================
+//function : GetMediumPos
+//purpose  : Return index and type of the shape  (EDGE or FACE only) to
+//          set a medium node on
+//=======================================================================
+
+std::pair<int, TopAbs_ShapeEnum> SMESH_MesherHelper::GetMediumPos(const SMDS_MeshNode* n1,
+                                                                  const SMDS_MeshNode* n2)
+{
+  TopAbs_ShapeEnum shapeType = TopAbs_SHAPE;
+  int              shapeID = -1;
+  TopoDS_Shape     shape;
+
+  if (( myShapeID == n1->getshapeId() || myShapeID == n2->getshapeId() ) && myShapeID > 0 )
+  {
+    shapeType = myShape.ShapeType();
+    shapeID   = myShapeID;
+  }
+  else if ( n1->getshapeId() == n2->getshapeId() )
+  {
+    shapeID = n2->getshapeId();
+    shape = GetSubShapeByNode( n1, GetMeshDS() );
+  }
+  else
+  {
+    const SMDS_TypeOfPosition Pos1 = n1->GetPosition()->GetTypeOfPosition();
+    const SMDS_TypeOfPosition Pos2 = n2->GetPosition()->GetTypeOfPosition();
+
+    if ( Pos1 == SMDS_TOP_3DSPACE || Pos2 == SMDS_TOP_3DSPACE )
+    {
+    }
+    else if ( Pos1 == SMDS_TOP_FACE || Pos2 == SMDS_TOP_FACE )
+    {
+      if ( Pos1 != SMDS_TOP_FACE || Pos2 != SMDS_TOP_FACE )
+      {
+        if ( Pos1 != SMDS_TOP_FACE ) std::swap( n1,n2 );
+        TopoDS_Shape F = GetSubShapeByNode( n1, GetMeshDS() );
+        TopoDS_Shape S = GetSubShapeByNode( n2, GetMeshDS() );
+        if ( IsSubShape( S, F ))
+        {
+          shapeType = TopAbs_FACE;
+          shapeID   = n1->getshapeId();
+        }
+      }
+    }
+    else if ( Pos1 == SMDS_TOP_EDGE && Pos2 == SMDS_TOP_EDGE )
+    {
+      TopoDS_Shape E1 = GetSubShapeByNode( n1, GetMeshDS() );
+      TopoDS_Shape E2 = GetSubShapeByNode( n2, GetMeshDS() );
+      shape = GetCommonAncestor( E1, E2, *myMesh, TopAbs_FACE );
+    }
+    else if ( Pos1 == SMDS_TOP_VERTEX && Pos2 == SMDS_TOP_VERTEX )
+    {
+      TopoDS_Shape V1 = GetSubShapeByNode( n1, GetMeshDS() );
+      TopoDS_Shape V2 = GetSubShapeByNode( n2, GetMeshDS() );
+      shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_FACE );
+      if ( shape.IsNull() ) shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_EDGE );
+    }
+    else // VERTEX and EDGE
+    {
+      if ( Pos1 != SMDS_TOP_VERTEX ) std::swap( n1,n2 );
+      TopoDS_Shape V = GetSubShapeByNode( n1, GetMeshDS() );
+      TopoDS_Shape E = GetSubShapeByNode( n2, GetMeshDS() );
+      if ( IsSubShape( V, E ))
+        shape = E;
+      else
+        shape = GetCommonAncestor( V, E, *myMesh, TopAbs_FACE );
+    }
+  }
+
+  if ( !shape.IsNull() )
+  {
+    if ( shapeID < 1 )
+      shapeID = GetMeshDS()->ShapeToIndex( shape );
+    shapeType = shape.ShapeType();
+  }
+  return make_pair( shapeID, shapeType );
+}
+
 //=======================================================================
 //function : GetMediumNode
 //purpose  : Return existing or create new medium nodes between given ones
@@ -920,58 +1000,32 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
 
   // get type of shape for the new medium node
   int faceID = -1, edgeID = -1;
-  const SMDS_PositionPtr Pos1 = n1->GetPosition();
-  const SMDS_PositionPtr Pos2 = n2->GetPosition();
-
   TopoDS_Edge E; double u [2];
   TopoDS_Face F; gp_XY  uv[2];
   bool uvOK[2] = { false, false };
 
-  if( myShape.IsNull() )
-  {
-    if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE ) {
-      faceID = n1->getshapeId();
-    }
-    else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE ) {
-      faceID = n2->getshapeId();
-    }
+  pair<int, TopAbs_ShapeEnum> pos = GetMediumPos( n1, n2 );
 
-    if( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
-      edgeID = n1->getshapeId();
-    }
-    if( Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE ) {
-      edgeID = n2->getshapeId();
-    }
-  }
   // get positions of the given nodes on shapes
-  TopAbs_ShapeEnum shapeType = myShape.IsNull() ? TopAbs_SHAPE : myShape.ShapeType();
-  if ( faceID>0 || shapeType == TopAbs_FACE)
+  if ( pos.second == TopAbs_FACE )
   {
-    if( myShape.IsNull() )
-      F = TopoDS::Face(meshDS->IndexToShape(faceID));
-    else {
-      F = TopoDS::Face(myShape);
-      faceID = myShapeID;
-    }
+    F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first ));
     uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]);
     uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]);
   }
-  else if (edgeID>0 || shapeType == TopAbs_EDGE)
+  else if ( pos.second == TopAbs_EDGE )
   {
-    if ( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE &&
-         Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE &&
-         n1->getshapeId() != n2->getshapeId() ) // issue 0021006
-    return getMediumNodeOnComposedWire(n1,n2,force3d);
-
-    if( myShape.IsNull() )
-      E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
-    else {
-      E = TopoDS::Edge(myShape);
-      edgeID = myShapeID;
-    }
+    if ( n1->GetPosition()->GetTypeOfPosition()==SMDS_TOP_EDGE &&
+         n2->GetPosition()->GetTypeOfPosition()==SMDS_TOP_EDGE &&
+         n1->getshapeId() != n2->getshapeId() )
+      // issue 0021006
+      return getMediumNodeOnComposedWire(n1,n2,force3d);
+    
+    E = TopoDS::Edge(meshDS->IndexToShape( edgeID = pos.first ));
     u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]);
     u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]);
   }
+
   if(!force3d)
   {
     // we try to create medium node using UV parameters of
@@ -1575,6 +1629,26 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
                                          SMESHDS_Mesh*      theMesh,
                                          SMESH_ProxyMesh*   theProxyMesh)
 {
+  return LoadNodeColumns(theParam2ColumnMap,
+                         theFace,
+                         std::list<TopoDS_Edge>(1,theBaseEdge),
+                         theMesh,
+                         theProxyMesh);
+}
+
+//=======================================================================
+//function : LoadNodeColumns
+//purpose  : Load nodes bound to face into a map of node columns
+//=======================================================================
+
+bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2ColumnMap,
+                                         const TopoDS_Face&            theFace,
+                                         const std::list<TopoDS_Edge>& theBaseSide,
+                                         SMESHDS_Mesh*                 theMesh,
+                                         SMESH_ProxyMesh*              theProxyMesh)
+{
+  // get a right submesh of theFace
+
   const SMESHDS_SubMesh* faceSubMesh = 0;
   if ( theProxyMesh )
   {
@@ -1593,22 +1667,41 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
   if ( !faceSubMesh || faceSubMesh->NbElements() == 0 )
     return false;
 
-  // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them
+  // get data of edges for normalization of params
 
-  map< double, const SMDS_MeshNode*> sortedBaseNodes;
-  if ( !SMESH_Algo::GetSortedNodesOnEdge( theMesh, theBaseEdge,/*noMedium=*/true, sortedBaseNodes)
-       || sortedBaseNodes.size() < 2 )
-    return false;
+  vector< double > length;
+  double fullLen = 0;
+  list<TopoDS_Edge>::const_iterator edge;
+  {
+    for ( edge = theBaseSide.begin(); edge != theBaseSide.end(); ++edge )
+    {
+      double len = std::max( 1e-10, SMESH_Algo::EdgeLength( *edge ));
+      fullLen += len;
+      length.push_back( len );
+    }
+  }
 
-  int nbRows = faceSubMesh->NbElements() / ( sortedBaseNodes.size()-1 ) + 1;
-  map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin();
-  double f = u_n->first, range = sortedBaseNodes.rbegin()->first - f;
-  for ( ; u_n != sortedBaseNodes.end(); u_n++ )
+  // get nodes on theBaseEdge sorted by param on edge and initialize theParam2ColumnMap with them
+  edge = theBaseSide.begin();
+  for ( int iE = 0; edge != theBaseSide.end(); ++edge, ++iE )
   {
-    double par = ( u_n->first - f ) / range;
-    vector<const SMDS_MeshNode*>& nCol = theParam2ColumnMap[ par ];
-    nCol.resize( nbRows );
-    nCol[0] = u_n->second;
+    map< double, const SMDS_MeshNode*> sortedBaseNodes;
+    SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNodes);
+    if ( sortedBaseNodes.empty() ) continue;
+
+    double f, l;
+    BRep_Tool::Range( *edge, f, l );
+    if ( edge->Orientation() == TopAbs_REVERSED ) std::swap( f, l );
+    const double coeff = 1. / ( l - f ) / length[iE] / fullLen;
+    const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first;
+    map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin();
+    for ( ; u_n != sortedBaseNodes.end(); u_n++ )
+    {
+      double par = prevPar + coeff * ( u_n->first - f );
+      TParam2ColumnMap::iterator u2nn =
+        theParam2ColumnMap.insert( theParam2ColumnMap.end(), make_pair( par, TNodeColumn()));
+      u2nn->second.push_back( u_n->second );
+    }
   }
   TParam2ColumnMap::iterator par_nVec_2, par_nVec_1 = theParam2ColumnMap.begin();
   if ( theProxyMesh )
@@ -1620,6 +1713,8 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
     }
   }
 
+  int nbRows = 1 + faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 );
+
   // fill theParam2ColumnMap column by column by passing from nodes on
   // theBaseEdge up via mesh faces on theFace
 
@@ -1630,6 +1725,8 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
   {
     vector<const SMDS_MeshNode*>& nCol1 = par_nVec_1->second;
     vector<const SMDS_MeshNode*>& nCol2 = par_nVec_2->second;
+    nCol1.resize( nbRows );
+    nCol2.resize( nbRows );
 
     int i1, i2, iRow = 0;
     const SMDS_MeshNode *n1 = nCol1[0], *n2 = nCol2[0];
@@ -1652,8 +1749,9 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
       }
       avoidSet.insert( face );
     }
-    if ( iRow + 1 < nbRows ) // compact if necessary
-      nCol1.resize( iRow + 1 ), nCol2.resize( iRow + 1 );
+    // set a real height
+    nCol1.resize( iRow + 1 );
+    nCol2.resize( iRow + 1 );
   }
   return theParam2ColumnMap.size() > 1 && theParam2ColumnMap.begin()->second.size() > 1;
 }
@@ -1888,6 +1986,30 @@ PShapeIteratorPtr SMESH_MesherHelper::GetAncestors(const TopoDS_Shape& shape,
   return PShapeIteratorPtr( new TAncestorsIterator( mesh.GetAncestors(shape), ancestorType));
 }
 
+//=======================================================================
+//function : GetCommonAncestor
+//purpose  : Find a common ancestors of two shapes of the given type
+//=======================================================================
+
+TopoDS_Shape SMESH_MesherHelper::GetCommonAncestor(const TopoDS_Shape& shape1,
+                                                   const TopoDS_Shape& shape2,
+                                                   const SMESH_Mesh&   mesh,
+                                                   TopAbs_ShapeEnum    ancestorType)
+{
+  TopoDS_Shape commonAnc;
+  if ( !shape1.IsNull() && !shape2.IsNull() )
+  {
+    PShapeIteratorPtr ancIt = GetAncestors( shape1, mesh, ancestorType );
+    while ( const TopoDS_Shape* anc = ancIt->next() )
+      if ( IsSubShape( shape2, *anc ))
+      {
+        commonAnc = *anc;
+        break;
+      }
+  }
+  return commonAnc;
+}
+
 //#include <Perf_Meter.hxx>
 
 //=======================================================================
@@ -2572,7 +2694,7 @@ namespace { // Structures used by FixQuadraticElements()
       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;
+      isStraight = _qlink->_nodeMove.SquareMagnitude() < 1/10./10. * faceSize2;
     }
     return isStraight;
   }
@@ -2945,6 +3067,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
       const SMDS_MeshElement* vol = elemIt->next();
       if ( !vol->IsQuadratic() || !volTool.Set( vol ))
         return;
+      double volMinSize2 = -1.;
       for ( int iF = 0; iF < volTool.NbFaces(); ++iF ) // loop on faces of volume
       {
         int nbN = volTool.NbFaceNodes( iF );
@@ -2957,10 +3080,17 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
           QLink link( faceNodes[iN], faceNodes[iN+2], faceNodes[iN+1] );
           pLink = links.insert( link ).first;
           faceLinks[ iN/2 ] = & *pLink;
-          if ( !isCurved )
-            isCurved = !link.IsStraight();
-          if ( link.MediumPos() == SMDS_TOP_3DSPACE && !link.IsStraight() )
-            return; // already fixed
+
+          if ( link.MediumPos() == SMDS_TOP_3DSPACE )
+          {
+            if ( !link.IsStraight() )
+              return; // already fixed
+          }
+          else if ( !isCurved )
+          {
+            if ( volMinSize2 < 0 ) volMinSize2 = volTool.MinLinearSize2();
+            isCurved = !isStraightLink( volMinSize2, link._nodeMove.SquareMagnitude() );
+          }
         }
         // store QFace
         pFace = faces.insert( QFace( faceLinks )).first;