Salome HOME
#16648 [CEA] RadialQuadrangle algorithm hypothesis change requires a Clear Mesh Data...
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index e03128dc84d92994c16584069d702820e70982b1..eceb92e3657040b6f37d0ec777c370f534dce5bf 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2019  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
@@ -36,6 +36,7 @@
 #include "SMESH_HypoFilter.hxx"
 #include "SMESH_Mesh.hxx"
 #include "SMESH_MeshAlgos.hxx"
+#include "SMESH_MeshEditor.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
 
@@ -73,7 +74,7 @@ using namespace std;
 
 namespace {
 
-  inline SMESH_TNodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_TNodeXYZ(n); }
+  inline SMESH_NodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_NodeXYZ(n); }
 
   enum { U_periodic = 1, V_periodic = 2 };
 }
@@ -643,22 +644,22 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
 {
   gp_Pnt2d uv( Precision::Infinite(), Precision::Infinite() );
 
-  const SMDS_PositionPtr Pos = n->GetPosition();
+  SMDS_PositionPtr pos = n->GetPosition();
   bool uvOK = false;
-  if ( Pos->GetTypeOfPosition() == SMDS_TOP_FACE )
+  if ( pos->GetTypeOfPosition() == SMDS_TOP_FACE )
   {
     // node has position on face
-    const SMDS_FacePosition* fpos = static_cast<const SMDS_FacePosition*>( Pos );
+    SMDS_FacePositionPtr fpos = pos;
     uv.SetCoord( fpos->GetUParameter(), fpos->GetVParameter() );
     if ( check )
       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 2.*getFaceMaxTol( F )); // 2. from 22830
   }
-  else if ( Pos->GetTypeOfPosition() == SMDS_TOP_EDGE )
+  else if ( pos->GetTypeOfPosition() == SMDS_TOP_EDGE )
   {
     // node has position on EDGE => it is needed to find
     // corresponding EDGE from FACE, get pcurve for this
     // EDGE and retrieve value from this pcurve
-    const SMDS_EdgePosition* epos = static_cast<const SMDS_EdgePosition*>( Pos );
+    SMDS_EdgePositionPtr epos = pos;
     const int              edgeID = n->getshapeId();
     const TopoDS_Edge& E = TopoDS::Edge( GetMeshDS()->IndexToShape( edgeID ));
     double f, l, u = epos->GetUParameter();
@@ -699,7 +700,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
       uv = newUV;
     }
   }
-  else if ( Pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
+  else if ( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
   {
     if ( int vertexID = n->getshapeId() ) {
       const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
@@ -804,33 +805,20 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
     // check that uv is correct
     TopLoc_Location loc;
     Handle(Geom_Surface) surface = BRep_Tool::Surface( F,loc );
-    gp_Pnt nodePnt = XYZ( n ), surfPnt(0,0,0);
+    SMESH_NodeXYZ nXYZ( n );
+    gp_Pnt nodePnt = nXYZ, surfPnt(0,0,0);
     double dist = 0;
     if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
     if ( infinit ||
          (dist = nodePnt.Distance( surfPnt = surface->Value( uv.X(), uv.Y() ))) > tol )
     {
       setPosOnShapeValidity( shapeID, false );
-      if ( !infinit && distXYZ ) {
-        surfPnt.Transform( loc );
-        distXYZ[0] = dist;
-        distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
-      }
       // uv incorrect, project the node to surface
-      GeomAPI_ProjectPointOnSurf& projector = GetProjector( F, loc, tol );
-      projector.Perform( nodePnt );
-      if ( !projector.IsDone() || projector.NbPoints() < 1 )
-      {
-        MESSAGE( "SMESH_MesherHelper::CheckNodeUV() failed to project" );
-        return false;
-      }
-      Quantity_Parameter U,V;
-      projector.LowerDistanceParameters(U,V);
-      uv.SetCoord( U,V );
-      surfPnt = surface->Value( U, V );
-      dist = nodePnt.Distance( surfPnt );
+      Handle(ShapeAnalysis_Surface) sprojector = GetSurface( F );
+      uv = sprojector->ValueOfUV( nXYZ, tol ).XY();
+      surfPnt = sprojector->Value( uv );
+      dist = surfPnt.Distance( nXYZ );
       if ( distXYZ ) {
-        surfPnt.Transform( loc );
         distXYZ[0] = dist;
         distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
       }
@@ -842,7 +830,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
       // store the fixed UV on the face
       if ( myShape.IsSame(F) && shapeID == myShapeID && myFixNodeParameters )
         const_cast<SMDS_MeshNode*>(n)->SetPosition
-          ( SMDS_PositionPtr( new SMDS_FacePosition( U, V )));
+          ( SMDS_PositionPtr( new SMDS_FacePosition( uv.X(), uv.Y() )));
     }
     else if ( myShape.IsSame(F) && uv.Modulus() > numeric_limits<double>::min() )
     {
@@ -854,7 +842,7 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
 
 //=======================================================================
 //function : GetProjector
-//purpose  : Return projector intitialized by given face without location, which is returned
+//purpose  : Return projector initialized by given face without location, which is returned
 //=======================================================================
 
 GeomAPI_ProjectPointOnSurf& SMESH_MesherHelper::GetProjector(const TopoDS_Face& F,
@@ -1031,8 +1019,7 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
   const SMDS_PositionPtr pos = n->GetPosition();
   if ( pos->GetTypeOfPosition()==SMDS_TOP_EDGE )
   {
-    const SMDS_EdgePosition* epos = static_cast<const SMDS_EdgePosition*>( pos );
-    param =  epos->GetUParameter();
+    param = pos->GetParameters()[0];
   }
   else if( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX )
   {
@@ -1049,6 +1036,16 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
       int vertexID = n->getshapeId();
       const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
       param =  BRep_Tool::Parameter( V, E );
+
+      if ( inEdgeNode )
+      {
+        BRepAdaptor_Curve curve( E );
+        if ( curve.IsPeriodic() )
+        {
+          double uInEdge = GetNodeU( E, inEdgeNode );
+          param += ShapeAnalysis::AdjustByPeriod( param, uInEdge, curve.Period() );
+        }
+      }
     }
   }
   if ( check )
@@ -1130,7 +1127,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
           MESSAGE( "SMESH_MesherHelper::CheckNodeU() failed to project" );
           return false;
         }
-        Quantity_Parameter U = projector->LowerDistanceParameter();
+        Standard_Real U = projector->LowerDistanceParameter();
         u = double( U );
         curvPnt = curve->Value( u );
         dist = nodePnt.Distance( curvPnt );
@@ -2179,13 +2176,30 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
     const SMDS_MeshNode* n14 = GetMediumNode( n1, n4, force3d, TopAbs_SOLID );
     const SMDS_MeshNode* n25 = GetMediumNode( n2, n5, force3d, TopAbs_SOLID );
     const SMDS_MeshNode* n36 = GetMediumNode( n3, n6, force3d, TopAbs_SOLID );
+    if ( myCreateBiQuadratic )
+    {
+      const SMDS_MeshNode* n1245 = GetCentralNode( n1,n2,n4,n5,n12,n25,n45,n14,force3d );
+      const SMDS_MeshNode* n1346 = GetCentralNode( n1,n3,n4,n6,n31,n36,n64,n14,force3d );
+      const SMDS_MeshNode* n2356 = GetCentralNode( n2,n3,n6,n5,n23,n36,n56,n25,force3d );
 
-    if(id)
-      elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
-                                     n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
+      if(id)
+        elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
+                                       n12, n23, n31, n45, n56, n64, n14, n25, n36,
+                                       n1245, n2356, n1346, id);
+      else
+        elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
+                                 n12, n23, n31, n45, n56, n64, n14, n25, n36,
+                                 n1245, n2356, n1346);
+    }
     else
-      elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
-                               n12, n23, n31, n45, n56, n64, n14, n25, n36);
+      {
+        if(id)
+          elem = meshDS->AddVolumeWithID(n1, n2, n3, n4, n5, n6,
+                                         n12, n23, n31, n45, n56, n64, n14, n25, n36, id);
+        else
+          elem = meshDS->AddVolume(n1, n2, n3, n4, n5, n6,
+                                   n12, n23, n31, n45, n56, n64, n14, n25, n36);
+      }
   }
   if ( mySetElemOnShape && myShapeID > 0 )
     meshDS->SetMeshElementOnShape( elem, myShapeID );
@@ -2669,7 +2683,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
   }
 
   // nb rows of nodes
-  size_t prevNbRows     = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here
+  size_t prevNbRows   = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here
   size_t expectNbRows = faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); // to be added
 
   // fill theParam2ColumnMap column by column by passing from nodes on
@@ -2887,7 +2901,7 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh*      faceSM,
     for ( size_t i = 0; i < nodes.size(); ++n, ++i )
       nodes[ i ] = *n;
 
-    // avoid elems on degenarate shapes as UV on them can be wrong
+    // avoid elems on degenerate shapes as UV on them can be wrong
     if ( helper.HasDegeneratedEdges() )
     {
       bool isOnDegen = false;
@@ -2946,7 +2960,7 @@ bool SMESH_MesherHelper::IsDistorted2D( SMESH_subMesh*      faceSM,
  * \brief Find out elements orientation on a geometrical face
  * \param theFace - The face correctly oriented in the shape being meshed
  * \retval bool - true if the face normal and the normal of first element
- *                in the correspoding submesh point in different directions
+ *                in the corresponding submesh point in different directions
  */
 //================================================================================
 
@@ -3443,6 +3457,24 @@ double SMESH_MesherHelper::GetOtherParam(const double param) const
   return fabs(param-myPar1[i]) < fabs(param-myPar2[i]) ? myPar2[i] : myPar1[i];
 }
 
+//=======================================================================
+//function : NbRealSeam
+//purpose  : Return a number of real seam edges in the shape set through
+//           IsQuadraticSubMesh() or SetSubShape(). A real seam edge encounters twice in a wire
+//=======================================================================
+
+size_t SMESH_MesherHelper::NbRealSeam() const
+{
+  size_t nb = 0;
+
+  std::set< int >::const_iterator id = mySeamShapeIds.begin();
+  for ( ; id != mySeamShapeIds.end(); ++id )
+    if ( *id < 0 ) ++nb;
+    else break;
+
+  return nb;
+}
+
 //=======================================================================
 //function : IsOnSeam
 //purpose  : Check if UV is on seam. Return 0 if not, 1 for U seam, 2 for V seam
@@ -4108,9 +4140,9 @@ namespace { // Structures used by FixQuadraticElements()
   {
     // code is valid for convex faces only
     gp_XYZ gc(0,0,0);
-    for ( TIDSortedNodeSet::const_iterator n = begin(); n!=end(); ++n)
-      gc += XYZ( *n ) / size();
-    for (unsigned i = 0; i < _sides.size(); ++i )
+    for ( TIDSortedNodeSet::const_iterator n = begin(); n != end(); ++n )
+      gc += XYZ( *n ) / double( size() );
+    for ( size_t i = 0; i < _sides.size(); ++i )
     {
       if ( _sides[i] == bentLink ) continue;
       gp_Vec linkNorm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
@@ -4126,12 +4158,12 @@ namespace { // Structures used by FixQuadraticElements()
         return true;
     }
     return false;
-    
+
   }
 
   //================================================================================
   /*!
-   * \brief Find pairs of continues faces 
+   * \brief Find pairs of continues faces
    */
   //================================================================================
 
@@ -4142,7 +4174,7 @@ namespace { // Structures used by FixQuadraticElements()
     //       |          Between _faces of link x2 two vertical faces are continues
     // x1----x2-----x3  and two horizontal faces are continues. We set vertical faces
     //       |          to _faces[0] and _faces[1] and horizontal faces to
-    //   v2  |   v3     _faces[2] and _faces[3] (or vise versa).
+    //   v2  |   v3     _faces[2] and _faces[3] (or vice versa).
     //       x4
 
     if ( _faces.empty() )
@@ -4343,7 +4375,7 @@ namespace { // Structures used by FixQuadraticElements()
                                              vector< TChain> &   resultChains,
                                              SMDS_TypeOfPosition pos )
   {
-    // put links in the set and evalute number of result chains by number of boundary links
+    // put links in the set and evaluate number of result chains by number of boundary links
     TLinkSet linkSet;
     size_t nbBndLinks = 0;
     for ( TChain::iterator lnk = allLinks.begin(); lnk != allLinks.end(); ++lnk ) {
@@ -4527,8 +4559,16 @@ namespace { // Structures used by FixQuadraticElements()
     TopoDS_Shape  shape = theHelper.GetSubShape().Oriented( TopAbs_FORWARD );
     if ( shape.IsNull() ) return;
 
-    if ( !theError ) theError = SMESH_ComputeError::New();
-
+    if ( !dynamic_cast<SMESH_BadInputElements*>( theError.get() ))
+    {
+      if ( !theError )
+        theError.reset( new SMESH_BadInputElements( meshDS ));
+      else
+        theError.reset( new SMESH_BadInputElements( meshDS,
+                                                    theError->myName,
+                                                    theError->myComment,
+                                                    theError->myAlgo));
+    }
     gp_XYZ faceNorm;
 
     if ( shape.ShapeType() == TopAbs_FACE ) // 2D
@@ -4646,13 +4686,13 @@ namespace { // Structures used by FixQuadraticElements()
                 gp_XYZ pMid3D = 0.5 * ( pN0 + SMESH_TNodeXYZ( nOnEdge[1] ));
                 meshDS->MoveNode( n, pMid3D.X(), pMid3D.Y(), pMid3D.Z() );
                 MSG( "move OUT of face " << n );
-                theError->myBadElements.push_back( f );
+                static_cast<SMESH_BadInputElements*>( theError.get() )->add( f );
               }
             }
           }
         }
       }
-      if ( !theError->myBadElements.empty() )
+      if ( theError->HasBadElems() )
         theError->myName = EDITERR_NO_MEDIUM_ON_GEOM;
       return;
 
@@ -4840,13 +4880,13 @@ namespace { // Structures used by FixQuadraticElements()
                   MSG( "move OUT of solid " << nMedium );
                 }
               }
-              theError->myBadElements.push_back( vol );
+              static_cast<SMESH_BadInputElements*>( theError.get() )->add( vol );
             }
           } // loop on volumes sharing a node on FACE
         } // loop on nodes on FACE
       }  // loop on FACEs of a SOLID
 
-      if ( !theError->myBadElements.empty() )
+      if ( theError->HasBadElems() )
         theError->myName = EDITERR_NO_MEDIUM_ON_GEOM;
     } // 3D case
   }
@@ -4858,7 +4898,7 @@ namespace { // Structures used by FixQuadraticElements()
  * \brief Move medium nodes of faces and volumes to fix distorted elements
  * \param error - container of fixed distorted elements
  * \param volumeOnly - to fix nodes on faces or not, if the shape is solid
- * 
+ *
  * Issue 0020307: EDF 992 SMESH : Linea/Quadratic with Medium Node on Geometry
  */
 //=======================================================================
@@ -4866,6 +4906,7 @@ namespace { // Structures used by FixQuadraticElements()
 void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
                                               bool                   volumeOnly)
 {
+  //MESSAGE("FixQuadraticElements " << volumeOnly);
   // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion
   if ( getenv("NO_FixQuadraticElements") )
     return;
@@ -4907,6 +4948,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
     nbfaces = faces.Extent(); /*avoid "unused varianbles": */ nbfaces++, nbfaces--; 
 #endif
     for ( TopTools_MapIteratorOfMapOfShape fIt( faces ); fIt.More(); fIt.Next() ) {
+      MESSAGE("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
       MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key()));
       SMESH_MesherHelper h(*myMesh);
       h.SetSubShape( fIt.Key() );
@@ -5107,7 +5149,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
             MSG("Internal chain - ignore");
             continue;
           }
-          // mesure chain length and compute link position along the chain
+          // measure chain length and compute link position along the chain
           double chainLen = 0;
           vector< double > linkPos;
           TChain savedChain; // backup
@@ -5231,8 +5273,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
               gp_XY newUV   = ApplyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added );
               gp_Pnt newPnt = s->Value( newUV.X(), newUV.Y());
               move = gp_Vec( XYZ((*link1)->_mediumNode), newPnt.Transformed(loc) );
-              if ( SMDS_FacePosition* nPos =
-                   dynamic_cast< SMDS_FacePosition* >((*link1)->_mediumNode->GetPosition()))
+              if ( SMDS_FacePositionPtr nPos = (*link1)->_mediumNode->GetPosition())
                 nPos->SetParameters( newUV.X(), newUV.Y() );
 #ifdef _DEBUG_
               if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() <
@@ -5263,10 +5304,11 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
   // 4. Move nodes
   // -------------
 
-  TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa;
+  TIDSortedElemSet biQuadQuas, biQuadTris, triQuadHexa, biQuadPenta;
   const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() +
                                    myMesh->NbBiQuadTriangles() +
-                                   myMesh->NbTriQuadraticHexas() );
+                                   myMesh->NbTriQuadraticHexas() +
+                                   myMesh->NbBiQuadPrisms());
   double distXYZ[4];
   faceHlp.ToFixNodeParameters( true );
 
@@ -5302,6 +5344,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
           case SMDSEntity_BiQuad_Quadrangle: biQuadQuas.insert( e ); break;
           case SMDSEntity_BiQuad_Triangle:   biQuadTris.insert( e ); break;
           case SMDSEntity_TriQuad_Hexa:      triQuadHexa.insert( e ); break;
+          case SMDSEntity_BiQuad_Penta:      biQuadPenta.insert( e ); break;
           default:;
           }
         }
@@ -5452,6 +5495,16 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
                              nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z());
     }
   }
+  // treat tri-quadratic hexahedra
+  {
+    SMDS_VolumeTool volExp;
+    TIDSortedElemSet::iterator pentIt = biQuadPenta.begin();
+    for ( ; pentIt != biQuadPenta.end(); ++pentIt )
+    {
+      MESSAGE("---");
+      volExp.Set( *pentIt, /*ignoreCentralNodes=*/false );
+    }
+  }
 #ifdef _DEBUG_
   // avoid warning: defined but not used operator<<()
   SMESH_Comment() << *links.begin() << *faces.begin();