X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Regular_1D.cxx;h=d3d92f9ea20d8a66a327c83c04ee3208f5623fbc;hb=da425c6a4aa2f29deb92e1ebbebbd0855f373b16;hp=be3e69f1d9fcc6549ec165fc5308af920f886309;hpb=8672ad3e7621ac25fffa8517599afa84ffea509a;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_Regular_1D.cxx b/src/StdMeshers/StdMeshers_Regular_1D.cxx index be3e69f1d..d3d92f9ea 100644 --- a/src/StdMeshers/StdMeshers_Regular_1D.cxx +++ b/src/StdMeshers/StdMeshers_Regular_1D.cxx @@ -30,19 +30,24 @@ 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 #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" @@ -84,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"); } //============================================================================= @@ -159,7 +165,7 @@ bool StdMeshers_Regular_1D::CheckHypothesis } if (_ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_TabFunc || _ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_ExprFunc) - _ivalue[ EXP_MODE_IND ] = (int) hyp->IsExponentMode(); + _ivalue[ CONV_MODE_IND ] = hyp->ConversionMode(); _hypType = NB_SEGMENTS; aStatus = SMESH_Hypothesis::HYP_OK; } @@ -198,6 +204,17 @@ bool StdMeshers_Regular_1D::CheckHypothesis _hypType = DEFLECTION; aStatus = SMESH_Hypothesis::HYP_OK; } + + else if (hypName == "AutomaticLength") + { + StdMeshers_AutomaticLength * hyp = const_cast + (dynamic_cast (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; @@ -247,240 +264,33 @@ static void compensateError(double a1, double an, } } -/*! - * \brief This class provides interface for a density function - */ -class Function -{ -public: - Function(bool expMode) : _expMode(expMode) {} - double operator() (double t) const; - virtual bool IsReady() const = 0; -protected: - virtual double compute(double t) const = 0; -private: - bool _expMode; -}; - -/*! - * \brief This class provides computation of density function given by a table - */ -class TabFunction: public Function -{ -public: - TabFunction(const vector& table, bool expMode); - virtual bool IsReady() const; -protected: - virtual double compute(double t) const; -private: - const vector& _table; -}; - -/*! - * \brief This class provides computation of density function given by an expression - */ -class ExprFunction: public Function -{ -public: - ExprFunction(const char* expr, bool expMode); - virtual bool IsReady() const; -protected: - virtual double compute(double t) const; -private: - Handle(Expr_GeneralExpression) _expression; - Expr_Array1OfNamedUnknown _var; - mutable TColStd_Array1OfReal _val; -}; - -double Function::operator() (double t) const -{ - double res = compute(t); - if (_expMode) - res = pow(10, res); - return res; -} - -TabFunction::TabFunction(const vector& table, bool expMode) - : Function(expMode), - _table(table) -{ -} - -bool TabFunction::IsReady() const -{ - return true; -} - -double TabFunction::compute (double t) const -{ - //find place of in table - int i; - for (i=0; i < _table.size()/2; i++) - if (_table[i*2] > t) - break; - if (i >= _table.size()/2) - i = _table.size()/2 - 1; - - if (i == 0) - return _table[1]; - - // interpolate function value on found interval - // (t - x[i-1]) / (x[i] - x[i-1]) = (y - f[i-1]) / (f[i] - f[i-1]) - // => y = f[i-1] + (f[i] - f[i-1]) * (t - x[i-1]) / (x[i] - x[i-1]) - double x1 = _table[(i-1)*2]; - double x2 = _table[i*2]; - double y1 = _table[(i-1)*2+1]; - double y2 = _table[i*2+1]; - if (x2 - x1 < Precision::Confusion()) - throw SALOME_Exception("TabFunction::compute : confused points"); - return y1 + (y2 - y1) * ((t - x1) / (x2 - x1)); -} - -ExprFunction::ExprFunction(const char* expr, bool expMode) - : Function(expMode), - _var(1,1), - _val(1,1) -{ - Handle( ExprIntrp_GenExp ) gen = ExprIntrp_GenExp::Create(); - gen->Process(TCollection_AsciiString((char*)expr)); - if (gen->IsDone()) - { - _expression = gen->Expression(); - _var(1) = new Expr_NamedUnknown("t"); - } -} - -bool ExprFunction::IsReady() const -{ - return !_expression.IsNull(); -} - -double ExprFunction::compute (double t) const -{ - ASSERT(!_expression.IsNull()); - _val(1) = t; - return _expression->Evaluate(_var, _val); -} - -//================================================================================ -/*! - * \brief Compute next abscissa when two previous ones are given - * \param sm2 - before previous abscissa - * \param sm1 - previous abscissa - * \param func - function of density - * \param reverse - the direction of next abscissa, increase (0) or decrease (1) - * \retval double - the new abscissa - * - * The abscissa s is given by the formulae - * - * ....|--------|----------------|..... - * sm2 sm1 s - * - * func(sm2) / func(sm1) = (sm1-sm2) / (s-sm1) - * => (s-sm1) * func(sm2) = (sm1-sm2) * func(sm1) - * => s = sm1 + (sm1-sm2) * func(sm1) / func(sm2) - */ -//================================================================================ - -static double nextAbscissa(double sm2, double sm1, const Function& func, int reverse) -{ - if (reverse) - { - sm1 = 1.0 - sm1; - sm2 = 1.0 - sm2; - } - return sm1 + (sm1-sm2) * func(sm1) / func(sm2); -} - -//================================================================================ -/*! - * \brief Compute distribution of points on a curve following the law of a function - * \param C3d - the curve to discretize - * \param first - the first parameter on the curve - * \param last - the last parameter on the curve - * \param theReverse - flag indicating that the curve must be reversed - * \param nbSeg - number of output segments - * \param func - the function f(t) - * \param theParams - output points - * \retval bool - true if success - */ -//================================================================================ - static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last, double length, bool theReverse, - int nbSeg, const Function& func, + int nbSeg, Function& func, list& theParams) { - if (!func.IsReady()) - return false; - vector xxx[2]; - int nbPnt = 1 + nbSeg; - int rev, i; - for (rev=0; rev < 2; rev++) - { - // curv abscisses initialisation - vector x(nbPnt, 0.); - // the first abscissa is 0.0 - - // The aim of the algorithm is to find a second abscisse x[1] such as the last - // one x[nbSeg] is very close to 1.0 with the epsilon precision + OSD::SetSignal( true ); - double x1_too_small = 0.0; - double x1_too_large = RealLast(); - double x1 = 1.0/nbSeg; - while (1) - { - x[1] = x1; - - // Check if the abscissa of the point 2 to N-1 - // are in the segment ... - - bool ok = true; - for (i=2; i <= nbSeg; i++) - { - x[i] = nextAbscissa(x[i-2], x[i-1], func, rev); - if (x[i] - 1.0 > Precision::Confusion()) - { - x[nbSeg] = x[i]; - ok = false; - break; - } - } - if (!ok) - { - // The segments are to large - // Decrease x1 ... - x1_too_large = x1; - x1 = (x1_too_small+x1_too_large)/2; - continue; - } - - // Look at the abscissa of the point N - // which is to be close to 1.0 - - // break condition --> algo converged !! - - if (1.0 - x[nbSeg] < Precision::Confusion()) - break; + if( nbSeg<=0 ) + return false; - // not ok ... + MESSAGE( "computeParamByFunc" ); - x1_too_small = x1; + int nbPnt = 1 + nbSeg; + vector x(nbPnt, 0.); - // Modify x1 value + if( !buildDistribution( func, 0.0, 1.0, nbSeg, x, 1E-4 ) ) + return false; - if (x1_too_large > 1e100) - x1 = 2*x1; - else - x1 = (x1_too_small+x1_too_large)/2; - } - xxx[rev] = x; + MESSAGE( "Points:\n" ); + char buf[1024]; + for( int i=0; i<=nbSeg; i++ ) + { + sprintf( buf, "%f\n", float(x[i] ) ); + MESSAGE( buf ); } + - // average - vector x(nbPnt, 0.); - for (i=0; i < nbPnt; i++) - x[i] = (xxx[0][i] + (1.0 - xxx[1][nbPnt-i])) / 2; // apply parameters in range [0,1] to the space of the curve double prevU = first; @@ -490,7 +300,7 @@ static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last, prevU = last; sign = -1.; } - for (i = 1; i < nbSeg; i++) + for( int i = 1; i < nbSeg; i++ ) { double curvLength = length * (x[i] - x[i-1]) * sign; GCPnts_AbscissaPoint Discret( C3d, curvLength, prevU ); @@ -503,7 +313,7 @@ static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last, return false; prevU = U; } - return false; + return true; } //============================================================================= @@ -561,7 +371,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge break; case StdMeshers_NumberOfSegments::DT_TabFunc: { - TabFunction func(_vvalue[ TAB_FUNC_IND ], (bool)_ivalue[ EXP_MODE_IND ]); + FunctionTable func(_vvalue[ TAB_FUNC_IND ], _ivalue[ CONV_MODE_IND ]); return computeParamByFunc(C3d, f, l, length, theReverse, _ivalue[ NB_SEGMENTS_IND ], func, theParams); @@ -569,7 +379,7 @@ bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge break; case StdMeshers_NumberOfSegments::DT_ExprFunc: { - ExprFunction func(_svalue[ EXPR_FUNC_IND ].c_str(), (bool)_ivalue[ EXP_MODE_IND ]); + 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); @@ -593,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; } @@ -697,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); @@ -750,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(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 { @@ -781,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(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; }