Salome HOME
0021559: EDF 2175 SMESH: Hexa/Tetra mixed meshes
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index eba6db261caa420e17a105739efdee79334505a4..25a0b6434a9c81ebfcc22a173964332c1e2f831b 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2011  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -78,7 +78,8 @@ namespace {
 //================================================================================
 
 SMESH_MesherHelper::SMESH_MesherHelper(SMESH_Mesh& theMesh)
-  : myParIndex(0), myMesh(&theMesh), myShapeID(0), myCreateQuadratic(false)
+  : myParIndex(0), myMesh(&theMesh), myShapeID(0), myCreateQuadratic(false),
+    myFixNodeParameters(false)
 {
   myPar1[0] = myPar2[0] = myPar1[1] = myPar2[1] = 0;
   mySetElemOnShape = ( ! myMesh->HasShapeToMesh() );
@@ -114,14 +115,21 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
 {
   SMESHDS_Mesh* meshDS = GetMeshDS();
   // we can create quadratic elements only if all elements
-  // created on subshapes of given shape are quadratic
+  // created on sub-shapes of given shape are quadratic
   // also we have to fill myTLinkNodeMap
   myCreateQuadratic = true;
   mySeamShapeIds.clear();
   myDegenShapeIds.clear();
   TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
+  if ( aSh.ShapeType()==TopAbs_COMPOUND )
+  {
+    TopoDS_Iterator subIt( aSh );
+    if ( subIt.More() )
+      subType = ( subIt.Value().ShapeType()==TopAbs_FACE ) ? TopAbs_EDGE : TopAbs_FACE;
+  }
   SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
 
+
   int nbOldLinks = myTLinkNodeMap.size();
 
   if ( !myMesh->HasShapeToMesh() )
@@ -279,7 +287,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
 //function : GetNodeUVneedInFaceNode
 //purpose  : Check if inFaceNode argument is necessary for call GetNodeUV(F,..)
 //           Return true if the face is periodic.
-//           If F is Null, answer about subshape set through IsQuadraticSubMesh() or
+//           If F is Null, answer about sub-shape set through IsQuadraticSubMesh() or
 //           * SetSubShape()
 //=======================================================================
 
@@ -433,6 +441,20 @@ void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const
   ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok));
 }
 
+//=======================================================================
+//function : ToFixNodeParameters
+//purpose  : Enables fixing node parameters on EDGEs and FACEs in 
+//           GetNodeU(...,check=true), GetNodeUV(...,check=true), CheckNodeUV() and
+//           CheckNodeU() in case if a node lies on a shape set via SetSubShape().
+//           Default is False
+//=======================================================================
+
+void SMESH_MesherHelper::ToFixNodeParameters(bool toFix)
+{
+  myFixNodeParameters = toFix;
+}
+
+
 //=======================================================================
 //function : GetUVOnSeam
 //purpose  : Select UV on either of 2 pcurves of a seam edge, closest to the given UV
@@ -643,7 +665,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
         return false;
       }
       // store the fixed UV on the face
-      if ( myShape.IsSame(F) && shapeID == myShapeID )
+      if ( myShape.IsSame(F) && shapeID == myShapeID && myFixNodeParameters )
         const_cast<SMDS_MeshNode*>(n)->SetPosition
           ( SMDS_PositionPtr( new SMDS_FacePosition( U, V )));
     }
@@ -866,7 +888,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
           return false;
         }
         // store the fixed U on the edge
-        if ( myShape.IsSame(E) && shapeID == myShapeID )
+        if ( myShape.IsSame(E) && shapeID == myShapeID && myFixNodeParameters )
           const_cast<SMDS_MeshNode*>(n)->SetPosition
             ( SMDS_PositionPtr( new SMDS_EdgePosition( U )));
       }
@@ -895,7 +917,8 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
 
 //=======================================================================
 //function : GetMediumPos
-//purpose  : Return index and type of the shape to set a medium node on
+//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,
@@ -910,51 +933,49 @@ std::pair<int, TopAbs_ShapeEnum> SMESH_MesherHelper::GetMediumPos(const SMDS_Mes
     shapeType = myShape.ShapeType();
     shapeID   = myShapeID;
   }
+  else if ( n1->getshapeId() == n2->getshapeId() )
+  {
+    shapeID = n2->getshapeId();
+    shape = GetSubShapeByNode( n1, GetMeshDS() );
+  }
   else
   {
-    const SMDS_PositionPtr Pos1 = n1->GetPosition();
-    const SMDS_PositionPtr Pos2 = n2->GetPosition();
+    const SMDS_TypeOfPosition Pos1 = n1->GetPosition()->GetTypeOfPosition();
+    const SMDS_TypeOfPosition Pos2 = n2->GetPosition()->GetTypeOfPosition();
 
-    if( Pos1->GetTypeOfPosition()==SMDS_TOP_FACE )
-    {
-      shapeType = TopAbs_FACE;
-      shapeID   = n1->getshapeId();
-    }
-    else if( Pos2->GetTypeOfPosition()==SMDS_TOP_FACE )
+    if ( Pos1 == SMDS_TOP_3DSPACE || Pos2 == SMDS_TOP_3DSPACE )
     {
-      shapeType = TopAbs_FACE;
-      shapeID   = n2->getshapeId();
     }
-    else if (Pos1->GetTypeOfPosition()==SMDS_TOP_3DSPACE ||
-             Pos2->GetTypeOfPosition()==SMDS_TOP_3DSPACE )
+    else if ( Pos1 == SMDS_TOP_FACE || Pos2 == SMDS_TOP_FACE )
     {
-    }
-    else if ( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE &&
-              Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE )
-    {
-      if ( n1->getshapeId() == n2->getshapeId() )
-      {
-        shapeType = TopAbs_EDGE;
-        shapeID   = n1->getshapeId();
-      }
-      else
+      if ( Pos1 != SMDS_TOP_FACE || Pos2 != SMDS_TOP_FACE )
       {
-        TopoDS_Shape E1 = GetSubShapeByNode( n1, GetMeshDS() );
-        TopoDS_Shape E2 = GetSubShapeByNode( n2, GetMeshDS() );
-        shape = GetCommonAncestor( E1, E2, *myMesh, TopAbs_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->GetTypeOfPosition()==SMDS_TOP_VERTEX &&
-              Pos2->GetTypeOfPosition()==SMDS_TOP_VERTEX )
+    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 );
+      shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_EDGE );
+      if ( shape.IsNull() ) shape = GetCommonAncestor( V1, V2, *myMesh, TopAbs_FACE );
     }
     else // VERTEX and EDGE
     {
-      if ( Pos1->GetTypeOfPosition()!=SMDS_TOP_VERTEX ) std::swap( n1,n2 );
+      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 ))
@@ -963,9 +984,11 @@ std::pair<int, TopAbs_ShapeEnum> SMESH_MesherHelper::GetMediumPos(const SMDS_Mes
         shape = GetCommonAncestor( V, E, *myMesh, TopAbs_FACE );
     }
   }
+
   if ( !shape.IsNull() )
   {
-    shapeID = GetMeshDS()->ShapeToIndex( shape );
+    if ( shapeID < 1 )
+      shapeID = GetMeshDS()->ShapeToIndex( shape );
     shapeType = shape.ShapeType();
   }
   return make_pair( shapeID, shapeType );
@@ -999,9 +1022,6 @@ 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 };
@@ -1017,23 +1037,27 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   }
   else if ( pos.second == TopAbs_EDGE )
   {
+    const SMDS_PositionPtr Pos1 = n1->GetPosition();
+    const SMDS_PositionPtr Pos2 = n2->GetPosition();
     if ( Pos1->GetTypeOfPosition()==SMDS_TOP_EDGE &&
          Pos2->GetTypeOfPosition()==SMDS_TOP_EDGE &&
-         n1->getshapeId() != n2->getshapeId() ) // issue 0021006
-    return getMediumNodeOnComposedWire(n1,n2,force3d);
-    
+         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)
+  if ( !force3d & uvOK[0] && uvOK[1] )
   {
     // we try to create medium node using UV parameters of
     // nodes, else - medium between corresponding 3d points
     if( ! F.IsNull() )
     {
-      if ( uvOK[0] && uvOK[1] )
+      //if ( uvOK[0] && uvOK[1] )
       {
         if ( IsDegenShape( n1->getshapeId() )) {
           if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 ));
@@ -1184,7 +1208,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_
 //purpose  : Creates a node
 //=======================================================================
 
-SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
+SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID,
+                                           double u, double v)
 {
   SMESHDS_Mesh * meshDS = GetMeshDS();
   SMDS_MeshNode* node = 0;
@@ -1194,11 +1219,11 @@ SMDS_MeshNode* SMESH_MesherHelper::AddNode(double x, double y, double z, int ID)
     node = meshDS->AddNode( x, y, z );
   if ( mySetElemOnShape && myShapeID > 0 ) {
     switch ( myShape.ShapeType() ) {
-    case TopAbs_SOLID:  meshDS->SetNodeInVolume( node, myShapeID); break;
-    case TopAbs_SHELL:  meshDS->SetNodeInVolume( node, myShapeID); break;
-    case TopAbs_FACE:   meshDS->SetNodeOnFace(   node, myShapeID); break;
-    case TopAbs_EDGE:   meshDS->SetNodeOnEdge(   node, myShapeID); break;
-    case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID); break;
+    case TopAbs_SOLID:  meshDS->SetNodeInVolume( node, myShapeID);       break;
+    case TopAbs_SHELL:  meshDS->SetNodeInVolume( node, myShapeID);       break;
+    case TopAbs_FACE:   meshDS->SetNodeOnFace(   node, myShapeID, u, v); break;
+    case TopAbs_EDGE:   meshDS->SetNodeOnEdge(   node, myShapeID, u);    break;
+    case TopAbs_VERTEX: meshDS->SetNodeOnVertex( node, myShapeID);       break;
     default: ;
     }
   }
@@ -1357,7 +1382,7 @@ SMDS_MeshFace* SMESH_MesherHelper::AddPolygonalFace (const vector<const SMDS_Mes
     for ( int i = 0; i < nodes.size(); ++i )
     {
       const SMDS_MeshNode* n1 = nodes[i];
-      const SMDS_MeshNode* n2 = nodes[(i+1)/nodes.size()];
+      const SMDS_MeshNode* n2 = nodes[(i+1)%nodes.size()];
       const SMDS_MeshNode* n12 = GetMediumNode(n1,n2,force3d);
       newNodes.push_back( n1 );
       newNodes.push_back( n12 );
@@ -1563,6 +1588,37 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
   return elem;
 }
 
+//=======================================================================
+//function : AddVolume
+//purpose  : Creates LINEAR!!!!!!!!! octahedron
+//=======================================================================
+
+SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
+                                               const SMDS_MeshNode* n2,
+                                               const SMDS_MeshNode* n3,
+                                               const SMDS_MeshNode* n4,
+                                               const SMDS_MeshNode* n5,
+                                               const SMDS_MeshNode* n6,
+                                               const SMDS_MeshNode* n7,
+                                               const SMDS_MeshNode* n8,
+                                               const SMDS_MeshNode* n9,
+                                               const SMDS_MeshNode* n10,
+                                               const SMDS_MeshNode* n11,
+                                               const SMDS_MeshNode* n12,
+                                               const int id, 
+                                               bool force3d)
+{
+  SMESHDS_Mesh * meshDS = GetMeshDS();
+  SMDS_MeshVolume* elem = 0;
+  if(id)
+    elem = meshDS->AddVolumeWithID(n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,id);
+  else
+    elem = meshDS->AddVolume(n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12);
+  if ( mySetElemOnShape && myShapeID > 0 )
+    meshDS->SetMeshElementOnShape( elem, myShapeID );
+  return elem;
+}
+
 //=======================================================================
 //function : AddPolyhedralVolume
 //purpose  : Creates polyhedron. In quadratic mesh, adds medium nodes
@@ -1619,6 +1675,24 @@ SMESH_MesherHelper::AddPolyhedralVolume (const std::vector<const SMDS_MeshNode*>
   return elem;
 }
 
+namespace
+{
+  //================================================================================
+  /*!
+   * \brief Check if a node belongs to any face of sub-mesh
+   */
+  //================================================================================
+
+  bool isNodeInSubMesh( const SMDS_MeshNode* n, const SMESHDS_SubMesh* sm )
+  {
+    SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
+    while ( fIt->more() )
+      if ( sm->Contains( fIt->next() ))
+        return true;
+    return false;
+  }
+}
+
 //=======================================================================
 //function : LoadNodeColumns
 //purpose  : Load nodes bound to face into a map of node columns
@@ -1690,13 +1764,36 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
     SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNodes);
     if ( sortedBaseNodes.empty() ) continue;
 
+    map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin();
+    if ( theProxyMesh ) // from sortedBaseNodes remove nodes not shared by faces of faceSubMesh
+    {
+      const SMDS_MeshNode* n1 = sortedBaseNodes.begin()->second;
+      const SMDS_MeshNode* n2 = sortedBaseNodes.rbegin()->second;
+      bool allNodesAreProxy = ( n1 != theProxyMesh->GetProxyNode( n1 ) &&
+                                n2 != theProxyMesh->GetProxyNode( n2 ));
+      if ( allNodesAreProxy )
+        for ( u_n = sortedBaseNodes.begin(); u_n != sortedBaseNodes.end(); u_n++ )
+          u_n->second = theProxyMesh->GetProxyNode( u_n->second );
+
+      if ( u_n = sortedBaseNodes.begin(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+      {
+        while ( ++u_n != sortedBaseNodes.end() && !isNodeInSubMesh( u_n->second, faceSubMesh ));
+        sortedBaseNodes.erase( sortedBaseNodes.begin(), u_n );
+      }
+      else if ( u_n = --sortedBaseNodes.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+      {
+        while ( u_n != sortedBaseNodes.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
+        sortedBaseNodes.erase( ++u_n, sortedBaseNodes.end() );
+      }
+      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 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++ )
+    for ( u_n = sortedBaseNodes.begin(); u_n != sortedBaseNodes.end(); u_n++ )
     {
       double par = prevPar + coeff * ( u_n->first - f );
       TParam2ColumnMap::iterator u2nn =
@@ -1704,21 +1801,16 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
       u2nn->second.push_back( u_n->second );
     }
   }
-  TParam2ColumnMap::iterator par_nVec_2, par_nVec_1 = theParam2ColumnMap.begin();
-  if ( theProxyMesh )
-  {
-    for ( ; par_nVec_1 != theParam2ColumnMap.end(); ++par_nVec_1 )
-    {
-      const SMDS_MeshNode* & n = par_nVec_1->second[0];
-      n = theProxyMesh->GetProxyNode( n );
-    }
-  }
+  if ( theParam2ColumnMap.empty() )
+    return false;
+
 
   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
 
+  TParam2ColumnMap::iterator par_nVec_1, par_nVec_2;
   par_nVec_2 = theParam2ColumnMap.begin();
   par_nVec_1 = par_nVec_2++;
   TIDSortedElemSet emptySet, avoidSet;
@@ -2982,7 +3074,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
   if ( getenv("NO_FixQuadraticElements") )
     return;
 
-  // 0. Apply algorithm to solids or geom faces
+  // 0. Apply algorithm to SOLIDs or FACEs
   // ----------------------------------------------
   if ( myShape.IsNull() ) {
     if ( !myMesh->HasShapeToMesh() ) return;
@@ -3021,6 +3113,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly)
       SMESH_MesherHelper h(*myMesh);
       h.SetSubShape( fIt.Key() );
       h.FixQuadraticElements(true);
+      h.ToFixNodeParameters(true);
     }
     //perf_print_all_meters(1);
     return;