1 // SMESH SMESH : implementaion of SMESH idl descriptions
3 // Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
24 // File : StdMeshers_Regular_1D.cxx
25 // Moved here from SMESH_Regular_1D.cxx
26 // Author : Paul RASCLE, EDF
32 #include "StdMeshers_Regular_1D.hxx"
33 #include "StdMeshers_Distribution.hxx"
34 #include "SMESH_Gen.hxx"
35 #include "SMESH_Mesh.hxx"
36 #include "SMESH_HypoFilter.hxx"
37 #include "SMESH_subMesh.hxx"
41 #include "StdMeshers_LocalLength.hxx"
42 #include "StdMeshers_NumberOfSegments.hxx"
43 #include "StdMeshers_Arithmetic1D.hxx"
44 #include "StdMeshers_StartEndLength.hxx"
45 #include "StdMeshers_Deflection1D.hxx"
46 #include "StdMeshers_AutomaticLength.hxx"
48 #include "SMDS_MeshElement.hxx"
49 #include "SMDS_MeshNode.hxx"
50 #include "SMDS_EdgePosition.hxx"
52 #include "Utils_SALOME_Exception.hxx"
53 #include "utilities.h"
55 #include <BRep_Tool.hxx>
56 #include <TopoDS_Edge.hxx>
57 #include <TopoDS_Shape.hxx>
58 #include <TopTools_ListIteratorOfListOfShape.hxx>
59 #include <GeomAdaptor_Curve.hxx>
60 #include <GCPnts_AbscissaPoint.hxx>
61 #include <GCPnts_UniformAbscissa.hxx>
62 #include <GCPnts_UniformDeflection.hxx>
63 #include <Standard_ErrorHandler.hxx>
64 #include <Precision.hxx>
65 #include <Expr_GeneralExpression.hxx>
66 #include <Expr_NamedUnknown.hxx>
67 #include <Expr_Array1OfNamedUnknown.hxx>
68 #include <TColStd_Array1OfReal.hxx>
69 #include <ExprIntrp_GenExp.hxx>
74 //=============================================================================
78 //=============================================================================
80 StdMeshers_Regular_1D::StdMeshers_Regular_1D(int hypId, int studyId,
81 SMESH_Gen * gen):SMESH_1D_Algo(hypId, studyId, gen)
83 MESSAGE("StdMeshers_Regular_1D::StdMeshers_Regular_1D");
85 _shapeType = (1 << TopAbs_EDGE);
87 _compatibleHypothesis.push_back("LocalLength");
88 _compatibleHypothesis.push_back("NumberOfSegments");
89 _compatibleHypothesis.push_back("StartEndLength");
90 _compatibleHypothesis.push_back("Deflection1D");
91 _compatibleHypothesis.push_back("Arithmetic1D");
92 _compatibleHypothesis.push_back("AutomaticLength");
95 //=============================================================================
99 //=============================================================================
101 StdMeshers_Regular_1D::~StdMeshers_Regular_1D()
105 //=============================================================================
109 //=============================================================================
111 bool StdMeshers_Regular_1D::CheckHypothesis
113 const TopoDS_Shape& aShape,
114 SMESH_Hypothesis::Hypothesis_Status& aStatus)
118 const list <const SMESHDS_Hypothesis * >&hyps = GetUsedHypothesis(aMesh, aShape);
119 if (hyps.size() == 0)
121 aStatus = SMESH_Hypothesis::HYP_MISSING;
122 return false; // can't work without a hypothesis
125 // use only the first hypothesis
126 const SMESHDS_Hypothesis *theHyp = hyps.front();
128 string hypName = theHyp->GetName();
130 if (hypName == "LocalLength")
132 const StdMeshers_LocalLength * hyp =
133 dynamic_cast <const StdMeshers_LocalLength * >(theHyp);
135 _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength();
136 ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
137 _hypType = LOCAL_LENGTH;
138 aStatus = SMESH_Hypothesis::HYP_OK;
141 else if (hypName == "NumberOfSegments")
143 const StdMeshers_NumberOfSegments * hyp =
144 dynamic_cast <const StdMeshers_NumberOfSegments * >(theHyp);
146 _ivalue[ NB_SEGMENTS_IND ] = hyp->GetNumberOfSegments();
147 ASSERT( _ivalue[ NB_SEGMENTS_IND ] > 0 );
148 _ivalue[ DISTR_TYPE_IND ] = (int) hyp->GetDistrType();
149 switch (_ivalue[ DISTR_TYPE_IND ])
151 case StdMeshers_NumberOfSegments::DT_Scale:
152 _value[ SCALE_FACTOR_IND ] = hyp->GetScaleFactor();
154 case StdMeshers_NumberOfSegments::DT_TabFunc:
155 _vvalue[ TAB_FUNC_IND ] = hyp->GetTableFunction();
157 case StdMeshers_NumberOfSegments::DT_ExprFunc:
158 _svalue[ EXPR_FUNC_IND ] = hyp->GetExpressionFunction();
160 case StdMeshers_NumberOfSegments::DT_Regular:
166 if (_ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_TabFunc ||
167 _ivalue[ DISTR_TYPE_IND ] == StdMeshers_NumberOfSegments::DT_ExprFunc)
168 _ivalue[ CONV_MODE_IND ] = hyp->ConversionMode();
169 _hypType = NB_SEGMENTS;
170 aStatus = SMESH_Hypothesis::HYP_OK;
173 else if (hypName == "Arithmetic1D")
175 const StdMeshers_Arithmetic1D * hyp =
176 dynamic_cast <const StdMeshers_Arithmetic1D * >(theHyp);
178 _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
179 _value[ END_LENGTH_IND ] = hyp->GetLength( false );
180 ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
181 _hypType = ARITHMETIC_1D;
182 aStatus = SMESH_Hypothesis::HYP_OK;
185 else if (hypName == "StartEndLength")
187 const StdMeshers_StartEndLength * hyp =
188 dynamic_cast <const StdMeshers_StartEndLength * >(theHyp);
190 _value[ BEG_LENGTH_IND ] = hyp->GetLength( true );
191 _value[ END_LENGTH_IND ] = hyp->GetLength( false );
192 ASSERT( _value[ BEG_LENGTH_IND ] > 0 && _value[ END_LENGTH_IND ] > 0 );
193 _hypType = BEG_END_LENGTH;
194 aStatus = SMESH_Hypothesis::HYP_OK;
197 else if (hypName == "Deflection1D")
199 const StdMeshers_Deflection1D * hyp =
200 dynamic_cast <const StdMeshers_Deflection1D * >(theHyp);
202 _value[ DEFLECTION_IND ] = hyp->GetDeflection();
203 ASSERT( _value[ DEFLECTION_IND ] > 0 );
204 _hypType = DEFLECTION;
205 aStatus = SMESH_Hypothesis::HYP_OK;
208 else if (hypName == "AutomaticLength")
210 StdMeshers_AutomaticLength * hyp = const_cast<StdMeshers_AutomaticLength *>
211 (dynamic_cast <const StdMeshers_AutomaticLength * >(theHyp));
213 _value[ BEG_LENGTH_IND ] = _value[ END_LENGTH_IND ] = hyp->GetLength( &aMesh, aShape );
214 ASSERT( _value[ BEG_LENGTH_IND ] > 0 );
215 _hypType = LOCAL_LENGTH;
216 aStatus = SMESH_Hypothesis::HYP_OK;
219 aStatus = SMESH_Hypothesis::HYP_INCOMPATIBLE;
221 return ( _hypType != NONE );
224 //=======================================================================
225 //function : compensateError
226 //purpose : adjust theParams so that the last segment length == an
227 //=======================================================================
229 static void compensateError(double a1, double an,
230 double U1, double Un,
232 GeomAdaptor_Curve& C3d,
233 list<double> & theParams)
235 int i, nPar = theParams.size();
236 if ( a1 + an < length && nPar > 1 )
238 list<double>::reverse_iterator itU = theParams.rbegin();
240 // dist from the last point to the edge end <Un>, it should be equal <an>
241 double Ln = GCPnts_AbscissaPoint::Length( C3d, Ul, Un );
242 double dLn = an - Ln; // error of <an>
243 if ( Abs( dLn ) <= Precision::Confusion() )
245 double dU = Abs( Ul - *itU ); // parametric length of the last but one segment
246 double dUn = dLn * Abs( Un - U1 ) / length; // parametric error of <an>
247 if ( dUn < 0.5 * dU ) { // last segment is a bit shorter than it should
248 dUn = -dUn; // move the last parameter to the edge beginning
250 else { // last segment is much shorter than it should -> remove the last param and
251 theParams.pop_back(); nPar--; // move the rest points toward the edge end
252 Ln = GCPnts_AbscissaPoint::Length( C3d, theParams.back(), Un );
253 dUn = ( an - Ln ) * Abs( Un - U1 ) / length;
254 if ( dUn < 0.5 * dU )
259 double q = dUn / ( nPar - 1 );
260 for ( itU = theParams.rbegin(), i = 1; i < nPar; itU++, i++ ) {
267 static bool computeParamByFunc(Adaptor3d_Curve& C3d, double first, double last,
268 double length, bool theReverse,
269 int nbSeg, Function& func,
270 list<double>& theParams)
272 OSD::SetSignal( true );
277 MESSAGE( "computeParamByFunc" );
279 int nbPnt = 1 + nbSeg;
280 vector<double> x(nbPnt, 0.);
282 if( !buildDistribution( func, 0.0, 1.0, nbSeg, x, 1E-4 ) )
285 MESSAGE( "Points:\n" );
287 for( int i=0; i<=nbSeg; i++ )
289 sprintf( buf, "%f\n", float(x[i] ) );
295 // apply parameters in range [0,1] to the space of the curve
296 double prevU = first;
303 for( int i = 1; i < nbSeg; i++ )
305 double curvLength = length * (x[i] - x[i-1]) * sign;
306 GCPnts_AbscissaPoint Discret( C3d, curvLength, prevU );
307 if ( !Discret.IsDone() )
309 double U = Discret.Parameter();
310 if ( U > first && U < last )
311 theParams.push_back( U );
319 //=============================================================================
323 //=============================================================================
324 bool StdMeshers_Regular_1D::computeInternalParameters(const TopoDS_Edge& theEdge,
325 list<double> & theParams,
326 const bool theReverse) const
331 Handle(Geom_Curve) Curve = BRep_Tool::Curve(theEdge, f, l);
332 GeomAdaptor_Curve C3d(Curve);
334 double length = EdgeLength(theEdge);
342 if ( _hypType == LOCAL_LENGTH )
344 // Local Length hypothesis
345 double nbseg = ceil(length / _value[ BEG_LENGTH_IND ]); // integer sup
347 nbseg = 1; // degenerated edge
348 eltSize = length / nbseg;
352 // Number Of Segments hypothesis
353 switch (_ivalue[ DISTR_TYPE_IND ])
355 case StdMeshers_NumberOfSegments::DT_Scale:
357 double scale = _value[ SCALE_FACTOR_IND ];
360 double alpha = pow( scale , 1.0 / (_ivalue[ NB_SEGMENTS_IND ] - 1));
361 double factor = (l - f) / (1 - pow( alpha,_ivalue[ NB_SEGMENTS_IND ]));
363 int i, NbPoints = 1 + _ivalue[ NB_SEGMENTS_IND ];
364 for ( i = 2; i < NbPoints; i++ )
366 double param = f + factor * (1 - pow(alpha, i - 1));
367 theParams.push_back( param );
372 case StdMeshers_NumberOfSegments::DT_TabFunc:
374 FunctionTable func(_vvalue[ TAB_FUNC_IND ], _ivalue[ CONV_MODE_IND ]);
375 return computeParamByFunc(C3d, f, l, length, theReverse,
376 _ivalue[ NB_SEGMENTS_IND ], func,
380 case StdMeshers_NumberOfSegments::DT_ExprFunc:
382 FunctionExpr func(_svalue[ EXPR_FUNC_IND ].c_str(), _ivalue[ CONV_MODE_IND ]);
383 return computeParamByFunc(C3d, f, l, length, theReverse,
384 _ivalue[ NB_SEGMENTS_IND ], func,
388 case StdMeshers_NumberOfSegments::DT_Regular:
389 eltSize = length / _ivalue[ NB_SEGMENTS_IND ];
396 GCPnts_UniformAbscissa Discret(C3d, eltSize, f, l);
397 if ( !Discret.IsDone() )
400 int NbPoints = Discret.NbPoints();
401 for ( int i = 2; i < NbPoints; i++ )
403 double param = Discret.Parameter(i);
404 theParams.push_back( param );
406 compensateError( eltSize, eltSize, f, l, length, C3d, theParams ); // for PAL9899
410 case BEG_END_LENGTH: {
412 // geometric progression: SUM(n) = ( a1 - an * q ) / ( 1 - q ) = length
414 double a1 = _value[ BEG_LENGTH_IND ];
415 double an = _value[ END_LENGTH_IND ];
416 double q = ( length - a1 ) / ( length - an );
418 double U1 = theReverse ? l : f;
419 double Un = theReverse ? f : l;
421 double eltSize = theReverse ? -a1 : a1;
423 // computes a point on a curve <C3d> at the distance <eltSize>
424 // from the point of parameter <param>.
425 GCPnts_AbscissaPoint Discret( C3d, eltSize, param );
426 if ( !Discret.IsDone() ) break;
427 param = Discret.Parameter();
428 if ( param > f && param < l )
429 theParams.push_back( param );
434 compensateError( a1, an, U1, Un, length, C3d, theParams );
438 case ARITHMETIC_1D: {
440 // arithmetic progression: SUM(n) = ( an - a1 + q ) * ( a1 + an ) / ( 2 * q ) = length
442 double a1 = _value[ BEG_LENGTH_IND ];
443 double an = _value[ END_LENGTH_IND ];
445 double q = ( an - a1 ) / ( 2 *length/( a1 + an ) - 1 );
446 int n = int( 1 + ( an - a1 ) / q );
448 double U1 = theReverse ? l : f;
449 double Un = theReverse ? f : l;
456 while ( n-- > 0 && eltSize * ( Un - U1 ) > 0 ) {
457 // computes a point on a curve <C3d> at the distance <eltSize>
458 // from the point of parameter <param>.
459 GCPnts_AbscissaPoint Discret( C3d, eltSize, param );
460 if ( !Discret.IsDone() ) break;
461 param = Discret.Parameter();
462 if ( param > f && param < l )
463 theParams.push_back( param );
468 compensateError( a1, an, U1, Un, length, C3d, theParams );
475 GCPnts_UniformDeflection Discret(C3d, _value[ DEFLECTION_IND ], true);
476 if ( !Discret.IsDone() )
479 int NbPoints = Discret.NbPoints();
480 for ( int i = 2; i < NbPoints; i++ )
482 double param = Discret.Parameter(i);
483 theParams.push_back( param );
495 //=============================================================================
499 //=============================================================================
501 bool StdMeshers_Regular_1D::Compute(SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
503 MESSAGE("StdMeshers_Regular_1D::Compute");
505 if ( _hypType == NONE )
508 SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
509 aMesh.GetSubMesh(aShape);
511 // quardatic mesh required?
512 SMESH_HypoFilter filter( SMESH_HypoFilter::HasName( "QuadraticMesh" ));
513 bool isQuadraticMesh = aMesh->GetHypothesis( aShape, filter, true );
515 const TopoDS_Edge & EE = TopoDS::Edge(aShape);
516 TopoDS_Edge E = TopoDS::Edge(EE.Oriented(TopAbs_FORWARD));
517 int shapeID = meshDS->ShapeToIndex( E );
520 Handle(Geom_Curve) Curve = BRep_Tool::Curve(E, f, l);
522 TopoDS_Vertex VFirst, VLast;
523 TopExp::Vertices(E, VFirst, VLast); // Vfirst corresponds to f and Vlast to l
525 ASSERT(!VFirst.IsNull());
526 SMDS_NodeIteratorPtr lid= aMesh.GetSubMesh(VFirst)->GetSubMeshDS()->GetNodes();
529 MESSAGE (" NO NODE BUILT ON VERTEX ");
532 const SMDS_MeshNode * idFirst = lid->next();
534 ASSERT(!VLast.IsNull());
535 lid=aMesh.GetSubMesh(VLast)->GetSubMeshDS()->GetNodes();
538 MESSAGE (" NO NODE BUILT ON VERTEX ");
541 const SMDS_MeshNode * idLast = lid->next();
545 list< double > params;
546 bool reversed = false;
547 if ( !_mainEdge.IsNull() )
548 reversed = aMesh.IsReversedInChain( EE, _mainEdge );
550 if ( ! computeInternalParameters( E, params, reversed ))
553 catch ( Standard_Failure ) {
557 // edge extrema (indexes : 1 & NbPoints) already in SMDS (TopoDS_Vertex)
558 // only internal nodes receive an edge position with param on curve
560 const SMDS_MeshNode * idPrev = idFirst;
562 for (list<double>::iterator itU = params.begin(); itU != params.end(); itU++)
565 gp_Pnt P = Curve->Value(param);
567 //Add the Node in the DataStructure
568 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
569 meshDS->SetNodeOnEdge(node, shapeID, param);
571 SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
572 meshDS->SetMeshElementOnShape(edge, shapeID);
575 SMDS_MeshEdge* edge = meshDS->AddEdge(idPrev, idLast);
576 meshDS->SetMeshElementOnShape(edge, shapeID);
580 // Edge is a degenerated Edge : We put n = 5 points on the edge.
582 BRep_Tool::Range(E, f, l);
583 double du = (l - f) / (NbPoints - 1);
584 //MESSAGE("************* Degenerated edge! *****************");
586 TopoDS_Vertex V1, V2;
587 TopExp::Vertices(E, V1, V2);
588 gp_Pnt P = BRep_Tool::Pnt(V1);
590 const SMDS_MeshNode * idPrev = idFirst;
591 for (int i = 2; i < NbPoints; i++)
593 double param = f + (i - 1) * du;
594 SMDS_MeshNode * node = meshDS->AddNode(P.X(), P.Y(), P.Z());
595 meshDS->SetNodeOnEdge(node, shapeID, param);
597 SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, node);
598 meshDS->SetMeshElementOnShape(edge, shapeID);
601 SMDS_MeshEdge * edge = meshDS->AddEdge(idPrev, idLast);
602 meshDS->SetMeshElementOnShape(edge, shapeID);
607 //=============================================================================
609 * See comments in SMESH_Algo.cxx
611 //=============================================================================
613 const list <const SMESHDS_Hypothesis *> & StdMeshers_Regular_1D::GetUsedHypothesis(
614 SMESH_Mesh & aMesh, const TopoDS_Shape & aShape)
616 _usedHypList.clear();
617 _usedHypList = GetAppliedHypothesis(aMesh, aShape); // copy
618 int nbHyp = _usedHypList.size();
622 // Check, if propagated from some other edge
623 if (aShape.ShapeType() == TopAbs_EDGE &&
624 aMesh.IsPropagatedHypothesis(aShape, _mainEdge))
626 // Propagation of 1D hypothesis from <aMainEdge> on this edge
627 //_usedHypList = GetAppliedHypothesis(aMesh, _mainEdge); // copy
628 // use a general method in order not to nullify _mainEdge
629 _usedHypList = SMESH_Algo::GetUsedHypothesis(aMesh, _mainEdge); // copy
630 nbHyp = _usedHypList.size();
635 TopTools_ListIteratorOfListOfShape ancIt( aMesh.GetAncestors( aShape ));
636 for (; ancIt.More(); ancIt.Next())
638 const TopoDS_Shape& ancestor = ancIt.Value();
639 _usedHypList = GetAppliedHypothesis(aMesh, ancestor); // copy
640 nbHyp = _usedHypList.size();
646 _usedHypList.clear(); //only one compatible hypothesis allowed
650 //=============================================================================
654 //=============================================================================
656 ostream & StdMeshers_Regular_1D::SaveTo(ostream & save)
661 //=============================================================================
665 //=============================================================================
667 istream & StdMeshers_Regular_1D::LoadFrom(istream & load)
672 //=============================================================================
676 //=============================================================================
678 ostream & operator <<(ostream & save, StdMeshers_Regular_1D & hyp)
680 return hyp.SaveTo( save );
683 //=============================================================================
687 //=============================================================================
689 istream & operator >>(istream & load, StdMeshers_Regular_1D & hyp)
691 return hyp.LoadFrom( load );