Salome HOME
Merge from V6_main_20120808 08Aug12
[modules/smesh.git] / src / StdMeshers / StdMeshers_Regular_1D.cxx
index ee4aa86e78a1e992e97c3c2a599a9ef13b4105d3..8d900955d284778f5dadb9fddb78416b11f5311e 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;
 
@@ -73,22 +75,24 @@ using namespace std;
 //=============================================================================
 
 StdMeshers_Regular_1D::StdMeshers_Regular_1D(int hypId, int studyId,
-       SMESH_Gen * gen):SMESH_1D_Algo(hypId, studyId, gen)
+        SMESH_Gen * gen):SMESH_1D_Algo(hypId, studyId, gen)
 {
-       MESSAGE("StdMeshers_Regular_1D::StdMeshers_Regular_1D");
-       _name = "Regular_1D";
-       _shapeType = (1 << TopAbs_EDGE);
-
-       _compatibleHypothesis.push_back("LocalLength");
-       _compatibleHypothesis.push_back("MaxLength");
-       _compatibleHypothesis.push_back("NumberOfSegments");
-       _compatibleHypothesis.push_back("StartEndLength");
-       _compatibleHypothesis.push_back("Deflection1D");
-       _compatibleHypothesis.push_back("Arithmetic1D");
-       _compatibleHypothesis.push_back("AutomaticLength");
-
-       _compatibleHypothesis.push_back("QuadraticMesh"); // auxiliary !!!
-       _compatibleHypothesis.push_back("Propagation"); // auxiliary !!!
+        MESSAGE("StdMeshers_Regular_1D::StdMeshers_Regular_1D");
+        _name = "Regular_1D";
+        _shapeType = (1 << TopAbs_EDGE);
+        _fpHyp = 0;
+
+        _compatibleHypothesis.push_back("LocalLength");
+        _compatibleHypothesis.push_back("MaxLength");
+        _compatibleHypothesis.push_back("NumberOfSegments");
+        _compatibleHypothesis.push_back("StartEndLength");
+        _compatibleHypothesis.push_back("Deflection1D");
+        _compatibleHypothesis.push_back("Arithmetic1D");
+        _compatibleHypothesis.push_back("FixedPoints1D");
+        _compatibleHypothesis.push_back("AutomaticLength");
+
+        _compatibleHypothesis.push_back("QuadraticMesh"); // auxiliary !!!
+        _compatibleHypothesis.push_back("Propagation"); // auxiliary !!!
 }
 
 //=============================================================================
@@ -115,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;
@@ -180,12 +183,15 @@ bool StdMeshers_Regular_1D::CheckHypothesis
     {
     case StdMeshers_NumberOfSegments::DT_Scale:
       _value[ SCALE_FACTOR_IND ] = hyp->GetScaleFactor();
+      _revEdgesIDs = hyp->GetReversedEdges();
       break;
     case StdMeshers_NumberOfSegments::DT_TabFunc:
       _vvalue[ TAB_FUNC_IND ] = hyp->GetTableFunction();
+      _revEdgesIDs = hyp->GetReversedEdges();
       break;
     case StdMeshers_NumberOfSegments::DT_ExprFunc:
       _svalue[ EXPR_FUNC_IND ] = hyp->GetExpressionFunction();
+      _revEdgesIDs = hyp->GetReversedEdges();
       break;
     case StdMeshers_NumberOfSegments::DT_Regular:
       break;
@@ -209,6 +215,19 @@ bool StdMeshers_Regular_1D::CheckHypothesis
     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
     _hypType = ARITHMETIC_1D;
+
+    _revEdgesIDs = hyp->GetReversedEdges();
+
+    aStatus = SMESH_Hypothesis::HYP_OK;
+  }
+
+  else if (hypName == "FixedPoints1D") {
+    _fpHyp = dynamic_cast <const StdMeshers_FixedPoints1D*>(theHyp);
+    ASSERT(_fpHyp);
+    _hypType = FIXED_POINTS_1D;
+
+    _revEdgesIDs = _fpHyp->GetReversedEdges();
+
     aStatus = SMESH_Hypothesis::HYP_OK;
   }
 
@@ -221,6 +240,9 @@ bool StdMeshers_Regular_1D::CheckHypothesis
     _value[ END_LENGTH_IND ] = hyp->GetLength( false );
     ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
     _hypType = BEG_END_LENGTH;
+
+    _revEdgesIDs = hyp->GetReversedEdges();
+
     aStatus = SMESH_Hypothesis::HYP_OK;
   }
 
@@ -303,6 +325,8 @@ static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
       return false;
     prevU = U;
   }
+  if ( theReverse )
+    theParams.reverse();
   return true;
 }
 
@@ -352,10 +376,14 @@ static void compensateError(double a1, double an,
     }
 
     double q  = dUn / ( nPar - 1 );
-    if ( !adjustNeighbors2an ) {
-      for ( itU = theParams.rbegin(), i = 1; i < nPar; itU++, i++ ) {
+    if ( !adjustNeighbors2an )
+    {
+      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 {
@@ -602,7 +630,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();
             }
           }
@@ -659,6 +687,14 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
               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;
@@ -706,6 +742,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;
@@ -734,9 +773,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;
@@ -764,6 +806,107 @@ bool StdMeshers_Regular_1D::computeInternalParameters(SMESH_Mesh &     theMesh,
     return true;
   }
 
+  case FIXED_POINTS_1D: {
+    const std::vector<double>& aPnts = _fpHyp->GetPoints();
+    const std::vector<int>& nbsegs = _fpHyp->GetNbSegments();
+    int i = 0;
+    TColStd_SequenceOfReal Params;
+    for(; i<aPnts.size(); i++) {
+      if( aPnts[i]<0.0001 || aPnts[i]>0.9999 ) continue;
+      int j=1;
+      bool IsExist = false;
+      for(; j<=Params.Length(); j++) {
+        if( fabs(aPnts[i]-Params.Value(j)) < 1e-4 ) {
+          IsExist = true;
+          break;
+        }
+        if( aPnts[i]<Params.Value(j) ) break;
+      }
+      if(!IsExist) Params.InsertBefore(j,aPnts[i]);
+    }
+    double par2, par1, lp;
+    par1 = f;
+    lp = l;
+    double sign = 1.0;
+    if(theReverse) {
+      par1 = l;
+      lp = f;
+      sign = -1.0;
+    }
+    double eltSize, segmentSize = 0.;
+    double currAbscissa = 0;
+    for(i=0; i<Params.Length(); i++) {
+      int nbseg = ( i > nbsegs.size()-1 ) ? nbsegs[0] : nbsegs[i];
+      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];
+    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
+    }
+    return true;
+  }
+
   case DEFLECTION: {
 
     GCPnts_UniformDeflection Discret(theC3d, _value[ DEFLECTION_IND ], f, l, true);
@@ -815,12 +958,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 ( !_mainEdge.IsNull() )
+    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;
+    }
+    // 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 );
@@ -845,7 +1019,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);
@@ -930,6 +1103,81 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & theMesh, const TopoDS_Shape & t
   return true;
 }
 
+
+//=============================================================================
+/*!
+ *  
+ */
+//=============================================================================
+
+bool StdMeshers_Regular_1D::Evaluate(SMESH_Mesh & theMesh,
+                                     const TopoDS_Shape & theShape,
+                                     MapShapeNbElems& aResMap)
+{
+  if ( _hypType == NONE )
+    return false;
+
+  //SMESHDS_Mesh * meshDS = theMesh.GetMeshDS();
+
+  const TopoDS_Edge & EE = TopoDS::Edge(theShape);
+  TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
+  //  int shapeID = meshDS->ShapeToIndex( E );
+
+  double f, l;
+  Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
+
+  TopoDS_Vertex VFirst, VLast;
+  TopExp::Vertices(E, VFirst, VLast);   // Vfirst corresponds to f and Vlast to l
+
+  ASSERT(!VFirst.IsNull());
+  ASSERT(!VLast.IsNull());
+
+  std::vector<int> aVec(SMDSEntity_Last,0);
+
+  if (!Curve.IsNull()) {
+    list< double > params;
+
+    BRepAdaptor_Curve C3d( E );
+    double length = EdgeLength( E );
+    if ( ! computeInternalParameters( theMesh, C3d, length, f, l, params, false, true )) {
+      SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
+      aResMap.insert(std::make_pair(sm,aVec));
+      SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
+      smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
+      return false;
+    }
+    redistributeNearVertices( theMesh, C3d, length, params, VFirst, VLast );
+
+    if(_quadraticMesh) {
+      aVec[SMDSEntity_Node] = 2*params.size() + 1;
+      aVec[SMDSEntity_Quad_Edge] = params.size() + 1;
+    }
+    else {
+      aVec[SMDSEntity_Node] = params.size();
+      aVec[SMDSEntity_Edge] = params.size() + 1;
+    }
+    
+  }
+  else {
+    //MESSAGE("************* Degenerated edge! *****************");
+    // Edge is a degenerated Edge : We put n = 5 points on the edge.
+    if(_quadraticMesh) {
+      aVec[SMDSEntity_Node] = 11;
+      aVec[SMDSEntity_Quad_Edge] = 6;
+    }
+    else {
+      aVec[SMDSEntity_Node] = 5;
+      aVec[SMDSEntity_Edge] = 6;
+    }
+  }
+
+  SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
+  aResMap.insert(std::make_pair(sm,aVec));
+
+  return true;
+}
+
+
 //=============================================================================
 /*!
  *  See comments in SMESH_Algo.cxx
@@ -946,10 +1194,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)