Salome HOME
Try to use length of edges for association in case of a 2-edges wire
[modules/smesh.git] / src / StdMeshers / StdMeshers_Regular_1D.cxx
index 8e50b059ba88c416dbc24ed7be91199d7e6b5f79..6ea520001c195e66985a5b9e329d9c9ff0e0c099 100644 (file)
@@ -1,30 +1,30 @@
-//  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
+// 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_Regular_1D.cxx
 //           Moved here from SMESH_Regular_1D.cxx
 //  Author : Paul RASCLE, EDF
 //  Module : SMESH
-
+//
 #include "StdMeshers_Regular_1D.hxx"
 #include "StdMeshers_Distribution.hxx"
 
 #include <TopExp_Explorer.hxx>
 #include <TopoDS.hxx>
 #include <TopoDS_Edge.hxx>
+#include <TopoDS_Vertex.hxx>
 
 #include <string>
+#include <limits>
 
 using namespace std;
 
@@ -117,9 +119,8 @@ bool StdMeshers_Regular_1D::CheckHypothesis
   _hypType = NONE;
   _quadraticMesh = false;
 
-  const bool ignoreAuxiliaryHyps = false;
   const list <const SMESHDS_Hypothesis * > & hyps =
-    GetUsedHypothesis(aMesh, aShape, ignoreAuxiliaryHyps);
+    GetUsedHypothesis(aMesh, aShape, /*ignoreAuxiliaryHyps=*/false);
 
   // find non-auxiliary hypothesis
   const SMESHDS_Hypothesis *theHyp = 0;
@@ -324,6 +325,8 @@ static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
       return false;
     prevU = U;
   }
+  if ( theReverse )
+    theParams.reverse();
   return true;
 }
 
@@ -351,7 +354,7 @@ static void compensateError(double a1, double an,
                             bool              adjustNeighbors2an = false)
 {
   int i, nPar = theParams.size();
-  if ( a1 + an < length && nPar > 1 )
+  if ( a1 + an <= length && nPar > 1 )
   {
     bool reverse = ( U1 > Un );
     GCPnts_AbscissaPoint Discret(C3d, reverse ? an : -an, Un);
@@ -372,14 +375,23 @@ static void compensateError(double a1, double an,
       dUn = Utgt - theParams.back();
     }
 
-    double q  = dUn / ( nPar - 1 );
-    if ( !adjustNeighbors2an ) {
-      for ( itU = theParams.rbegin(), i = 1; i < nPar; itU++, i++ ) {
+    if ( !adjustNeighbors2an )
+    {
+      double q = dUn / ( Utgt - Un ); // (signed) factor of segment length change
+      for ( itU = theParams.rbegin(), i = 1; i < nPar; i++ ) {
+        double prevU = *itU;
         (*itU) += dUn;
-        dUn -= q;
+        ++itU;
+        dUn = q * (*itU - prevU) * (prevU-U1)/(Un-U1);
       }
     }
-    else {
+    else if ( nPar == 1 )
+    {
+      theParams.back() += dUn;
+    }
+    else
+    {
+      double q  = dUn / ( nPar - 1 );
       theParams.back() += dUn;
       double sign = reverse ? -1 : 1;
       double prevU = theParams.back();
@@ -601,12 +613,14 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
   case NB_SEGMENTS: {
 
     double eltSize = 1;
+    int nbSegments;
     if ( _hypType == MAX_LENGTH )
     {
       double nbseg = ceil(theLength / _value[ BEG_LENGTH_IND ]); // integer sup
       if (nbseg <= 0)
         nbseg = 1;                        // degenerated edge
       eltSize = theLength / nbseg;
+      nbSegments = (int) nbseg;
     }
     else if ( _hypType == LOCAL_LENGTH )
     {
@@ -623,7 +637,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
           bool computed = sm->IsMeshComputed();
           if (!computed) {
             if (sm->GetComputeState() == SMESH_subMesh::READY_TO_COMPUTE) {
-              sm->ComputeStateEngine(SMESH_subMesh::COMPUTE);
+              _gen->Compute( theMesh, _mainEdge, /*anUpward=*/true);
               computed = sm->IsMeshComputed();
             }
           }
@@ -647,13 +661,14 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
       if (nbseg <= 0)
         nbseg = 1;                        // degenerated edge
       eltSize = theLength / nbseg;
+      nbSegments = (int) nbseg;
     }
     else
     {
       // Number Of Segments hypothesis
-      int NbSegm = _ivalue[ NB_SEGMENTS_IND ];
-      if ( NbSegm < 1 )  return false;
-      if ( NbSegm == 1 ) return true;
+      nbSegments = _ivalue[ NB_SEGMENTS_IND ];
+      if ( nbSegments < 1 )  return false;
+      if ( nbSegments == 1 ) return true;
 
       switch (_ivalue[ DISTR_TYPE_IND ])
       {
@@ -663,8 +678,8 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
 
           if (fabs(scale - 1.0) < Precision::Confusion()) {
             // special case to avoid division by zero
-            for (int i = 1; i < NbSegm; i++) {
-              double param = f + (l - f) * i / NbSegm;
+            for (int i = 1; i < nbSegments; i++) {
+              double param = f + (l - f) * i / nbSegments;
               theParams.push_back( param );
             }
           } else {
@@ -672,14 +687,22 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
             if ( theReverse )
               scale = 1.0 / scale;
 
-            double alpha = pow(scale, 1.0 / (NbSegm - 1));
-            double factor = (l - f) / (1.0 - pow(alpha, NbSegm));
+            double alpha = pow(scale, 1.0 / (nbSegments - 1));
+            double factor = (l - f) / (1.0 - pow(alpha, nbSegments));
 
-            for (int i = 1; i < NbSegm; i++) {
+            for (int i = 1; i < nbSegments; i++) {
               double param = f + factor * (1.0 - pow(alpha, i));
               theParams.push_back( param );
             }
           }
+          const double lenFactor = theLength/(l-f);
+          list<double>::iterator u = theParams.begin(), uEnd = theParams.end();
+          for ( ; u != uEnd; ++u )
+          {
+            GCPnts_AbscissaPoint Discret( theC3d, ((*u)-f) * lenFactor, f );
+            if ( Discret.IsDone() )
+              *u = Discret.Parameter();
+          }
           return true;
         }
         break;
@@ -700,7 +723,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
         }
         break;
       case StdMeshers_NumberOfSegments::DT_Regular:
-        eltSize = theLength / _ivalue[ NB_SEGMENTS_IND ];
+        eltSize = theLength / nbSegments;
         break;
       default:
         return false;
@@ -710,13 +733,13 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
     if ( !Discret.IsDone() )
       return error( "GCPnts_UniformAbscissa failed");
 
-    int NbPoints = Discret.NbPoints();
-    for ( int i = 2; i < NbPoints; i++ )
+    int NbPoints = Min( Discret.NbPoints(), nbSegments + 1 );
+    for ( int i = 2; i < NbPoints; i++ ) // skip 1st and last points
     {
       double param = Discret.Parameter(i);
       theParams.push_back( param );
     }
-    compensateError( eltSize, eltSize, f, l, theLength, theC3d, theParams ); // for PAL9899
+    compensateError( eltSize, eltSize, f, l, theLength, theC3d, theParams, true ); // for PAL9899
     return true;
   }
 
@@ -727,6 +750,9 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
     double a1 = _value[ BEG_LENGTH_IND ];
     double an = _value[ END_LENGTH_IND ];
     double q  = ( theLength - a1 ) / ( theLength - an );
+    if ( q < theLength/1e6 || 1.01*theLength < a1 + an)
+      return error ( SMESH_Comment("Invalid segment lengths (")<<a1<<" and "<<an<<") "<<
+                     "for an edge of length "<<theLength);
 
     double U1 = theReverse ? l : f;
     double Un = theReverse ? f : l;
@@ -755,9 +781,12 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
 
     double a1 = _value[ BEG_LENGTH_IND ];
     double an = _value[ END_LENGTH_IND ];
+    if ( 1.01*theLength < a1 + an)
+      return error ( SMESH_Comment("Invalid segment lengths (")<<a1<<" and "<<an<<") "<<
+                     "for an edge of length "<<theLength);
 
     double  q = ( an - a1 ) / ( 2 *theLength/( a1 + an ) - 1 );
-    int     n = int( 1 + ( an - a1 ) / q );
+    int n = int(fabs(q) > numeric_limits<double>::min() ? ( 1+( an-a1 )/q ) : ( 1+theLength/a1 ));
 
     double U1 = theReverse ? l : f;
     double Un = theReverse ? f : l;
@@ -803,38 +832,83 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
       }
       if(!IsExist) Params.InsertBefore(j,aPnts[i]);
     }
-    double pf, pl, par2, par1, psize;
-    if (theReverse) {
-      pf = l;
-      pl = f;
-    }
-    else {
-      pf = f;
-      pl = l;
+    double par2, par1, lp;
+    par1 = f;
+    lp = l;
+    double sign = 1.0;
+    if(theReverse) {
+      par1 = l;
+      lp = f;
+      sign = -1.0;
     }
-    psize = pl - pf;
-    par1 = pf;
-    //cout<<"aPnts.size() = "<<aPnts.size()<<"  Params.Length() = "
-    //    <<Params.Length()<<"   nbsegs.size() = "<<nbsegs.size()<<endl;
+    double eltSize, segmentSize = 0.;
+    double currAbscissa = 0;
     for(i=0; i<Params.Length(); i++) {
-      par2 = pf + Params.Value(i+1)*psize;
       int nbseg = ( i > nbsegs.size()-1 ) ? nbsegs[0] : nbsegs[i];
-      double dp = (par2-par1)/nbseg;
-      int j = 1;
-      for(; j<=nbseg; j++) {
-        double param = par1 + dp*j;
-        theParams.push_back( param );
+      segmentSize = Params.Value(i+1)*theLength - currAbscissa;
+      currAbscissa += segmentSize;
+      GCPnts_AbscissaPoint APnt(theC3d, sign*segmentSize, par1);
+      if( !APnt.IsDone() )
+        return error( "GCPnts_AbscissaPoint failed");
+      par2 = APnt.Parameter();
+      eltSize = segmentSize/nbseg;
+      GCPnts_UniformAbscissa Discret(theC3d, eltSize, par1, par2);
+      if(theReverse)
+        Discret.Initialize(theC3d, eltSize, par2, par1);
+      else
+        Discret.Initialize(theC3d, eltSize, par1, par2);
+      if ( !Discret.IsDone() )
+        return error( "GCPnts_UniformAbscissa failed");
+      int NbPoints = Discret.NbPoints();
+      list<double> tmpParams;
+      for(int i=2; i<NbPoints; i++) {
+        double param = Discret.Parameter(i);
+        tmpParams.push_back( param );
       }
+      if (theReverse) {
+        compensateError( eltSize, eltSize, par2, par1, segmentSize, theC3d, tmpParams );
+        tmpParams.reverse();
+      }
+      else {
+        compensateError( eltSize, eltSize, par1, par2, segmentSize, theC3d, tmpParams );
+      }
+      list<double>::iterator itP = tmpParams.begin();
+      for(; itP != tmpParams.end(); itP++) {
+        theParams.push_back( *(itP) );
+      }
+      theParams.push_back( par2 );
+
       par1 = par2;
     }
     // add for last
     int nbseg = ( nbsegs.size() > Params.Length() ) ? nbsegs[Params.Length()] : nbsegs[0];
-    double dp = (pl-par1)/nbseg;
-    int j = 1;
-    for(; j<nbseg; j++) {
-      double param = par1 + dp*j;
-      theParams.push_back( param );
+    segmentSize = theLength - currAbscissa;
+    eltSize = segmentSize/nbseg;
+    GCPnts_UniformAbscissa Discret;
+    if(theReverse)
+      Discret.Initialize(theC3d, eltSize, par1, lp);
+    else
+      Discret.Initialize(theC3d, eltSize, lp, par1);
+    if ( !Discret.IsDone() )
+      return error( "GCPnts_UniformAbscissa failed");
+    int NbPoints = Discret.NbPoints();
+    list<double> tmpParams;
+    for(int i=2; i<NbPoints; i++) {
+      double param = Discret.Parameter(i);
+      tmpParams.push_back( param );
+    }
+    if (theReverse) {
+      compensateError( eltSize, eltSize, lp, par1, segmentSize, theC3d, tmpParams );
+      tmpParams.reverse();
+    }
+    else {
+      compensateError( eltSize, eltSize, par1, lp, segmentSize, theC3d, tmpParams );
+    }
+    list<double>::iterator itP = tmpParams.begin();
+    for(; itP != tmpParams.end(); itP++) {
+      theParams.push_back( *(itP) );
     }
+
     if (theReverse) {
       theParams.reverse(); // NPAL18025
     }
@@ -892,23 +966,43 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & t
   if (!idFirst || !idLast)
     return error( COMPERR_BAD_INPUT_MESH, "No node on vertex");
 
+  // remove elements created by e.g. patern mapping (PAL21999)
+  // CLEAN event is incorrectly ptopagated seemingly due to Propagation hyp
+  // so TEMPORARY solution is to clean the submesh manually
+  //theMesh.GetSubMesh(theShape)->ComputeStateEngine( SMESH_subMesh::CLEAN );
+  if (SMESHDS_SubMesh * subMeshDS = meshDS->MeshElements(theShape))
+  {
+    SMDS_ElemIteratorPtr ite = subMeshDS->GetElements();
+    while (ite->more())
+      meshDS->RemoveFreeElement(ite->next(), subMeshDS);
+    SMDS_NodeIteratorPtr itn = subMeshDS->GetNodes();
+    while (itn->more()) {
+      const SMDS_MeshNode * node = itn->next();
+      if ( node->NbInverseElements() == 0 )
+        meshDS->RemoveFreeNode(node, subMeshDS);
+      else
+        meshDS->RemoveNode(node);
+    }
+  }
+
   if (!Curve.IsNull())
   {
     list< double > params;
     bool reversed = false;
     if ( theMesh.GetShapeToMesh().ShapeType() >= TopAbs_WIRE ) {
+      // if the shape to mesh is WIRE or EDGE
       reversed = ( EE.Orientation() == TopAbs_REVERSED );
     }
     if ( !_mainEdge.IsNull() ) {
+      // take into account reversing the edge the hypothesis is propagated from
       reversed = ( _mainEdge.Orientation() == TopAbs_REVERSED );
+      int mainID = meshDS->ShapeToIndex(_mainEdge);
+      if ( std::find( _revEdgesIDs.begin(), _revEdgesIDs.end(), mainID) != _revEdgesIDs.end())
+        reversed = !reversed;
     }
-    else if ( _revEdgesIDs.size() > 0 ) {
-      for ( int i = 0; i < _revEdgesIDs.size(); i++) {
-        if ( _revEdgesIDs[i] == shapeID ) {
-          reversed = !reversed;
-        }
-      }
-    }
+    // take into account this edge reversing
+    if ( std::find( _revEdgesIDs.begin(), _revEdgesIDs.end(), shapeID) != _revEdgesIDs.end())
+      reversed = !reversed;
 
     BRepAdaptor_Curve C3d( E );
     double length = EdgeLength( E );
@@ -933,7 +1027,6 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & t
       parLast = f;
     }
     */
-
     for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++) {
       double param = *itU;
       gp_Pnt P = Curve->Value(param);
@@ -1032,7 +1125,7 @@ bool StdMeshers_Regular_1D::Evaluate(SMESH_Mesh & theMesh,
   if ( _hypType == NONE )
     return false;
 
-  SMESHDS_Mesh * meshDS = theMesh.GetMeshDS();
+  //SMESHDS_Mesh * meshDS = theMesh.GetMeshDS();
 
   const TopoDS_Edge & EE = TopoDS::Edge(theShape);
   TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
@@ -1047,8 +1140,7 @@ bool StdMeshers_Regular_1D::Evaluate(SMESH_Mesh & theMesh,
   ASSERT(!VFirst.IsNull());
   ASSERT(!VLast.IsNull());
 
-  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 (!Curve.IsNull()) {
     list< double > params;
@@ -1110,10 +1202,9 @@ StdMeshers_Regular_1D::GetUsedHypothesis(SMESH_Mesh &         aMesh,
 
   SMESH_HypoFilter auxiliaryFilter, compatibleFilter;
   auxiliaryFilter.Init( SMESH_HypoFilter::IsAuxiliary() );
-  const bool ignoreAux = true;
-  InitCompatibleHypoFilter( compatibleFilter, ignoreAux );
+  InitCompatibleHypoFilter( compatibleFilter, /*ignoreAux=*/true );
 
-  // get non-auxiliary assigned to aShape
+  // get non-auxiliary assigned directly to aShape
   int nbHyp = aMesh.GetHypotheses( aShape, compatibleFilter, _usedHypList, false );
 
   if (nbHyp == 0 && aShape.ShapeType() == TopAbs_EDGE)