Salome HOME
IPAL54529: Hexahedron(ijk) fails on a block with composite sides if Viscous Layers...
[modules/smesh.git] / src / StdMeshers / StdMeshers_Quadrangle_2D.cxx
index d3d071863c0df229e424762daaab1ccbc735d7b4..6173866f4b9ebd7d8c2ea74a803edc9c57494321 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
@@ -80,9 +80,9 @@ using namespace std;
  */
 //=============================================================================
 
-StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId,
+StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId,
                                                     SMESH_Gen* gen)
-  : SMESH_2D_Algo(hypId, studyId, gen),
+  : SMESH_2D_Algo(hypId, gen),
     myQuadranglePreference(false),
     myTrianglePreference(false),
     myTriaVertexID(-1),
@@ -127,6 +127,7 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis
   myTrianglePreference   = false;
   myHelper               = (SMESH_MesherHelper*)NULL;
   myParams               = NULL;
+  myProxyMesh.reset();
   myQuadList.clear();
 
   aStatus = SMESH_Hypothesis::HYP_OK;
@@ -251,7 +252,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh&         aMesh,
     int nfull = n1+n2+n3+n4;
     if ((nfull % 2) == 0 && ((n1 != n3) || (n2 != n4)))
     {
-      // special path genarating only quandrangle faces
+      // special path generating only quandrangle faces
       res = computeQuadPref( aMesh, F, quad );
     }
   }
@@ -643,9 +644,11 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
               SMESH_ComputeErrorPtr& err = aMesh.GetSubMesh( aFace )->GetComputeError();
               if ( !err || err->IsOK() || err->myName < COMPERR_WARNING )
               {
-                err.reset( new SMESH_ComputeError( COMPERR_WARNING,
-                                                   "Bad quality quad created"));
-                err->myBadElements.push_back( face );
+                SMESH_BadInputElements* badElems =
+                  new SMESH_BadInputElements( meshDS, COMPERR_WARNING,
+                                              "Bad quality quad created");
+                badElems->add( face );
+                err.reset( badElems );
               }
             }
             --i;
@@ -818,9 +821,11 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
               SMESH_ComputeErrorPtr& err = aMesh.GetSubMesh( aFace )->GetComputeError();
               if ( !err || err->IsOK() || err->myName < COMPERR_WARNING )
               {
-                err.reset( new SMESH_ComputeError( COMPERR_WARNING,
-                                                   "Bad quality quad created"));
-                err->myBadElements.push_back( face );
+                SMESH_BadInputElements* badElems =
+                  new SMESH_BadInputElements( meshDS, COMPERR_WARNING,
+                                              "Bad quality quad created");
+                badElems->add( face );
+                err.reset( badElems );
               }
             }
             --i;
@@ -960,8 +965,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);
-  for (int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aVec[i] = 0;
+  std::vector<int> aVec(SMDSEntity_Last,0);
   if (IsQuadratic) {
     aVec[SMDSEntity_Quad_Triangle] = nbFaces3;
     aVec[SMDSEntity_Quad_Quadrangle] = nbFaces4;
@@ -1044,6 +1048,74 @@ namespace
     return true;
   }
 
+  //================================================================================
+  /*!
+   * \brief Return angle between mesh segments of given EDGEs meeting at theVertexNode
+   */
+  //================================================================================
+
+  double getAngleByNodes( const int                  theE1Index,
+                          const int                  theE2Index,
+                          const SMDS_MeshNode*       theVertexNode,
+                          const StdMeshers_FaceSide& theFaceSide,
+                          const gp_Vec&              theFaceNormal)
+  {
+    int eID1 = theFaceSide.EdgeID( theE1Index );
+    int eID2 = theFaceSide.EdgeID( theE2Index );
+
+    const SMDS_MeshNode *n1 = 0, *n2 = 0;
+    bool is1st;
+    SMDS_ElemIteratorPtr segIt = theVertexNode->GetInverseElementIterator( SMDSAbs_Edge );
+    while ( segIt->more() )
+    {
+      const SMDS_MeshElement* seg = segIt->next();
+      int shapeID = seg->GetShapeID();
+      if      ( shapeID == eID1 )
+        is1st = true;
+      else if ( shapeID == eID2 )
+        is1st = false;
+      else
+        continue;
+      ( is1st ? n1 : n2 ) = seg->GetNodeWrap( 1 + seg->GetNodeIndex( theVertexNode ));
+    }
+
+    if ( !n1 || !n2 )
+    {
+      std::vector<const SMDS_MeshNode*> nodes;
+      for ( int is2nd = 0; is2nd < 2; ++is2nd )
+      {
+        const SMDS_MeshNode* & n = is2nd ? n2 : n1;
+        if ( n ) continue;
+        nodes.clear();
+        if ( is2nd ) theFaceSide.GetEdgeNodes( theE2Index, nodes );
+        else         theFaceSide.GetEdgeNodes( theE1Index, nodes );
+        if ( nodes.size() >= 2 )
+        {
+          if ( nodes[0] == theVertexNode )
+            n = nodes[1];
+          else
+            n = nodes[ nodes.size() - 2 ];
+        }
+      }
+    }
+    double angle = -2 * M_PI;
+    if ( n1 && n2 )
+    {
+      SMESH_NodeXYZ p1 = n1, p2 = theVertexNode, p3 = n2;
+      gp_Vec v1( p1, p2 ), v2( p2, p3 );
+      try
+      {
+        angle = v1.AngleWithRef( v2, theFaceNormal );
+      }
+      catch(...)
+      {
+      }
+      if ( std::isnan( angle ))
+        angle = -2 * M_PI;
+    }
+    return angle;
+  }
+
   //--------------------------------------------------------------------------------
   /*!
    * \brief EDGE of a FACE
@@ -1053,6 +1125,7 @@ namespace
     TopoDS_Edge   myEdge;
     TopoDS_Vertex my1stVertex;
     int           myIndex;
+    bool          myIsCorner;   // is fixed corner
     double        myAngle;      // angle at my1stVertex
     int           myNbSegments; // discretization
     Edge*         myPrev;       // preceding EDGE
@@ -1081,6 +1154,7 @@ namespace
 
     // quality criteria to minimize
     int    myOppDiff;
+    int    myIsFixedCorner;
     double myQuartDiff;
     double mySumAngle;
 
@@ -1089,7 +1163,8 @@ namespace
     void AddSelf( QuadQuality::set& theVariants )
     {
       if ( myCornerE[2] == myCornerE[1] || // exclude invalid variants
-           myCornerE[2] == myCornerE[3] )
+           myCornerE[2] == myCornerE[3] ||
+           myCornerE[0] == myCornerE[3] )
         return;
 
       // count nb segments between corners
@@ -1107,9 +1182,14 @@ namespace
       myOppDiff = ( Abs( myNbSeg[0] - myNbSeg[2] ) +
                     Abs( myNbSeg[1] - myNbSeg[3] ));
 
+      myIsFixedCorner = - totNbSeg * ( myCornerE[0]->myIsCorner +
+                                       myCornerE[1]->myIsCorner +
+                                       myCornerE[2]->myIsCorner +
+                                       myCornerE[3]->myIsCorner );
+
       double nbSideIdeal = totNbSeg / 4.;
       myQuartDiff = -( Min( Min( myNbSeg[0], myNbSeg[1] ),
-                            Min( myNbSeg[1], myNbSeg[2] )) / nbSideIdeal );
+                            Min( myNbSeg[2], myNbSeg[3] )) / nbSideIdeal );
 
       theVariants.insert( *this );
 
@@ -1120,7 +1200,7 @@ namespace
     };
 
     // first criterion - equality of nbSeg of opposite sides
-    int    crit1() const { return myOppDiff; }
+    int    crit1() const { return myOppDiff + myIsFixedCorner; }
 
     // second criterion - equality of nbSeg of adjacent sides and sharpness of angles
     double crit2() const { return myQuartDiff + mySumAngle; }
@@ -1138,7 +1218,7 @@ namespace
   //================================================================================
   /*!
    * \brief Unite EDGEs to get a required number of sides
-   *  \param [in] theNbCorners - the required number of sides
+   *  \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 [out] theVertices - the found corner vertices
@@ -1148,11 +1228,10 @@ namespace
   void uniteEdges( const int                   theNbCorners,
                    const bool                  theConsiderMesh,
                    const StdMeshers_FaceSide&  theFaceSide,
-                   const TopoDS_Shape&         theBaseVertex,
-                   std::vector<TopoDS_Vertex>& theVertices )
+                   const TopTools_MapOfShape&  theFixedVertices,
+                   std::vector<TopoDS_Vertex>& theVertices,
+                   bool&                       theHaveConcaveVertices)
   {
-    theVertices.clear();
-
     // form a circular list of EDGEs
     std::vector< Edge > edges( theFaceSide.NbEdges() );
     boost::intrusive::circular_list_algorithms< Edge > circularList;
@@ -1179,33 +1258,59 @@ namespace
 
     // sort edges by angle
     std::multimap< double, Edge* > edgeByAngle;
-    int i, iBase = -1, nbConvexAngles = 0;
+    int i, nbConvexAngles = 0, nbSharpAngles = 0;
+    const SMDS_MeshNode* vertNode = 0;
+    gp_Vec faceNormal;
+    const double angTol     = 5. / 180 * M_PI;
+    const double sharpAngle = 0.5 * M_PI - angTol;
     Edge* e = edge0;
     for ( i = 0; i < nbEdges; ++i, e = e->myNext )
     {
       e->my1stVertex = SMESH_MesherHelper::IthVertex( 0, e->myEdge );
-      if ( e->my1stVertex.IsSame( theBaseVertex ))
-        iBase = e->myIndex;
+      e->myIsCorner = theFixedVertices.Contains( e->my1stVertex );
 
       e->myAngle = -2 * M_PI;
-      if ( !theConsiderMesh || theFaceSide.VertexNode( e->myIndex ))
+      if ( !theConsiderMesh || ( vertNode = theFaceSide.VertexNode( e->myIndex )))
       {
         e->myAngle = SMESH_MesherHelper::GetAngle( e->myPrev->myEdge, e->myEdge,
-                                                   theFaceSide.Face(), e->my1stVertex );
+                                                   theFaceSide.Face(), e->my1stVertex,
+                                                   &faceNormal );
         if ( e->myAngle > 2 * M_PI ) // GetAngle() failed
           e->myAngle *= -1.;
+        else if ( vertNode && ( 0. <= e->myAngle ) && ( e->myAngle <= angTol ))
+          e->myAngle = getAngleByNodes( e->myPrev->myIndex, e->myIndex,
+                                        vertNode, theFaceSide, faceNormal );
       }
       edgeByAngle.insert( std::make_pair( e->myAngle, e ));
-      nbConvexAngles += ( e->myAngle > 0 );
+      nbConvexAngles += ( e->myAngle > angTol );
+      nbSharpAngles  += ( e->myAngle > sharpAngle );
     }
 
-    if ( !theConsiderMesh || theNbCorners < 4 || nbConvexAngles <= theNbCorners )
+    theHaveConcaveVertices = ( nbConvexAngles < nbEdges );
+
+    if ((int) theVertices.size() == theNbCorners )
+      return;
+
+    theVertices.clear();
+
+    if ( !theConsiderMesh || theNbCorners < 4 ||
+         nbConvexAngles <= theNbCorners ||
+         nbSharpAngles  == theNbCorners )
     {
+      if ( nbEdges == theNbCorners ) // return all vertices
+      {
+        for ( e = edge0; (int) theVertices.size() < theNbCorners; e = e->myNext )
+          theVertices.push_back( e->my1stVertex );
+        return;
+      }
+
       // return corners with maximal angles
 
       std::set< int > cornerIndices;
-      if ( iBase != -1 )
-        cornerIndices.insert( iBase );
+      if ( !theFixedVertices.IsEmpty() )
+        for ( i = 0, e = edge0; i < nbEdges; ++i, e = e->myNext )
+          if ( e->myIsCorner )
+            cornerIndices.insert( e->myIndex );
 
       std::multimap< double, Edge* >::reverse_iterator a2e = edgeByAngle.rbegin();
       for (; (int) cornerIndices.size() < theNbCorners; ++a2e )
@@ -1224,9 +1329,16 @@ namespace
     for ( i = 0, e = edge0; i < nbEdges; ++i, e = e->myNext )
     {
       nodes.clear();
-      theFaceSide.GetEdgeNodes( e->myIndex, nodes, /*addVertex=*/false, false );
-      e->myNbSegments += nodes.size() + 1;
-      totNbSeg += nodes.size() + 1;
+      theFaceSide.GetEdgeNodes( e->myIndex, nodes, /*addVertex=*/true, true );
+      if ( nodes.size() == 2 && nodes[0] == nodes[1] ) // all nodes merged
+      {
+        e->myAngle = -1; // to remove
+      }
+      else
+      {
+        e->myNbSegments += nodes.size() - 1;
+        totNbSeg        += nodes.size() - 1;
+      }
 
       // join with the previous edge those edges with concave angles
       if ( e->myAngle <= 0 )
@@ -1376,7 +1488,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
   // check nb of available EDGEs
   if ( faceSide.NbEdges() < nbCorners )
     return error(COMPERR_BAD_SHAPE,
-                 TComm("Face must have 4 sides and not ") << faceSide.NbEdges() );
+                 TComm("Face must have 4 sides but not ") << faceSide.NbEdges() );
 
   if ( theConsiderMesh )
   {
@@ -1389,7 +1501,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
   {
     if ( theVertices.size() < 3 )
       return error(COMPERR_BAD_SHAPE,
-                   TComm("Face must have 3 meshed sides and not ") << theVertices.size() );
+                   TComm("Face must have 3 meshed sides but not ") << theVertices.size() );
   }
   else // triaVertex not defined or invalid
   {
@@ -1408,14 +1520,26 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
     }
     if ( theVertices.size() + theNbDegenEdges < 4 )
       return error(COMPERR_BAD_SHAPE,
-                   TComm("Face must have 4 meshed sides and not ") << theVertices.size() );
+                   TComm("Face must have 4 meshed sides but not ") << theVertices.size() );
   }
 
-  if ((int) theVertices.size() > nbCorners )
+  myCheckOri = false;
+  if ( theVertices.size() > 3 )
   {
-    // there are more EDGEs than required nb of sides;
-    // unite some EDGEs to fix the nb of sides
-    uniteEdges( nbCorners, theConsiderMesh, faceSide, triaVertex, theVertices );
+    TopTools_MapOfShape fixedVertices;
+    if ( !triaVertex.IsNull() )
+      fixedVertices.Add( triaVertex );
+    if ( myParams )
+    {
+      const std::vector< int >& vIDs = myParams->GetCorners();
+      for ( size_t i = 0; i < vIDs.size(); ++i )
+      {
+        const TopoDS_Shape& vertex = helper.GetMeshDS()->IndexToShape( vIDs[ i ]);
+        if ( !vertex.IsNull() )
+          fixedVertices.Add( vertex );
+      }
+    }
+    uniteEdges( nbCorners, theConsiderMesh, faceSide, fixedVertices, theVertices, myCheckOri );
   }
 
   if ( nbCorners == 3 && !triaVertex.IsSame( theVertices[0] ))
@@ -1688,6 +1812,8 @@ bool StdMeshers_Quadrangle_2D::checkNbEdgesForEvaluate(SMESH_Mesh&          aMes
         }
       }
       list<TopoDS_Edge>::iterator ite = sideEdges.begin();
+      if ( nbSides >= (int)aNbNodes.size() )
+        return false;
       aNbNodes[nbSides] = 1;
       for (; ite!=sideEdges.end(); ite++) {
         SMESH_subMesh * sm = aMesh.GetSubMesh(*ite);
@@ -4305,6 +4431,7 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad)
       SMESH_TNodeXYZ a2( quad->UVPt( nbhoriz-1, nbvertic-1 ).node );
       SMESH_TNodeXYZ a3( quad->UVPt( 0,         nbvertic-1 ).node );
 
+      // compute TFI
       for (int i = 1; i < nbhoriz-1; i++)
       {
         SMESH_TNodeXYZ p0( quad->UVPt( i, 0          ).node );
@@ -4316,13 +4443,39 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad)
 
           UVPtStruct& uvp = quad->UVPt( i, j );
 
-          gp_Pnt    p = myHelper->calcTFI(uvp.x,uvp.y, a0,a1,a2,a3, p0,p1,p2,p3);
+          gp_Pnt pnew = myHelper->calcTFI(uvp.x,uvp.y, a0,a1,a2,a3, p0,p1,p2,p3);
+          meshDS->MoveNode( uvp.node, pnew.X(), pnew.Y(), pnew.Z() );
+        }
+      }
+      // project to surface
+      double cellSize;
+      for (int i = 1; i < nbhoriz-1; i++)
+      {
+        for (int j = 1; j < nbvertic-1; j++)
+        {
+          UVPtStruct& uvp = quad->UVPt( i, j );
+          SMESH_NodeXYZ p = uvp.node;
+
+          cellSize = Max( p.SquareDistance( quad->UVPt( i+1, j ).node ),
+                          p.SquareDistance( quad->UVPt( i-1, j ).node ));
+          cellSize = Max( p.SquareDistance( quad->UVPt( i, j+1 ).node ), cellSize );
+          cellSize = Max( p.SquareDistance( quad->UVPt( i, j-1 ).node ), cellSize );
+
           gp_Pnt2d uv = surface->NextValueOfUV( uvp.UV(), p, 10*tol );
           gp_Pnt pnew = surface->Value( uv );
-
-          meshDS->MoveNode( uvp.node, pnew.X(), pnew.Y(), pnew.Z() );
-          uvp.u = uv.X();
-          uvp.v = uv.Y();
+          bool     ok = ( pnew.SquareDistance( p ) < 2 * cellSize );
+          if ( !ok )
+          {
+            uv   = surface->ValueOfUV( p, 10*tol );
+            pnew = surface->Value( uv );
+            ok   = ( pnew.SquareDistance( p ) < 2 * cellSize );
+          }
+          if ( ok )
+          {
+            meshDS->MoveNode( uvp.node, pnew.X(), pnew.Y(), pnew.Z() );
+            uvp.u = uv.X();
+            uvp.v = uv.Y();
+          }
         }
       }
     }
@@ -4643,9 +4796,11 @@ bool StdMeshers_Quadrangle_2D::check()
   {
     SMESH_subMesh* fSM = myHelper->GetMesh()->GetSubMesh( geomFace );
     SMESH_ComputeErrorPtr& err = fSM->GetComputeError();
-    err.reset ( new SMESH_ComputeError( COMPERR_ALGO_FAILED,
-                                        "Inverted elements generated"));
-    err->myBadElements.swap( badFaces );
+    SMESH_BadInputElements* badElems =
+      new SMESH_BadInputElements( meshDS, COMPERR_ALGO_FAILED,
+                                  "Inverted elements generated");
+    badElems->myBadElements.swap( badFaces );
+    err.reset( badElems );
 
     return !isOK;
   }