fix bug PAL10974: add IsReversedSubMesh()

src/SMESH/SMESH_Algo.cxx

@@ -31,16 +31,28 @@ using namespace std;
 #include "SMESH_Gen.hxx"
 #include "SMESH_Mesh.hxx"
 #include "SMESH_HypoFilter.hxx"
+#include "SMDS_FacePosition.hxx"
+#include "SMDS_MeshElement.hxx"
+#include "SMDS_MeshNode.hxx"
+#include "SMESHDS_Mesh.hxx"
+#include "SMESHDS_SubMesh.hxx"
 
-#include <GeomAdaptor_Curve.hxx>
 #include <BRep_Tool.hxx>
 #include <GCPnts_AbscissaPoint.hxx>
+#include <GeomAdaptor_Curve.hxx>
+#include <Geom_Surface.hxx>
+#include <TopLoc_Location.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
+#include <TopTools_ListOfShape.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Face.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Pnt2d.hxx>
+#include <gp_Vec.hxx>
 
 #include "utilities.h"
 
 #include <algorithm>
-#include <TopTools_ListOfShape.hxx>
-#include <TopTools_ListIteratorOfListOfShape.hxx>
 
 //=============================================================================
 /*!
@@ -147,3 +159,111 @@ double SMESH_Algo::EdgeLength(const TopoDS_Edge & E)
 	double length = gabs.Length(AdaptCurve, UMin, UMax);
 	return length;
 }
+
+//================================================================================
+/*!
+ * \brief Find out elements orientation on a geometrical face
+  * \param theFace - The face correctly oriented in the shape being meshed
+  * \param theMeshDS - The mesh data structure
+  * \retval bool - true if the face normal and the normal of first element
+  *                in the correspoding submesh point in different directions
+ */
+//================================================================================
+
+bool SMESH_Algo::IsReversedSubMesh (const TopoDS_Face&  theFace,
+                                    SMESHDS_Mesh*       theMeshDS)
+{
+  if ( theFace.IsNull() || !theMeshDS )
+    return false;
+
+  // find out orientation of a meshed face
+  int faceID = theMeshDS->ShapeToIndex( theFace );
+  TopoDS_Shape aMeshedFace = theMeshDS->IndexToShape( faceID );
+  bool isReversed = ( theFace.Orientation() != aMeshedFace.Orientation() );
+
+  const SMESHDS_SubMesh * aSubMeshDSFace = theMeshDS->MeshElements( faceID );
+  if ( !aSubMeshDSFace )
+    return isReversed;
+
+  // find element with node located on face and get its normal
+  const SMDS_FacePosition* facePos = 0;
+  int vertexID = 0;
+  gp_Pnt nPnt[3];
+  gp_Vec Ne;
+  bool normalOK = false;
+  SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements();
+  while ( iteratorElem->more() ) // loop on elements on theFace
+  {
+    const SMDS_MeshElement* elem = iteratorElem->next();
+    if ( elem && elem->NbNodes() > 2 ) {
+      SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator();
+      const SMDS_FacePosition* fPos = 0;
+      int i = 0, vID = 0;
+      while ( nodesIt->more() ) { // loop on nodes
+        const SMDS_MeshNode* node
+          = static_cast<const SMDS_MeshNode *>(nodesIt->next());
+        if ( i == 3 ) i = 2;
+        nPnt[ i++ ].SetCoord( node->X(), node->Y(), node->Z() );
+        // check position
+        const SMDS_PositionPtr& pos = node->GetPosition();
+        if ( !pos ) continue;
+        if ( pos->GetTypeOfPosition() == SMDS_TOP_FACE ) {
+          fPos = dynamic_cast< const SMDS_FacePosition* >( pos.get() );
+        }
+        else if ( pos->GetTypeOfPosition() == SMDS_TOP_VERTEX ) {
+          vID = pos->GetShapeId();
+        }
+      }
+      if ( fPos || ( !normalOK && vID )) {
+        // compute normal
+        gp_Vec v01( nPnt[0], nPnt[1] ), v02( nPnt[0], nPnt[2] );
+        if ( v01.SquareMagnitude() > RealSmall() &&
+             v02.SquareMagnitude() > RealSmall() )
+        {
+          Ne = v01 ^ v02;
+          normalOK = ( Ne.SquareMagnitude() > RealSmall() );
+        }
+        // we need position on theFace or at least on vertex
+        if ( normalOK ) {
+          vertexID = vID;
+          if ((facePos = fPos))
+            break;
+        }
+      }
+    }
+  }
+  if ( !normalOK )
+    return isReversed;
+
+  // node position on face
+  double u,v;
+  if ( facePos ) {
+    u = facePos->GetUParameter();
+    v = facePos->GetVParameter();
+  }
+  else if ( vertexID ) {
+    TopoDS_Shape V = theMeshDS->IndexToShape( vertexID );
+    if ( V.IsNull() || V.ShapeType() != TopAbs_VERTEX )
+      return isReversed;
+    gp_Pnt2d uv = BRep_Tool::Parameters( TopoDS::Vertex( V ), theFace );
+    u = uv.X();
+    v = uv.Y();
+  }
+  else
+  {
+    return isReversed;
+  }
+
+  // face normal at node position
+  TopLoc_Location loc;
+  Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace, loc );
+  if ( surf.IsNull() || surf->Continuity() < GeomAbs_C1 ) return isReversed;
+  gp_Vec d1u, d1v;
+  surf->D1( u, v, nPnt[0], d1u, d1v );
+  gp_Vec Nf = (d1u ^ d1v).Transformed( loc );
+
+  if ( theFace.Orientation() == TopAbs_REVERSED )
+    Nf.Reverse();
+
+  return Ne * Nf < 0.;
+}

src/SMESH/SMESH_Algo.hxx

@@ -40,6 +40,8 @@
 
 class SMESH_Gen;
 class SMESH_Mesh;
+class TopoDS_Face;
+class SMESHDS_Mesh;
 
 class SMESH_Algo:public SMESH_Hypothesis
 {
@@ -62,6 +64,16 @@ class SMESH_Algo:public SMESH_Hypothesis
 
 	static double EdgeLength(const TopoDS_Edge & E);
 
+  /*!
+   * \brief Find out elements orientation on a geometrical face
+   * \param theFace - The face correctly oriented in the shape being meshed
+   * \param theMeshDS - The mesh data structure
+   * \retval bool - true if the face normal and the normal of first element
+   *                in the correspoding submesh point in different directions
+   */
+  static bool IsReversedSubMesh (const TopoDS_Face&  theFace,
+                                 SMESHDS_Mesh*       theMeshDS);
+
  public:
         // algo features