Salome HOME
add QuadraticMesh hypothesis
[modules/smesh.git] / src / StdMeshers / StdMeshers_Regular_1D.cxx
index 017af4544b42d5b05d5f3e79ecf9c7b4e67bcd66..d3d92f9ea20d8a66a327c83c04ee3208f5623fbc 100644 (file)
 using namespace std;
 
 #include "StdMeshers_Regular_1D.hxx"
+#include "StdMeshers_Distribution.hxx"
 #include "SMESH_Gen.hxx"
 #include "SMESH_Mesh.hxx"
+#include "SMESH_HypoFilter.hxx"
+#include "SMESH_subMesh.hxx"
+
+#include <OSD.hxx>
 
 #include "StdMeshers_LocalLength.hxx"
 #include "StdMeshers_NumberOfSegments.hxx"
 #include "StdMeshers_Arithmetic1D.hxx"
 #include "StdMeshers_StartEndLength.hxx"
 #include "StdMeshers_Deflection1D.hxx"
+#include "StdMeshers_AutomaticLength.hxx"
 
 #include "SMDS_MeshElement.hxx"
 #include "SMDS_MeshNode.hxx"
 #include "SMDS_EdgePosition.hxx"
-#include "SMESH_subMesh.hxx"
 
+#include "Utils_SALOME_Exception.hxx"
 #include "utilities.h"
 
 #include <BRep_Tool.hxx>
@@ -56,9 +62,14 @@ using namespace std;
 #include <GCPnts_UniformDeflection.hxx>
 #include <Standard_ErrorHandler.hxx>
 #include <Precision.hxx>
+#include <Expr_GeneralExpression.hxx>
+#include <Expr_NamedUnknown.hxx>
+#include <Expr_Array1OfNamedUnknown.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <ExprIntrp_GenExp.hxx>
 
 #include <string>
-//#include <algorithm>
+#include <math.h>
 
 //=============================================================================
 /*!
@@ -78,6 +89,7 @@ StdMeshers_Regular_1D::StdMeshers_Regular_1D(int hypId, int studyId,
        _compatibleHypothesis.push_back("StartEndLength");
        _compatibleHypothesis.push_back("Deflection1D");
        _compatibleHypothesis.push_back("Arithmetic1D");
+       _compatibleHypothesis.push_back("AutomaticLength");
 }
 
 //=============================================================================
@@ -131,9 +143,29 @@ bool StdMeshers_Regular_1D::CheckHypothesis
     const StdMeshers_NumberOfSegments * hyp =
       dynamic_cast <const StdMeshers_NumberOfSegments * >(theHyp);
     ASSERT(hyp);
-    _value[ NB_SEGMENTS_IND  ] = hyp->GetNumberOfSegments();
-    _value[ SCALE_FACTOR_IND ] = hyp->GetScaleFactor();
-    ASSERT( _value[ NB_SEGMENTS_IND ] > 0 );
+    _ivalue[ NB_SEGMENTS_IND  ] = hyp->GetNumberOfSegments();
+    ASSERT( _ivalue[ NB_SEGMENTS_IND ] > 0 );
+    _ivalue[ DISTR_TYPE_IND ] = (int) hyp->GetDistrType();
+    switch (_ivalue[ DISTR_TYPE_IND ])
+    {
+    case StdMeshers_NumberOfSegments::DT_Scale:
+      _value[ SCALE_FACTOR_IND ] = hyp->GetScaleFactor();
+      break;
+    case StdMeshers_NumberOfSegments::DT_TabFunc:
+      _vvalue[ TAB_FUNC_IND ] = hyp->GetTableFunction();
+      break;
+    case StdMeshers_NumberOfSegments::DT_ExprFunc:
+      _svalue[ EXPR_FUNC_IND ] = hyp->GetExpressionFunction();
+      break;
+    case StdMeshers_NumberOfSegments::DT_Regular:
+      break;
+    default:
+      ASSERT(0);
+      break;
+    }
+    if (_ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_TabFunc ||
+        _ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_ExprFunc)
+        _ivalue[ CONV_MODE_IND ] = hyp->ConversionMode();
     _hypType = NB_SEGMENTS;
     aStatus = SMESH_Hypothesis::HYP_OK;
   }
@@ -172,6 +204,17 @@ bool StdMeshers_Regular_1D::CheckHypothesis
     _hypType = DEFLECTION;
     aStatus = SMESH_Hypothesis::HYP_OK;
   }
+
+  else if (hypName == "AutomaticLength")
+  {
+    StdMeshers_AutomaticLength * hyp = const_cast<StdMeshers_AutomaticLength *>
+      (dynamic_cast <const StdMeshers_AutomaticLength * >(theHyp));
+    ASSERT(hyp);
+    _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength( &aMesh, aShape );
+    ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
+    _hypType = LOCAL_LENGTH;
+    aStatus = SMESH_Hypothesis::HYP_OK;
+  }
   else
     aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
 
@@ -221,6 +264,58 @@ static void compensateError(double a1, double an,
   }
 }
 
+static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
+                               double length, bool theReverse, 
+                               int nbSeg, Function& func,
+                               list<double>& theParams)
+{
+  OSD::SetSignal( true );
+
+  if( nbSeg<=0 )
+    return false;
+
+  MESSAGE( "computeParamByFunc" );
+
+  int nbPnt = 1 + nbSeg;
+  vector<double> x(nbPnt, 0.);
+
+  if( !buildDistribution( func, 0.0, 1.0, nbSeg, x, 1E-4 ) )
+     return false;
+
+  MESSAGE( "Points:\n" );
+  char buf[1024];
+  for( int i=0; i<=nbSeg; i++ )
+  {
+    sprintf(  buf, "%f\n", float(x[i] ) );
+    MESSAGE( buf );
+  }
+    
+
+
+  // apply parameters in range [0,1] to the space of the curve
+  double prevU = first;
+  double sign = 1.;
+  if (theReverse)
+  {
+    prevU = last;
+    sign = -1.;
+  }
+  for( int i = 1; i < nbSeg; i++ )
+  {
+    double curvLength = length * (x[i] - x[i-1]) * sign;
+    GCPnts_AbscissaPoint Discret( C3d, curvLength, prevU );
+    if ( !Discret.IsDone() )
+      return false;
+    double U = Discret.Parameter();
+    if ( U > first && U < last )
+      theParams.push_back( U );
+    else
+      return false;
+    prevU = U;
+  }
+  return true;
+}
+
 //=============================================================================
 /*!
  *  
@@ -246,6 +341,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge
     double eltSize = 1;
     if ( _hypType == LOCAL_LENGTH )
     {
+      // Local Length hypothesis
       double nbseg = ceil(length / _value[ BEG_LENGTH_IND ]); // integer sup
       if (nbseg <= 0)
         nbseg = 1;                        // degenerated edge
@@ -253,26 +349,47 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge
     }
     else
     {
-      double epsilon = 0.001;
-      if (fabs(_value[ SCALE_FACTOR_IND ] - 1.0) > epsilon)
+      // Number Of Segments hypothesis
+      switch (_ivalue[ DISTR_TYPE_IND ])
       {
-        double scale = _value[ SCALE_FACTOR_IND ];
-        if ( theReverse )
-          scale = 1. / scale;
-        double alpha = pow( scale , 1.0 / (_value[ NB_SEGMENTS_IND ] - 1));
-        double factor = (l - f) / (1 - pow( alpha,_value[ NB_SEGMENTS_IND ]));
-
-        int i, NbPoints = 1 + (int) _value[ NB_SEGMENTS_IND ];
-        for ( i = 2; i < NbPoints; i++ )
+      case StdMeshers_NumberOfSegments::DT_Scale:
         {
-          double param = f + factor * (1 - pow(alpha, i - 1));
-          theParams.push_back( param );
+          double scale = _value[ SCALE_FACTOR_IND ];
+          if ( theReverse )
+            scale = 1. / scale;
+          double alpha = pow( scale , 1.0 / (_ivalue[ NB_SEGMENTS_IND ] - 1));
+          double factor = (l - f) / (1 - pow( alpha,_ivalue[ NB_SEGMENTS_IND ]));
+
+          int i, NbPoints = 1 + _ivalue[ NB_SEGMENTS_IND ];
+          for ( i = 2; i < NbPoints; i++ )
+          {
+            double param = f + factor * (1 - pow(alpha, i - 1));
+            theParams.push_back( param );
+          }
+          return true;
         }
-        return true;
-      }
-      else
-      {
-        eltSize = length / _value[ NB_SEGMENTS_IND ];
+        break;
+      case StdMeshers_NumberOfSegments::DT_TabFunc:
+        {
+          FunctionTable func(_vvalue[ TAB_FUNC_IND ], _ivalue[ CONV_MODE_IND ]);
+          return computeParamByFunc(C3d, f, l, length, theReverse,
+                                    _ivalue[ NB_SEGMENTS_IND ], func,
+                                    theParams);
+        }
+        break;
+      case StdMeshers_NumberOfSegments::DT_ExprFunc:
+        {
+          FunctionExpr func(_svalue[ EXPR_FUNC_IND ].c_str(), _ivalue[ CONV_MODE_IND ]);
+          return computeParamByFunc(C3d, f, l, length, theReverse,
+                                    _ivalue[ NB_SEGMENTS_IND ], func,
+                                    theParams);
+        }
+        break;
+      case StdMeshers_NumberOfSegments::DT_Regular:
+        eltSize = length / _ivalue[ NB_SEGMENTS_IND ];
+        break;
+      default:
+        return false;
       }
     }
 
@@ -286,6 +403,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge
       double param = Discret.Parameter(i);
       theParams.push_back( param );
     }
+    compensateError( eltSize, eltSize, f, l, length, C3d, theParams ); // for PAL9899
     return true;
   }
 
@@ -390,8 +508,13 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
   aMesh.GetSubMesh(aShape);
 
+  // quardatic mesh required?
+  SMESH_HypoFilter filter( SMESH_HypoFilter::HasName( "QuadraticMesh" ));
+  bool isQuadraticMesh = aMesh->GetHypothesis( aShape, filter, true );
+
   const TopoDS_Edge & EE = TopoDS::Edge(aShape);
   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);
@@ -443,19 +566,14 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
 
       //Add the Node in the DataStructure
       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
-      meshDS->SetNodeOnEdge(node, E);
-
-      // **** edgePosition associe au point = param. 
-      SMDS_EdgePosition* epos =
-        dynamic_cast<SMDS_EdgePosition *>(node->GetPosition().get());
-      epos->SetUParameter(param);
+      meshDS->SetNodeOnEdge(node, shapeID, param);
 
       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
-      meshDS->SetMeshElementOnShape(edge, E);
+      meshDS->SetMeshElementOnShape(edge, shapeID);
       idPrev = node;
     }
     SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast);
-    meshDS->SetMeshElementOnShape(edge, E);
+    meshDS->SetMeshElementOnShape(edge, shapeID);
   }
   else
   {
@@ -474,18 +592,14 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh
     {
       double param = f + (i - 1) * du;
       SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
-      meshDS->SetNodeOnEdge(node, E);
-
-      SMDS_EdgePosition* epos =
-        dynamic_cast<SMDS_EdgePosition*>(node->GetPosition().get());
-      epos->SetUParameter(param);
+      meshDS->SetNodeOnEdge(node, shapeID, param);
 
       SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
-      meshDS->SetMeshElementOnShape(edge, E);
+      meshDS->SetMeshElementOnShape(edge, shapeID);
       idPrev = node;
     }
     SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast);
-    meshDS->SetMeshElementOnShape(edge, E);
+    meshDS->SetMeshElementOnShape(edge, shapeID);
   }
   return true;
 }