Salome HOME
0020832: EDF 1359 SMESH : Automatic meshing of boundary layers
authoreap <eap@opencascade.com>
Tue, 18 Jan 2011 12:09:32 +0000 (12:09 +0000)
committereap <eap@opencascade.com>
Tue, 18 Jan 2011 12:09:32 +0000 (12:09 +0000)
   1) make CheckNodeUV() and CheckNodeU() optionally return XYZ of node
      projection to shape
   2) prevent ancestors iterator from returning duplicates
   3) move IsClosedEdge() to SMESH_MesherHelper from StdMeshers_ProjectionUtils
   4) for hexa 3D
   static bool LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
                               const TopoDS_Face& theFace,
                               const TopoDS_Edge& theBaseEdge,
                               SMESHDS_Mesh*      theMesh,
+                              SMESH_ProxyMesh*   theProxyMesh=0);

src/SMESH/SMESH_MesherHelper.cxx
src/SMESH/SMESH_MesherHelper.hxx

index 3c343354857ef31003a273fbba3cc4115a90b77e..66e630f0a70093a4280a8f474d9538749f031b97 100644 (file)
@@ -30,6 +30,7 @@
 #include "SMDS_EdgePosition.hxx"
 #include "SMDS_VolumeTool.hxx"
 #include "SMESH_subMesh.hxx"
+#include "SMESH_ProxyMesh.hxx"
 
 #include <BRepAdaptor_Surface.hxx>
 #include <BRepTools.hxx>
@@ -291,7 +292,7 @@ bool SMESH_MesherHelper::IsMedium(const SMDS_MeshNode*      node,
 //=======================================================================
 
 TopoDS_Shape SMESH_MesherHelper::GetSubShapeByNode(const SMDS_MeshNode* node,
-                                                   SMESHDS_Mesh*        meshDS)
+                                                   const SMESHDS_Mesh*  meshDS)
 {
   int shapeID = node->getshapeId();
   if ( 0 < shapeID && shapeID <= meshDS->MaxShapeIndex() )
@@ -500,7 +501,8 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
                                      const SMDS_MeshNode* n,
                                      gp_XY&               uv,
                                      const double         tol,
-                                     const bool           force) const
+                                     const bool           force,
+                                     double               distXYZ[4]) const
 {
   int shapeID = n->getshapeId();
   if ( force || toCheckPosOnShape( shapeID ))
@@ -511,13 +513,19 @@ 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 );
+    gp_Pnt nodePnt = XYZ( n ), surfPnt(0,0,0);
+    double dist = 0;
     if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
-    if ( Precision::IsInfinite( uv.X() ) ||
-         Precision::IsInfinite( uv.Y() ) ||
-         nodePnt.Distance( surface->Value( uv.X(), uv.Y() )) > toldis )
+    bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() ));
+    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 );
@@ -529,7 +537,14 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
       Quantity_Parameter U,V;
       projector.LowerDistanceParameters(U,V);
       uv.SetCoord( U,V );
-      if ( nodePnt.Distance( surface->Value( U, V )) > toldis )
+      surfPnt = surface->Value( U, V );
+      dist = nodePnt.Distance( surfPnt );
+      if ( distXYZ ) {
+        surfPnt.Transform( loc );
+        distXYZ[0] = dist;
+        distXYZ[1] = surfPnt.X(); distXYZ[2] = surfPnt.Y(); distXYZ[3]=surfPnt.Z();
+      }
+      if ( dist > tol )
       {
         MESSAGE( "SMESH_MesherHelper::CheckNodeUV(), invalid projection" );
         return false;
@@ -689,14 +704,14 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
                                     double&              u,
                                     const double         tol,
                                     const bool           force,
-                                    double*              distance) const
+                                    double               distXYZ[4]) const
 {
   int shapeID = n->getshapeId();
   if ( force || toCheckPosOnShape( shapeID ))
   {
-    double toldis = tol;
-    double tolmin = 1.e-7*myMesh->GetMeshDS()->getMaxDim(); // nodes coordinates are stored in float format
-    if (toldis < tolmin) toldis = tolmin;
+    //double toldis = tol;
+    //double tolmin = 1.e-7*myMesh->GetMeshDS()->getMaxDim(); // nodes coordinates are stored in float format
+    //if (toldis < tolmin) toldis = tolmin;
     // check that u is correct
     TopLoc_Location loc; double f,l;
     Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l );
@@ -712,9 +727,14 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
     {
       gp_Pnt nodePnt = SMESH_MeshEditor::TNodeXYZ( n );
       if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
-      double dist = nodePnt.Distance( curve->Value( u ));
-      if ( distance ) *distance = dist;
-      if ( dist > toldis )
+      gp_Pnt curvPnt = curve->Value( u );
+      double dist = nodePnt.Distance( curvPnt );
+      if ( distXYZ ) {
+        curvPnt.Transform( loc );
+        distXYZ[0] = dist;
+        distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
+      }
+      if ( dist > tol /*toldis*/ )
       {
         setPosOnShapeValidity( shapeID, false );
         // u incorrect, project the node to the curve
@@ -736,12 +756,17 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
         }
         Quantity_Parameter U = projector->LowerDistanceParameter();
         u = double( U );
-        dist = nodePnt.Distance( curve->Value( U ));
-        if ( distance ) *distance = dist;
-        if ( dist > toldis )
+        curvPnt = curve->Value( u );
+        dist = nodePnt.Distance( curvPnt );
+        if ( distXYZ ) {
+          curvPnt.Transform( loc );
+          distXYZ[0] = dist;
+          distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
+        }
+        if ( dist > tol /*toldis*/)
         {
           MESSAGE( "SMESH_MesherHelper::CheckNodeU(), invalid projection" );
-          MESSAGE("distance " << nodePnt.Distance(curve->Value( U )) << " " << toldis);
+          MESSAGE("distance " << dist << " " << tol/*dis*/);
           return false;
         }
         // store the fixed U on the edge
@@ -948,7 +973,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_
   // To find position on edge and 3D position for n12,
   // project <middle> to 2 edges and select projection most close to <middle>
 
-  double u = 0, distMiddleProj = Precision::Infinite();
+  double u = 0, distMiddleProj = Precision::Infinite(), distXYZ[4];
   int iOkEdge = 0;
   TopoDS_Edge edges[2];
   for ( int is2nd = 0; is2nd < 2; ++is2nd )
@@ -962,11 +987,11 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_
     // project to get U of projection and distance from middle to projection
     TopoDS_Edge edge = edges[ is2nd ] = TopoDS::Edge( shape );
     double node2MiddleDist = middle.Distance( XYZ(n) );
-    double foundU = GetNodeU( edge, n ), foundDist = node2MiddleDist;
-    CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, &foundDist );
-    if ( foundDist < node2MiddleDist )
+    double foundU = GetNodeU( edge, n );
+    CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, distXYZ );
+    if ( distXYZ[0] < node2MiddleDist )
     {
-      distMiddleProj = foundDist;
+      distMiddleProj = distXYZ[0];
       u = foundU;
       iOkEdge = is2nd;
     }
@@ -1452,9 +1477,24 @@ SMESH_MesherHelper::AddPolyhedralVolume (const std::vector<const SMDS_MeshNode*>
 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
                                          const TopoDS_Face& theFace,
                                          const TopoDS_Edge& theBaseEdge,
-                                         SMESHDS_Mesh*      theMesh)
+                                         SMESHDS_Mesh*      theMesh,
+                                         SMESH_ProxyMesh*   theProxyMesh)
 {
-  SMESHDS_SubMesh* faceSubMesh = theMesh->MeshElements( theFace );
+  const SMESHDS_SubMesh* faceSubMesh = 0;
+  if ( theProxyMesh )
+  {
+    faceSubMesh = theProxyMesh->GetSubMesh( theFace );
+    if ( !faceSubMesh ||
+         faceSubMesh->NbElements() == 0 ||
+         theProxyMesh->IsTemporary( faceSubMesh->GetElements()->next() ))
+    {
+      // can use a proxy sub-mesh with not temporary elements only
+      faceSubMesh = 0;
+      theProxyMesh = 0;
+    }
+  }
+  if ( !faceSubMesh )
+    faceSubMesh = theMesh->MeshElements( theFace );
   if ( !faceSubMesh || faceSubMesh->NbElements() == 0 )
     return false;
 
@@ -1475,12 +1515,21 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
     nCol.resize( nbRows );
     nCol[0] = 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 );
+    }
+  }
 
   // fill theParam2ColumnMap column by column by passing from nodes on
   // theBaseEdge up via mesh faces on theFace
 
-  TParam2ColumnMap::iterator par_nVec_2 = theParam2ColumnMap.begin();
-  TParam2ColumnMap::iterator par_nVec_1 = par_nVec_2++;
+  par_nVec_2 = theParam2ColumnMap.begin();
+  par_nVec_1 = par_nVec_2++;
   TIDSortedElemSet emptySet, avoidSet;
   for ( ; par_nVec_2 != theParam2ColumnMap.end(); ++par_nVec_1, ++par_nVec_2 )
   {
@@ -1511,7 +1560,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
     if ( iRow + 1 < nbRows ) // compact if necessary
       nCol1.resize( iRow + 1 ), nCol2.resize( iRow + 1 );
   }
-  return true;
+  return theParam2ColumnMap.size() > 1 && theParam2ColumnMap.begin()->second.size() > 1;
 }
 
 //=======================================================================
@@ -1611,6 +1660,21 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
   return tol;
 }
 
+//================================================================================
+/*!
+ * \brief Check if the first and last vertices of an edge are the same
+ * \param anEdge - the edge to check
+ * \retval bool - true if same
+ */
+//================================================================================
+
+bool SMESH_MesherHelper::IsClosedEdge( const TopoDS_Edge& anEdge )
+{
+  if ( anEdge.Orientation() >= TopAbs_INTERNAL )
+    return IsClosedEdge( TopoDS::Edge( anEdge.Oriented( TopAbs_FORWARD )));
+  return TopExp::FirstVertex( anEdge ).IsSame( TopExp::LastVertex( anEdge ));
+}
+
 //=======================================================================
 //function : IsQuadraticMesh
 //purpose  : Check mesh without geometry for: if all elements on this shape are quadratic,
@@ -2932,10 +2996,14 @@ struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
 {
   TopTools_ListIteratorOfListOfShape _ancIter;
   TopAbs_ShapeEnum                   _type;
+  TopTools_MapOfShape                _encountered;
   TAncestorsIterator( const TopTools_ListOfShape& ancestors, TopAbs_ShapeEnum type)
     : _ancIter( ancestors ), _type( type )
   {
-    if ( _ancIter.More() && _ancIter.Value().ShapeType() != _type ) next();
+    if ( _ancIter.More() ) {
+      if ( _ancIter.Value().ShapeType() != _type ) next();
+      else _encountered.Add( _ancIter.Value() );
+    }
   }
   virtual bool more()
   {
@@ -2946,7 +3014,7 @@ struct TAncestorsIterator : public SMDS_Iterator<const TopoDS_Shape*>
     const TopoDS_Shape* s = _ancIter.More() ? & _ancIter.Value() : 0;
     if ( _ancIter.More() )
       for ( _ancIter.Next();  _ancIter.More(); _ancIter.Next())
-        if ( _ancIter.Value().ShapeType() == _type )
+        if ( _ancIter.Value().ShapeType() == _type && _encountered.Add( _ancIter.Value() ))
           break;
     return s;
   }
index d904f95ded3491c0ec6ef67954b8cadf812ed355..7a53c71eb064edc04ee9757588d08618d5095189 100644 (file)
@@ -43,6 +43,7 @@
 
 class GeomAPI_ProjectPointOnSurf;
 class GeomAPI_ProjectPointOnCurve;
+class SMESH_ProxyMesh;
 
 typedef std::map<SMESH_TLink, const SMDS_MeshNode*>           TLinkNodeMap;
 typedef std::map<SMESH_TLink, const SMDS_MeshNode*>::iterator ItTLinkNode;
@@ -96,7 +97,8 @@ public:
   static bool LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
                               const TopoDS_Face& theFace,
                               const TopoDS_Edge& theBaseEdge,
-                              SMESHDS_Mesh*      theMesh);
+                              SMESHDS_Mesh*      theMesh,
+                              SMESH_ProxyMesh*   theProxyMesh=0);
   /*!
    * \brief Return support shape of a node
    * \param node - the node
@@ -104,7 +106,7 @@ public:
    * \retval TopoDS_Shape - found support shape
    */
   static TopoDS_Shape GetSubShapeByNode(const SMDS_MeshNode* node,
-                                        SMESHDS_Mesh*        meshDS);
+                                        const SMESHDS_Mesh*  meshDS);
 
   /*!
    * \brief Return a valid node index, fixing the given one if necessary
@@ -143,6 +145,8 @@ public:
 
   static double MaxTolerance( const TopoDS_Shape& shape );
 
+  static bool IsClosedEdge( const TopoDS_Edge& anEdge );
+
 
 public:
   // ---------- PUBLIC INSTANCE METHODS ----------
@@ -301,16 +305,19 @@ public:
   /*!
    * \brief Check and fix node UV on a face
    *  \param force - check even if checks of other nodes on this face passed OK
+   *  \param distXYZ - returns result distance and point coordinates
    *  \retval bool - false if UV is bad and could not be fixed
    */
   bool CheckNodeUV(const TopoDS_Face&   F,
                    const SMDS_MeshNode* n,
                    gp_XY&               uv,
                    const double         tol,
-                   const bool           force=false) const;
+                   const bool           force=false,
+                   double               distXYZ[4]=0) const;
   /*!
    * \brief Check and fix node U on an edge
    *  \param force - check even if checks of other nodes on this edge passed OK
+   *  \param distXYZ - returns result distance and point coordinates
    *  \retval bool - false if U is bad and could not be fixed
    */
   bool CheckNodeU(const TopoDS_Edge&   E,
@@ -318,7 +325,7 @@ public:
                   double&              u,
                   const double         tol,
                   const bool           force=false,
-                  double*              distance=0) const;
+                  double               distXYZ[4]=0) const;
   /*!
    * \brief Return middle UV taking in account surface period
    */