1 // Copyright (c) 1996-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
15 // Jeannine PANCIATICI le 06/06/96
16 // Igor FEOKTISTOV 14/12/98 - correction of Approximate() and Init().
17 // Approximation d une MultiLine de points decrite par le tool MLineTool.
18 // avec criteres variationnels
20 #include <AppDef_MultiLine.hxx>
21 #include <AppDef_SmoothCriterion.hxx>
22 #include <AppDef_Variational.hxx>
23 #include <AppParCurves_MultiBSpCurve.hxx>
24 #include <FEmTool_Assembly.hxx>
25 #include <FEmTool_Curve.hxx>
26 #include <gp_VectorWithNullMagnitude.hxx>
27 #include <math_Matrix.hxx>
28 #include <PLib_Base.hxx>
29 #include <Standard_ConstructionError.hxx>
30 #include <Standard_DimensionError.hxx>
31 #include <Standard_DomainError.hxx>
32 #include <Standard_OutOfRange.hxx>
33 #include <StdFail_NotDone.hxx>
35 #define No_Standard_RangeError
36 #define No_Standard_OutOfRange
37 #define No_Standard_DimensionError
38 #define No_Standard_ConstructionError
40 #include <Standard_Stream.hxx>
42 #include <AppParCurves.hxx>
43 #include <AppParCurves_Constraint.hxx>
44 #include <AppParCurves_HArray1OfConstraintCouple.hxx>
45 #include <AppParCurves_Array1OfMultiPoint.hxx>
46 #include <AppParCurves_MultiPoint.hxx>
47 #include <AppParCurves_MultiBSpCurve.hxx>
48 #include <AppDef_LinearCriteria.hxx>
49 #include <Convert_CompPolynomialToPoles.hxx>
51 #include <gp_Pnt2d.hxx>
53 #include <gp_Vec2d.hxx>
54 #include <TColgp_Array1OfPnt.hxx>
55 #include <TColgp_Array1OfPnt2d.hxx>
56 #include <TColgp_Array1OfVec.hxx>
57 #include <TColgp_Array1OfVec2d.hxx>
58 #include <TColStd_Array1OfInteger.hxx>
59 #include <TColStd_HArray1OfInteger.hxx>
60 #include <TColStd_Array1OfReal.hxx>
61 #include <TColStd_HArray1OfReal.hxx>
62 #include <TColStd_HArray2OfReal.hxx>
63 #include <StdFail_NotDone.hxx>
64 #include <Standard_SStream.hxx>
65 #include <Standard_NoSuchObject.hxx>
66 #include <Precision.hxx>
67 #include <AppDef_MyLineTool.hxx>
69 #include <TColStd_HArray2OfInteger.hxx>
70 #include <TColStd_Array2OfInteger.hxx>
71 #include <TColStd_Array2OfReal.hxx>
72 #include <FEmTool_Assembly.hxx>
73 #include <FEmTool_AssemblyTable.hxx>
74 #include <FEmTool_Curve.hxx>
75 #include <math_Matrix.hxx>
76 #include <math_Vector.hxx>
77 #include <PLib_Base.hxx>
78 #include <PLib_JacobiPolynomial.hxx>
79 #include <PLib_HermitJacobi.hxx>
80 #include <FEmTool_HAssemblyTable.hxx>
91 //=======================================================================
92 //function : AppDef_Variational
93 //purpose : Initialization of the fields.
94 //=======================================================================
96 AppDef_Variational::AppDef_Variational(const AppDef_MultiLine& SSP,
97 const Standard_Integer FirstPoint,
98 const Standard_Integer LastPoint,
99 const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
100 const Standard_Integer MaxDegree,
101 const Standard_Integer MaxSegment,
102 const GeomAbs_Shape Continuity,
103 const Standard_Boolean WithMinMax,
104 const Standard_Boolean WithCutting,
105 const Standard_Real Tolerance,
106 const Standard_Integer NbIterations):
108 myFirstPoint(FirstPoint),
109 myLastPoint(LastPoint),
110 myConstraints(TheConstraints),
111 myMaxDegree(MaxDegree),
112 myMaxSegment(MaxSegment),
113 myNbIterations(NbIterations),
114 myTolerance(Tolerance),
115 myContinuity(Continuity),
116 myWithMinMax(WithMinMax),
117 myWithCutting(WithCutting)
120 if (myMaxDegree < 1) Standard_DomainError::Raise();
121 myMaxDegree = Min (30, myMaxDegree);
123 if (myMaxSegment < 1) Standard_DomainError::Raise();
125 if (myWithMinMax != 0 && myWithMinMax !=1 ) Standard_DomainError::Raise();
126 if (myWithCutting != 0 && myWithCutting !=1 ) Standard_DomainError::Raise();
128 myIsOverConstr = Standard_False;
129 myIsCreated = Standard_False;
130 myIsDone = Standard_False;
131 switch (myContinuity) {
142 Standard_ConstructionError::Raise();
145 myNbP2d = AppDef_MyLineTool::NbP2d(SSP);
146 myNbP3d = AppDef_MyLineTool::NbP3d(SSP);
147 myDimension = 2 * myNbP2d + 3* myNbP3d ;
152 myKnots= new TColStd_HArray1OfReal(1,2);
153 myKnots->SetValue(1,0.);
154 myKnots->SetValue(2,1.);
158 mySmoothCriterion = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint);
159 myParameters = new TColStd_HArray1OfReal(myFirstPoint, myLastPoint);
160 myNbPoints=myLastPoint-myFirstPoint+1;
161 if (myNbPoints <= 0) Standard_ConstructionError::Raise();
163 myTabPoints= new TColStd_HArray1OfReal(1,myDimension*myNbPoints);
165 // Table of Points initialization
167 Standard_Integer ipoint,jp2d,jp3d,index;
168 TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
169 TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
174 for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
177 if(myNbP2d !=0 && myNbP3d ==0 ) {
178 AppDef_MyLineTool::Value(mySSP,ipoint,TabP2d);
180 for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
181 { P2d = TabP2d.Value(jp2d);
183 myTabPoints->SetValue(index++,P2d.X());
184 myTabPoints->SetValue(index++,P2d.Y());
187 if(myNbP3d !=0 && myNbP2d == 0 ) {
188 AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d);
190 for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
192 { P3d=TabP3d.Value(jp3d);
194 myTabPoints->SetValue(index++,P3d.X());
195 myTabPoints->SetValue(index++,P3d.Y());
196 myTabPoints->SetValue(index++,P3d.Z());
200 if(myNbP3d !=0 && myNbP2d != 0 ) {
201 AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d,TabP2d);
203 for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
205 { P3d=TabP3d.Value(jp3d);
207 myTabPoints->SetValue(index++,P3d.X());
208 myTabPoints->SetValue(index++,P3d.Y());
209 myTabPoints->SetValue(index++,P3d.Z());
212 for ( jp2d = 1 ; jp2d <= myNbP2d ; jp2d++)
214 { P2d=TabP2d.Value(jp2d);
216 myTabPoints->SetValue(index++,P2d.X());
217 myTabPoints->SetValue(index++,P2d.Y());
224 //=======================================================================
226 //purpose : Initializes the tables of constraints
227 // and verifies if the problem is not over-constrained.
228 // This method is used in the Create and the method SetConstraint.
229 //=======================================================================
231 void AppDef_Variational::Init()
234 Standard_Integer ipoint,jp2d,jp3d,index,jndex;
235 Standard_Integer CurMultyPoint;
236 TColgp_Array1OfVec TabV3d(1, Max(1,myNbP3d));
237 TColgp_Array1OfVec2d TabV2d(1, Max(1,myNbP2d));
238 TColgp_Array1OfVec TabV3dcurv(1, Max(1,myNbP3d));
239 TColgp_Array1OfVec2d TabV2dcurv(1, Max(1,myNbP2d));
244 myNbConstraints=myConstraints->Length();
245 if (myNbConstraints < 0) Standard_ConstructionError::Raise();
247 myTypConstraints = new TColStd_HArray1OfInteger(1,Max(1,2*myNbConstraints));
248 myTabConstraints = new TColStd_HArray1OfReal(1,Max(1,2*myDimension*myNbConstraints));
249 myTtheta = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));
250 myTfthet = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints));
254 // Table of types initialization
255 Standard_Integer iconstr;
262 AppParCurves_Constraint valcontr;
264 for ( iconstr = myConstraints->Lower() ; iconstr <= myConstraints->Upper() ; iconstr++)
266 ipoint=(myConstraints->Value(iconstr)).Index();
268 valcontr=(myConstraints->Value(iconstr)).Constraint();
270 case AppParCurves_NoConstraint:
271 CurMultyPoint -= myNbP3d * 6 + myNbP2d * 2;
273 case AppParCurves_PassPoint:
274 myTypConstraints->SetValue(index++,ipoint);
275 myTypConstraints->SetValue(index++,0);
277 if(myNbP2d !=0 ) jndex=jndex+4*myNbP2d;
278 if(myNbP3d !=0 ) jndex=jndex+6*myNbP3d;
280 case AppParCurves_TangencyPoint:
281 myTypConstraints->SetValue(index++,ipoint);
282 myTypConstraints->SetValue(index++,1);
284 if(myNbP2d !=0 && myNbP3d == 0 )
286 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False)
287 Standard_ConstructionError::Raise();
288 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
290 Vt2d=TabV2d.Value(jp2d);
292 myTabConstraints->SetValue(jndex++,Vt2d.X());
293 myTabConstraints->SetValue(jndex++,Vt2d.Y());
295 InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
298 if(myNbP3d !=0 && myNbP2d == 0)
300 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False)
301 Standard_ConstructionError::Raise();
302 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
304 Vt3d=TabV3d.Value(jp3d);
306 myTabConstraints->SetValue(jndex++,Vt3d.X());
308 myTabConstraints->SetValue(jndex++,Vt3d.Y());
310 myTabConstraints->SetValue(jndex++,Vt3d.Z());
312 InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
315 if(myNbP3d !=0 && myNbP2d != 0)
317 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False)
318 Standard_ConstructionError::Raise();
319 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
321 Vt3d=TabV3d.Value(jp3d);
323 myTabConstraints->SetValue(jndex++,Vt3d.X());
324 myTabConstraints->SetValue(jndex++,Vt3d.Y());
325 myTabConstraints->SetValue(jndex++,Vt3d.Z());
327 InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
330 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
332 Vt2d=TabV2d.Value(jp2d);
334 myTabConstraints->SetValue(jndex++,Vt2d.X());
335 myTabConstraints->SetValue(jndex++,Vt2d.Y());
337 InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
344 case AppParCurves_CurvaturePoint:
345 myTypConstraints->SetValue(index++,ipoint);
346 myTypConstraints->SetValue(index++,2);
348 if(myNbP2d !=0 && myNbP3d == 0)
350 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False )
351 Standard_ConstructionError::Raise();
352 if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV2dcurv) == Standard_False)
353 Standard_ConstructionError::Raise();
354 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
356 Vt2d=TabV2d.Value(jp2d);
358 Vc2d=TabV2dcurv.Value(jp2d);
359 if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
360 Standard_ConstructionError::Raise();
361 myTabConstraints->SetValue(jndex++,Vt2d.X());
362 myTabConstraints->SetValue(jndex++,Vt2d.Y());
363 myTabConstraints->SetValue(jndex++,Vc2d.X());
364 myTabConstraints->SetValue(jndex++,Vc2d.Y());
365 InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
369 if(myNbP3d !=0 && myNbP2d == 0 )
371 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False )
372 Standard_ConstructionError::Raise();
373 if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv) == Standard_False)
374 Standard_ConstructionError::Raise();
375 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
377 Vt3d=TabV3d.Value(jp3d);
379 Vc3d=TabV3dcurv.Value(jp3d);
380 if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False)
381 Standard_ConstructionError::Raise();
382 myTabConstraints->SetValue(jndex++,Vt3d.X());
383 myTabConstraints->SetValue(jndex++,Vt3d.Y());
384 myTabConstraints->SetValue(jndex++,Vt3d.Z());
385 myTabConstraints->SetValue(jndex++,Vc3d.X());
386 myTabConstraints->SetValue(jndex++,Vc3d.Y());
387 myTabConstraints->SetValue(jndex++,Vc3d.Z());
388 InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
391 if(myNbP3d !=0 && myNbP2d != 0 )
393 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False )
394 Standard_ConstructionError::Raise();
395 if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv,TabV2dcurv) == Standard_False)
396 Standard_ConstructionError::Raise();
397 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
399 Vt3d=TabV3d.Value(jp3d);
401 Vc3d=TabV3dcurv.Value(jp3d);
402 if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False)
403 Standard_ConstructionError::Raise();
404 myTabConstraints->SetValue(jndex++,Vt3d.X());
405 myTabConstraints->SetValue(jndex++,Vt3d.Y());
406 myTabConstraints->SetValue(jndex++,Vt3d.Z());
407 myTabConstraints->SetValue(jndex++,Vc3d.X());
408 myTabConstraints->SetValue(jndex++,Vc3d.Y());
409 myTabConstraints->SetValue(jndex++,Vc3d.Z());
410 InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
412 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
414 Vt2d=TabV2d.Value(jp2d);
416 Vc2d=TabV2dcurv.Value(jp2d);
417 if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
418 Standard_ConstructionError::Raise();
419 myTabConstraints->SetValue(jndex++,Vt2d.X());
420 myTabConstraints->SetValue(jndex++,Vt2d.Y());
421 myTabConstraints->SetValue(jndex++,Vc2d.X());
422 myTabConstraints->SetValue(jndex++,Vc2d.Y());
423 InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
429 Standard_ConstructionError::Raise();
431 CurMultyPoint += myNbP3d * 6 + myNbP2d * 2;
433 // OverConstraint Detection
434 Standard_Integer MaxSeg;
435 if(myWithCutting == Standard_True) MaxSeg = myMaxSegment ;
437 if (((myMaxDegree-myNivCont)*MaxSeg-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
439 myIsOverConstr =Standard_True;
440 myIsCreated = Standard_False;
443 InitSmoothCriterion();
444 myIsCreated = Standard_True;
450 //=======================================================================
451 //function : Approximate
452 //purpose : Makes the approximation with the current fields.
453 //=======================================================================
455 void AppDef_Variational::Approximate()
458 if (myIsCreated == Standard_False ) StdFail_NotDone:: Raise();
461 Standard_Real WQuadratic, WQuality;
463 TColStd_Array1OfReal Ecarts(myFirstPoint, myLastPoint);
465 mySmoothCriterion->GetWeight(WQuadratic, WQuality);
467 Handle(FEmTool_Curve) TheCurve;
469 mySmoothCriterion->GetCurve(TheCurve);
471 //---------------------------------------------------------------------
473 TheMotor(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
476 if(myWithMinMax && myTolerance < myMaxError)
477 Adjusting(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
479 //---------------------------------------------------------------------
481 Standard_Integer jp2d,jp3d,index,ipole,
482 NbElem = TheCurve->NbElements();
484 TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
485 TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
486 Standard_Real debfin[2] = {-1., 1};
494 Handle(TColStd_HArray2OfReal) PolynomialIntervalsPtr =
495 new TColStd_HArray2OfReal(1,NbElem,1,2) ;
497 Handle(TColStd_HArray1OfInteger) NbCoeffPtr =
498 new TColStd_HArray1OfInteger(1, myMaxSegment);
500 Standard_Integer size = myMaxSegment * (myMaxDegree + 1) * myDimension ;
501 Handle(TColStd_HArray1OfReal) CoeffPtr = new TColStd_HArray1OfReal(1,size);
505 Handle(TColStd_HArray1OfReal) IntervallesPtr =
506 new TColStd_HArray1OfReal(1, NbElem + 1);
508 IntervallesPtr->ChangeArray1() = TheCurve->Knots();
510 TheCurve->GetPolynom(CoeffPtr->ChangeArray1());
514 for(ii = 1; ii <= NbElem; ii++)
515 NbCoeffPtr->SetValue(ii, TheCurve->Degree(ii)+1);
518 for (ii = PolynomialIntervalsPtr->LowerRow() ;
519 ii <= PolynomialIntervalsPtr->UpperRow() ;ii++)
521 PolynomialIntervalsPtr->SetValue(ii,1,debfin[0]) ;
522 PolynomialIntervalsPtr->SetValue(ii,2,debfin[1]) ;
525 printf("\n =========== Parameters for converting\n");
526 printf("nb_courbes : %d \n", NbElem);
527 printf("tab_option[4] : %d \n", myNivCont);
528 printf("myDimension : %d \n", myDimension);
529 printf("myMaxDegree : %d \n", myMaxDegree);
530 printf("\n NbcoeffPtr :\n");
531 for(ii = 1; ii <= NbElem; ii++) printf("NbCoeffPtr(%d) = %d \n", ii, NbCoeffPtr->Value(ii));
532 size = NbElem*(myMaxDegree + 1) * myDimension;
533 printf("\n CoeffPtr :\n");
534 for(ii = 1; ii <= size; ii++) printf("CoeffPtr(%d) = %.8e \n", ii, CoeffPtr->Value(ii));
535 printf("\n PolinimialIntervalsPtr :\n");
536 for (ii = PolynomialIntervalsPtr->LowerRow() ;
537 ii <= PolynomialIntervalsPtr->UpperRow() ;ii++)
539 printf(" %d : [%f, %f] \n", ii, PolynomialIntervalsPtr->Value(ii,1),
540 PolynomialIntervalsPtr->Value(ii,2)) ;
542 printf("\n IntervallesPtr :\n");
543 for (ii = IntervallesPtr->Lower();
544 ii <= IntervallesPtr->Upper() - 1; ii++)
546 printf(" %d : [%f, %f] \n", ii, IntervallesPtr->Value(ii),
547 IntervallesPtr->Value(ii+1)) ;
550 Convert_CompPolynomialToPoles AConverter(NbElem,
556 PolynomialIntervalsPtr,
558 if (AConverter.IsDone())
560 Handle(TColStd_HArray2OfReal) PolesPtr ;
561 Handle(TColStd_HArray1OfInteger) Mults;
562 Standard_Integer NbPoles=AConverter.NbPoles();
563 // Standard_Integer Deg=AConverter.Degree();
564 AppParCurves_Array1OfMultiPoint TabMU(1,NbPoles);
565 AConverter.Poles(PolesPtr) ;
566 AConverter.Knots(myKnots) ;
567 AConverter.Multiplicities(Mults) ;
569 for (ipole=PolesPtr->LowerRow();ipole<=PolesPtr->UpperRow();ipole++)
571 index=PolesPtr->LowerCol();
574 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
576 P2d.SetX(PolesPtr->Value(ipole,index++));
577 P2d.SetY(PolesPtr->Value(ipole,index++));
578 TabP2d.SetValue(jp2d,P2d);
583 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
585 // cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
586 P3d.SetX(PolesPtr->Value(ipole,index++));
587 // cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
588 P3d.SetY(PolesPtr->Value(ipole,index++));
589 // cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
590 P3d.SetZ(PolesPtr->Value(ipole,index++));
591 TabP3d.SetValue(jp3d,P3d);
596 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
598 P2d.SetX(PolesPtr->Value(ipole,index++));
599 P2d.SetY(PolesPtr->Value(ipole,index++));
600 TabP2d.SetValue(jp2d,P2d);
603 if(myNbP2d !=0 && myNbP3d !=0)
605 AppParCurves_MultiPoint aMultiPoint(TabP3d,TabP2d);
606 TabMU.SetValue(ipole,aMultiPoint);
608 else if (myNbP2d !=0)
610 AppParCurves_MultiPoint aMultiPoint(TabP2d);
611 TabMU.SetValue(ipole,aMultiPoint);
616 AppParCurves_MultiPoint aMultiPoint(TabP3d);
617 TabMU.SetValue(ipole,aMultiPoint);
621 AppParCurves_MultiBSpCurve aCurve(TabMU,myKnots->Array1(),Mults->Array1());
623 myIsDone = Standard_True;
629 //=======================================================================
630 //function : IsCreated
631 //purpose : returns True if the creation is done
632 //=======================================================================
634 Standard_Boolean AppDef_Variational::IsCreated() const
639 //=======================================================================
641 //purpose : returns True if the approximation is ok
642 //=======================================================================
644 Standard_Boolean AppDef_Variational::IsDone() const
649 //=======================================================================
650 //function : IsOverConstrained
651 //purpose : returns True if the problem is overconstrained
652 // in this case, approximation cannot be done.
653 //=======================================================================
655 Standard_Boolean AppDef_Variational::IsOverConstrained() const
657 return myIsOverConstr;
660 //=======================================================================
662 //purpose : returns all the BSpline curves approximating the
663 // MultiLine SSP after minimization of the parameter.
665 //=======================================================================
667 AppParCurves_MultiBSpCurve AppDef_Variational::Value() const
669 if (myIsDone == Standard_False) StdFail_NotDone::Raise();
674 //=======================================================================
675 //function : MaxError
676 //purpose : returns the maximum of the distances between
677 // the points of the multiline and the approximation
679 //=======================================================================
681 Standard_Real AppDef_Variational::MaxError() const
683 if (myIsDone == Standard_False) StdFail_NotDone::Raise();
687 //=======================================================================
688 //function : MaxErrorIndex
689 //purpose : returns the index of the MultiPoint of ErrorMax
690 //=======================================================================
692 Standard_Integer AppDef_Variational::MaxErrorIndex() const
694 if (myIsDone == Standard_False) StdFail_NotDone::Raise();
695 return myMaxErrorIndex;
698 //=======================================================================
699 //function : QuadraticError
700 //purpose : returns the quadratic average of the distances between
701 // the points of the multiline and the approximation
703 //=======================================================================
705 Standard_Real AppDef_Variational::QuadraticError() const
707 if (myIsDone == Standard_False) StdFail_NotDone::Raise();
708 return myCriterium[0];
711 //=======================================================================
712 //function : Distance
713 //purpose : returns the distances between the points of the
714 // multiline and the approximation curves.
715 //=======================================================================
717 void AppDef_Variational::Distance(math_Matrix& mat)
720 if (myIsDone == Standard_False) StdFail_NotDone::Raise();
721 Standard_Integer ipoint,jp2d,jp3d,index;
722 TColgp_Array1OfPnt TabP3d(1,Max(1,myNbP3d));
723 TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
724 Standard_Integer j0 = mat.LowerCol() - myFirstPoint;
733 for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
737 AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d);
739 for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
741 { P3d=TabP3d.Value(jp3d);
742 myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt3d);
743 mat(index++, j0 + ipoint)=P3d.Distance(Pt3d);
748 if(myNbP3d == 0 ) AppDef_MyLineTool::Value(mySSP,ipoint,TabP2d);
749 else AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d,TabP2d);
750 for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
752 { P2d = TabP2d.Value(jp2d);
753 myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt2d);
754 mat(index++, j0 + ipoint)=P2d.Distance(Pt2d);
761 //=======================================================================
762 //function : AverageError
763 //purpose : returns the average error between
764 // the MultiLine and the approximation.
765 //=======================================================================
767 Standard_Real AppDef_Variational::AverageError() const
769 if (myIsDone == Standard_False) StdFail_NotDone::Raise();
770 return myAverageError;
773 //=======================================================================
774 //function : Parameters
775 //purpose : returns the parameters uses to the approximations
776 //=======================================================================
778 const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Parameters() const
780 if (myIsDone == Standard_False) StdFail_NotDone::Raise();
784 //=======================================================================
786 //purpose : returns the knots uses to the approximations
787 //=======================================================================
789 const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Knots() const
791 if (myIsDone == Standard_False) StdFail_NotDone::Raise();
795 //=======================================================================
796 //function : Criterium
797 //purpose : returns the values of the quality criterium.
798 //=======================================================================
800 void AppDef_Variational::Criterium(Standard_Real& VFirstOrder, Standard_Real& VSecondOrder, Standard_Real& VThirdOrder) const
802 if (myIsDone == Standard_False) StdFail_NotDone::Raise();
803 VFirstOrder=myCriterium[1] ;
804 VSecondOrder=myCriterium[2];
805 VThirdOrder=myCriterium[3];
808 //=======================================================================
809 //function : CriteriumWeight
810 //purpose : returns the Weights (as percent) associed to the criterium used in
812 //=======================================================================
814 void AppDef_Variational::CriteriumWeight(Standard_Real& Percent1, Standard_Real& Percent2, Standard_Real& Percent3) const
816 Percent1 = myPercent[0];
817 Percent2 = myPercent[1];
818 Percent3 = myPercent[2] ;
821 //=======================================================================
822 //function : MaxDegree
823 //purpose : returns the Maximum Degree used in the approximation
824 //=======================================================================
826 Standard_Integer AppDef_Variational::MaxDegree() const
831 //=======================================================================
832 //function : MaxSegment
833 //purpose : returns the Maximum of segment used in the approximation
834 //=======================================================================
836 Standard_Integer AppDef_Variational::MaxSegment() const
841 //=======================================================================
842 //function : Continuity
843 //purpose : returns the Continuity used in the approximation
844 //=======================================================================
846 GeomAbs_Shape AppDef_Variational::Continuity() const
851 //=======================================================================
852 //function : WithMinMax
853 //purpose : returns if the approximation search to minimize the
854 // maximum Error or not.
855 //=======================================================================
857 Standard_Boolean AppDef_Variational::WithMinMax() const
862 //=======================================================================
863 //function : WithCutting
864 //purpose : returns if the approximation can insert new Knots or not.
865 //=======================================================================
867 Standard_Boolean AppDef_Variational::WithCutting() const
869 return myWithCutting;
872 //=======================================================================
873 //function : Tolerance
874 //purpose : returns the tolerance used in the approximation.
875 //=======================================================================
877 Standard_Real AppDef_Variational::Tolerance() const
882 //=======================================================================
883 //function : NbIterations
884 //purpose : returns the number of iterations used in the approximation.
885 //=======================================================================
887 Standard_Integer AppDef_Variational::NbIterations() const
889 return myNbIterations;
892 //=======================================================================
894 //purpose : Prints on the stream o information on the current state
896 //=======================================================================
898 void AppDef_Variational::Dump(Standard_OStream& o) const
900 o << " \nVariational Smoothing " << endl;
901 o << " Number of multipoints " << myNbPoints << endl;
902 o << " Number of 2d par multipoint " << myNbP2d << endl;
903 o << " Nombre of 3d par multipoint " << myNbP3d << endl;
904 o << " Number of PassagePoint " << myNbPassPoints << endl;
905 o << " Number of TangencyPoints " << myNbTangPoints << endl;
906 o << " Number of CurvaturePoints " << myNbCurvPoints << endl;
907 o << " \nTolerance " << o.setf(ios::scientific) << setprecision(3) << setw(9) << myTolerance;
908 if ( WithMinMax()) { o << " as Max Error." << endl;}
909 else { o << " as size Error." << endl;}
910 o << "CriteriumWeights : " << myPercent[0] << " , "
911 << myPercent[1] << " , " << myPercent[2] << endl;
914 o << " MaxError " << setprecision(3) << setw(9) << myMaxError << endl;
915 o << " Index of MaxError " << myMaxErrorIndex << endl;
916 o << " Average Error " << setprecision(3) << setw(9) << myAverageError << endl;
917 o << " Quadratic Error " << setprecision(3) << setw(9) << myCriterium[0] << endl;
918 o << " Tension Criterium " << setprecision(3) << setw(9) << myCriterium[1] << endl;
919 o << " Flexion Criterium " << setprecision(3) << setw(9) << myCriterium[2] << endl;
920 o << " Jerk Criterium " << setprecision(3) << setw(9) << myCriterium[3] << endl;
921 o << " NbSegments " << myKnots->Length()-1 << endl;
924 { if (myIsOverConstr) o << "The probleme is overconstraint " << endl;
925 else o << " Erreur dans l''approximation" << endl;
929 //=======================================================================
930 //function : SetConstraints
931 //purpose : Define the constraints to approximate
932 // If this value is incompatible with the others fields
933 // this method modify nothing and returns false
934 //=======================================================================
936 Standard_Boolean AppDef_Variational::SetConstraints(const Handle(AppParCurves_HArray1OfConstraintCouple)& aConstraint)
939 myConstraints=aConstraint;
941 if (myIsOverConstr ) return Standard_False;
942 else return Standard_True;
945 //=======================================================================
946 //function : SetParameters
947 //purpose : Defines the parameters used by the approximations.
948 //=======================================================================
950 void AppDef_Variational::SetParameters(const Handle(TColStd_HArray1OfReal)& param)
952 myParameters->ChangeArray1() = param->Array1();
955 //=======================================================================
956 //function : SetKnots
957 //purpose : Defines the knots used by the approximations
958 // -- If this value is incompatible with the others fields
959 // -- this method modify nothing and returns false
960 //=======================================================================
962 Standard_Boolean AppDef_Variational::SetKnots(const Handle(TColStd_HArray1OfReal)& knots)
964 myKnots->ChangeArray1() = knots->Array1();
965 return Standard_True;
968 //=======================================================================
969 //function : SetMaxDegree
970 //purpose : Define the Maximum Degree used in the approximation
971 // If this value is incompatible with the others fields
972 // this method modify nothing and returns false
973 //=======================================================================
975 Standard_Boolean AppDef_Variational::SetMaxDegree(const Standard_Integer Degree)
977 if (((Degree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
978 return Standard_False;
983 InitSmoothCriterion();
985 return Standard_True;
990 //=======================================================================
991 //function : SetMaxSegment
992 //purpose : Define the maximum number of segments used in the approximation
993 // If this value is incompatible with the others fields
994 // this method modify nothing and returns false
995 //=======================================================================
997 Standard_Boolean AppDef_Variational::SetMaxSegment(const Standard_Integer NbSegment)
999 if ( myWithCutting == Standard_True &&
1000 ((myMaxDegree-myNivCont)*NbSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1001 return Standard_False;
1004 myMaxSegment=NbSegment;
1005 return Standard_True;
1009 //=======================================================================
1010 //function : SetContinuity
1011 //purpose : Define the Continuity used in the approximation
1012 // If this value is incompatible with the others fields
1013 // this method modify nothing and returns false
1014 //=======================================================================
1016 Standard_Boolean AppDef_Variational::SetContinuity(const GeomAbs_Shape C)
1018 Standard_Integer NivCont=0;
1030 Standard_ConstructionError::Raise();
1032 if (((myMaxDegree-NivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1033 return Standard_False;
1039 InitSmoothCriterion();
1040 return Standard_True;
1044 //=======================================================================
1045 //function : SetWithMinMax
1046 //purpose : Define if the approximation search to minimize the
1047 // maximum Error or not.
1048 //=======================================================================
1050 void AppDef_Variational::SetWithMinMax(const Standard_Boolean MinMax)
1052 myWithMinMax=MinMax;
1054 InitSmoothCriterion();
1057 //=======================================================================
1058 //function : SetWithCutting
1059 //purpose : Define if the approximation can insert new Knots or not.
1060 // If this value is incompatible with the others fields
1061 // this method modify nothing and returns false
1062 //=======================================================================
1064 Standard_Boolean AppDef_Variational::SetWithCutting(const Standard_Boolean Cutting)
1066 if (Cutting == Standard_False)
1068 if (((myMaxDegree-myNivCont)*myKnots->Length()-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1069 return Standard_False;
1073 myWithCutting=Cutting;
1074 InitSmoothCriterion();
1075 return Standard_True;
1080 if (((myMaxDegree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1081 return Standard_False;
1085 myWithCutting=Cutting;
1086 InitSmoothCriterion();
1087 return Standard_True;
1092 //=======================================================================
1093 //function : SetCriteriumWeight
1094 //purpose : define the Weights (as percent) associed to the criterium used in
1095 // the optimization.
1096 //=======================================================================
1098 void AppDef_Variational::SetCriteriumWeight(const Standard_Real Percent1, const Standard_Real Percent2, const Standard_Real Percent3)
1100 if (Percent1 < 0 || Percent2 < 0 || Percent3 < 0 ) Standard_DomainError::Raise();
1101 Standard_Real Total = Percent1 + Percent2 + Percent3;
1102 myPercent[0] = Percent1/Total;
1103 myPercent[1] = Percent2/Total;
1104 myPercent[2] = Percent3/Total;
1106 InitSmoothCriterion();
1109 //=======================================================================
1110 //function : SetCriteriumWeight
1111 //purpose : define the Weight (as percent) associed to the
1112 // criterium Order used in the optimization : Others
1113 // weights are updated.
1114 //=======================================================================
1116 void AppDef_Variational::SetCriteriumWeight(const Standard_Integer Order, const Standard_Real Percent)
1118 if ( Percent < 0 ) Standard_DomainError::Raise();
1119 if ( Order < 1 || Order > 3 ) Standard_ConstructionError::Raise();
1120 myPercent[Order-1] = Percent;
1121 Standard_Real Total = myPercent[0] + myPercent[1] + myPercent[2];
1122 myPercent[0] = myPercent[0] / Total;
1123 myPercent[1] = myPercent[1] / Total;
1124 myPercent[2] = myPercent[2] / Total;
1126 InitSmoothCriterion();
1130 //=======================================================================
1131 //function : SetTolerance
1132 //purpose : define the tolerance used in the approximation.
1133 //=======================================================================
1135 void AppDef_Variational::SetTolerance(const Standard_Real Tol)
1138 InitSmoothCriterion();
1141 //=======================================================================
1142 //function : SetNbIterations
1143 //purpose : define the number of iterations used in the approximation.
1144 //=======================================================================
1146 void AppDef_Variational::SetNbIterations(const Standard_Integer Iter)
1148 myNbIterations=Iter;
1151 //====================== Private Methodes =============================//
1152 //=======================================================================
1153 //function : TheMotor
1154 //purpose : Smoothing's motor like STRIM routine "MOTLIS"
1155 //=======================================================================
1156 void AppDef_Variational::TheMotor(
1157 Handle(AppDef_SmoothCriterion)& J,
1158 // const Standard_Real WQuadratic,
1159 const Standard_Real ,
1160 const Standard_Real WQuality,
1161 Handle(FEmTool_Curve)& TheCurve,
1162 TColStd_Array1OfReal& Ecarts)
1166 const Standard_Real BigValue = 1.e37, SmallValue = 1.e-6, SmallestValue = 1.e-9;
1168 Handle(TColStd_HArray1OfReal) CurrentTi, NewTi, OldTi;
1169 Handle(TColStd_HArray2OfInteger) Dependence;
1170 Standard_Boolean lestim, lconst, ToOptim, iscut;
1171 Standard_Boolean isnear = Standard_False, again = Standard_True;
1172 Standard_Integer NbEst, ICDANA, NumPnt, Iter;
1173 Standard_Integer MaxNbEst =5;
1174 Standard_Real VOCRI[3] = {BigValue, BigValue, BigValue}, EROLD = BigValue,
1175 VALCRI[3], ERRMAX = BigValue, ERRMOY, ERRQUA;
1176 Standard_Real CBLONG, LNOLD;
1177 Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
1178 Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1179 Handle(FEmTool_Curve) CCurrent, COld, CNew;
1180 Standard_Real EpsLength = SmallValue;
1181 Standard_Real EpsDeg;
1183 Standard_Real e1, e2, e3;
1184 Standard_Real J1min, J2min, J3min;
1185 Standard_Integer iprog;
1189 J->GetEstimation(e1, e2, e3);
1190 J1min = 1.e-8; J2min = J3min = (e1 + 1.e-8) * 1.e-6;
1192 if(e1 < J1min) e1 = J1min;// Like in
1193 if(e2 < J2min) e2 = J2min;// MOTLIS
1194 if(e3 < J3min) e3 = J3min;
1197 J->SetEstimation(e1, e2, e3);
1199 CCurrent = TheCurve;
1200 CurrentTi = new TColStd_HArray1OfReal(1, myParameters->Length());
1201 CurrentTi->ChangeArray1() = myParameters->Array1();
1202 OldTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1203 OldTi->ChangeArray1() = CurrentTi->Array1();
1205 LNOLD = CBLONG = J->EstLength();
1206 Dependence = J->DependenceTable();
1208 J->SetCurve(CCurrent);
1209 FEmTool_Assembly * TheAssembly =
1210 new FEmTool_Assembly (Dependence->Array2(), J->AssemblyTable());
1212 //============ Optimization ============================
1213 // Standard_Integer inagain = 0;
1216 // (1) Loop Optimization / Estimation
1217 lestim = Standard_True;
1218 lconst = Standard_True;
1221 J->SetCurve(CCurrent);
1225 // (1.1) Curve's Optimization.
1226 EpsLength = SmallValue * CBLONG / NbrPnt;
1227 CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
1228 CCurrent->NbElements(), CCurrent->Base(), EpsLength);
1229 CNew->Knots() = CCurrent->Knots();
1231 J->SetParameters(CurrentTi);
1232 EpsDeg = Min(WQuality * .1, CBLONG * .001);
1234 Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
1236 lconst = Standard_False;
1238 // (1.2) calcul des criteres de qualites et amelioration
1240 ICDANA = J->QualityValues(J1min, J2min, J3min,
1241 VALCRI[0], VALCRI[1], VALCRI[2]);
1243 if(ICDANA > 0) lconst = Standard_True;
1245 J->ErrorValues(ERRMAX, ERRQUA, ERRMOY);
1247 isnear = ((Sqrt(ERRQUA / NbrPnt) < 2*WQuality) &&
1248 (myNbIterations > 1));
1250 // (1.3) Optimisation des ti par proj orthogonale
1251 // et calcul de l'erreur aux points.
1254 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1255 Project(CNew, CurrentTi->Array1(),
1256 NewTi->ChangeArray1(),
1258 ERRMAX, ERRQUA, ERRMOY, 2);
1260 else NewTi = CurrentTi;
1263 // (1.4) Progression's test
1265 if ((EROLD > WQuality) && (ERRMAX < 0.95*EROLD)) iprog++;
1266 if ((EROLD > WQuality) && (ERRMAX < 0.8*EROLD)) iprog++;
1267 if ((EROLD > WQuality) && (ERRMAX < WQuality)) iprog++;
1268 if ((EROLD > WQuality) && (ERRMAX < 0.99*EROLD)
1269 && (ERRMAX < 1.1*WQuality)) iprog++;
1270 if ( VALCRI[0] < 0.975 * VOCRI[0]) iprog++;
1271 if ( VALCRI[0] < 0.9 * VOCRI[0]) iprog++;
1272 if ( VALCRI[1] < 0.95 * VOCRI[1]) iprog++;
1273 if ( VALCRI[1] < 0.8 * VOCRI[1]) iprog++;
1274 if ( VALCRI[2] < 0.95 * VOCRI[2]) iprog++;
1275 if ( VALCRI[2] < 0.8 * VOCRI[2]) iprog++;
1276 if ((VOCRI[1]>SmallestValue)&&(VOCRI[2]>SmallestValue)) {
1277 if ((VALCRI[1]/VOCRI[1] + 2*VALCRI[2]/VOCRI[2]) < 2.8) iprog++;
1280 if (iprog < 2 && NbEst == 0 ) {
1281 // (1.5) Invalidation of new knots.
1282 VALCRI[0] = VOCRI[0];
1283 VALCRI[1] = VOCRI[1];
1284 VALCRI[2] = VOCRI[2];
1293 VOCRI[0] = VALCRI[0];
1294 VOCRI[1] = VALCRI[1];
1295 VOCRI[2] = VALCRI[2];
1302 // (1.6) Test if the Estimations seems OK, else repeat
1304 lestim = ( (NbEst<MaxNbEst) && (ICDANA == 2)&& (iprog > 0) );
1306 if (lestim && isnear) {
1307 // (1.7) Optimization of ti by ACR.
1309 std::stable_sort (CurrentTi->begin(), CurrentTi->end());
1311 Standard_Integer Decima = 4;
1313 CCurrent->Length(0., 1., CBLONG);
1314 J->EstLength() = CBLONG;
1316 ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
1317 lconst = Standard_True;
1323 // (2) loop of parametric / geometric optimization
1326 ToOptim = ((Iter < myNbIterations) && (isnear));
1330 // (2.1) Save curent results
1331 VOCRI[0] = VALCRI[0];
1332 VOCRI[1] = VALCRI[1];
1333 VOCRI[2] = VALCRI[2];
1337 OldTi->ChangeArray1() = CurrentTi->Array1();
1339 // (2.2) Optimization des ti by ACR.
1341 std::stable_sort (CurrentTi->begin(), CurrentTi->end());
1343 Standard_Integer Decima = 4;
1345 CCurrent->Length(0., 1., CBLONG);
1346 J->EstLength() = CBLONG;
1348 ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
1349 lconst = Standard_True;
1351 // (2.3) Optimisation des courbes
1352 EpsLength = SmallValue * CBLONG / NbrPnt;
1354 CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
1355 CCurrent->NbElements(), CCurrent->Base(), EpsLength);
1356 CNew->Knots() = CCurrent->Knots();
1358 J->SetParameters(CurrentTi);
1360 EpsDeg = Min(WQuality * .1, CBLONG * .001);
1361 Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
1365 // (2.4) calcul des criteres de qualites et amelioration
1368 ICDANA = J->QualityValues(J1min, J2min, J3min, VALCRI[0], VALCRI[1], VALCRI[2]);
1369 if(ICDANA > 0) lconst = Standard_True;
1371 J->GetEstimation(e1, e2, e3);
1372 // (2.5) Optimisation des ti par proj orthogonale
1374 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1375 Project(CCurrent, CurrentTi->Array1(),
1376 NewTi->ChangeArray1(),
1378 ERRMAX, ERRQUA, ERRMOY, 2);
1380 // (2.6) Test de non regression
1382 Standard_Integer iregre = 0;
1383 if (NbrConstraint < NbrPnt) {
1384 if ( (ERRMAX > WQuality) && (ERRMAX > 1.05*EROLD)) iregre++;
1385 if ( (ERRMAX > WQuality) && (ERRMAX > 2*EROLD)) iregre++;
1386 if ( (EROLD > WQuality) && (ERRMAX <= 0.5*EROLD)) iregre--;
1388 Standard_Real E1, E2, E3;
1389 J->GetEstimation(E1, E2, E3);
1390 if ( (VALCRI[0] > E1) && (VALCRI[0] > 1.1*VOCRI[0])) iregre++;
1391 if ( (VALCRI[1] > E2) && (VALCRI[1] > 1.1*VOCRI[1])) iregre++;
1392 if ( (VALCRI[2] > E3) && (VALCRI[2] > 1.1*VOCRI[2])) iregre++;
1396 // if (iregre >= 1) {
1397 // (2.7) on restaure l'iteration precedente
1398 VALCRI[0] = VOCRI[0];
1399 VALCRI[1] = VOCRI[1];
1400 VALCRI[2] = VOCRI[2];
1404 CurrentTi->ChangeArray1() = OldTi->Array1();
1405 ToOptim = Standard_False;
1407 else { // Iteration is Ok.
1411 if (Iter >= myNbIterations) ToOptim = Standard_False;
1414 // (3) Decoupe eventuelle
1416 if ( (CCurrent->NbElements() < myMaxSegment) && myWithCutting ) {
1418 // (3.1) Sauvgarde de l'etat precedent
1419 VOCRI[0] = VALCRI[0];
1420 VOCRI[1] = VALCRI[1];
1421 VOCRI[2] = VALCRI[2];
1424 OldTi->ChangeArray1() = CurrentTi->Array1();
1426 // (3.2) On arrange les ti : Trie + recadrage sur (0,1)
1427 // ---> On trie, afin d'assurer l'ordre par la suite.
1429 std::stable_sort (CurrentTi->begin(), CurrentTi->end());
1431 if ((CurrentTi->Value(1)!= 0.) ||
1432 (CurrentTi->Value(NbrPnt)!= 1.)) {
1433 Standard_Real t, DelatT =
1434 1.0 /(CurrentTi->Value(NbrPnt)-CurrentTi->Value(1));
1435 for (Standard_Integer ii=2; ii<NbrPnt; ii++) {
1436 t = (CurrentTi->Value(ii)-CurrentTi->Value(1))*DelatT;
1437 CurrentTi->SetValue(ii, t);
1439 CurrentTi->SetValue(1, 0.);
1440 CurrentTi->SetValue(NbrPnt, 1.);
1443 // (3.3) Insert new Knots
1445 SplitCurve(CCurrent, CurrentTi->Array1(), EpsLength, CNew, iscut);
1446 if (!iscut) again = Standard_False;
1449 // New Knots => New Assembly.
1452 TheAssembly = new FEmTool_Assembly (Dependence->Array2(),
1453 J->AssemblyTable());
1456 else { again = Standard_False;}
1459 // ================ Great loop end ===================
1463 // (4) Compute the best Error.
1464 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1465 Project(CCurrent, CurrentTi->Array1(),
1466 NewTi->ChangeArray1(),
1468 ERRMAX, ERRQUA, ERRMOY, 10);
1470 // (5) field's update
1472 TheCurve = CCurrent;
1473 J->EstLength() = CBLONG;
1474 myParameters->ChangeArray1() = NewTi->Array1();
1475 myCriterium[0] = ERRQUA;
1476 myCriterium[1] = Sqrt(VALCRI[0]);
1477 myCriterium[2] = Sqrt(VALCRI[1]);
1478 myCriterium[3] = Sqrt(VALCRI[2]);
1479 myMaxError = ERRMAX;
1480 myMaxErrorIndex = NumPnt;
1481 if(NbrPnt > NbrConstraint)
1482 myAverageError = ERRMOY / (NbrPnt - NbrConstraint);
1484 myAverageError = ERRMOY / NbrConstraint;
1489 //=======================================================================
1490 //function : Optimization
1491 //purpose : (like FORTRAN subroutine MINIMI)
1492 //=======================================================================
1493 void AppDef_Variational::Optimization(Handle(AppDef_SmoothCriterion)& J,
1494 FEmTool_Assembly& A,
1495 const Standard_Boolean ToAssemble,
1496 const Standard_Real EpsDeg,
1497 Handle(FEmTool_Curve)& Curve,
1498 const TColStd_Array1OfReal& Parameters) const
1500 Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
1501 NbElm = Curve->NbElements(),
1502 NbDim = Curve->Dimension();
1504 Handle(FEmTool_HAssemblyTable) AssTable;
1506 math_Matrix H(0, MxDeg, 0, MxDeg);
1507 math_Vector G(0, MxDeg), Sol(1, A.NbGlobVar());
1509 Standard_Integer el, dim;
1511 A.GetAssemblyTable(AssTable);
1512 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1514 Standard_Real CBLONG = J->EstLength();
1516 // Updating Assembly
1522 for (el = 1; el <= NbElm; el++) {
1524 J->Hessian(el, 1, 1, H);
1525 for(dim = 1; dim <= NbDim; dim++)
1526 A.AddMatrix(el, dim, dim, H);
1529 for(dim = 1; dim <= NbDim; dim++) {
1530 J->Gradient(el, dim, G);
1531 A.AddVector(el, dim, G);
1536 // Solution of system
1538 if(NbConstr != 0) { //Treatment of constraints
1539 AssemblingConstraints(Curve, Parameters, CBLONG, A);
1548 J->InputVector(Sol, AssTable);
1550 // Updating Curve and reduction of degree
1552 Standard_Integer Newdeg;
1553 Standard_Real MaxError;
1556 for(el = 1; el <= NbElm; el++) {
1557 Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
1562 TColStd_Array1OfReal& TabInt = Curve->Knots();
1563 Standard_Integer Icnt = 1, p0 = Parameters.Lower() - myFirstPoint, point;
1564 for(el = 1; el <= NbElm; el++) {
1565 while((Icnt < NbConstr) &&
1566 (Parameters(p0 + myTypConstraints->Value(2 * Icnt - 1)) <= TabInt(el))) Icnt++;
1567 point = p0 + myTypConstraints->Value(2 * Icnt - 1);
1568 if(Parameters(point) <= TabInt(el) || Parameters(point) >= TabInt(el + 1))
1569 Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
1571 if(Curve->Degree(el) < MxDeg) Curve->SetDegree(el, MxDeg);
1576 void AppDef_Variational::Project(const Handle(FEmTool_Curve)& C,
1577 const TColStd_Array1OfReal& Ti,
1578 TColStd_Array1OfReal& ProjTi,
1579 TColStd_Array1OfReal& Distance,
1580 Standard_Integer& NumPoints,
1581 Standard_Real& MaxErr,
1582 Standard_Real& QuaErr,
1583 Standard_Real& AveErr,
1584 const Standard_Integer NbIterations) const
1589 const Standard_Real Seuil = 1.e-9, Eps = 1.e-12;
1591 MaxErr = QuaErr = AveErr = 0.;
1593 Standard_Integer Ipnt, NItCv, Iter, i, i0 = -myDimension, d0 = Distance.Lower() - 1;
1595 Standard_Real TNew, Dist, T0, Dist0, F1, F2, Aux, DF, Ecart;
1597 Standard_Boolean EnCour;
1599 TColStd_Array1OfReal ValOfC(1, myDimension), FirstDerOfC(1, myDimension),
1600 SecndDerOfC(1, myDimension);
1602 for(Ipnt = 1; Ipnt <= ProjTi.Length(); Ipnt++) {
1608 EnCour = Standard_True;
1611 C->D0(TNew, ValOfC);
1614 for(i = 1; i <= myDimension; i++) {
1615 Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
1620 // ------- Newton's method for solving (C'(t),C(t) - P) = 0
1628 C->D2(TNew, SecndDerOfC);
1629 C->D1(TNew, FirstDerOfC);
1632 for(i = 1; i <= myDimension; i++) {
1633 Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
1634 DF = FirstDerOfC(i);
1635 F1 += Aux*DF; // (C'(t),C(t) - P)
1636 F2 += DF*DF + Aux * SecndDerOfC(i); // ((C'(t),C(t) - P))'
1640 EnCour = Standard_False;
1642 // Formula of Newton x(k+1) = x(k) - F(x(k))/F'(x(k))
1644 if(TNew < 0.) TNew = 0.;
1645 if(TNew > 1.) TNew = 1.;
1648 // Analysis of result
1650 C->D0(TNew, ValOfC);
1653 for(i = 1; i <= myDimension; i++) {
1654 Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
1659 Ecart = Dist0 - Dist;
1661 if(Ecart <= -Seuil) {
1662 // Pas d'amelioration on s'arrete
1663 EnCour = Standard_False;
1667 else if(Ecart <= Seuil)
1673 if((NItCv >= 2) || (Iter >= NbIterations)) EnCour = Standard_False;
1679 ProjTi(Ipnt) = TNew;
1680 Distance(d0 + Ipnt) = Dist;
1685 QuaErr += Dist * Dist;
1689 NumPoints = NumPoints + myFirstPoint - 1;// Setting NumPoints to interval [myFirstPoint, myLastPoint]
1693 void AppDef_Variational::ACR(Handle(FEmTool_Curve)& Curve,
1694 TColStd_Array1OfReal& Ti,
1695 const Standard_Integer Decima) const
1698 const Standard_Real Eps = 1.e-8;
1700 TColStd_Array1OfReal& Knots = Curve->Knots();
1701 Standard_Integer NbrPnt = Ti.Length(), TiFirst = Ti.Lower(), TiLast = Ti.Upper(),
1702 KFirst = Knots.Lower(), KLast = Knots.Upper();
1704 Standard_Real CbLong, DeltaT, VTest, UNew, UOld, DU, TPara, TOld, DTInv, Ratio;
1705 Standard_Integer ipnt, ii, IElm, IOld, POld, PCnt, ICnt=0;
1706 Standard_Integer NbCntr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1708 // (1) Calcul de la longueur de courbe
1710 Curve->Length(Ti(TiFirst), Ti(TiLast), CbLong);
1712 // (2) Mise de l'acr dans Ti
1716 // (2.0) Initialisation
1717 DeltaT = (Ti(TiLast) - Ti(TiFirst)) / Decima;
1718 VTest = Ti(TiFirst) + DeltaT;
1721 PCnt = myTypConstraints->Value(1) - myFirstPoint + TiFirst;
1737 for(ipnt = TiFirst + 1; ipnt <= TiLast; ipnt++) {
1739 while((ICnt <= NbCntr) && (PCnt < ipnt)) {
1741 PCnt = myTypConstraints->Value(2*ICnt-1) - myFirstPoint + TiFirst;
1746 if(TPara >= VTest || PCnt == ipnt) {
1748 if ( Ti(TiLast) - TPara <= 1.e-2*DeltaT) {
1752 // (2.2), (2.3) Cacul la longueur de courbe
1753 Curve->Length(Ti(TiFirst), TPara, UNew);
1757 while(Knots(IElm + 1) < TPara && IElm < KLast - 1) IElm++;
1759 // (2.4) Mise a jours des parametres de decoupe
1760 DTInv = 1. / (TPara - TOld);
1763 for(ii = IOld+1; ii <= IElm; ii++) {
1764 Ratio = (Knots(ii) - TOld) * DTInv;
1765 Knots(ii) = UOld + Ratio * DU;
1768 // (2.5) Mise a jours des parametres de points.
1770 //Very strange loop, because it never works (POld+1 > ipnt-1)
1771 for(ii = POld+1; ii <= ipnt-1; ii++) {
1772 Ratio = ( Ti(ii) - TOld ) * DTInv;
1773 Ti(ii) = UOld + Ratio * DU;
1784 // --> Nouveau seuil parametrique pour le decimage
1786 if(TPara >= VTest) {
1787 // ii = RealToInt((TPara - VTest + Eps) / DeltaT);
1788 // VTest += (ii + 1) * DeltaT;
1789 VTest += Ceiling((TPara - VTest + Eps) / DeltaT) * DeltaT;
1790 if(VTest > 1. - Eps) VTest = 1.;
1795 // --- On ajuste les valeurs extremes
1800 while ( Ti(ii) > Knots(KLast) ) {
1809 //----------------------------------------------------------//
1810 // Standard_Integer NearIndex //
1811 // Purpose: searching nearest index of TabPar corresponding //
1812 // given value T (is similar to fortran routine MSRRE2) //
1813 //----------------------------------------------------------//
1815 static Standard_Integer NearIndex(const Standard_Real T,
1816 const TColStd_Array1OfReal& TabPar,
1817 const Standard_Real Eps, Standard_Integer& Flag)
1819 Standard_Integer Loi = TabPar.Lower(), Upi = TabPar.Upper();
1823 if(T < TabPar(Loi)) { Flag = -1; return Loi;}
1824 if(T > TabPar(Upi)) { Flag = 1; return Upi;}
1826 Standard_Integer Ibeg = Loi, Ifin = Upi, Imidl;
1828 while(Ibeg + 1 != Ifin) {
1829 Imidl = (Ibeg + Ifin) / 2;
1830 if((T >= TabPar(Ibeg)) && (T <= TabPar(Imidl)))
1836 if(Abs(T - TabPar(Ifin)) < Eps) return Ifin;
1842 //----------------------------------------------------------//
1843 // void GettingKnots //
1844 // Purpose: calculating values of new Knots for elements //
1845 // with degree that is equal Deg //
1846 //----------------------------------------------------------//
1848 static void GettingKnots(const TColStd_Array1OfReal& TabPar,
1849 const Handle(FEmTool_Curve)& InCurve,
1850 const Standard_Integer Deg,
1851 Standard_Integer& NbElm,
1852 TColStd_Array1OfReal& NewKnots)
1856 const Standard_Real Eps = 1.e-12;
1858 TColStd_Array1OfReal& OldKnots = InCurve->Knots();
1859 Standard_Integer NbMaxOld = InCurve->NbElements();
1860 Standard_Integer NbMax = NewKnots.Upper(), Ipt, Ipt1, Ipt2;
1861 Standard_Integer el = 0, i1 = OldKnots.Lower(), i0 = i1 - 1, Flag;
1864 while((NbElm < NbMax) && (el < NbMaxOld)) {
1866 el++; i0++; i1++; // i0, i1 are indexes of left and right knots of element el
1868 if(InCurve->Degree(el) == Deg) {
1872 Ipt1 = NearIndex(OldKnots(i0), TabPar, Eps, Flag);
1873 if(Flag != 0) Ipt1 = TabPar.Lower();
1874 Ipt2 = NearIndex(OldKnots(i1), TabPar, Eps, Flag);
1875 if(Flag != 0) Ipt2 = TabPar.Upper();
1877 if(Ipt2 - Ipt1 >= 1) {
1879 Ipt = (Ipt1 + Ipt2) / 2;
1880 if(2 * Ipt == Ipt1 + Ipt2)
1881 TPar = 2. * TabPar(Ipt);
1883 TPar = TabPar(Ipt) + TabPar(Ipt + 1);
1885 NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1) + TPar) / 4.;
1888 NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1)) / 2.;
1893 void AppDef_Variational::SplitCurve(const Handle(FEmTool_Curve)& InCurve,
1894 const TColStd_Array1OfReal& Ti,
1895 const Standard_Real CurveTol,
1896 Handle(FEmTool_Curve)& OutCurve,
1897 Standard_Boolean& iscut) const
1899 Standard_Integer NbElmOld = InCurve->NbElements();
1901 if(NbElmOld >= myMaxSegment) {iscut = Standard_False; return;}
1903 Standard_Integer MaxDegree =
1905 InCurve->Base()->WorkDegree();
1906 Standard_Integer NbElm = NbElmOld;
1907 TColStd_Array1OfReal NewKnots(NbElm + 1, myMaxSegment);
1909 GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree(), NbElm, NewKnots);
1910 GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree() - 1, NbElm, NewKnots);
1912 GettingKnots(Ti, InCurve, MaxDegree, NbElm, NewKnots);
1913 GettingKnots(Ti, InCurve, MaxDegree - 1, NbElm, NewKnots);
1916 if(NbElm > NbElmOld) {
1918 iscut = Standard_True;
1920 OutCurve = new FEmTool_Curve(InCurve->Dimension(), NbElm, InCurve->Base(), CurveTol);
1921 TColStd_Array1OfReal& OutKnots = OutCurve->Knots();
1922 TColStd_Array1OfReal& InKnots = InCurve->Knots();
1924 Standard_Integer i, i0 = OutKnots.Lower();
1925 for(i = InKnots.Lower(); i <= InKnots.Upper(); i++) OutKnots(i) = InKnots(i);
1926 for(i = NbElmOld + 1; i <= NbElm; i++) OutKnots(i + i0) = NewKnots(i);
1928 std::sort (OutKnots.begin(), OutKnots.end());
1931 iscut = Standard_False;
1935 //=======================================================================
1936 //function : InitSmoothCriterion
1937 //purpose : Initializes the SmoothCriterion
1938 //=======================================================================
1939 void AppDef_Variational::InitSmoothCriterion()
1942 const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9;
1943 // const Standard_Real J1 = .01, J2 = .001, J3 = .001;
1947 Standard_Real Length;
1949 InitParameters(Length);
1951 mySmoothCriterion->SetParameters(myParameters);
1953 Standard_Real E1, E2, E3;
1955 InitCriterionEstimations(Length, E1, E2, E3);
1957 J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6;
1959 if(E1 < J1) E1 = J1;
1960 if(E2 < J2) E2 = J2;
1961 if(E3 < J3) E3 = J3;
1963 mySmoothCriterion->EstLength() = Length;
1964 mySmoothCriterion->SetEstimation(E1, E2, E3);
1966 Standard_Real WQuadratic, WQuality;
1968 if(!myWithMinMax && myTolerance != 0.)
1969 WQuality = myTolerance;
1970 else if (myTolerance == 0.)
1973 WQuality = Max(myTolerance, Eps2 * Length);
1975 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1976 WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality;
1977 if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic;
1979 if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.);
1982 mySmoothCriterion->SetWeight(WQuadratic, WQuality,
1983 myPercent[0], myPercent[1], myPercent[2]);
1986 Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity);
1987 Handle(FEmTool_Curve) TheCurve;
1988 Standard_Integer NbElem;
1989 Standard_Real CurvTol = Eps2 * Length / myNbPoints;
1991 // Decoupe de l'intervalle en fonction des contraintes
1992 if(myWithCutting == Standard_True && NbConstr != 0) {
1994 InitCutting(TheBase, CurvTol, TheCurve);
2000 TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol);
2001 TheCurve->Knots().SetValue(TheCurve->Knots().Lower(),
2002 myParameters->Value(myFirstPoint));
2003 TheCurve->Knots().SetValue(TheCurve->Knots().Upper(),
2004 myParameters->Value(myLastPoint));
2008 mySmoothCriterion->SetCurve(TheCurve);
2015 //=======================================================================
2016 //function : InitParameters
2017 //purpose : Calculation of initial estimation of the multicurve's length
2018 // and parameters for multipoints.
2019 //=======================================================================
2021 void AppDef_Variational::InitParameters(Standard_Real& Length)
2024 const Standard_Real Eps1 = Precision::Confusion() * .01;
2026 Standard_Real aux, dist;
2027 Standard_Integer i, i0, i1 = 0, ipoint;
2031 myParameters->SetValue(myFirstPoint, Length);
2033 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) {
2034 i0 = i1; i1 += myDimension;
2036 for(i = 1; i <= myDimension; i++) {
2037 aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i);
2040 Length += Sqrt(dist);
2041 myParameters->SetValue(ipoint, Length);
2046 Standard_ConstructionError::Raise("AppDef_Variational::InitParameters");
2049 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++)
2050 myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length);
2052 myParameters->SetValue(myLastPoint, 1.);
2054 // Avec peu de point il y a sans doute sous estimation ...
2055 if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1));
2059 //=======================================================================
2060 //function : InitCriterionEstimations
2061 //purpose : Calculation of initial estimation of the linear criteria
2063 //=======================================================================
2065 void AppDef_Variational::InitCriterionEstimations(const Standard_Real Length,
2068 Standard_Real& E3) const
2070 E1 = Length * Length;
2072 const Standard_Real Eps1 = Precision::Confusion() * .01;
2074 math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension),
2075 VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension);
2077 // ========== Treatment of first point =================
2079 Standard_Integer ipnt = myFirstPoint;
2081 EstTangent(ipnt, VTang1);
2083 EstTangent(ipnt, VTang2);
2085 EstTangent(ipnt, VTang3);
2087 ipnt = myFirstPoint;
2088 EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1);
2090 EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2);
2092 // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 Begin
2093 // Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt));
2094 Standard_Integer anInd = ipnt;
2095 Standard_Real Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt));
2096 // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 End
2098 if(Delta <= Eps1) Delta = 1.;
2100 E2 = VScnd1.Norm2() * Delta;
2102 E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
2103 // ========== Treatment of internal points =================
2105 Standard_Integer CurrPoint = 2;
2107 for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) {
2109 Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1));
2111 if(CurrPoint == 1) {
2112 if(ipnt + 1 != myLastPoint) {
2113 EstTangent(ipnt + 2, VTang3);
2114 EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2);
2117 EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2);
2119 E2 += VScnd1.Norm2() * Delta;
2120 E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.;
2123 else if(CurrPoint == 2) {
2124 if(ipnt + 1 != myLastPoint) {
2125 EstTangent(ipnt + 2, VTang1);
2126 EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3);
2129 EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3);
2131 E2 += VScnd2.Norm2() * Delta;
2132 E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
2136 if(ipnt + 1 != myLastPoint) {
2137 EstTangent(ipnt + 2, VTang2);
2138 EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1);
2141 EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1);
2143 E2 += VScnd3.Norm2() * Delta;
2144 E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.;
2148 CurrPoint++; if(CurrPoint == 4) CurrPoint = 1;
2151 // ========== Treatment of last point =================
2153 Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1));
2154 if(Delta <= Eps1) Delta = 1.;
2158 if(CurrPoint == 1) {
2160 E2 += VScnd1.Norm2() * Delta;
2161 aux = VScnd1.Subtracted(VScnd3).Norm2();
2162 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
2165 else if(CurrPoint == 2) {
2167 E2 += VScnd2.Norm2() * Delta;
2168 aux = VScnd2.Subtracted(VScnd1).Norm2();
2169 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
2174 E2 += VScnd3.Norm2() * Delta;
2175 aux = VScnd3.Subtracted(VScnd2).Norm2();
2176 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
2180 aux = Length * Length;
2182 E2 *= aux; E3 *= aux;
2188 //=======================================================================
2189 //function : EstTangent
2190 //purpose : Calculation of estimation of the Tangent
2191 // (see fortran routine MMLIPRI)
2192 //=======================================================================
2195 void AppDef_Variational::EstTangent(const Standard_Integer ipnt,
2196 math_Vector& VTang) const
2199 Standard_Integer i ;
2200 const Standard_Real Eps1 = Precision::Confusion() * .01;
2201 const Standard_Real EpsNorm = 1.e-9;
2203 Standard_Real Wpnt = 1.;
2206 if(ipnt == myFirstPoint) {
2207 // Estimation at first point
2212 Standard_Integer adr1 = 1, adr2 = adr1 + myDimension,
2213 adr3 = adr2 + myDimension;
2215 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
2216 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
2217 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
2219 // Parabolic interpolation
2220 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
2221 // first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d))
2222 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d
2223 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
2224 Standard_Real V2 = 0.;
2225 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
2227 Standard_Real d = V1 / (V1 + V2), d1;
2228 d1 = 1. / (d * (1 - d)); d *= d;
2229 VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1;
2232 // Simple 2-point estimation
2234 VTang = Pnt2 - Pnt1;
2238 else if (ipnt == myLastPoint) {
2239 // Estimation at last point
2244 Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1,
2245 adr2 = adr1 + myDimension,
2246 adr3 = adr2 + myDimension;
2248 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
2249 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
2250 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
2252 // Parabolic interpolation
2253 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
2254 // first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d))
2255 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2)
2256 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
2257 Standard_Real V2 = 0.;
2258 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
2260 Standard_Real d = V1 / (V1 + V2), d1;
2261 d1 = 1. / (d * (1 - d)); d *= d - 2;
2262 VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1;
2265 // Simple 2-point estimation
2267 VTang = Pnt3 - Pnt2;
2273 Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1,
2274 adr2 = adr1 + 2 * myDimension;
2276 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
2277 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
2279 VTang = Pnt2 - Pnt1;
2283 Standard_Real Vnorm = VTang.Norm();
2285 if(Vnorm <= EpsNorm)
2290 // Estimation with constraints
2292 Standard_Real Wcnt = 0.;
2293 Standard_Integer IdCnt = 1;
2295 // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
2297 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2298 math_Vector VCnt(1, myDimension, 0.);
2302 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
2303 IdCnt <= NbConstr) IdCnt++;
2304 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
2305 (myTypConstraints->Value(2 * IdCnt) >= 1)) {
2307 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0;
2308 for( i = 1; i <= myNbP3d; i++) {
2309 for(Standard_Integer j = 1; j <= 3; j++)
2310 VCnt(++k) = myTabConstraints->Value(++i0);
2313 for(i = 1; i <= myNbP2d; i++) {
2314 for(Standard_Integer j = 1; j <= 2; j++)
2315 VCnt(++k) = myTabConstraints->Value(++i0);
2321 // Averaging of estimation
2323 Standard_Real Denom = Wpnt + Wcnt;
2324 if(Denom == 0.) Denom = 1.;
2325 else Denom = 1. / Denom;
2327 VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom;
2329 Vnorm = VTang.Norm();
2331 if(Vnorm <= EpsNorm)
2340 //=======================================================================
2341 //function : EstSecnd
2342 //purpose : Calculation of estimation of the second derivative
2343 // (see fortran routine MLIMSCN)
2344 //=======================================================================
2346 void AppDef_Variational::EstSecnd(const Standard_Integer ipnt,
2347 const math_Vector& VTang1,
2348 const math_Vector& VTang2,
2349 const Standard_Real Length,
2350 math_Vector& VScnd) const
2352 Standard_Integer i ;
2354 const Standard_Real Eps = 1.e-9;
2356 Standard_Real Wpnt = 1.;
2360 if(ipnt == myFirstPoint)
2361 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt);
2362 else if(ipnt == myLastPoint)
2363 aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1);
2365 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1);
2372 VScnd = (VTang2 - VTang1) * aux;
2374 // Estimation with constraints
2376 Standard_Real Wcnt = 0.;
2377 Standard_Integer IdCnt = 1;
2379 // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
2381 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2382 math_Vector VCnt(1, myDimension, 0.);
2386 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
2387 IdCnt <= NbConstr) IdCnt++;
2389 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
2390 (myTypConstraints->Value(2 * IdCnt) >= 2)) {
2392 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0;
2393 for( i = 1; i <= myNbP3d; i++) {
2394 for(Standard_Integer j = 1; j <= 3; j++)
2395 VCnt(++k) = myTabConstraints->Value(++i0);
2399 for(i = 1; i <= myNbP2d; i++) {
2400 for(Standard_Integer j = 1; j <= 2; j++)
2401 VCnt(++k) = myTabConstraints->Value(++i0);
2407 // Averaging of estimation
2409 Standard_Real Denom = Wpnt + Wcnt;
2410 if(Denom == 0.) Denom = 1.;
2411 else Denom = 1. / Denom;
2413 VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom;
2419 //=======================================================================
2420 //function : InitCutting
2421 //purpose : Realisation of curve's cutting a priory accordingly to
2422 // constraints (see fortran routine MLICUT)
2423 //=======================================================================
2425 void AppDef_Variational::InitCutting(const Handle(PLib_Base)& aBase,
2426 const Standard_Real CurvTol,
2427 Handle(FEmTool_Curve)& aCurve) const
2430 // Definition of number of elements
2431 Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem;
2432 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2434 for(i = 1; i <= NbConstr; i++) {
2435 kk = Abs(myTypConstraints->Value(2 * i)) + 1;
2436 ORCMx = Max(ORCMx, kk);
2440 if(ORCMx > myMaxDegree - myNivCont)
2441 Standard_ConstructionError::Raise("AppDef_Variational::InitCutting");
2443 Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4,
2446 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
2448 while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) {
2451 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
2456 if(NbElem > myMaxSegment)
2457 Standard_ConstructionError::Raise("AppDef_Variational::InitCutting");
2460 aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol);
2462 Standard_Integer NCnt = (NCont - 1) / NbElem + 1;
2463 Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont);
2465 TColStd_Array1OfReal& Knot = aCurve->Knots();
2467 Standard_Integer IDeb = 0, IFin = NbConstr + 1,
2469 IndEl = Knot.Lower(), IUpper = Knot.Upper(),
2473 Knot(IndEl) = myParameters->Value(myFirstPoint);
2474 Knot(IUpper) = myParameters->Value(myLastPoint);
2476 while(NbElem - NbEl > 1) {
2479 if(NPlus == 0) NCnt--;
2481 while(NDeb < NCnt && IDeb < IFin) {
2483 NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1;
2489 myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1))
2491 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
2493 Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) +
2494 myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2;
2499 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
2503 if(NPlus == 0) NCnt--;
2505 if(NbElem - NbEl == 1) break;
2509 while(NFin < NCnt && IDeb < IFin) {
2511 NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1;
2516 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
2517 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
2521 if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1))
2522 Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1));
2524 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
2525 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
2530 //=======================================================================
2531 //function : Adjusting
2532 //purpose : Smoothing's adjusting like STRIM routine "MAJLIS"
2533 //=======================================================================
2534 void AppDef_Variational::Adjusting(
2535 Handle(AppDef_SmoothCriterion)& J,
2536 Standard_Real& WQuadratic,
2537 Standard_Real& WQuality,
2538 Handle(FEmTool_Curve)& TheCurve,
2539 TColStd_Array1OfReal& Ecarts)
2542 // cout << "=========== Adjusting =============" << endl;
2544 /* Initialized data */
2546 const Standard_Integer mxiter = 2;
2547 const Standard_Real eps1 = 1e-6;
2548 Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
2549 Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2550 Standard_Real CurvTol = eps1 * J->EstLength() / NbrPnt;
2553 /* Local variables */
2554 Standard_Integer iter, ipnt;
2555 Standard_Real ecart, emold, erold, tpara;
2556 Standard_Real vocri[4], j1cibl, vtest, vseuil;
2557 Standard_Integer i, numint, flag;
2558 TColStd_Array1OfReal tbpoid(myFirstPoint, myLastPoint);
2559 Standard_Boolean loptim, lrejet;
2560 Handle(AppDef_SmoothCriterion) JNew;
2561 Handle(FEmTool_Curve) CNew;
2562 Standard_Real E1, E2, E3;
2565 /* (0.b) Initialisations */
2567 loptim = Standard_True;
2572 /* ============ boucle sur le moteur de lissage ============== */
2574 vtest = WQuality * .9;
2575 j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
2581 /* (1) Sauvegarde de l'etat precedents */
2583 vocri[0] = myCriterium[0];
2584 vocri[1] = myCriterium[1];
2585 vocri[2] = myCriterium[2];
2586 vocri[3] = myCriterium[3];
2588 emold = myAverageError;
2590 /* (2) Augmentation du poids des moindre carre */
2592 if (j1cibl > vtest) {
2593 WQuadratic = j1cibl / vtest * WQuadratic;
2596 /* (3) Augmentation du poid associe aux points a problemes */
2598 vseuil = WQuality * .88;
2600 for (ipnt = myFirstPoint; ipnt <= myLastPoint; ++ipnt) {
2601 if (Ecarts(ipnt) > vtest) {
2602 ecart = (Ecarts(ipnt) - vseuil) / WQuality;
2603 tbpoid(ipnt) = (ecart * 3 + 1.) * tbpoid(ipnt);
2607 /* (4) Decoupe force */
2609 if (TheCurve->NbElements() < myMaxSegment && myWithCutting) {
2611 numint = NearIndex(myParameters->Value(myMaxErrorIndex), TheCurve->Knots(), 0, flag);
2613 tpara = (TheCurve->Knots()(numint) + TheCurve->Knots()(numint + 1) +
2614 myParameters->Value(myMaxErrorIndex) * 2) / 4;
2616 CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements() + 1,
2617 TheCurve->Base(), CurvTol);
2619 for(i = 1; i <= numint; i++) CNew->Knots()(i) = TheCurve->Knots()(i);
2620 for(i = numint + 1; i <= TheCurve->Knots().Length(); i++)
2621 CNew->Knots()(i + 1) = TheCurve->Knots()(i);
2623 CNew->Knots()(numint + 1) = tpara;
2627 CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements(), TheCurve->Base(), CurvTol);
2629 CNew->Knots() = TheCurve->Knots();
2633 JNew = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint);
2635 JNew->EstLength() = J->EstLength();
2637 J->GetEstimation(E1, E2, E3);
2639 JNew->SetEstimation(E1, E2, E3);
2641 JNew->SetParameters(myParameters);
2643 JNew->SetWeight(WQuadratic, WQuality, myPercent[0], myPercent[1], myPercent[2]);
2645 JNew->SetWeight(tbpoid);
2647 JNew->SetCurve(CNew);
2651 TheMotor(JNew, WQuadratic, WQuality, CNew, Ecarts);
2653 /* (6) Tests de rejet */
2655 j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
2656 vseuil = Sqrt(vocri[1]) + (erold - myMaxError) * 4;
2658 lrejet = ((myMaxError > WQuality) && (myMaxError > erold * 1.01)) || (Sqrt(myCriterium[1]) > vseuil * 1.05);
2661 myCriterium[0] = vocri[0];
2662 myCriterium[1] = vocri[1];
2663 myCriterium[2] = vocri[2];
2664 myCriterium[3] = vocri[3];
2666 myAverageError = emold;
2668 loptim = Standard_False;
2673 J->SetCurve(TheCurve);
2676 /* (7) Test de convergence */
2678 if (((iter >= mxiter) && (myMaxSegment == CNew->NbElements())) || myMaxError < WQuality) {
2679 loptim = Standard_False;
2684 static Standard_Boolean NotParallel(gp_Vec& T, gp_Vec& V)
2688 if (V.CrossMagnitude(T) > 1.e-12)
2689 return Standard_True;
2691 if (V.CrossMagnitude(T) > 1.e-12)
2692 return Standard_True;
2694 if (V.CrossMagnitude(T) > 1.e-12)
2695 return Standard_True;
2696 return Standard_False;
2699 void AppDef_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve,
2700 const TColStd_Array1OfReal& Parameters,
2701 const Standard_Real CBLONG,
2702 FEmTool_Assembly& A ) const
2705 Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
2706 NbElm = Curve->NbElements(),
2707 NbDim = Curve->Dimension();
2709 TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg);
2710 math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg),
2711 V1((Standard_Real*)&G1(0), 0, MxDeg),
2712 V2((Standard_Real*)&G2(0), 0, MxDeg);
2714 Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1,
2717 p0 = Parameters.Lower() - myFirstPoint,
2718 curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim,
2719 jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d;
2720 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2722 // Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints;
2723 // Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints;
2724 Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints;
2725 Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints;
2726 NBeg2d = Ng3d * myNbP3d;
2727 // NgPC1 = NbConstr + myNbCurvPoints;
2728 NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints;
2730 NTang3d = 3 * NgPC1;
2731 NTang2d = 2 * NgPC1;
2733 TColStd_Array1OfReal& Intervals = Curve->Knots();
2735 Standard_Real t, R1, R2;
2737 Handle(PLib_Base) myBase = Curve->Base();
2738 Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase);
2739 Standard_Integer Order = myHermitJacobi->NivConstr() + 1;
2741 Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1;
2743 A.NullifyConstraint();
2747 for(i = 1; i <= NbConstr; i++) {
2749 ipnt += 2; ityp += 2;
2751 Point = myTypConstraints->Value(ipnt);
2752 TypOfConstr = myTypConstraints->Value(ityp);
2754 t = Parameters(p0 + Point);
2756 for(el = curel; el <= NbElm; ) {
2757 if( t <= Intervals(++el) ) {
2764 UFirst = Intervals(curel); ULast = Intervals(curel + 1);
2765 coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.;
2767 t = (t - c0) / coeff;
2769 if(TypOfConstr == 0) {
2771 for(k = 1; k < Order; k++) {
2772 mfact = Pow(coeff, k);
2774 G0(k + Order) *= mfact;
2777 else if(TypOfConstr == 1) {
2778 myBase->D1(t, G0, G1);
2779 for(k = 1; k < Order; k++) {
2780 mfact = Pow(coeff, k);
2782 G0(k + Order) *= mfact;
2784 G1(k + Order) *= mfact;
2787 for(k = 0; k <= MxDeg; k++) {
2792 myBase->D2(t, G0, G1, G2);
2793 for(k = 1; k < Order; k++) {
2794 mfact = Pow(coeff, k);
2796 G0(k + Order) *= mfact;
2798 G1(k + Order) *= mfact;
2800 G2(k + Order) *= mfact;
2803 mfact1 = mfact / coeff;
2804 for(k = 0; k <= MxDeg; k++) {
2812 j = NbDim * (Point - myFirstPoint);
2813 Standard_Integer n0 = NPass;
2815 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2816 IndexOfConstraint = n0;
2817 for(k = 1; k <= 3; k++) {
2819 A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
2820 IndexOfConstraint += NgPC1;
2826 n0 = NPass + NBeg2d;
2827 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2828 IndexOfConstraint = n0;
2829 for(k = 1; k <= 2; k++) {
2831 A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
2832 IndexOfConstraint += NgPC1;
2838 /* if(TypOfConstr == 1) {
2840 IndexOfConstraint = NTang3d + 1;
2841 jt = Ntheta * (i - 1);
2842 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2843 for(k = 1; k <= 3; k++) {
2844 A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
2845 A.AddConstraint(IndexOfConstraint + 1, curel, k, myTtheta->Value(jt + 3 + k) * V1, 0.);
2847 IndexOfConstraint += Ng3d;
2851 IndexOfConstraint = NBeg2d + NTang2d + 1;
2852 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2853 for(k = 1; k <= 2; k++) {
2854 A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
2856 IndexOfConstraint += Ng2d;
2863 if(TypOfConstr == 1) {
2867 j = 2 * NbDim * (i - 1);
2869 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2870 IndexOfConstraint = n0;
2871 for(k = 1; k <= 3; k++) {
2873 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2874 IndexOfConstraint += NgPC1;
2880 n0 = NPass + NBeg2d;
2881 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2882 IndexOfConstraint = n0;
2883 for(k = 1; k <= 2; k++) {
2885 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2886 IndexOfConstraint += NgPC1;
2892 if(TypOfConstr == 2) {
2896 j = 2 * NbDim * (i - 1);
2898 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2899 IndexOfConstraint = n0;
2900 for(k = 1; k <= 3; k++) {
2902 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2903 IndexOfConstraint += NgPC1;
2909 n0 = NPass + NBeg2d;
2910 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2911 IndexOfConstraint = n0;
2912 for(k = 1; k <= 2; k++) {
2914 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2915 IndexOfConstraint += NgPC1;
2921 j = 2 * NbDim * (i - 1) + 3;
2922 jt = Ntheta * (i - 1);
2923 IndexOfConstraint = NTang3d + 1;
2925 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2927 for(k = 1; k <= 3; k++) {
2928 R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
2929 R2 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + 3 + k);
2931 R1 *= CBLONG * CBLONG;
2932 R2 *= CBLONG * CBLONG;
2933 for(k = 1; k <= 3; k++) {
2935 if(k > 1) R1 = R2 = 0.;
2936 A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
2937 A.AddConstraint(IndexOfConstraint + 1, curel, curdim, myTfthet->Value(jt + 3 + k) * V2, R2);
2939 IndexOfConstraint += Ng3d;
2945 IndexOfConstraint = NBeg2d + NTang2d + 1;
2946 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2948 for(k = 1; k <= 2; k++) {
2949 R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
2951 R1 *= CBLONG * CBLONG;
2952 for(k = 1; k <= 2; k++) {
2955 A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
2957 IndexOfConstraint += Ng2d;
2970 Standard_Boolean AppDef_Variational::InitTthetaF(const Standard_Integer ndimen,
2971 const AppParCurves_Constraint typcon,
2972 const Standard_Integer begin,
2973 const Standard_Integer jndex)
2975 if ((ndimen < 2)||(ndimen >3))
2976 return Standard_False;
2978 gp_Vec theta1, theta2;
2980 Standard_Real XX, XY, YY, XZ, YZ, ZZ;
2982 if ((typcon == AppParCurves_TangencyPoint)||(typcon == AppParCurves_CurvaturePoint))
2984 T.SetX(myTabConstraints->Value(jndex));
2985 T.SetY(myTabConstraints->Value(jndex + 1));
2987 T.SetZ(myTabConstraints->Value(jndex + 2));
2997 if (!NotParallel(T, V))
2998 return Standard_False;
3001 myTtheta->SetValue(begin, theta1.X());
3002 myTtheta->SetValue(begin + 1, theta1.Y());
3005 theta2 = T ^ theta1;
3007 myTtheta->SetValue(begin + 2, theta1.Z());
3008 myTtheta->SetValue(begin + 3, theta2.X());
3009 myTtheta->SetValue(begin + 4, theta2.Y());
3010 myTtheta->SetValue(begin + 5, theta2.Z());
3013 // Calculation of myTfthet
3014 if (typcon == AppParCurves_CurvaturePoint)
3021 F.SetX(YY * theta1.X() - XY * theta1.Y());
3022 F.SetY(XX * theta1.Y() - XY * theta1.X());
3023 myTfthet->SetValue(begin, F.X());
3024 myTfthet->SetValue(begin + 1, F.Y());
3032 F.SetX((ZZ + YY) * theta1.X() - XY * theta1.Y() - XZ * theta1.Z());
3033 F.SetY((XX + ZZ) * theta1.Y() - XY * theta1.X() - YZ * theta1.Z());
3034 F.SetZ((XX + YY) * theta1.Z() - XZ * theta1.X() - YZ * theta1.Y());
3035 myTfthet->SetValue(begin, F.X());
3036 myTfthet->SetValue(begin + 1, F.Y());
3037 myTfthet->SetValue(begin + 2, F.Z());
3038 F.SetX((ZZ + YY) * theta2.X() - XY * theta2.Y() - XZ * theta2.Z());
3039 F.SetY((XX + ZZ) * theta2.Y() - XY * theta2.X() - YZ * theta2.Z());
3040 F.SetZ((XX + YY) * theta2.Z() - XZ * theta2.X() - YZ * theta2.Y());
3041 myTfthet->SetValue(begin + 3, F.X());
3042 myTfthet->SetValue(begin + 4, F.Y());
3043 myTfthet->SetValue(begin + 5, F.Z());
3047 return Standard_True;