X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Regular_1D.cxx;h=78519c8617299661b7242bb2da885c1c6c3537db;hp=51df8669dd44badccf07f0fe71fc074784688a2c;hb=57b43b4d010e2d0a1529d3c131bbb9d416e63258;hpb=2ad752b10cb247a162b27593121a7ae13173258d diff --git a/src/StdMeshers/StdMeshers_Regular_1D.cxx b/src/StdMeshers/StdMeshers_Regular_1D.cxx index 51df8669d..78519c861 100644 --- a/src/StdMeshers/StdMeshers_Regular_1D.cxx +++ b/src/StdMeshers/StdMeshers_Regular_1D.cxx @@ -30,20 +30,22 @@ 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 "StdMeshers_LocalLength.hxx" #include "StdMeshers_NumberOfSegments.hxx" #include "StdMeshers_Arithmetic1D.hxx" #include "StdMeshers_StartEndLength.hxx" #include "StdMeshers_Deflection1D.hxx" -#include +#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" @@ -63,16 +65,13 @@ using namespace std; #include #include #include - -#include -#include -#include #include -#include #include #include +using namespace std; + //============================================================================= /*! * @@ -92,6 +91,8 @@ StdMeshers_Regular_1D::StdMeshers_Regular_1D(int hypId, int studyId, _compatibleHypothesis.push_back("Deflection1D"); _compatibleHypothesis.push_back("Arithmetic1D"); _compatibleHypothesis.push_back("AutomaticLength"); + + _compatibleHypothesis.push_back("QuadraticMesh"); // auxiliary !!! } //============================================================================= @@ -111,22 +112,37 @@ StdMeshers_Regular_1D::~StdMeshers_Regular_1D() //============================================================================= bool StdMeshers_Regular_1D::CheckHypothesis - (SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, + (SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, SMESH_Hypothesis::Hypothesis_Status& aStatus) { _hypType = NONE; + _quadraticMesh = false; + + const bool ignoreAuxiliaryHyps = false; + const list & hyps = + GetUsedHypothesis(aMesh, aShape, ignoreAuxiliaryHyps); + + // find non-auxiliary hypothesis + const SMESHDS_Hypothesis *theHyp = 0; + list ::const_iterator h = hyps.begin(); + for ( ; h != hyps.end(); ++h ) { + if ( static_cast(*h)->IsAuxiliary() ) { + if ( strcmp( "QuadraticMesh", (*h)->GetName() ) == 0 ) + _quadraticMesh = true; + } + else { + if ( !theHyp ) + theHyp = *h; // use only the first non-auxiliary hypothesis + } + } - const list &hyps = GetUsedHypothesis(aMesh, aShape); - if (hyps.size() == 0) + if ( !theHyp ) { aStatus = SMESH_Hypothesis::HYP_MISSING; return false; // can't work without a hypothesis } - // use only the first hypothesis - const SMESHDS_Hypothesis *theHyp = hyps.front(); - string hypName = theHyp->GetName(); if (hypName == "LocalLength") @@ -167,7 +183,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; } @@ -266,304 +282,13 @@ static void compensateError(double a1, double an, } } -class Function -{ -public: - Function( const bool exp ) - : myExp( exp ) - { - } - - virtual ~Function() - { - } - - virtual bool value( const double, double& f ) - { - if( myExp ) - f = pow( 10, f ); - return true; - } - virtual double integral( const double, const double ) = 0; - -private: - bool myExp; -}; - -class FunctionIntegral : public Function -{ -public: - FunctionIntegral( Function*, const double ); - virtual ~FunctionIntegral(); - virtual bool value( const double, double& ); - virtual double integral( const double, const double ); - -private: - Function* myFunc; - double myStart; -}; - -FunctionIntegral::FunctionIntegral( Function* f, const double st ) -: Function( false ) -{ - myFunc = f; - myStart = st; -} - -FunctionIntegral::~FunctionIntegral() -{ -} - -bool FunctionIntegral::value( const double t, double& f ) -{ - f = myFunc ? myFunc->integral( myStart, t ) : 0; - return myFunc!=0 && Function::value( t, f ); -} - -double FunctionIntegral::integral( const double, const double ) -{ - return 0; -} - -class FunctionTable : public Function -{ -public: - FunctionTable( const std::vector&, const bool ); - virtual ~FunctionTable(); - virtual bool value( const double, double& ); - virtual double integral( const double, const double ); - -private: - bool findBounds( const double, int&, int& ) const; - - //integral from x[i] to x[i+1] - double integral( const int i ); - - //integral from x[i] to x[i]+d - //warning: function is presented as linear on interaval from x[i] to x[i]+d, - // for correct result d must be >=0 and <=x[i+1]-x[i] - double integral( const int i, const double d ); - -private: - std::vector myData; -}; - -FunctionTable::FunctionTable( const std::vector& data, const bool exp ) -: Function( exp ) -{ - myData = data; -} - -FunctionTable::~FunctionTable() -{ -} - -bool FunctionTable::value( const double t, double& f ) -{ - int i1, i2; - if( !findBounds( t, i1, i2 ) ) - return false; - - double - x1 = myData[2*i1], y1 = myData[2*i1+1], - x2 = myData[2*i2], y2 = myData[2*i2+1]; - - Function::value( x1, y1 ); - Function::value( x2, y2 ); - - f = y1 + ( y2-y1 ) * ( t-x1 ) / ( x2-x1 ); - return true; -} - -double FunctionTable::integral( const int i ) -{ - if( i>=0 && iProcess( ( Standard_CString )str ); - if( !myExpr->IsDone() ) - myExpr.Nullify(); - - myVars.ChangeValue( 1 ) = new Expr_NamedUnknown( "t" ); -} - -FunctionExpr::~FunctionExpr() -{ -} - -Standard_Boolean FunctionExpr::Value( Standard_Real T, Standard_Real& F ) -{ - double f; - Standard_Boolean res = value( T, f ); - F = f; - return res; -} - -bool FunctionExpr::value( const double t, double& f ) -{ - if( myExpr.IsNull() ) - return false; - - CASCatch_CatchSignals aCatchSignals; - aCatchSignals.Activate(); - - myValues.ChangeValue( 1 ) = t; - bool ok = true; - CASCatch_TRY { - f = myExpr->Expression()->Evaluate( myVars, myValues ); - } - CASCatch_CATCH(CASCatch_Failure) { - aCatchSignals.Deactivate(); - Handle(CASCatch_Failure) aFail = CASCatch_Failure::Caught(); - f = 0.0; - } - - aCatchSignals.Deactivate(); - ok = Function::value( t, f ) && ok; - return ok; -} - -double FunctionExpr::integral( const double a, const double b ) -{ - double res = 0.0; - CASCatch_TRY - { - math_GaussSingleIntegration _int( *this, a, b, 20 ); - if( _int.IsDone() ) - res = _int.Value(); - } - CASCatch_CATCH(CASCatch_Failure) - { - res = 0.0; - MESSAGE( "Exception in integral calculating" ); - } - return res; -} - - - - - - - -double dihotomySolve( Function& f, const double val, const double _start, const double _fin, const double eps, bool& ok ) -{ - double start = _start, fin = _fin, start_val, fin_val; bool ok1, ok2; - ok1 = f.value( start, start_val ); - ok2 = f.value( fin, fin_val ); - - if( !ok1 || !ok2 ) - { - ok = false; - return 0.0; - } - - bool start_pos = start_val>=val, fin_pos = fin_val>=val; - ok = true; - - while( fin-start>eps ) - { - double mid = ( start+fin )/2.0, mid_val; - ok = f.value( mid, mid_val ); - if( !ok ) - return 0.0; - -// char buf[1024]; -// sprintf( buf, "start=%f\nfin=%f\nmid_val=%f\n", float( start ), float( fin ), float( mid_val ) ); -// MESSAGE( buf ); - - bool mid_pos = mid_val>=val; - if( start_pos!=mid_pos ) - { - fin_pos = mid_pos; - fin = mid; - } - else if( fin_pos!=mid_pos ) - { - start_pos = mid_pos; - start = mid; - } - else - break; - } - return (start+fin)/2.0; -} - static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last, double length, bool theReverse, int nbSeg, Function& func, list& theParams) { OSD::SetSignal( true ); + if( nbSeg<=0 ) return false; @@ -572,18 +297,9 @@ static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last, int nbPnt = 1 + nbSeg; vector x(nbPnt, 0.); - x[0] = 0.0; - double J = func.integral( 0.0, 1.0 ) / nbSeg; - bool ok; - for( int i=1; iGetSubMeshDS()->GetNodes(); - if (!lid->more()) - { + if (!lid->more()) { MESSAGE (" NO NODE BUILT ON VERTEX "); return false; } const SMDS_MeshNode * idLast = lid->next(); - if (!Curve.IsNull()) - { + if (!Curve.IsNull()) { list< double > params; bool reversed = false; if ( !_mainEdge.IsNull() ) reversed = aMesh.IsReversedInChain( EE, _mainEdge ); try { - if ( ! computeInternalParameters( E, params, reversed )) + if ( ! computeInternalParameters( E, params, reversed )) { + //cout << "computeInternalParameters() failed" <::iterator itU = params.begin(); itU != params.end(); itU++) - { + for (list::iterator itU = params.begin(); itU != params.end(); itU++) { double param = *itU; gp_Pnt P = Curve->Value(param); @@ -865,17 +600,39 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z()); meshDS->SetNodeOnEdge(node, shapeID, param); - SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node); - meshDS->SetMeshElementOnShape(edge, shapeID); + if(_quadraticMesh) { + // create medium node + double prm = ( parPrev + param )/2; + gp_Pnt PM = Curve->Value(prm); + SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z()); + meshDS->SetNodeOnEdge(NM, shapeID, prm); + SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM); + meshDS->SetMeshElementOnShape(edge, shapeID); + } + else { + SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node); + meshDS->SetMeshElementOnShape(edge, shapeID); + } + idPrev = node; + parPrev = param; + } + if(_quadraticMesh) { + double prm = ( parPrev + parLast )/2; + gp_Pnt PM = Curve->Value(prm); + SMDS_MeshNode * NM = meshDS->AddNode(PM.X(), PM.Y(), PM.Z()); + meshDS->SetNodeOnEdge(NM, shapeID, prm); + SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM); + meshDS->SetMeshElementOnShape(edge, shapeID); + } + else { + SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast); + meshDS->SetMeshElementOnShape(edge, shapeID); } - SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast); - meshDS->SetMeshElementOnShape(edge, shapeID); } - else - { + else { // Edge is a degenerated Edge : We put n = 5 points on the edge. - int NbPoints = 5; + const int NbPoints = 5; BRep_Tool::Range(E, f, l); double du = (l - f) / (NbPoints - 1); //MESSAGE("************* Degenerated edge! *****************"); @@ -885,18 +642,36 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh gp_Pnt P = BRep_Tool::Pnt(V1); const SMDS_MeshNode * idPrev = idFirst; - for (int i = 2; i < NbPoints; i++) - { + for (int i = 2; i < NbPoints; i++) { double param = f + (i - 1) * du; SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z()); + if(_quadraticMesh) { + // create medium node + double prm = param - du/2.; + SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnEdge(NM, shapeID, prm); + SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node, NM); + meshDS->SetMeshElementOnShape(edge, shapeID); + } + else { + SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node); + meshDS->SetMeshElementOnShape(edge, shapeID); + } meshDS->SetNodeOnEdge(node, shapeID, param); - - SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node); - meshDS->SetMeshElementOnShape(edge, shapeID); idPrev = node; } - SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast); - meshDS->SetMeshElementOnShape(edge, shapeID); + if(_quadraticMesh) { + // create medium node + double prm = l - du/2.; + SMDS_MeshNode * NM = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnEdge(NM, shapeID, prm); + SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast, NM); + meshDS->SetMeshElementOnShape(edge, shapeID); + } + else { + SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast); + meshDS->SetMeshElementOnShape(edge, shapeID); + } } return true; } @@ -907,40 +682,47 @@ bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aSh */ //============================================================================= -const list & StdMeshers_Regular_1D::GetUsedHypothesis( - SMESH_Mesh & aMesh, const TopoDS_Shape & aShape) +const list & +StdMeshers_Regular_1D::GetUsedHypothesis(SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape, + const bool ignoreAuxiliary) { _usedHypList.clear(); - _usedHypList = GetAppliedHypothesis(aMesh, aShape); // copy - int nbHyp = _usedHypList.size(); _mainEdge.Nullify(); + + SMESH_HypoFilter auxiliaryFilter, compatibleFilter; + auxiliaryFilter.Init( SMESH_HypoFilter::IsAuxiliary() ); + const bool ignoreAux = true; + InitCompatibleHypoFilter( compatibleFilter, ignoreAux ); + + // get non-auxiliary assigned to aShape + int nbHyp = aMesh.GetHypotheses( aShape, compatibleFilter, _usedHypList, false ); + if (nbHyp == 0) { // Check, if propagated from some other edge if (aShape.ShapeType() == TopAbs_EDGE && aMesh.IsPropagatedHypothesis(aShape, _mainEdge)) { - // Propagation of 1D hypothesis from on this edge - //_usedHypList = GetAppliedHypothesis(aMesh, _mainEdge); // copy - // use a general method in order not to nullify _mainEdge - _usedHypList = SMESH_Algo::GetUsedHypothesis(aMesh, _mainEdge); // copy - nbHyp = _usedHypList.size(); + // Propagation of 1D hypothesis from on this edge; + // get non-auxiliary assigned to _mainEdge + nbHyp = aMesh.GetHypotheses( _mainEdge, compatibleFilter, _usedHypList, true ); } } - if (nbHyp == 0) + + if (nbHyp == 0) // nothing propagated nor assigned to aShape { - TopTools_ListIteratorOfListOfShape ancIt( aMesh.GetAncestors( aShape )); - for (; ancIt.More(); ancIt.Next()) - { - const TopoDS_Shape& ancestor = ancIt.Value(); - _usedHypList = GetAppliedHypothesis(aMesh, ancestor); // copy - nbHyp = _usedHypList.size(); - if (nbHyp == 1) - break; - } + SMESH_Algo::GetUsedHypothesis( aMesh, aShape, ignoreAuxiliary ); + nbHyp = _usedHypList.size(); } - if (nbHyp > 1) - _usedHypList.clear(); //only one compatible hypothesis allowed + else + { + // get auxiliary hyps from aShape + aMesh.GetHypotheses( aShape, auxiliaryFilter, _usedHypList, true ); + } + if ( nbHyp > 1 && ignoreAuxiliary ) + _usedHypList.clear(); //only one compatible non-auxiliary hypothesis allowed + return _usedHypList; }