Salome HOME
Merge from V6_main 12/11/2012
[modules/smesh.git] / src / StdMeshers / StdMeshers_RadialQuadrangle_1D2D.cxx
index a94600a30891408b3617cc088db003668cea7ce0..3aeee6062afbbcd9380a9e6d52b016dc91d0f519 100644 (file)
@@ -1,30 +1,26 @@
-//  Copyright (C) 2007-2008  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012  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
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
 //
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
-//
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
+
 //  SMESH SMESH : implementaion of SMESH idl descriptions
 // File      : StdMeshers_RadialQuadrangle_1D2D.cxx
 // Module    : SMESH
-// Created   : Fri Oct 20 11:37:07 2006
-// Author    : Edward AGAPOV (eap)
-//
+
 #include "StdMeshers_RadialQuadrangle_1D2D.hxx"
 
 #include "StdMeshers_NumberOfLayers.hxx"
@@ -52,6 +48,7 @@
 #include <Geom_TrimmedCurve.hxx>
 #include <TColgp_SequenceOfPnt.hxx>
 #include <TColgp_SequenceOfPnt2d.hxx>
+#include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopoDS.hxx>
@@ -80,7 +77,7 @@ StdMeshers_RadialQuadrangle_1D2D::StdMeshers_RadialQuadrangle_1D2D(int hypId,
   _compatibleHypothesis.push_back("NumberOfLayers2D");
   myNbLayerHypo = 0;
   myDistributionHypo = 0;
-  _requireDescretBoundary = false;
+  _requireDiscreteBoundary = false;
   _supportSubmeshes = true;
 }
 
@@ -148,7 +145,8 @@ namespace
    */
   class TEdgeMarker : public SMESH_subMeshEventListener
   {
-    TEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false) {}
+    TEdgeMarker(): SMESH_subMeshEventListener(/*isDeletable=*/false,
+                                              "StdMeshers_RadialQuadrangle_1D2D::TEdgeMarker") {}
   public:
     //!<  Return static listener
     static SMESH_subMeshEventListener* getListener()
@@ -428,7 +426,6 @@ void StdMeshers_RadialQuadrangle_1D2D::SubmeshRestored(SMESH_subMesh* faceSubMes
 bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
                                                const TopoDS_Shape& aShape)
 {
-  TopExp_Explorer exp;
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 
   myHelper = new SMESH_MesherHelper( aMesh );
@@ -480,8 +477,8 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
       for(; itn != theNodes.end(); itn++ ) {
         CNodes.push_back( (*itn).second );
         double ang = (*itn).first - fang;
-        if( ang>PI ) ang = ang - 2*PI;
-        if( ang<-PI ) ang = ang + 2*PI;
+        if( ang>M_PI ) ang = ang - 2.*M_PI;
+        if( ang<-M_PI ) ang = ang + 2.*M_PI;
         Angles.Append( ang ); 
       }
     }
@@ -491,8 +488,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
     if ( !computeLayerPositions(P0,P1))
       return false;
 
-    exp.Init( CircEdge, TopAbs_VERTEX );
-    TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
+    TopoDS_Vertex V1 = myHelper->IthVertex(0, CircEdge );
     gp_Pnt2d p2dV = BRep_Tool::Parameters( V1, TopoDS::Face(aShape) );
 
     NC = meshDS->AddNode(P0.X(), P0.Y(), P0.Z());
@@ -529,7 +525,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
     // a segment of line
     double fp, lp;
     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
-    if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
+    if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
       // not half of circle
       return error(COMPERR_BAD_SHAPE);
     }
@@ -544,23 +540,21 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
     map< double, const SMDS_MeshNode* > theNodes;
     if ( !GetSortedNodesOnEdge(aMesh.GetMeshDS(),CircEdge,true,theNodes) )
       return error("Circular edge is incorrectly meshed");
-    if (theNodes.size()%2 == 0 )
-      return error("Circular edge is incorrectly meshed, number of segments must be even");
 
-    CNodes.clear();
     map< double, const SMDS_MeshNode* >::iterator itn = theNodes.begin();
+    CNodes.clear();
+    CNodes.push_back( itn->second );
     double fang = (*itn).first;
     itn++;
     for(; itn != theNodes.end(); itn++ ) {
       CNodes.push_back( (*itn).second );
       double ang = (*itn).first - fang;
-      if( ang>PI ) ang = ang - 2*PI;
-      if( ang<-PI ) ang = ang + 2*PI;
+      if( ang>M_PI ) ang = ang - 2.*M_PI;
+      if( ang<-M_PI ) ang = ang + 2.*M_PI;
       Angles.Append( ang );
     }
     const SMDS_MeshNode* NF = theNodes.begin()->second;
     const SMDS_MeshNode* NL = theNodes.rbegin()->second;
-    CNodes.push_back( NF );
     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
     P0 = aCirc->Location();
@@ -574,6 +568,8 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
       if (!GetSortedNodesOnEdge(aMesh.GetMeshDS(),LinEdge1,true,theNodes))
         return error("Invalid mesh on a straight edge");
 
+      Nodes1.resize( myLayerPositions.size()+1 );
+      Nodes2.resize( myLayerPositions.size()+1 );
       vector< const SMDS_MeshNode* > *pNodes1 = &Nodes1, *pNodes2 = &Nodes2;
       bool nodesFromP0ToP1 = ( theNodes.rbegin()->second == NF );
       if ( !nodesFromP0ToP1 ) std::swap( pNodes1, pNodes2 );
@@ -584,8 +580,8 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
       {
         (*pNodes1)[i] = ritn->second;
         (*pNodes2)[i] =  itn->second;
-        Points.Append( gpXYZ( Nodes1[i]));
-        Pnts2d1.Append( myHelper->GetNodeUV( F, Nodes1[i]));
+        Points.Prepend( gpXYZ( Nodes1[i]));
+        Pnts2d1.Prepend( myHelper->GetNodeUV( F, Nodes1[i]));
       }
       NC = const_cast<SMDS_MeshNode*>( itn->second );
       Points.Remove( Nodes1.size() );
@@ -691,14 +687,23 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
     for(; itn != theNodes.end(); itn++ ) {
       CNodes.push_back( (*itn).second );
       double ang = (*itn).first - fang;
-      if( ang>PI ) ang = ang - 2*PI;
-      if( ang<-PI ) ang = ang + 2*PI;
+      if( ang>M_PI ) ang = ang - 2.*M_PI;
+      if( ang<-M_PI ) ang = ang + 2.*M_PI;
       Angles.Append( ang );
     }
     P1 = gp_Pnt( NF->X(), NF->Y(), NF->Z() );
     gp_Pnt P2( NL->X(), NL->Y(), NL->Z() );
     P0 = aCirc->Location();
 
+    // make P1 belong to LinEdge1
+    TopoDS_Vertex V1 = myHelper->IthVertex( 0, LinEdge1 );
+    TopoDS_Vertex V2 = myHelper->IthVertex( 1, LinEdge1 );
+    gp_Pnt PE1 = BRep_Tool::Pnt(V1);
+    gp_Pnt PE2 = BRep_Tool::Pnt(V2);
+    if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
+        ( P1.Distance(PE2) > Precision::Confusion() ) )
+      std::swap( LinEdge1, LinEdge2 );
+
     bool linEdge1Computed, linEdge2Computed;
     if ( !computeLayerPositions(P0,P1,LinEdge1,&linEdge1Computed))
       return false;
@@ -707,23 +712,12 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
     Nodes2.resize( myLayerPositions.size()+1 );
 
     // check that both linear edges have same hypotheses
-    if ( !computeLayerPositions(P0,P1,LinEdge2, &linEdge2Computed))
+    if ( !computeLayerPositions(P0,P2,LinEdge2, &linEdge2Computed))
          return false;
     if ( Nodes1.size() != myLayerPositions.size()+1 )
       return error("Different hypotheses apply to radial edges");
 
-    exp.Init( LinEdge1, TopAbs_VERTEX );
-    TopoDS_Vertex V1 = TopoDS::Vertex( exp.Current() );
-    exp.Next();
-    TopoDS_Vertex V2 = TopoDS::Vertex( exp.Current() );
-    gp_Pnt PE1 = BRep_Tool::Pnt(V1);
-    gp_Pnt PE2 = BRep_Tool::Pnt(V2);
-    if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
-        ( P1.Distance(PE2) > Precision::Confusion() ) )
-    {
-      std::swap( LinEdge1, LinEdge2 );
-      std::swap( linEdge1Computed, linEdge2Computed );
-    }
+    // find the central vertex
     TopoDS_Vertex VC = V2;
     if( ( P1.Distance(PE1) > Precision::Confusion() ) &&
         ( P2.Distance(PE1) > Precision::Confusion() ) )
@@ -882,6 +876,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
 
   // orientation
   bool IsForward = ( CircEdge.Orientation()==TopAbs_FORWARD );
+  const double angleSign = ( F.Orientation() == TopAbs_REVERSED ? -1.0 : 1.0 );
 
   // create nodes and mesh elements on face
   // find axis of rotation
@@ -900,7 +895,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Compute(SMESH_Mesh&         aMesh,
     gp_Ax1 theAxis(P0,gp_Dir(Axis));
     aTrsf.SetRotation( theAxis, Angles.Value(i) );
     gp_Trsf2d aTrsf2d;
-    aTrsf2d.SetRotation( PC, Angles.Value(i) );
+    aTrsf2d.SetRotation( PC, Angles.Value(i) * angleSign );
     // create nodes
     int j = 1;
     for(; j<=Points.Length(); j++) {
@@ -1001,13 +996,14 @@ bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt&
   }
   if ( hyp1D ) // try to compute with hyp1D
   {
-    if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions,p1,p2,*mesh,hyp1D ))
+    if ( !TNodeDistributor::GetDistributor(*mesh)->Compute( myLayerPositions,p1,p2,*mesh,hyp1D )) {
       if ( myDistributionHypo ) { // bad hyp assigned 
         return error( TNodeDistributor::GetDistributor(*mesh)->GetComputeError() );
       }
       else {
         // bad hyp found, its Ok, lets try with default nb of segnents
       }
+    }
   }
   
   if ( myLayerPositions.empty() ) // try to use nb of layers
@@ -1034,9 +1030,25 @@ bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt&
     vector< double > nodeParams;
     GetNodeParamOnEdge( mesh->GetMeshDS(), linEdge, nodeParams );
 
+    // nb of present nodes must be different in cases of 1 and 2 straight edges
+
+    TopoDS_Vertex VV[2];
+    TopExp::Vertices( linEdge, VV[0], VV[1]);
+    const gp_Pnt* points[] = { &p1, &p2 };
+    gp_Pnt       vPoints[] = { BRep_Tool::Pnt(VV[0]), BRep_Tool::Pnt(VV[1]) };
+    const double     tol[] = { BRep_Tool::Tolerance(VV[0]), BRep_Tool::Tolerance(VV[1]) };
+    bool pointsAreOnVertices = true;
+    for ( int iP = 0; iP < 2 && pointsAreOnVertices; ++iP )
+      pointsAreOnVertices = ( points[iP]->Distance( vPoints[0] ) < tol[0] ||
+                              points[iP]->Distance( vPoints[1] ) < tol[1] );
+
+    int nbNodes = nodeParams.size() - 2; // 2 straight edges
+    if ( !pointsAreOnVertices )
+      nbNodes = ( nodeParams.size() - 3 ) / 2; // 1 straight edge
+
     if ( myLayerPositions.empty() )
     {
-      myLayerPositions.resize( nodeParams.size() - 2 );
+      myLayerPositions.resize( nbNodes );
     }
     else if ( myDistributionHypo || myNbLayerHypo )
     {
@@ -1052,8 +1064,10 @@ bool StdMeshers_RadialQuadrangle_1D2D::computeLayerPositions(const gp_Pnt&
         mesh->GetSubMesh( linEdge )->ComputeStateEngine( SMESH_subMesh::CLEAN );
         if ( linEdgeComputed ) *linEdgeComputed = false;
       }
-      else if ( myLayerPositions.size() != nodeParams.size()-2 ) {
-        return error("Radial edge is meshed by other algorithm");
+      else {
+        
+        if ( myLayerPositions.size() != nbNodes )
+          return error("Radial edge is meshed by other algorithm");
       }
     }
   }
@@ -1130,7 +1144,7 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
     // a segment of line
     double fp, lp;
     Handle(Geom_Circle) aCirc = Handle(Geom_Circle)::DownCast( getCurve( CircEdge, &fp, &lp ));
-    if( fabs(fabs(lp-fp)-PI) > Precision::Confusion() ) {
+    if( fabs(fabs(lp-fp)-M_PI) > Precision::Confusion() ) {
       // not half of circle
       return error(COMPERR_BAD_SHAPE);
     }
@@ -1268,4 +1282,3 @@ bool StdMeshers_RadialQuadrangle_1D2D::Evaluate(SMESH_Mesh& aMesh,
   return false;
 
 }
-