Salome HOME
bos #20256: [CEA 18523] Porting SMESH to int 64 bits
[modules/smesh.git] / src / StdMeshers / StdMeshers_Quadrangle_2D.cxx
index 6173866f4b9ebd7d8c2ea74a803edc9c57494321..914eebff3837e353de40c316cdce4ecd0c76acdb 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2019  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2021  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
@@ -51,6 +51,7 @@
 #include <Geom_Surface.hxx>
 #include <NCollection_DefineArray2.hxx>
 #include <Precision.hxx>
+#include <ShapeAnalysis.hxx>
 #include <TColStd_SequenceOfInteger.hxx>
 #include <TColStd_SequenceOfReal.hxx>
 #include <TColgp_SequenceOfXY.hxx>
@@ -203,12 +204,12 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis
 
 //=============================================================================
 /*!
- *
+ * Compute the mesh on the given shape
  */
 //=============================================================================
 
-bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh&         aMesh,
-                                        const TopoDS_Shape& aShape)
+bool StdMeshers_Quadrangle_2D::ComputeSMESH_Mesh&         aMesh,
+                                        const TopoDS_Shape& aShape )
 {
   const TopoDS_Face& F = TopoDS::Face(aShape);
   aMesh.GetSubMesh( F );
@@ -916,7 +917,7 @@ bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh&         aMesh,
   std::vector<int> aNbNodes(4);
   bool IsQuadratic = false;
   if (!checkNbEdgesForEvaluate(aMesh, aFace, aResMap, aNbNodes, IsQuadratic)) {
-    std::vector<int> aResVec(SMDSEntity_Last);
+    std::vector<smIdType> aResVec(SMDSEntity_Last);
     for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
     SMESH_subMesh * sm = aMesh.GetSubMesh(aFace);
     aResMap.insert(std::make_pair(sm,aResVec));
@@ -965,7 +966,7 @@ bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh&         aMesh,
   //int nbFaces4 = (nbhoriz-1-kdh)*(nbvertic-1-kdv);
   int nbFaces4 = (nbhoriz-1)*(nbvertic-1);
 
-  std::vector<int> aVec(SMDSEntity_Last,0);
+  std::vector<smIdType> aVec(SMDSEntity_Last,0);
   if (IsQuadratic) {
     aVec[SMDSEntity_Quad_Triangle] = nbFaces3;
     aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4;
@@ -1221,7 +1222,9 @@ namespace
    *  \param [in] theNbCorners - the required number of sides, 3 or 4
    *  \param [in] theConsiderMesh - to considered only meshed VERTEXes
    *  \param [in] theFaceSide - the FACE EDGEs
+   *  \param [in] theFixedVertices - VERTEXes to be used as corners
    *  \param [out] theVertices - the found corner vertices
+   *  \param [out] theHaveConcaveVertices - return if there are concave vertices
    */
   //================================================================================
 
@@ -1432,6 +1435,39 @@ namespace
     return;
   }
 
+  //================================================================================
+  /*!
+   * \brief Remove a seam and degenerated edge from a wire if the shape is
+   *        a quadrangle with a seam inside.
+   */
+  //================================================================================
+
+  bool removeInternalSeam( std::list<TopoDS_Edge>& theWire,
+                           SMESH_MesherHelper&     theHelper)
+  {
+    if ( !theHelper.HasRealSeam() ||
+         theHelper.NbDegeneratedEdges() != 2 ) // 1 EDGE + 1 VERTEX
+      return false;
+
+    typedef std::list<TopoDS_Edge>::iterator TEdgeIter;
+    std::vector< TEdgeIter > edgesToRemove;
+    edgesToRemove.reserve( 5 );
+    for ( TEdgeIter eIt = theWire.begin(); eIt != theWire.end(); ++eIt )
+    {
+      int eID = theHelper.ShapeToIndex( *eIt );
+      if ( theHelper.IsRealSeam( eID ) || theHelper.IsDegenShape( eID ))
+        edgesToRemove.push_back( eIt );
+    }
+
+    if ( theWire.size() - edgesToRemove.size() < 4 )
+      return false; // cone e.g.
+
+    for ( size_t i = 0; i < edgesToRemove.size(); ++i )
+      theWire.erase( edgesToRemove[ i ]);
+
+    return true;
+  }
+
 } // namespace
 
 //================================================================================
@@ -1462,6 +1498,9 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
   if ( myHelper )
     helper.CopySubShapeInfo( *myHelper );
 
+  if ( removeInternalSeam( theWire, helper ))
+    theNbDegenEdges = 1;
+
   StdMeshers_FaceSide faceSide( theFace, theWire, &theMesh,
                                 /*isFwd=*/true, /*skipMedium=*/true, &helper );
 
@@ -1492,7 +1531,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
 
   if ( theConsiderMesh )
   {
-    const int nbSegments = Max( faceSide.NbPoints()-1, faceSide.NbSegments() );
+    const smIdType nbSegments = std::max( faceSide.NbPoints()-1, faceSide.NbSegments() );
     if ( nbSegments < nbCorners )
       return error(COMPERR_BAD_INPUT_MESH, TComm("Too few boundary nodes: ") << nbSegments);
   }
@@ -1567,7 +1606,8 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
 
 //=============================================================================
 /*!
- *
+ * Return FaceQuadStruct where sides ordered CCW, top and left sides
+ *        reversed to be co-directed with bottom and right sides
  */
 //=============================================================================
 
@@ -1728,7 +1768,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh&          aMes
   if (anIt==aResMap.end()) {
     return false;
   }
-  std::vector<int> aVec = (*anIt).second;
+  std::vector<smIdType> aVec = (*anIt).second;
   IsQuadratic = (aVec[SMDSEntity_Quad_Edge] > aVec[SMDSEntity_Edge]);
   if (nbEdgesInWire.front() == 3) { // exactly 3 edges
     if (myTriaVertexID>0) {
@@ -1750,7 +1790,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh&          aMes
         SMESH_subMesh * sm = aMesh.GetSubMesh(E1);
         MapShapeNbElemsItr anIt = aResMap.find(sm);
         if (anIt==aResMap.end()) return false;
-        std::vector<int> aVec = (*anIt).second;
+        std::vector<smIdType> aVec = (*anIt).second;
         if (IsQuadratic)
           aNbNodes[0] = (aVec[SMDSEntity_Node]-1)/2 + 2;
         else
@@ -1784,7 +1824,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh&          aMes
       if (anIt==aResMap.end()) {
         return false;
       }
-      std::vector<int> aVec = (*anIt).second;
+      std::vector<smIdType> aVec = (*anIt).second;
       if (IsQuadratic)
         aNbNodes[nbSides] = (aVec[SMDSEntity_Node]-1)/2 + 2;
       else
@@ -1821,7 +1861,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh&          aMes
         if (anIt==aResMap.end()) {
           return false;
         }
-        std::vector<int> aVec = (*anIt).second;
+        std::vector<smIdType> aVec = (*anIt).second;
         if (IsQuadratic)
           aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
         else
@@ -1862,7 +1902,7 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh&          aMes
           if (anIt==aResMap.end()) {
             return false;
           }
-          std::vector<int> aVec = (*anIt).second;
+          std::vector<smIdType> aVec = (*anIt).second;
           if (IsQuadratic)
             aNbNodes[nbSides] += (aVec[SMDSEntity_Node]-1)/2 + 1;
           else
@@ -3139,7 +3179,7 @@ bool StdMeshers_Quadrangle_2D::evaluateQuadPref(SMESH_Mesh &        aMesh,
     nbFaces += (drl+addv)*(nb-1) + (nt-1);
   } // end new version implementation
 
-  std::vector<int> aVec(SMDSEntity_Last);
+  std::vector<smIdType> aVec(SMDSEntity_Last);
   for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
   if (IsQuadratic) {
     aVec[SMDSEntity_Quad_Quadrangle] = nbFaces;
@@ -3169,8 +3209,8 @@ bool StdMeshers_Quadrangle_2D::evaluateQuadPref(SMESH_Mesh &        aMesh,
  */
 //=============================================================================
 
-void StdMeshers_Quadrangle_2D::splitQuadFace(SMESHDS_Mesh *       theMeshDS,
-                                             int                  theFaceID,
+void StdMeshers_Quadrangle_2D::splitQuadFace(SMESHDS_Mesh *       /*theMeshDS*/,
+                                             int                  /*theFaceID*/,
                                              const SMDS_MeshNode* theNode1,
                                              const SMDS_MeshNode* theNode2,
                                              const SMDS_MeshNode* theNode3,
@@ -4736,6 +4776,7 @@ bool StdMeshers_Quadrangle_2D::check()
 
     const SMDS_MeshNode* nInFace = 0;
     if ( myHelper->HasSeam() )
+    {
       for ( int i = 0; i < nbN && !nInFace; ++i )
         if ( !myHelper->IsSeamShape( nn[i]->getshapeId() ))
         {
@@ -4744,6 +4785,33 @@ bool StdMeshers_Quadrangle_2D::check()
           if ( myHelper->IsOnSeam( uv ))
             nInFace = NULL;
         }
+    }
+    if ( myHelper->GetPeriodicIndex() && !nInFace )
+    {
+      for ( int i = 0; i < nbN && !nInFace; ++i )
+        if ( fSubMesh->Contains( nn[i] ))
+          nInFace = nn[i];
+      if ( !nInFace )
+        for ( int i = 0; i < nbN && !nInFace; ++i )
+        {
+          SMDS_ElemIteratorPtr fIt = nn[i]->GetInverseElementIterator( SMDSAbs_Face );
+          while ( fIt->more() && !nInFace )
+          {
+            const SMDS_MeshElement* face = fIt->next();
+            if ( !fSubMesh->Contains( face ))
+              continue;
+            for ( int iN = 0, nN = face->NbCornerNodes(); iN < nN; ++iN )
+            {
+              const SMDS_MeshNode* n = face->GetNode( iN );
+              if ( fSubMesh->Contains( n ))
+              {
+                nInFace = n;
+                break;
+              }
+            }
+          }
+        }
+    }
 
     toCheckUV = true;
     for ( int i = 0; i < nbN; ++i )
@@ -4774,20 +4842,35 @@ bool StdMeshers_Quadrangle_2D::check()
     default:;
     }
 
-    // if ( isBad && myHelper->HasRealSeam() )
-    // {
-    //   // detect a case where a face intersects the seam
-    //   for ( int iPar = 1; iPar < 3; ++iPar )
-    //     if ( iPar & myHelper->GetPeriodicIndex() )
-    //     {
-    //       double min = uv[0].Coord( iPar ), max = uv[0].Coord( iPar );
-    //       for ( int i = 1; i < nbN; ++i )
-    //       {
-    //         min = Min( min, uv[i].Coord( iPar ));
-    //         max = Max( max, uv[i].Coord( iPar ));
-    //       }
-    //     }
-    // }
+    if ( isBad && myHelper->HasRealSeam() )
+    {
+      // fix uv for a case where a face intersects the seam
+      for ( int iPar = 1; iPar < 3; ++iPar )
+        if ( iPar & myHelper->GetPeriodicIndex() )
+        {
+          double max = uv[0].Coord( iPar );
+          for ( int i = 1; i < nbN; ++i )
+            max = Max( max, uv[i].Coord( iPar ));
+
+          for ( int i = 0; i < nbN; ++i )
+          {
+            double par   = uv[i].Coord( iPar );
+            double shift = ShapeAnalysis::AdjustByPeriod( par, max, myHelper->GetPeriod( iPar ));
+            uv[i].SetCoord( iPar, par + shift );
+          }
+        }
+      double sign1 = getArea( uv[0], uv[1], uv[2] );
+      double sign2 = getArea( uv[0], uv[2], uv[3] );
+      if ( sign1 * sign2 < 0 )
+      {
+        sign2 = getArea( uv[1], uv[2], uv[3] );
+        sign1 = getArea( uv[1], uv[3], uv[0] );
+        if ( sign1 * sign2 < 0 )
+          continue; // this should not happen
+      }
+      isBad = ( sign1 * okSign < 0 );
+    }
+
     if ( isBad )
       badFaces.push_back ( f );
   }