1 // Copyright (c) 1996-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
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.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14
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
19
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>
34
35 #define No_Standard_RangeError
36 #define No_Standard_OutOfRange
37 #define No_Standard_DimensionError
38 #define No_Standard_ConstructionError
39
40 #include <Standard_Stream.hxx>
41
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>
50 #include <gp_Pnt.hxx>
51 #include <gp_Pnt2d.hxx>
52 #include <gp_Vec.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>
68
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>
81
82 // Add this line:
83 #include <algorithm>
84
85 #if defined(_MSC_VER)
86 # include <stdio.h>
87 # include <stdarg.h>
88 #endif /* _MSC_VER */
89
90 //
91 //=======================================================================
92 //function : AppDef_Variational
93 //purpose : Initialization of the fields.
94 //=======================================================================
95 //
AppDef_Variational(const AppDef_MultiLine & SSP,const Standard_Integer FirstPoint,const Standard_Integer LastPoint,const Handle (AppParCurves_HArray1OfConstraintCouple)& TheConstraints,const Standard_Integer MaxDegree,const Standard_Integer MaxSegment,const GeomAbs_Shape Continuity,const Standard_Boolean WithMinMax,const Standard_Boolean WithCutting,const Standard_Real Tolerance,const Standard_Integer NbIterations)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):
107 mySSP(SSP),
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)
118 {
119 // Verifications:
120 if (myMaxDegree < 1) throw Standard_DomainError();
121 myMaxDegree = Min (30, myMaxDegree);
122 //
123 if (myMaxSegment < 1) throw Standard_DomainError();
124 //
125 if (myWithMinMax != 0 && myWithMinMax !=1 ) throw Standard_DomainError();
126 if (myWithCutting != 0 && myWithCutting !=1 ) throw Standard_DomainError();
127 //
128 myIsOverConstr = Standard_False;
129 myIsCreated = Standard_False;
130 myIsDone = Standard_False;
131 switch (myContinuity) {
132 case GeomAbs_C0:
133 myNivCont=0;
134 break ;
135 case GeomAbs_C1:
136 myNivCont=1;
137 break ;
138 case GeomAbs_C2:
139 myNivCont=2;
140 break ;
141 default:
142 throw Standard_ConstructionError();
143 }
144 //
145 myNbP2d = AppDef_MyLineTool::NbP2d(SSP);
146 myNbP3d = AppDef_MyLineTool::NbP3d(SSP);
147 myDimension = 2 * myNbP2d + 3* myNbP3d ;
148 //
149 myPercent[0]=0.4;
150 myPercent[1]=0.2;
151 myPercent[2]=0.4;
152 myKnots= new TColStd_HArray1OfReal(1,2);
153 myKnots->SetValue(1,0.);
154 myKnots->SetValue(2,1.);
155
156 // Declaration
157 //
158 mySmoothCriterion = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint);
159 myParameters = new TColStd_HArray1OfReal(myFirstPoint, myLastPoint);
160 myNbPoints=myLastPoint-myFirstPoint+1;
161 if (myNbPoints <= 0) throw Standard_ConstructionError();
162 //
163 myTabPoints= new TColStd_HArray1OfReal(1,myDimension*myNbPoints);
164 //
165 // Table of Points initialization
166 //
167 Standard_Integer ipoint,jp2d,jp3d,index;
168 TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
169 TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
170 gp_Pnt2d P2d;
171 gp_Pnt P3d;
172 index=1;
173
174 for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
175 {
176
177 if(myNbP2d !=0 && myNbP3d ==0 ) {
178 AppDef_MyLineTool::Value(mySSP,ipoint,TabP2d);
179
180 for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++)
181 { P2d = TabP2d.Value(jp2d);
182
183 myTabPoints->SetValue(index++,P2d.X());
184 myTabPoints->SetValue(index++,P2d.Y());
185 }
186 }
187 if(myNbP3d !=0 && myNbP2d == 0 ) {
188 AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d);
189
190 for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
191
192 { P3d=TabP3d.Value(jp3d);
193
194 myTabPoints->SetValue(index++,P3d.X());
195 myTabPoints->SetValue(index++,P3d.Y());
196 myTabPoints->SetValue(index++,P3d.Z());
197
198 }
199 }
200 if(myNbP3d !=0 && myNbP2d != 0 ) {
201 AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d,TabP2d);
202
203 for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
204
205 { P3d=TabP3d.Value(jp3d);
206
207 myTabPoints->SetValue(index++,P3d.X());
208 myTabPoints->SetValue(index++,P3d.Y());
209 myTabPoints->SetValue(index++,P3d.Z());
210
211 }
212 for ( jp2d = 1 ; jp2d <= myNbP2d ; jp2d++)
213
214 { P2d=TabP2d.Value(jp2d);
215
216 myTabPoints->SetValue(index++,P2d.X());
217 myTabPoints->SetValue(index++,P2d.Y());
218 }
219 }
220 }
221 Init();
222 }
223 //
224 //=======================================================================
225 //function : Init
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 //=======================================================================
230 //
Init()231 void AppDef_Variational::Init()
232 {
233
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));
240
241 gp_Vec Vt3d, Vc3d;
242 gp_Vec2d Vt2d, Vc2d;
243
244 myNbConstraints=myConstraints->Length();
245 if (myNbConstraints < 0) throw Standard_ConstructionError();
246
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));
251
252
253 //
254 // Table of types initialization
255 Standard_Integer iconstr;
256 index=1;
257 jndex=1;
258 CurMultyPoint = 1;
259 myNbPassPoints=0;
260 myNbTangPoints=0;
261 myNbCurvPoints=0;
262 AppParCurves_Constraint valcontr;
263
264 for ( iconstr = myConstraints->Lower() ; iconstr <= myConstraints->Upper() ; iconstr++)
265 {
266 ipoint=(myConstraints->Value(iconstr)).Index();
267
268 valcontr=(myConstraints->Value(iconstr)).Constraint();
269 switch (valcontr) {
270 case AppParCurves_NoConstraint:
271 CurMultyPoint -= myNbP3d * 6 + myNbP2d * 2;
272 break ;
273 case AppParCurves_PassPoint:
274 myTypConstraints->SetValue(index++,ipoint);
275 myTypConstraints->SetValue(index++,0);
276 myNbPassPoints++;
277 if(myNbP2d !=0 ) jndex=jndex+4*myNbP2d;
278 if(myNbP3d !=0 ) jndex=jndex+6*myNbP3d;
279 break ;
280 case AppParCurves_TangencyPoint:
281 myTypConstraints->SetValue(index++,ipoint);
282 myTypConstraints->SetValue(index++,1);
283 myNbTangPoints++;
284 if(myNbP2d !=0 && myNbP3d == 0 )
285 {
286 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False)
287 throw Standard_ConstructionError();
288 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
289 {
290 Vt2d=TabV2d.Value(jp2d);
291 Vt2d.Normalize();
292 myTabConstraints->SetValue(jndex++,Vt2d.X());
293 myTabConstraints->SetValue(jndex++,Vt2d.Y());
294 jndex=jndex+2;
295 InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4);
296 }
297 }
298 if(myNbP3d !=0 && myNbP2d == 0)
299 {
300 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False)
301 throw Standard_ConstructionError();
302 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
303 {
304 Vt3d=TabV3d.Value(jp3d);
305 Vt3d.Normalize();
306 myTabConstraints->SetValue(jndex++,Vt3d.X());
307
308 myTabConstraints->SetValue(jndex++,Vt3d.Y());
309
310 myTabConstraints->SetValue(jndex++,Vt3d.Z());
311 jndex=jndex+3;
312 InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
313 }
314 }
315 if(myNbP3d !=0 && myNbP2d != 0)
316 {
317 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False)
318 throw Standard_ConstructionError();
319 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
320 {
321 Vt3d=TabV3d.Value(jp3d);
322 Vt3d.Normalize();
323 myTabConstraints->SetValue(jndex++,Vt3d.X());
324 myTabConstraints->SetValue(jndex++,Vt3d.Y());
325 myTabConstraints->SetValue(jndex++,Vt3d.Z());
326 jndex=jndex+3;
327 InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6);
328 }
329
330 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
331 {
332 Vt2d=TabV2d.Value(jp2d);
333 Vt2d.Normalize();
334 myTabConstraints->SetValue(jndex++,Vt2d.X());
335 myTabConstraints->SetValue(jndex++,Vt2d.Y());
336 jndex=jndex+2;
337 InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4);
338 }
339 }
340
341
342 break ;
343
344 case AppParCurves_CurvaturePoint:
345 myTypConstraints->SetValue(index++,ipoint);
346 myTypConstraints->SetValue(index++,2);
347 myNbCurvPoints++;
348 if(myNbP2d !=0 && myNbP3d == 0)
349 {
350 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False )
351 throw Standard_ConstructionError();
352 if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV2dcurv) == Standard_False)
353 throw Standard_ConstructionError();
354 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
355 {
356 Vt2d=TabV2d.Value(jp2d);
357 Vt2d.Normalize();
358 Vc2d=TabV2dcurv.Value(jp2d);
359 if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
360 throw Standard_ConstructionError();
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);
366 }
367 }
368
369 if(myNbP3d !=0 && myNbP2d == 0 )
370 {
371 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False )
372 throw Standard_ConstructionError();
373 if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv) == Standard_False)
374 throw Standard_ConstructionError();
375 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
376 {
377 Vt3d=TabV3d.Value(jp3d);
378 Vt3d.Normalize();
379 Vc3d=TabV3dcurv.Value(jp3d);
380 if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False)
381 throw Standard_ConstructionError();
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);
389 }
390 }
391 if(myNbP3d !=0 && myNbP2d != 0 )
392 {
393 if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False )
394 throw Standard_ConstructionError();
395 if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv,TabV2dcurv) == Standard_False)
396 throw Standard_ConstructionError();
397 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
398 {
399 Vt3d=TabV3d.Value(jp3d);
400 Vt3d.Normalize();
401 Vc3d=TabV3dcurv.Value(jp3d);
402 if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False)
403 throw Standard_ConstructionError();
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);
411 }
412 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
413 {
414 Vt2d=TabV2d.Value(jp2d);
415 Vt2d.Normalize();
416 Vc2d=TabV2dcurv.Value(jp2d);
417 if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular())
418 throw Standard_ConstructionError();
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);
424 }
425
426 }
427 break ;
428 default:
429 throw Standard_ConstructionError();
430 }
431 CurMultyPoint += myNbP3d * 6 + myNbP2d * 2;
432 }
433 // OverConstraint Detection
434 Standard_Integer MaxSeg;
435 if(myWithCutting == Standard_True) MaxSeg = myMaxSegment ;
436 else MaxSeg = 1;
437 if (((myMaxDegree-myNivCont)*MaxSeg-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
438 {
439 myIsOverConstr =Standard_True;
440 myIsCreated = Standard_False;
441 }
442 else {
443 InitSmoothCriterion();
444 myIsCreated = Standard_True;
445 }
446
447
448 }
449 //
450 //=======================================================================
451 //function : Approximate
452 //purpose : Makes the approximation with the current fields.
453 //=======================================================================
454 //
Approximate()455 void AppDef_Variational::Approximate()
456
457 {
458 if (myIsCreated == Standard_False ) throw StdFail_NotDone();
459
460
461 Standard_Real WQuadratic, WQuality;
462
463 TColStd_Array1OfReal Ecarts(myFirstPoint, myLastPoint);
464
465 mySmoothCriterion->GetWeight(WQuadratic, WQuality);
466
467 Handle(FEmTool_Curve) TheCurve;
468
469 mySmoothCriterion->GetCurve(TheCurve);
470
471 //---------------------------------------------------------------------
472
473 TheMotor(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
474
475
476 if(myWithMinMax && myTolerance < myMaxError)
477 Adjusting(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts);
478
479 //---------------------------------------------------------------------
480
481 Standard_Integer jp2d,jp3d,index,ipole,
482 NbElem = TheCurve->NbElements();
483
484 TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
485 TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
486 Standard_Real debfin[2] = {-1., 1};
487
488 gp_Pnt2d P2d;
489 gp_Pnt P3d;
490
491 index=0;
492
493 {
494 Handle(TColStd_HArray2OfReal) PolynomialIntervalsPtr =
495 new TColStd_HArray2OfReal(1,NbElem,1,2) ;
496
497 Handle(TColStd_HArray1OfInteger) NbCoeffPtr =
498 new TColStd_HArray1OfInteger(1, myMaxSegment);
499
500 Standard_Integer size = myMaxSegment * (myMaxDegree + 1) * myDimension ;
501 Handle(TColStd_HArray1OfReal) CoeffPtr = new TColStd_HArray1OfReal(1,size);
502
503 CoeffPtr->Init(0.);
504
505 Handle(TColStd_HArray1OfReal) IntervallesPtr =
506 new TColStd_HArray1OfReal(1, NbElem + 1);
507
508 IntervallesPtr->ChangeArray1() = TheCurve->Knots();
509
510 TheCurve->GetPolynom(CoeffPtr->ChangeArray1());
511
512 Standard_Integer ii;
513
514 for(ii = 1; ii <= NbElem; ii++)
515 NbCoeffPtr->SetValue(ii, TheCurve->Degree(ii)+1);
516
517
518 for (ii = PolynomialIntervalsPtr->LowerRow() ;
519 ii <= PolynomialIntervalsPtr->UpperRow() ;ii++)
520 {
521 PolynomialIntervalsPtr->SetValue(ii,1,debfin[0]) ;
522 PolynomialIntervalsPtr->SetValue(ii,2,debfin[1]) ;
523 }
524 /*
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++)
538 {
539 printf(" %d : [%f, %f] \n", ii, PolynomialIntervalsPtr->Value(ii,1),
540 PolynomialIntervalsPtr->Value(ii,2)) ;
541 }
542 printf("\n IntervallesPtr :\n");
543 for (ii = IntervallesPtr->Lower();
544 ii <= IntervallesPtr->Upper() - 1; ii++)
545 {
546 printf(" %d : [%f, %f] \n", ii, IntervallesPtr->Value(ii),
547 IntervallesPtr->Value(ii+1)) ;
548 }
549 */
550 Convert_CompPolynomialToPoles AConverter(NbElem,
551 myNivCont,
552 myDimension,
553 myMaxDegree,
554 NbCoeffPtr,
555 CoeffPtr,
556 PolynomialIntervalsPtr,
557 IntervallesPtr) ;
558 if (AConverter.IsDone())
559 {
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) ;
568
569 for (ipole=PolesPtr->LowerRow();ipole<=PolesPtr->UpperRow();ipole++)
570 {
571 index=PolesPtr->LowerCol();
572 /* if(myNbP2d !=0 )
573 {
574 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
575 {
576 P2d.SetX(PolesPtr->Value(ipole,index++));
577 P2d.SetY(PolesPtr->Value(ipole,index++));
578 TabP2d.SetValue(jp2d,P2d);
579 }
580 }*/
581 if(myNbP3d !=0 )
582 {
583 for (jp3d=1;jp3d<=myNbP3d;jp3d++)
584 {
585 // std::cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
586 P3d.SetX(PolesPtr->Value(ipole,index++));
587 // std::cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
588 P3d.SetY(PolesPtr->Value(ipole,index++));
589 // std::cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
590 P3d.SetZ(PolesPtr->Value(ipole,index++));
591 TabP3d.SetValue(jp3d,P3d);
592 }
593 }
594 if(myNbP2d !=0 )
595 {
596 for (jp2d=1;jp2d<=myNbP2d;jp2d++)
597 {
598 P2d.SetX(PolesPtr->Value(ipole,index++));
599 P2d.SetY(PolesPtr->Value(ipole,index++));
600 TabP2d.SetValue(jp2d,P2d);
601 }
602 }
603 if(myNbP2d !=0 && myNbP3d !=0)
604 {
605 AppParCurves_MultiPoint aMultiPoint(TabP3d,TabP2d);
606 TabMU.SetValue(ipole,aMultiPoint);
607 }
608 else if (myNbP2d !=0)
609 {
610 AppParCurves_MultiPoint aMultiPoint(TabP2d);
611 TabMU.SetValue(ipole,aMultiPoint);
612 }
613 else
614 {
615
616 AppParCurves_MultiPoint aMultiPoint(TabP3d);
617 TabMU.SetValue(ipole,aMultiPoint);
618 }
619
620 }
621 AppParCurves_MultiBSpCurve aCurve(TabMU,myKnots->Array1(),Mults->Array1());
622 myMBSpCurve=aCurve;
623 myIsDone = Standard_True;
624 }
625 }
626
627 }
628 //
629 //=======================================================================
630 //function : IsCreated
631 //purpose : returns True if the creation is done
632 //=======================================================================
633 //
IsCreated() const634 Standard_Boolean AppDef_Variational::IsCreated() const
635 {
636 return myIsCreated;
637 }
638 //
639 //=======================================================================
640 //function : IsDone
641 //purpose : returns True if the approximation is ok
642 //=======================================================================
643 //
IsDone() const644 Standard_Boolean AppDef_Variational::IsDone() const
645 {
646 return myIsDone;
647 }
648 //
649 //=======================================================================
650 //function : IsOverConstrained
651 //purpose : returns True if the problem is overconstrained
652 // in this case, approximation cannot be done.
653 //=======================================================================
654 //
IsOverConstrained() const655 Standard_Boolean AppDef_Variational::IsOverConstrained() const
656 {
657 return myIsOverConstr;
658 }
659 //
660 //=======================================================================
661 //function : Value
662 //purpose : returns all the BSpline curves approximating the
663 // MultiLine SSP after minimization of the parameter.
664
665 //=======================================================================
666 //
Value() const667 AppParCurves_MultiBSpCurve AppDef_Variational::Value() const
668 {
669 if (myIsDone == Standard_False) throw StdFail_NotDone();
670 return myMBSpCurve;
671
672 }
673 //
674 //=======================================================================
675 //function : MaxError
676 //purpose : returns the maximum of the distances between
677 // the points of the multiline and the approximation
678 // curves.
679 //=======================================================================
680 //
MaxError() const681 Standard_Real AppDef_Variational::MaxError() const
682 {
683 if (myIsDone == Standard_False) throw StdFail_NotDone();
684 return myMaxError;
685 }
686 //
687 //=======================================================================
688 //function : MaxErrorIndex
689 //purpose : returns the index of the MultiPoint of ErrorMax
690 //=======================================================================
691 //
MaxErrorIndex() const692 Standard_Integer AppDef_Variational::MaxErrorIndex() const
693 {
694 if (myIsDone == Standard_False) throw StdFail_NotDone();
695 return myMaxErrorIndex;
696 }
697 //
698 //=======================================================================
699 //function : QuadraticError
700 //purpose : returns the quadratic average of the distances between
701 // the points of the multiline and the approximation
702 // curves.
703 //=======================================================================
704 //
QuadraticError() const705 Standard_Real AppDef_Variational::QuadraticError() const
706 {
707 if (myIsDone == Standard_False) throw StdFail_NotDone();
708 return myCriterium[0];
709 }
710 //
711 //=======================================================================
712 //function : Distance
713 //purpose : returns the distances between the points of the
714 // multiline and the approximation curves.
715 //=======================================================================
716 //
Distance(math_Matrix & mat)717 void AppDef_Variational::Distance(math_Matrix& mat)
718
719 {
720 if (myIsDone == Standard_False) throw StdFail_NotDone();
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;
725
726 gp_Pnt2d P2d;
727 gp_Pnt P3d;
728
729
730 gp_Pnt Pt3d;
731 gp_Pnt2d Pt2d;
732
733 for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++)
734 {
735 index=1;
736 if(myNbP3d !=0 ) {
737 AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d);
738
739 for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++)
740
741 { P3d=TabP3d.Value(jp3d);
742 myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt3d);
743 mat(index++, j0 + ipoint)=P3d.Distance(Pt3d);
744
745 }
746 }
747 if(myNbP2d !=0 ) {
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++)
751
752 { P2d = TabP2d.Value(jp2d);
753 myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt2d);
754 mat(index++, j0 + ipoint)=P2d.Distance(Pt2d);
755 }
756 }
757 }
758
759 }
760 //
761 //=======================================================================
762 //function : AverageError
763 //purpose : returns the average error between
764 // the MultiLine and the approximation.
765 //=======================================================================
766 //
AverageError() const767 Standard_Real AppDef_Variational::AverageError() const
768 {
769 if (myIsDone == Standard_False) throw StdFail_NotDone();
770 return myAverageError;
771 }
772 //
773 //=======================================================================
774 //function : Parameters
775 //purpose : returns the parameters uses to the approximations
776 //=======================================================================
777 //
Handle(TColStd_HArray1OfReal)778 const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Parameters() const
779 {
780 if (myIsDone == Standard_False) throw StdFail_NotDone();
781 return myParameters;
782 }
783 //
784 //=======================================================================
785 //function : Knots
786 //purpose : returns the knots uses to the approximations
787 //=======================================================================
788 //
Handle(TColStd_HArray1OfReal)789 const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Knots() const
790 {
791 if (myIsDone == Standard_False) throw StdFail_NotDone();
792 return myKnots;
793 }
794 //
795 //=======================================================================
796 //function : Criterium
797 //purpose : returns the values of the quality criterium.
798 //=======================================================================
799 //
Criterium(Standard_Real & VFirstOrder,Standard_Real & VSecondOrder,Standard_Real & VThirdOrder) const800 void AppDef_Variational::Criterium(Standard_Real& VFirstOrder, Standard_Real& VSecondOrder, Standard_Real& VThirdOrder) const
801 {
802 if (myIsDone == Standard_False) throw StdFail_NotDone();
803 VFirstOrder=myCriterium[1] ;
804 VSecondOrder=myCriterium[2];
805 VThirdOrder=myCriterium[3];
806 }
807 //
808 //=======================================================================
809 //function : CriteriumWeight
810 //purpose : returns the Weights (as percent) associed to the criterium used in
811 // the optimization.
812 //=======================================================================
813 //
CriteriumWeight(Standard_Real & Percent1,Standard_Real & Percent2,Standard_Real & Percent3) const814 void AppDef_Variational::CriteriumWeight(Standard_Real& Percent1, Standard_Real& Percent2, Standard_Real& Percent3) const
815 {
816 Percent1 = myPercent[0];
817 Percent2 = myPercent[1];
818 Percent3 = myPercent[2] ;
819 }
820 //
821 //=======================================================================
822 //function : MaxDegree
823 //purpose : returns the Maximum Degree used in the approximation
824 //=======================================================================
825 //
MaxDegree() const826 Standard_Integer AppDef_Variational::MaxDegree() const
827 {
828 return myMaxDegree;
829 }
830 //
831 //=======================================================================
832 //function : MaxSegment
833 //purpose : returns the Maximum of segment used in the approximation
834 //=======================================================================
835 //
MaxSegment() const836 Standard_Integer AppDef_Variational::MaxSegment() const
837 {
838 return myMaxSegment;
839 }
840 //
841 //=======================================================================
842 //function : Continuity
843 //purpose : returns the Continuity used in the approximation
844 //=======================================================================
845 //
Continuity() const846 GeomAbs_Shape AppDef_Variational::Continuity() const
847 {
848 return myContinuity;
849 }
850 //
851 //=======================================================================
852 //function : WithMinMax
853 //purpose : returns if the approximation search to minimize the
854 // maximum Error or not.
855 //=======================================================================
856 //
WithMinMax() const857 Standard_Boolean AppDef_Variational::WithMinMax() const
858 {
859 return myWithMinMax;
860 }
861 //
862 //=======================================================================
863 //function : WithCutting
864 //purpose : returns if the approximation can insert new Knots or not.
865 //=======================================================================
866 //
WithCutting() const867 Standard_Boolean AppDef_Variational::WithCutting() const
868 {
869 return myWithCutting;
870 }
871 //
872 //=======================================================================
873 //function : Tolerance
874 //purpose : returns the tolerance used in the approximation.
875 //=======================================================================
876 //
Tolerance() const877 Standard_Real AppDef_Variational::Tolerance() const
878 {
879 return myTolerance;
880 }
881 //
882 //=======================================================================
883 //function : NbIterations
884 //purpose : returns the number of iterations used in the approximation.
885 //=======================================================================
886 //
NbIterations() const887 Standard_Integer AppDef_Variational::NbIterations() const
888 {
889 return myNbIterations;
890 }
891 //
892 //=======================================================================
893 //function : Dump
894 //purpose : Prints on the stream o information on the current state
895 // of the object.
896 //=======================================================================
897 //
Dump(Standard_OStream & o) const898 void AppDef_Variational::Dump(Standard_OStream& o) const
899 {
900 o << " \nVariational Smoothing " << std::endl;
901 o << " Number of multipoints " << myNbPoints << std::endl;
902 o << " Number of 2d par multipoint " << myNbP2d << std::endl;
903 o << " Nombre of 3d par multipoint " << myNbP3d << std::endl;
904 o << " Number of PassagePoint " << myNbPassPoints << std::endl;
905 o << " Number of TangencyPoints " << myNbTangPoints << std::endl;
906 o << " Number of CurvaturePoints " << myNbCurvPoints << std::endl;
907 o << " \nTolerance " << o.setf(std::ios::scientific) << std::setprecision(3) << std::setw(9) << myTolerance;
908 if ( WithMinMax()) { o << " as Max Error." << std::endl;}
909 else { o << " as size Error." << std::endl;}
910 o << "CriteriumWeights : " << myPercent[0] << " , "
911 << myPercent[1] << " , " << myPercent[2] << std::endl;
912
913 if (myIsDone ) {
914 o << " MaxError " << std::setprecision(3) << std::setw(9) << myMaxError << std::endl;
915 o << " Index of MaxError " << myMaxErrorIndex << std::endl;
916 o << " Average Error " << std::setprecision(3) << std::setw(9) << myAverageError << std::endl;
917 o << " Quadratic Error " << std::setprecision(3) << std::setw(9) << myCriterium[0] << std::endl;
918 o << " Tension Criterium " << std::setprecision(3) << std::setw(9) << myCriterium[1] << std::endl;
919 o << " Flexion Criterium " << std::setprecision(3) << std::setw(9) << myCriterium[2] << std::endl;
920 o << " Jerk Criterium " << std::setprecision(3) << std::setw(9) << myCriterium[3] << std::endl;
921 o << " NbSegments " << myKnots->Length()-1 << std::endl;
922 }
923 else
924 {
925 o << (myIsOverConstr
926 ? " The problem is overconstraint"
927 : " Error in approximation") << std::endl;
928 }
929 }
930
931 //=======================================================================
932 //function : SetConstraints
933 //purpose : Define the constraints to approximate
934 // If this value is incompatible with the others fields
935 // this method modify nothing and returns false
936 //=======================================================================
937 //
SetConstraints(const Handle (AppParCurves_HArray1OfConstraintCouple)& aConstraint)938 Standard_Boolean AppDef_Variational::SetConstraints(const Handle(AppParCurves_HArray1OfConstraintCouple)& aConstraint)
939
940 {
941 myConstraints=aConstraint;
942 Init();
943 if (myIsOverConstr ) return Standard_False;
944 else return Standard_True;
945 }
946 //
947 //=======================================================================
948 //function : SetParameters
949 //purpose : Defines the parameters used by the approximations.
950 //=======================================================================
951 //
SetParameters(const Handle (TColStd_HArray1OfReal)& param)952 void AppDef_Variational::SetParameters(const Handle(TColStd_HArray1OfReal)& param)
953 {
954 myParameters->ChangeArray1() = param->Array1();
955 }
956 //
957 //=======================================================================
958 //function : SetKnots
959 //purpose : Defines the knots used by the approximations
960 // -- If this value is incompatible with the others fields
961 // -- this method modify nothing and returns false
962 //=======================================================================
963 //
SetKnots(const Handle (TColStd_HArray1OfReal)& knots)964 Standard_Boolean AppDef_Variational::SetKnots(const Handle(TColStd_HArray1OfReal)& knots)
965 {
966 myKnots->ChangeArray1() = knots->Array1();
967 return Standard_True;
968 }
969 //
970 //=======================================================================
971 //function : SetMaxDegree
972 //purpose : Define the Maximum Degree used in the approximation
973 // If this value is incompatible with the others fields
974 // this method modify nothing and returns false
975 //=======================================================================
976 //
SetMaxDegree(const Standard_Integer Degree)977 Standard_Boolean AppDef_Variational::SetMaxDegree(const Standard_Integer Degree)
978 {
979 if (((Degree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
980 return Standard_False;
981 else
982 {
983 myMaxDegree=Degree;
984
985 InitSmoothCriterion();
986
987 return Standard_True;
988 }
989
990 }
991 //
992 //=======================================================================
993 //function : SetMaxSegment
994 //purpose : Define the maximum number of segments used in the approximation
995 // If this value is incompatible with the others fields
996 // this method modify nothing and returns false
997 //=======================================================================
998 //
SetMaxSegment(const Standard_Integer NbSegment)999 Standard_Boolean AppDef_Variational::SetMaxSegment(const Standard_Integer NbSegment)
1000 {
1001 if ( myWithCutting == Standard_True &&
1002 ((myMaxDegree-myNivCont)*NbSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1003 return Standard_False;
1004 else
1005 {
1006 myMaxSegment=NbSegment;
1007 return Standard_True;
1008 }
1009 }
1010 //
1011 //=======================================================================
1012 //function : SetContinuity
1013 //purpose : Define the Continuity used in the approximation
1014 // If this value is incompatible with the others fields
1015 // this method modify nothing and returns false
1016 //=======================================================================
1017 //
SetContinuity(const GeomAbs_Shape C)1018 Standard_Boolean AppDef_Variational::SetContinuity(const GeomAbs_Shape C)
1019 {
1020 Standard_Integer NivCont=0;
1021 switch (C) {
1022 case GeomAbs_C0:
1023 NivCont=0;
1024 break ;
1025 case GeomAbs_C1:
1026 NivCont=1;
1027 break ;
1028 case GeomAbs_C2:
1029 NivCont=2;
1030 break ;
1031 default:
1032 throw Standard_ConstructionError();
1033 }
1034 if (((myMaxDegree-NivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1035 return Standard_False;
1036 else
1037 {
1038 myContinuity= C;
1039 myNivCont=NivCont;
1040
1041 InitSmoothCriterion();
1042 return Standard_True;
1043 }
1044 }
1045 //
1046 //=======================================================================
1047 //function : SetWithMinMax
1048 //purpose : Define if the approximation search to minimize the
1049 // maximum Error or not.
1050 //=======================================================================
1051 //
SetWithMinMax(const Standard_Boolean MinMax)1052 void AppDef_Variational::SetWithMinMax(const Standard_Boolean MinMax)
1053 {
1054 myWithMinMax=MinMax;
1055
1056 InitSmoothCriterion();
1057 }
1058 //
1059 //=======================================================================
1060 //function : SetWithCutting
1061 //purpose : Define if the approximation can insert new Knots or not.
1062 // If this value is incompatible with the others fields
1063 // this method modify nothing and returns false
1064 //=======================================================================
1065 //
SetWithCutting(const Standard_Boolean Cutting)1066 Standard_Boolean AppDef_Variational::SetWithCutting(const Standard_Boolean Cutting)
1067 {
1068 if (Cutting == Standard_False)
1069 {
1070 if (((myMaxDegree-myNivCont)*myKnots->Length()-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1071 return Standard_False;
1072
1073 else
1074 {
1075 myWithCutting=Cutting;
1076 InitSmoothCriterion();
1077 return Standard_True;
1078 }
1079 }
1080 else
1081 {
1082 if (((myMaxDegree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1083 return Standard_False;
1084
1085 else
1086 {
1087 myWithCutting=Cutting;
1088 InitSmoothCriterion();
1089 return Standard_True;
1090 }
1091 }
1092 }
1093 //
1094 //=======================================================================
1095 //function : SetCriteriumWeight
1096 //purpose : define the Weights (as percent) associed to the criterium used in
1097 // the optimization.
1098 //=======================================================================
1099 //
SetCriteriumWeight(const Standard_Real Percent1,const Standard_Real Percent2,const Standard_Real Percent3)1100 void AppDef_Variational::SetCriteriumWeight(const Standard_Real Percent1, const Standard_Real Percent2, const Standard_Real Percent3)
1101 {
1102 if (Percent1 < 0 || Percent2 < 0 || Percent3 < 0 ) throw Standard_DomainError();
1103 Standard_Real Total = Percent1 + Percent2 + Percent3;
1104 myPercent[0] = Percent1/Total;
1105 myPercent[1] = Percent2/Total;
1106 myPercent[2] = Percent3/Total;
1107
1108 InitSmoothCriterion();
1109 }
1110 //
1111 //=======================================================================
1112 //function : SetCriteriumWeight
1113 //purpose : define the Weight (as percent) associed to the
1114 // criterium Order used in the optimization : Others
1115 // weights are updated.
1116 //=======================================================================
1117 //
SetCriteriumWeight(const Standard_Integer Order,const Standard_Real Percent)1118 void AppDef_Variational::SetCriteriumWeight(const Standard_Integer Order, const Standard_Real Percent)
1119 {
1120 if ( Percent < 0 ) throw Standard_DomainError();
1121 if ( Order < 1 || Order > 3 ) throw Standard_ConstructionError();
1122 myPercent[Order-1] = Percent;
1123 Standard_Real Total = myPercent[0] + myPercent[1] + myPercent[2];
1124 myPercent[0] = myPercent[0] / Total;
1125 myPercent[1] = myPercent[1] / Total;
1126 myPercent[2] = myPercent[2] / Total;
1127
1128 InitSmoothCriterion();
1129
1130 }
1131 //
1132 //=======================================================================
1133 //function : SetTolerance
1134 //purpose : define the tolerance used in the approximation.
1135 //=======================================================================
1136 //
SetTolerance(const Standard_Real Tol)1137 void AppDef_Variational::SetTolerance(const Standard_Real Tol)
1138 {
1139 myTolerance=Tol;
1140 InitSmoothCriterion();
1141 }
1142 //
1143 //=======================================================================
1144 //function : SetNbIterations
1145 //purpose : define the number of iterations used in the approximation.
1146 //=======================================================================
1147 //
SetNbIterations(const Standard_Integer Iter)1148 void AppDef_Variational::SetNbIterations(const Standard_Integer Iter)
1149 {
1150 myNbIterations=Iter;
1151 }
1152
1153 //====================== Private Methodes =============================//
1154 //=======================================================================
1155 //function : TheMotor
1156 //purpose : Smoothing's motor like STRIM routine "MOTLIS"
1157 //=======================================================================
TheMotor(Handle (AppDef_SmoothCriterion)& J,const Standard_Real,const Standard_Real WQuality,Handle (FEmTool_Curve)& TheCurve,TColStd_Array1OfReal & Ecarts)1158 void AppDef_Variational::TheMotor(
1159 Handle(AppDef_SmoothCriterion)& J,
1160 // const Standard_Real WQuadratic,
1161 const Standard_Real ,
1162 const Standard_Real WQuality,
1163 Handle(FEmTool_Curve)& TheCurve,
1164 TColStd_Array1OfReal& Ecarts)
1165 {
1166 // ...
1167
1168 const Standard_Real BigValue = 1.e37, SmallValue = 1.e-6, SmallestValue = 1.e-9;
1169
1170 Handle(TColStd_HArray1OfReal) CurrentTi, NewTi, OldTi;
1171 Handle(TColStd_HArray2OfInteger) Dependence;
1172 Standard_Boolean lestim, lconst, ToOptim, iscut;
1173 Standard_Boolean isnear = Standard_False, again = Standard_True;
1174 Standard_Integer NbEst, ICDANA, NumPnt, Iter;
1175 Standard_Integer MaxNbEst =5;
1176 Standard_Real VOCRI[3] = {BigValue, BigValue, BigValue}, EROLD = BigValue,
1177 VALCRI[3], ERRMAX = BigValue, ERRMOY, ERRQUA;
1178 Standard_Real CBLONG, LNOLD;
1179 Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
1180 Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1181 Handle(FEmTool_Curve) CCurrent, COld, CNew;
1182 Standard_Real EpsLength = SmallValue;
1183 Standard_Real EpsDeg;
1184
1185 Standard_Real e1, e2, e3;
1186 Standard_Real J1min, J2min, J3min;
1187 Standard_Integer iprog;
1188
1189 // (0) Init
1190
1191 J->GetEstimation(e1, e2, e3);
1192 J1min = 1.e-8; J2min = J3min = (e1 + 1.e-8) * 1.e-6;
1193
1194 if(e1 < J1min) e1 = J1min;// Like in
1195 if(e2 < J2min) e2 = J2min;// MOTLIS
1196 if(e3 < J3min) e3 = J3min;
1197
1198
1199 J->SetEstimation(e1, e2, e3);
1200
1201 CCurrent = TheCurve;
1202 CurrentTi = new TColStd_HArray1OfReal(1, myParameters->Length());
1203 CurrentTi->ChangeArray1() = myParameters->Array1();
1204 OldTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1205 OldTi->ChangeArray1() = CurrentTi->Array1();
1206 COld = CCurrent;
1207 LNOLD = CBLONG = J->EstLength();
1208 Dependence = J->DependenceTable();
1209
1210 J->SetCurve(CCurrent);
1211 FEmTool_Assembly * TheAssembly =
1212 new FEmTool_Assembly (Dependence->Array2(), J->AssemblyTable());
1213
1214 //============ Optimization ============================
1215 // Standard_Integer inagain = 0;
1216 while (again) {
1217
1218 // (1) Loop Optimization / Estimation
1219 lestim = Standard_True;
1220 lconst = Standard_True;
1221 NbEst = 0;
1222
1223 J->SetCurve(CCurrent);
1224
1225 while(lestim) {
1226
1227 // (1.1) Curve's Optimization.
1228 EpsLength = SmallValue * CBLONG / NbrPnt;
1229 CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
1230 CCurrent->NbElements(), CCurrent->Base(), EpsLength);
1231 CNew->Knots() = CCurrent->Knots();
1232
1233 J->SetParameters(CurrentTi);
1234 EpsDeg = Min(WQuality * .1, CBLONG * .001);
1235
1236 Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
1237
1238 lconst = Standard_False;
1239
1240 // (1.2) calcul des criteres de qualites et amelioration
1241 // des estimation.
1242 ICDANA = J->QualityValues(J1min, J2min, J3min,
1243 VALCRI[0], VALCRI[1], VALCRI[2]);
1244
1245 if(ICDANA > 0) lconst = Standard_True;
1246
1247 J->ErrorValues(ERRMAX, ERRQUA, ERRMOY);
1248
1249 isnear = ((Sqrt(ERRQUA / NbrPnt) < 2*WQuality) &&
1250 (myNbIterations > 1));
1251
1252 // (1.3) Optimisation des ti par proj orthogonale
1253 // et calcul de l'erreur aux points.
1254
1255 if (isnear) {
1256 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1257 Project(CNew, CurrentTi->Array1(),
1258 NewTi->ChangeArray1(),
1259 Ecarts, NumPnt,
1260 ERRMAX, ERRQUA, ERRMOY, 2);
1261 }
1262 else NewTi = CurrentTi;
1263
1264
1265 // (1.4) Progression's test
1266 iprog = 0;
1267 if ((EROLD > WQuality) && (ERRMAX < 0.95*EROLD)) iprog++;
1268 if ((EROLD > WQuality) && (ERRMAX < 0.8*EROLD)) iprog++;
1269 if ((EROLD > WQuality) && (ERRMAX < WQuality)) iprog++;
1270 if ((EROLD > WQuality) && (ERRMAX < 0.99*EROLD)
1271 && (ERRMAX < 1.1*WQuality)) iprog++;
1272 if ( VALCRI[0] < 0.975 * VOCRI[0]) iprog++;
1273 if ( VALCRI[0] < 0.9 * VOCRI[0]) iprog++;
1274 if ( VALCRI[1] < 0.95 * VOCRI[1]) iprog++;
1275 if ( VALCRI[1] < 0.8 * VOCRI[1]) iprog++;
1276 if ( VALCRI[2] < 0.95 * VOCRI[2]) iprog++;
1277 if ( VALCRI[2] < 0.8 * VOCRI[2]) iprog++;
1278 if ((VOCRI[1]>SmallestValue)&&(VOCRI[2]>SmallestValue)) {
1279 if ((VALCRI[1]/VOCRI[1] + 2*VALCRI[2]/VOCRI[2]) < 2.8) iprog++;
1280 }
1281
1282 if (iprog < 2 && NbEst == 0 ) {
1283 // (1.5) Invalidation of new knots.
1284 VALCRI[0] = VOCRI[0];
1285 VALCRI[1] = VOCRI[1];
1286 VALCRI[2] = VOCRI[2];
1287 ERRMAX = EROLD;
1288 CBLONG = LNOLD;
1289 CCurrent = COld;
1290 CurrentTi = OldTi;
1291
1292 goto L8000; // exit
1293 }
1294
1295 VOCRI[0] = VALCRI[0];
1296 VOCRI[1] = VALCRI[1];
1297 VOCRI[2] = VALCRI[2];
1298 LNOLD = CBLONG;
1299 EROLD = ERRMAX;
1300
1301 CCurrent = CNew;
1302 CurrentTi = NewTi;
1303
1304 // (1.6) Test if the Estimations seems OK, else repeat
1305 NbEst++;
1306 lestim = ( (NbEst<MaxNbEst) && (ICDANA == 2)&& (iprog > 0) );
1307
1308 if (lestim && isnear) {
1309 // (1.7) Optimization of ti by ACR.
1310
1311 std::stable_sort (CurrentTi->begin(), CurrentTi->end());
1312
1313 Standard_Integer Decima = 4;
1314
1315 CCurrent->Length(0., 1., CBLONG);
1316 J->EstLength() = CBLONG;
1317
1318 ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
1319 lconst = Standard_True;
1320
1321 }
1322 }
1323
1324
1325 // (2) loop of parametric / geometric optimization
1326
1327 Iter = 1;
1328 ToOptim = ((Iter < myNbIterations) && (isnear));
1329
1330 while(ToOptim) {
1331 Iter++;
1332 // (2.1) Save current results
1333 VOCRI[0] = VALCRI[0];
1334 VOCRI[1] = VALCRI[1];
1335 VOCRI[2] = VALCRI[2];
1336 EROLD = ERRMAX;
1337 LNOLD = CBLONG;
1338 COld = CCurrent;
1339 OldTi->ChangeArray1() = CurrentTi->Array1();
1340
1341 // (2.2) Optimization des ti by ACR.
1342
1343 std::stable_sort (CurrentTi->begin(), CurrentTi->end());
1344
1345 Standard_Integer Decima = 4;
1346
1347 CCurrent->Length(0., 1., CBLONG);
1348 J->EstLength() = CBLONG;
1349
1350 ACR(CCurrent, CurrentTi->ChangeArray1(), Decima);
1351 lconst = Standard_True;
1352
1353 // (2.3) Optimisation des courbes
1354 EpsLength = SmallValue * CBLONG / NbrPnt;
1355
1356 CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
1357 CCurrent->NbElements(), CCurrent->Base(), EpsLength);
1358 CNew->Knots() = CCurrent->Knots();
1359
1360 J->SetParameters(CurrentTi);
1361
1362 EpsDeg = Min(WQuality * .1, CBLONG * .001);
1363 Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
1364
1365 CCurrent = CNew;
1366
1367 // (2.4) calcul des criteres de qualites et amelioration
1368 // des estimation.
1369
1370 ICDANA = J->QualityValues(J1min, J2min, J3min, VALCRI[0], VALCRI[1], VALCRI[2]);
1371 if(ICDANA > 0) lconst = Standard_True;
1372
1373 J->GetEstimation(e1, e2, e3);
1374 // (2.5) Optimisation des ti par proj orthogonale
1375
1376 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1377 Project(CCurrent, CurrentTi->Array1(),
1378 NewTi->ChangeArray1(),
1379 Ecarts, NumPnt,
1380 ERRMAX, ERRQUA, ERRMOY, 2);
1381
1382 // (2.6) Test de non regression
1383
1384 Standard_Integer iregre = 0;
1385 if (NbrConstraint < NbrPnt) {
1386 if ( (ERRMAX > WQuality) && (ERRMAX > 1.05*EROLD)) iregre++;
1387 if ( (ERRMAX > WQuality) && (ERRMAX > 2*EROLD)) iregre++;
1388 if ( (EROLD > WQuality) && (ERRMAX <= 0.5*EROLD)) iregre--;
1389 }
1390 Standard_Real E1, E2, E3;
1391 J->GetEstimation(E1, E2, E3);
1392 if ( (VALCRI[0] > E1) && (VALCRI[0] > 1.1*VOCRI[0])) iregre++;
1393 if ( (VALCRI[1] > E2) && (VALCRI[1] > 1.1*VOCRI[1])) iregre++;
1394 if ( (VALCRI[2] > E3) && (VALCRI[2] > 1.1*VOCRI[2])) iregre++;
1395
1396
1397 if (iregre >= 2) {
1398 // if (iregre >= 1) {
1399 // (2.7) on restaure l'iteration precedente
1400 VALCRI[0] = VOCRI[0];
1401 VALCRI[1] = VOCRI[1];
1402 VALCRI[2] = VOCRI[2];
1403 ERRMAX = EROLD;
1404 CBLONG = LNOLD;
1405 CCurrent = COld;
1406 CurrentTi->ChangeArray1() = OldTi->Array1();
1407 ToOptim = Standard_False;
1408 }
1409 else { // Iteration is Ok.
1410 CCurrent = CNew;
1411 CurrentTi = NewTi;
1412 }
1413 if (Iter >= myNbIterations) ToOptim = Standard_False;
1414 }
1415
1416 // (3) Decoupe eventuelle
1417
1418 if ( (CCurrent->NbElements() < myMaxSegment) && myWithCutting ) {
1419
1420 // (3.1) Sauvgarde de l'etat precedent
1421 VOCRI[0] = VALCRI[0];
1422 VOCRI[1] = VALCRI[1];
1423 VOCRI[2] = VALCRI[2];
1424 EROLD = ERRMAX;
1425 COld = CCurrent;
1426 OldTi->ChangeArray1() = CurrentTi->Array1();
1427
1428 // (3.2) On arrange les ti : Trie + recadrage sur (0,1)
1429 // ---> On trie, afin d'assurer l'ordre par la suite.
1430
1431 std::stable_sort (CurrentTi->begin(), CurrentTi->end());
1432
1433 if ((CurrentTi->Value(1)!= 0.) ||
1434 (CurrentTi->Value(NbrPnt)!= 1.)) {
1435 Standard_Real t, DelatT =
1436 1.0 /(CurrentTi->Value(NbrPnt)-CurrentTi->Value(1));
1437 for (Standard_Integer ii=2; ii<NbrPnt; ii++) {
1438 t = (CurrentTi->Value(ii)-CurrentTi->Value(1))*DelatT;
1439 CurrentTi->SetValue(ii, t);
1440 }
1441 CurrentTi->SetValue(1, 0.);
1442 CurrentTi->SetValue(NbrPnt, 1.);
1443 }
1444
1445 // (3.3) Insert new Knots
1446
1447 SplitCurve(CCurrent, CurrentTi->Array1(), EpsLength, CNew, iscut);
1448 if (!iscut) again = Standard_False;
1449 else {
1450 CCurrent = CNew;
1451 // New Knots => New Assembly.
1452 J->SetCurve(CNew);
1453 delete TheAssembly;
1454 TheAssembly = new FEmTool_Assembly (Dependence->Array2(),
1455 J->AssemblyTable());
1456 }
1457 }
1458 else { again = Standard_False;}
1459 }
1460
1461 // ================ Great loop end ===================
1462
1463 L8000:
1464
1465 // (4) Compute the best Error.
1466 NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1467 Project(CCurrent, CurrentTi->Array1(),
1468 NewTi->ChangeArray1(),
1469 Ecarts, NumPnt,
1470 ERRMAX, ERRQUA, ERRMOY, 10);
1471
1472 // (5) field's update
1473
1474 TheCurve = CCurrent;
1475 J->EstLength() = CBLONG;
1476 myParameters->ChangeArray1() = NewTi->Array1();
1477 myCriterium[0] = ERRQUA;
1478 myCriterium[1] = Sqrt(VALCRI[0]);
1479 myCriterium[2] = Sqrt(VALCRI[1]);
1480 myCriterium[3] = Sqrt(VALCRI[2]);
1481 myMaxError = ERRMAX;
1482 myMaxErrorIndex = NumPnt;
1483 if(NbrPnt > NbrConstraint)
1484 myAverageError = ERRMOY / (NbrPnt - NbrConstraint);
1485 else
1486 myAverageError = ERRMOY / NbrConstraint;
1487
1488 delete TheAssembly;
1489 }
1490
1491 //=======================================================================
1492 //function : Optimization
1493 //purpose : (like FORTRAN subroutine MINIMI)
1494 //=======================================================================
Optimization(Handle (AppDef_SmoothCriterion)& J,FEmTool_Assembly & A,const Standard_Boolean ToAssemble,const Standard_Real EpsDeg,Handle (FEmTool_Curve)& Curve,const TColStd_Array1OfReal & Parameters) const1495 void AppDef_Variational::Optimization(Handle(AppDef_SmoothCriterion)& J,
1496 FEmTool_Assembly& A,
1497 const Standard_Boolean ToAssemble,
1498 const Standard_Real EpsDeg,
1499 Handle(FEmTool_Curve)& Curve,
1500 const TColStd_Array1OfReal& Parameters) const
1501 {
1502 Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
1503 NbElm = Curve->NbElements(),
1504 NbDim = Curve->Dimension();
1505
1506 Handle(FEmTool_HAssemblyTable) AssTable;
1507
1508 math_Matrix H(0, MxDeg, 0, MxDeg);
1509 math_Vector G(0, MxDeg), Sol(1, A.NbGlobVar());
1510
1511 Standard_Integer el, dim;
1512
1513 A.GetAssemblyTable(AssTable);
1514 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1515
1516 Standard_Real CBLONG = J->EstLength();
1517
1518 // Updating Assembly
1519 if (ToAssemble)
1520 A.NullifyMatrix();
1521 A.NullifyVector();
1522
1523
1524 for (el = 1; el <= NbElm; el++) {
1525 if (ToAssemble) {
1526 J->Hessian(el, 1, 1, H);
1527 for(dim = 1; dim <= NbDim; dim++)
1528 A.AddMatrix(el, dim, dim, H);
1529 }
1530
1531 for(dim = 1; dim <= NbDim; dim++) {
1532 J->Gradient(el, dim, G);
1533 A.AddVector(el, dim, G);
1534 }
1535 }
1536
1537
1538 // Solution of system
1539 if (ToAssemble) {
1540 if(NbConstr != 0) { //Treatment of constraints
1541 AssemblingConstraints(Curve, Parameters, CBLONG, A);
1542 }
1543 A.Solve();
1544 }
1545 A.Solution(Sol);
1546
1547
1548 // Updating J
1549 J->SetCurve(Curve);
1550 J->InputVector(Sol, AssTable);
1551
1552 // Updating Curve and reduction of degree
1553
1554 Standard_Integer Newdeg;
1555 Standard_Real MaxError;
1556
1557 if(NbConstr == 0) {
1558 for(el = 1; el <= NbElm; el++) {
1559 Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
1560 }
1561 }
1562 else {
1563
1564 TColStd_Array1OfReal& TabInt = Curve->Knots();
1565 Standard_Integer Icnt = 1, p0 = Parameters.Lower() - myFirstPoint, point;
1566 for(el = 1; el <= NbElm; el++) {
1567 while((Icnt < NbConstr) &&
1568 (Parameters(p0 + myTypConstraints->Value(2 * Icnt - 1)) <= TabInt(el))) Icnt++;
1569 point = p0 + myTypConstraints->Value(2 * Icnt - 1);
1570 if(Parameters(point) <= TabInt(el) || Parameters(point) >= TabInt(el + 1))
1571 Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError);
1572 else
1573 if(Curve->Degree(el) < MxDeg) Curve->SetDegree(el, MxDeg);
1574 }
1575 }
1576 }
1577
Project(const Handle (FEmTool_Curve)& C,const TColStd_Array1OfReal & Ti,TColStd_Array1OfReal & ProjTi,TColStd_Array1OfReal & Distance,Standard_Integer & NumPoints,Standard_Real & MaxErr,Standard_Real & QuaErr,Standard_Real & AveErr,const Standard_Integer NbIterations) const1578 void AppDef_Variational::Project(const Handle(FEmTool_Curve)& C,
1579 const TColStd_Array1OfReal& Ti,
1580 TColStd_Array1OfReal& ProjTi,
1581 TColStd_Array1OfReal& Distance,
1582 Standard_Integer& NumPoints,
1583 Standard_Real& MaxErr,
1584 Standard_Real& QuaErr,
1585 Standard_Real& AveErr,
1586 const Standard_Integer NbIterations) const
1587
1588 {
1589 // Initialisation
1590
1591 const Standard_Real Seuil = 1.e-9, Eps = 1.e-12;
1592
1593 MaxErr = QuaErr = AveErr = 0.;
1594
1595 Standard_Integer Ipnt, NItCv, Iter, i, i0 = -myDimension, d0 = Distance.Lower() - 1;
1596
1597 Standard_Real TNew, Dist, T0, Dist0, F1, F2, Aux, DF, Ecart;
1598
1599 Standard_Boolean EnCour;
1600
1601 TColStd_Array1OfReal ValOfC(1, myDimension), FirstDerOfC(1, myDimension),
1602 SecndDerOfC(1, myDimension);
1603
1604 for(Ipnt = 1; Ipnt <= ProjTi.Length(); Ipnt++) {
1605
1606 i0 += myDimension;
1607
1608 TNew = Ti(Ipnt);
1609
1610 EnCour = Standard_True;
1611 NItCv = 0;
1612 Iter = 0;
1613 C->D0(TNew, ValOfC);
1614
1615 Dist = 0;
1616 for(i = 1; i <= myDimension; i++) {
1617 Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
1618 Dist += Aux * Aux;
1619 }
1620 Dist = Sqrt(Dist);
1621
1622 // ------- Newton's method for solving (C'(t),C(t) - P) = 0
1623
1624 while( EnCour ) {
1625
1626 Iter++;
1627 T0 = TNew;
1628 Dist0 = Dist;
1629
1630 C->D2(TNew, SecndDerOfC);
1631 C->D1(TNew, FirstDerOfC);
1632
1633 F1 = F2 = 0.;
1634 for(i = 1; i <= myDimension; i++) {
1635 Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
1636 DF = FirstDerOfC(i);
1637 F1 += Aux*DF; // (C'(t),C(t) - P)
1638 F2 += DF*DF + Aux * SecndDerOfC(i); // ((C'(t),C(t) - P))'
1639 }
1640
1641 if(Abs(F2) < Eps)
1642 EnCour = Standard_False;
1643 else {
1644 // Formula of Newton x(k+1) = x(k) - F(x(k))/F'(x(k))
1645 TNew -= F1 / F2;
1646 if(TNew < 0.) TNew = 0.;
1647 if(TNew > 1.) TNew = 1.;
1648
1649
1650 // Analysis of result
1651
1652 C->D0(TNew, ValOfC);
1653
1654 Dist = 0;
1655 for(i = 1; i <= myDimension; i++) {
1656 Aux = ValOfC(i) - myTabPoints->Value(i0 + i);
1657 Dist += Aux * Aux;
1658 }
1659 Dist = Sqrt(Dist);
1660
1661 Ecart = Dist0 - Dist;
1662
1663 if(Ecart <= -Seuil) {
1664 // Pas d'amelioration on s'arrete
1665 EnCour = Standard_False;
1666 TNew = T0;
1667 Dist = Dist0;
1668 }
1669 else if(Ecart <= Seuil)
1670 // Convergence
1671 NItCv++;
1672 else
1673 NItCv = 0;
1674
1675 if((NItCv >= 2) || (Iter >= NbIterations)) EnCour = Standard_False;
1676
1677 }
1678 }
1679
1680
1681 ProjTi(Ipnt) = TNew;
1682 Distance(d0 + Ipnt) = Dist;
1683 if(Dist > MaxErr) {
1684 MaxErr = Dist;
1685 NumPoints = Ipnt;
1686 }
1687 QuaErr += Dist * Dist;
1688 AveErr += Dist;
1689 }
1690
1691 NumPoints = NumPoints + myFirstPoint - 1;// Setting NumPoints to interval [myFirstPoint, myLastPoint]
1692
1693 }
1694
ACR(Handle (FEmTool_Curve)& Curve,TColStd_Array1OfReal & Ti,const Standard_Integer Decima) const1695 void AppDef_Variational::ACR(Handle(FEmTool_Curve)& Curve,
1696 TColStd_Array1OfReal& Ti,
1697 const Standard_Integer Decima) const
1698 {
1699
1700 const Standard_Real Eps = 1.e-8;
1701
1702 TColStd_Array1OfReal& Knots = Curve->Knots();
1703 Standard_Integer NbrPnt = Ti.Length(), TiFirst = Ti.Lower(), TiLast = Ti.Upper(),
1704 KFirst = Knots.Lower(), KLast = Knots.Upper();
1705
1706 Standard_Real CbLong, DeltaT, VTest, UNew, UOld, DU, TPara, TOld, DTInv, Ratio;
1707 Standard_Integer ipnt, ii, IElm, IOld, POld, PCnt, ICnt=0;
1708 Standard_Integer NbCntr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1709
1710 // (1) Calcul de la longueur de courbe
1711
1712 Curve->Length(Ti(TiFirst), Ti(TiLast), CbLong);
1713
1714 // (2) Mise de l'acr dans Ti
1715
1716 if(NbrPnt >= 2) {
1717
1718 // (2.0) Initialisation
1719 DeltaT = (Ti(TiLast) - Ti(TiFirst)) / Decima;
1720 VTest = Ti(TiFirst) + DeltaT;
1721
1722 if(NbCntr > 0) {
1723 PCnt = myTypConstraints->Value(1) - myFirstPoint + TiFirst;
1724 ICnt = 1;
1725 }
1726 else
1727 PCnt = TiLast + 1;
1728
1729 UOld = 0.;
1730
1731 TOld = Ti(TiFirst);
1732 POld = TiFirst;
1733
1734 IElm = KFirst;
1735 IOld = IElm;
1736
1737 Ti(TiFirst) = 0.;
1738
1739 for(ipnt = TiFirst + 1; ipnt <= TiLast; ipnt++) {
1740
1741 while((ICnt <= NbCntr) && (PCnt < ipnt)) {
1742 ICnt++;
1743 PCnt = myTypConstraints->Value(2*ICnt-1) - myFirstPoint + TiFirst;
1744 }
1745
1746 TPara = Ti(ipnt);
1747
1748 if(TPara >= VTest || PCnt == ipnt) {
1749
1750 if ( Ti(TiLast) - TPara <= 1.e-2*DeltaT) {
1751 ipnt = TiLast;
1752 TPara = Ti(ipnt);
1753 }
1754 // (2.2), (2.3) Cacul la longueur de courbe
1755 Curve->Length(Ti(TiFirst), TPara, UNew);
1756
1757 UNew /= CbLong;
1758
1759 while(Knots(IElm + 1) < TPara && IElm < KLast - 1) IElm++;
1760
1761 // (2.4) Mise a jours des parametres de decoupe
1762 DTInv = 1. / (TPara - TOld);
1763 DU = UNew - UOld;
1764
1765 for(ii = IOld+1; ii <= IElm; ii++) {
1766 Ratio = (Knots(ii) - TOld) * DTInv;
1767 Knots(ii) = UOld + Ratio * DU;
1768 }
1769
1770 // (2.5) Mise a jours des parametres de points.
1771
1772 //Very strange loop, because it never works (POld+1 > ipnt-1)
1773 for(ii = POld+1; ii <= ipnt-1; ii++) {
1774 Ratio = ( Ti(ii) - TOld ) * DTInv;
1775 Ti(ii) = UOld + Ratio * DU;
1776 }
1777
1778 Ti(ipnt) = UNew;
1779
1780 UOld = UNew;
1781 IOld = IElm;
1782 TOld = TPara;
1783 POld = ipnt;
1784
1785 }
1786 // --> Nouveau seuil parametrique pour le decimage
1787
1788 if(TPara >= VTest) {
1789 // ii = RealToInt((TPara - VTest + Eps) / DeltaT);
1790 // VTest += (ii + 1) * DeltaT;
1791 VTest += Ceiling((TPara - VTest + Eps) / DeltaT) * DeltaT;
1792 if(VTest > 1. - Eps) VTest = 1.;
1793 }
1794 }
1795 }
1796
1797 // --- On ajuste les valeurs extremes
1798
1799 Ti(TiFirst) = 0.;
1800 Ti(TiLast) = 1.;
1801 ii = TiLast - 1;
1802 while ( Ti(ii) > Knots(KLast) ) {
1803 Ti(ii) = 1.;
1804 --ii;
1805 }
1806 Knots(KFirst) = 0.;
1807 Knots(KLast) = 1.;
1808
1809 }
1810
1811 //----------------------------------------------------------//
1812 // Standard_Integer NearIndex //
1813 // Purpose: searching nearest index of TabPar corresponding //
1814 // given value T (is similar to fortran routine MSRRE2) //
1815 //----------------------------------------------------------//
1816
NearIndex(const Standard_Real T,const TColStd_Array1OfReal & TabPar,const Standard_Real Eps,Standard_Integer & Flag)1817 static Standard_Integer NearIndex(const Standard_Real T,
1818 const TColStd_Array1OfReal& TabPar,
1819 const Standard_Real Eps, Standard_Integer& Flag)
1820 {
1821 Standard_Integer Loi = TabPar.Lower(), Upi = TabPar.Upper();
1822
1823 Flag = 0;
1824
1825 if(T < TabPar(Loi)) { Flag = -1; return Loi;}
1826 if(T > TabPar(Upi)) { Flag = 1; return Upi;}
1827
1828 Standard_Integer Ibeg = Loi, Ifin = Upi, Imidl;
1829
1830 while(Ibeg + 1 != Ifin) {
1831 Imidl = (Ibeg + Ifin) / 2;
1832 if((T >= TabPar(Ibeg)) && (T <= TabPar(Imidl)))
1833 Ifin = Imidl;
1834 else
1835 Ibeg = Imidl;
1836 }
1837
1838 if(Abs(T - TabPar(Ifin)) < Eps) return Ifin;
1839
1840 return Ibeg;
1841 }
1842
1843
1844 //----------------------------------------------------------//
1845 // void GettingKnots //
1846 // Purpose: calculating values of new Knots for elements //
1847 // with degree that is equal Deg //
1848 //----------------------------------------------------------//
1849
GettingKnots(const TColStd_Array1OfReal & TabPar,const Handle (FEmTool_Curve)& InCurve,const Standard_Integer Deg,Standard_Integer & NbElm,TColStd_Array1OfReal & NewKnots)1850 static void GettingKnots(const TColStd_Array1OfReal& TabPar,
1851 const Handle(FEmTool_Curve)& InCurve,
1852 const Standard_Integer Deg,
1853 Standard_Integer& NbElm,
1854 TColStd_Array1OfReal& NewKnots)
1855
1856 {
1857
1858 const Standard_Real Eps = 1.e-12;
1859
1860 TColStd_Array1OfReal& OldKnots = InCurve->Knots();
1861 Standard_Integer NbMaxOld = InCurve->NbElements();
1862 Standard_Integer NbMax = NewKnots.Upper(), Ipt, Ipt1, Ipt2;
1863 Standard_Integer el = 0, i1 = OldKnots.Lower(), i0 = i1 - 1, Flag;
1864 Standard_Real TPar;
1865
1866 while((NbElm < NbMax) && (el < NbMaxOld)) {
1867
1868 el++; i0++; i1++; // i0, i1 are indexes of left and right knots of element el
1869
1870 if(InCurve->Degree(el) == Deg) {
1871
1872 NbElm++;
1873
1874 Ipt1 = NearIndex(OldKnots(i0), TabPar, Eps, Flag);
1875 if(Flag != 0) Ipt1 = TabPar.Lower();
1876 Ipt2 = NearIndex(OldKnots(i1), TabPar, Eps, Flag);
1877 if(Flag != 0) Ipt2 = TabPar.Upper();
1878
1879 if(Ipt2 - Ipt1 >= 1) {
1880
1881 Ipt = (Ipt1 + Ipt2) / 2;
1882 if(2 * Ipt == Ipt1 + Ipt2)
1883 TPar = 2. * TabPar(Ipt);
1884 else
1885 TPar = TabPar(Ipt) + TabPar(Ipt + 1);
1886
1887 NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1) + TPar) / 4.;
1888 }
1889 else
1890 NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1)) / 2.;
1891 }
1892 }
1893 }
1894
SplitCurve(const Handle (FEmTool_Curve)& InCurve,const TColStd_Array1OfReal & Ti,const Standard_Real CurveTol,Handle (FEmTool_Curve)& OutCurve,Standard_Boolean & iscut) const1895 void AppDef_Variational::SplitCurve(const Handle(FEmTool_Curve)& InCurve,
1896 const TColStd_Array1OfReal& Ti,
1897 const Standard_Real CurveTol,
1898 Handle(FEmTool_Curve)& OutCurve,
1899 Standard_Boolean& iscut) const
1900 {
1901 Standard_Integer NbElmOld = InCurve->NbElements();
1902
1903 if(NbElmOld >= myMaxSegment) {iscut = Standard_False; return;}
1904 #ifdef OCCT_DEBUG
1905 Standard_Integer MaxDegree =
1906 #endif
1907 InCurve->Base()->WorkDegree();
1908 Standard_Integer NbElm = NbElmOld;
1909 TColStd_Array1OfReal NewKnots(NbElm + 1, myMaxSegment);
1910 #ifndef OCCT_DEBUG
1911 GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree(), NbElm, NewKnots);
1912 GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree() - 1, NbElm, NewKnots);
1913 #else
1914 GettingKnots(Ti, InCurve, MaxDegree, NbElm, NewKnots);
1915 GettingKnots(Ti, InCurve, MaxDegree - 1, NbElm, NewKnots);
1916
1917 #endif
1918 if(NbElm > NbElmOld) {
1919
1920 iscut = Standard_True;
1921
1922 OutCurve = new FEmTool_Curve(InCurve->Dimension(), NbElm, InCurve->Base(), CurveTol);
1923 TColStd_Array1OfReal& OutKnots = OutCurve->Knots();
1924 TColStd_Array1OfReal& InKnots = InCurve->Knots();
1925
1926 Standard_Integer i, i0 = OutKnots.Lower();
1927 for(i = InKnots.Lower(); i <= InKnots.Upper(); i++) OutKnots(i) = InKnots(i);
1928 for(i = NbElmOld + 1; i <= NbElm; i++) OutKnots(i + i0) = NewKnots(i);
1929
1930 std::sort (OutKnots.begin(), OutKnots.end());
1931 }
1932 else
1933 iscut = Standard_False;
1934
1935 }
1936
1937 //=======================================================================
1938 //function : InitSmoothCriterion
1939 //purpose : Initializes the SmoothCriterion
1940 //=======================================================================
InitSmoothCriterion()1941 void AppDef_Variational::InitSmoothCriterion()
1942 {
1943
1944 const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9;
1945 // const Standard_Real J1 = .01, J2 = .001, J3 = .001;
1946
1947
1948
1949 Standard_Real Length;
1950
1951 InitParameters(Length);
1952
1953 mySmoothCriterion->SetParameters(myParameters);
1954
1955 Standard_Real E1, E2, E3;
1956
1957 InitCriterionEstimations(Length, E1, E2, E3);
1958 /*
1959 J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6;
1960
1961 if(E1 < J1) E1 = J1;
1962 if(E2 < J2) E2 = J2;
1963 if(E3 < J3) E3 = J3;
1964 */
1965 mySmoothCriterion->EstLength() = Length;
1966 mySmoothCriterion->SetEstimation(E1, E2, E3);
1967
1968 Standard_Real WQuadratic, WQuality;
1969
1970 if(!myWithMinMax && myTolerance != 0.)
1971 WQuality = myTolerance;
1972 else if (myTolerance == 0.)
1973 WQuality = 1.;
1974 else
1975 WQuality = Max(myTolerance, Eps2 * Length);
1976
1977 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1978 WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality;
1979 if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic;
1980
1981 if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.);
1982
1983
1984 mySmoothCriterion->SetWeight(WQuadratic, WQuality,
1985 myPercent[0], myPercent[1], myPercent[2]);
1986
1987
1988 Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity);
1989 Handle(FEmTool_Curve) TheCurve;
1990 Standard_Integer NbElem;
1991 Standard_Real CurvTol = Eps2 * Length / myNbPoints;
1992
1993 // Decoupe de l'intervalle en fonction des contraintes
1994 if(myWithCutting == Standard_True && NbConstr != 0) {
1995
1996 InitCutting(TheBase, CurvTol, TheCurve);
1997
1998 }
1999 else {
2000
2001 NbElem = 1;
2002 TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol);
2003 TheCurve->Knots().SetValue(TheCurve->Knots().Lower(),
2004 myParameters->Value(myFirstPoint));
2005 TheCurve->Knots().SetValue(TheCurve->Knots().Upper(),
2006 myParameters->Value(myLastPoint));
2007
2008 }
2009
2010 mySmoothCriterion->SetCurve(TheCurve);
2011
2012 return;
2013
2014 }
2015
2016 //
2017 //=======================================================================
2018 //function : InitParameters
2019 //purpose : Calculation of initial estimation of the multicurve's length
2020 // and parameters for multipoints.
2021 //=======================================================================
2022 //
InitParameters(Standard_Real & Length)2023 void AppDef_Variational::InitParameters(Standard_Real& Length)
2024 {
2025
2026 const Standard_Real Eps1 = Precision::Confusion() * .01;
2027
2028 Standard_Real aux, dist;
2029 Standard_Integer i, i0, i1 = 0, ipoint;
2030
2031
2032 Length = 0.;
2033 myParameters->SetValue(myFirstPoint, Length);
2034
2035 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) {
2036 i0 = i1; i1 += myDimension;
2037 dist = 0;
2038 for(i = 1; i <= myDimension; i++) {
2039 aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i);
2040 dist += aux * aux;
2041 }
2042 Length += Sqrt(dist);
2043 myParameters->SetValue(ipoint, Length);
2044 }
2045
2046
2047 if(Length <= Eps1)
2048 throw Standard_ConstructionError("AppDef_Variational::InitParameters");
2049
2050
2051 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++)
2052 myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length);
2053
2054 myParameters->SetValue(myLastPoint, 1.);
2055
2056 // Avec peu de point il y a sans doute sous estimation ...
2057 if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1));
2058 }
2059
2060 //
2061 //=======================================================================
2062 //function : InitCriterionEstimations
2063 //purpose : Calculation of initial estimation of the linear criteria
2064 //
2065 //=======================================================================
2066 //
InitCriterionEstimations(const Standard_Real Length,Standard_Real & E1,Standard_Real & E2,Standard_Real & E3) const2067 void AppDef_Variational::InitCriterionEstimations(const Standard_Real Length,
2068 Standard_Real& E1,
2069 Standard_Real& E2,
2070 Standard_Real& E3) const
2071 {
2072 E1 = Length * Length;
2073
2074 const Standard_Real Eps1 = Precision::Confusion() * .01;
2075
2076 math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension),
2077 VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension);
2078
2079 // ========== Treatment of first point =================
2080
2081 Standard_Integer ipnt = myFirstPoint;
2082
2083 EstTangent(ipnt, VTang1);
2084 ipnt++;
2085 EstTangent(ipnt, VTang2);
2086 ipnt++;
2087 EstTangent(ipnt, VTang3);
2088
2089 ipnt = myFirstPoint;
2090 EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1);
2091 ipnt++;
2092 EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2);
2093
2094 // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 Begin
2095 // Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt));
2096 Standard_Integer anInd = ipnt;
2097 Standard_Real Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt));
2098 // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 End
2099
2100 if(Delta <= Eps1) Delta = 1.;
2101
2102 E2 = VScnd1.Norm2() * Delta;
2103
2104 E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
2105 // ========== Treatment of internal points =================
2106
2107 Standard_Integer CurrPoint = 2;
2108
2109 for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) {
2110
2111 Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1));
2112
2113 if(CurrPoint == 1) {
2114 if(ipnt + 1 != myLastPoint) {
2115 EstTangent(ipnt + 2, VTang3);
2116 EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2);
2117 }
2118 else
2119 EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2);
2120
2121 E2 += VScnd1.Norm2() * Delta;
2122 E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.;
2123
2124 }
2125 else if(CurrPoint == 2) {
2126 if(ipnt + 1 != myLastPoint) {
2127 EstTangent(ipnt + 2, VTang1);
2128 EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3);
2129 }
2130 else
2131 EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3);
2132
2133 E2 += VScnd2.Norm2() * Delta;
2134 E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
2135
2136 }
2137 else {
2138 if(ipnt + 1 != myLastPoint) {
2139 EstTangent(ipnt + 2, VTang2);
2140 EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1);
2141 }
2142 else
2143 EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1);
2144
2145 E2 += VScnd3.Norm2() * Delta;
2146 E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.;
2147
2148 }
2149
2150 CurrPoint++; if(CurrPoint == 4) CurrPoint = 1;
2151 }
2152
2153 // ========== Treatment of last point =================
2154
2155 Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1));
2156 if(Delta <= Eps1) Delta = 1.;
2157
2158 Standard_Real aux;
2159
2160 if(CurrPoint == 1) {
2161
2162 E2 += VScnd1.Norm2() * Delta;
2163 aux = VScnd1.Subtracted(VScnd3).Norm2();
2164 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
2165
2166 }
2167 else if(CurrPoint == 2) {
2168
2169 E2 += VScnd2.Norm2() * Delta;
2170 aux = VScnd2.Subtracted(VScnd1).Norm2();
2171 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
2172
2173 }
2174 else {
2175
2176 E2 += VScnd3.Norm2() * Delta;
2177 aux = VScnd3.Subtracted(VScnd2).Norm2();
2178 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
2179
2180 }
2181
2182 aux = Length * Length;
2183
2184 E2 *= aux; E3 *= aux;
2185
2186
2187 }
2188
2189 //
2190 //=======================================================================
2191 //function : EstTangent
2192 //purpose : Calculation of estimation of the Tangent
2193 // (see fortran routine MMLIPRI)
2194 //=======================================================================
2195 //
2196
EstTangent(const Standard_Integer ipnt,math_Vector & VTang) const2197 void AppDef_Variational::EstTangent(const Standard_Integer ipnt,
2198 math_Vector& VTang) const
2199
2200 {
2201 Standard_Integer i ;
2202 const Standard_Real Eps1 = Precision::Confusion() * .01;
2203 const Standard_Real EpsNorm = 1.e-9;
2204
2205 Standard_Real Wpnt = 1.;
2206
2207
2208 if(ipnt == myFirstPoint) {
2209 // Estimation at first point
2210 if(myNbPoints < 3)
2211 Wpnt = 0.;
2212 else {
2213
2214 Standard_Integer adr1 = 1, adr2 = adr1 + myDimension,
2215 adr3 = adr2 + myDimension;
2216
2217 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
2218 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
2219 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
2220
2221 // Parabolic interpolation
2222 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
2223 // first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d))
2224 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d
2225 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
2226 Standard_Real V2 = 0.;
2227 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
2228 if(V2 > Eps1) {
2229 Standard_Real d = V1 / (V1 + V2), d1;
2230 d1 = 1. / (d * (1 - d)); d *= d;
2231 VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1;
2232 }
2233 else {
2234 // Simple 2-point estimation
2235
2236 VTang = Pnt2 - Pnt1;
2237 }
2238 }
2239 }
2240 else if (ipnt == myLastPoint) {
2241 // Estimation at last point
2242 if(myNbPoints < 3)
2243 Wpnt = 0.;
2244 else {
2245
2246 Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1,
2247 adr2 = adr1 + myDimension,
2248 adr3 = adr2 + myDimension;
2249
2250 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
2251 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
2252 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
2253
2254 // Parabolic interpolation
2255 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
2256 // first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d))
2257 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2)
2258 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
2259 Standard_Real V2 = 0.;
2260 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
2261 if(V2 > Eps1) {
2262 Standard_Real d = V1 / (V1 + V2), d1;
2263 d1 = 1. / (d * (1 - d)); d *= d - 2;
2264 VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1;
2265 }
2266 else {
2267 // Simple 2-point estimation
2268
2269 VTang = Pnt3 - Pnt2;
2270 }
2271 }
2272 }
2273 else {
2274
2275 Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1,
2276 adr2 = adr1 + 2 * myDimension;
2277
2278 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
2279 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
2280
2281 VTang = Pnt2 - Pnt1;
2282
2283 }
2284
2285 Standard_Real Vnorm = VTang.Norm();
2286
2287 if(Vnorm <= EpsNorm)
2288 VTang.Init(0.);
2289 else
2290 VTang /= Vnorm;
2291
2292 // Estimation with constraints
2293
2294 Standard_Real Wcnt = 0.;
2295 Standard_Integer IdCnt = 1;
2296
2297 // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
2298
2299 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2300 math_Vector VCnt(1, myDimension, 0.);
2301
2302 if(NbConstr > 0) {
2303
2304 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
2305 IdCnt <= NbConstr) IdCnt++;
2306 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
2307 (myTypConstraints->Value(2 * IdCnt) >= 1)) {
2308 Wcnt = 1.;
2309 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0;
2310 for( i = 1; i <= myNbP3d; i++) {
2311 for(Standard_Integer j = 1; j <= 3; j++)
2312 VCnt(++k) = myTabConstraints->Value(++i0);
2313 i0 += 3;
2314 }
2315 for(i = 1; i <= myNbP2d; i++) {
2316 for(Standard_Integer j = 1; j <= 2; j++)
2317 VCnt(++k) = myTabConstraints->Value(++i0);
2318 i0 += 2;
2319 }
2320 }
2321 }
2322
2323 // Averaging of estimation
2324
2325 Standard_Real Denom = Wpnt + Wcnt;
2326 if(Denom == 0.) Denom = 1.;
2327 else Denom = 1. / Denom;
2328
2329 VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom;
2330
2331 Vnorm = VTang.Norm();
2332
2333 if(Vnorm <= EpsNorm)
2334 VTang.Init(0.);
2335 else
2336 VTang /= Vnorm;
2337
2338
2339 }
2340
2341 //
2342 //=======================================================================
2343 //function : EstSecnd
2344 //purpose : Calculation of estimation of the second derivative
2345 // (see fortran routine MLIMSCN)
2346 //=======================================================================
2347 //
EstSecnd(const Standard_Integer ipnt,const math_Vector & VTang1,const math_Vector & VTang2,const Standard_Real Length,math_Vector & VScnd) const2348 void AppDef_Variational::EstSecnd(const Standard_Integer ipnt,
2349 const math_Vector& VTang1,
2350 const math_Vector& VTang2,
2351 const Standard_Real Length,
2352 math_Vector& VScnd) const
2353 {
2354 Standard_Integer i ;
2355
2356 const Standard_Real Eps = 1.e-9;
2357
2358 Standard_Real Wpnt = 1.;
2359
2360 Standard_Real aux;
2361
2362 if(ipnt == myFirstPoint)
2363 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt);
2364 else if(ipnt == myLastPoint)
2365 aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1);
2366 else
2367 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1);
2368
2369 if(aux <= Eps)
2370 aux = 1.;
2371 else
2372 aux = 1. / aux;
2373
2374 VScnd = (VTang2 - VTang1) * aux;
2375
2376 // Estimation with constraints
2377
2378 Standard_Real Wcnt = 0.;
2379 Standard_Integer IdCnt = 1;
2380
2381 // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
2382
2383 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2384 math_Vector VCnt(1, myDimension, 0.);
2385
2386 if(NbConstr > 0) {
2387
2388 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
2389 IdCnt <= NbConstr) IdCnt++;
2390
2391 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
2392 (myTypConstraints->Value(2 * IdCnt) >= 2)) {
2393 Wcnt = 1.;
2394 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0;
2395 for( i = 1; i <= myNbP3d; i++) {
2396 for(Standard_Integer j = 1; j <= 3; j++)
2397 VCnt(++k) = myTabConstraints->Value(++i0);
2398 i0 += 3;
2399 }
2400 i0--;
2401 for(i = 1; i <= myNbP2d; i++) {
2402 for(Standard_Integer j = 1; j <= 2; j++)
2403 VCnt(++k) = myTabConstraints->Value(++i0);
2404 i0 += 2;
2405 }
2406 }
2407 }
2408
2409 // Averaging of estimation
2410
2411 Standard_Real Denom = Wpnt + Wcnt;
2412 if(Denom == 0.) Denom = 1.;
2413 else Denom = 1. / Denom;
2414
2415 VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom;
2416
2417 }
2418
2419
2420 //
2421 //=======================================================================
2422 //function : InitCutting
2423 //purpose : Realisation of curve's cutting a priory accordingly to
2424 // constraints (see fortran routine MLICUT)
2425 //=======================================================================
2426 //
InitCutting(const Handle (PLib_Base)& aBase,const Standard_Real CurvTol,Handle (FEmTool_Curve)& aCurve) const2427 void AppDef_Variational::InitCutting(const Handle(PLib_Base)& aBase,
2428 const Standard_Real CurvTol,
2429 Handle(FEmTool_Curve)& aCurve) const
2430 {
2431
2432 // Definition of number of elements
2433 Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem;
2434 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2435
2436 for(i = 1; i <= NbConstr; i++) {
2437 kk = Abs(myTypConstraints->Value(2 * i)) + 1;
2438 ORCMx = Max(ORCMx, kk);
2439 NCont += kk;
2440 }
2441
2442 if(ORCMx > myMaxDegree - myNivCont)
2443 throw Standard_ConstructionError("AppDef_Variational::InitCutting");
2444
2445 Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4,
2446 myNivCont + 1);
2447
2448 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
2449
2450 while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) {
2451
2452 NLibre++;
2453 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
2454
2455 }
2456
2457
2458 if(NbElem > myMaxSegment)
2459 throw Standard_ConstructionError("AppDef_Variational::InitCutting");
2460
2461
2462 aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol);
2463
2464 Standard_Integer NCnt = (NCont - 1) / NbElem + 1;
2465 Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont);
2466
2467 TColStd_Array1OfReal& Knot = aCurve->Knots();
2468
2469 Standard_Integer IDeb = 0, IFin = NbConstr + 1,
2470 NDeb = 0, NFin = 0,
2471 IndEl = Knot.Lower(), IUpper = Knot.Upper(),
2472 NbEl = 0;
2473
2474
2475 Knot(IndEl) = myParameters->Value(myFirstPoint);
2476 Knot(IUpper) = myParameters->Value(myLastPoint);
2477
2478 while(NbElem - NbEl > 1) {
2479
2480 IndEl++; NbEl++;
2481 if(NPlus == 0) NCnt--;
2482
2483 while(NDeb < NCnt && IDeb < IFin) {
2484 IDeb++;
2485 NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1;
2486 }
2487
2488 if(NDeb == NCnt) {
2489 NDeb = 0;
2490 if(NPlus == 1 &&
2491 myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1))
2492
2493 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
2494 else
2495 Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) +
2496 myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2;
2497
2498 }
2499 else {
2500 NDeb -= NCnt;
2501 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
2502 }
2503
2504 NPlus--;
2505 if(NPlus == 0) NCnt--;
2506
2507 if(NbElem - NbEl == 1) break;
2508
2509 NbEl++;
2510
2511 while(NFin < NCnt && IDeb < IFin) {
2512 IFin--;
2513 NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1;
2514 }
2515
2516 if(NFin == NCnt) {
2517 NFin = 0;
2518 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
2519 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
2520 }
2521 else {
2522 NFin -= NCnt;
2523 if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1))
2524 Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1));
2525 else
2526 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
2527 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
2528 }
2529 }
2530 }
2531
2532 //=======================================================================
2533 //function : Adjusting
2534 //purpose : Smoothing's adjusting like STRIM routine "MAJLIS"
2535 //=======================================================================
Adjusting(Handle (AppDef_SmoothCriterion)& J,Standard_Real & WQuadratic,Standard_Real & WQuality,Handle (FEmTool_Curve)& TheCurve,TColStd_Array1OfReal & Ecarts)2536 void AppDef_Variational::Adjusting(
2537 Handle(AppDef_SmoothCriterion)& J,
2538 Standard_Real& WQuadratic,
2539 Standard_Real& WQuality,
2540 Handle(FEmTool_Curve)& TheCurve,
2541 TColStd_Array1OfReal& Ecarts)
2542 {
2543
2544 // std::cout << "=========== Adjusting =============" << std::endl;
2545
2546 /* Initialized data */
2547
2548 const Standard_Integer mxiter = 2;
2549 const Standard_Real eps1 = 1e-6;
2550 Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
2551 Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2552 Standard_Real CurvTol = eps1 * J->EstLength() / NbrPnt;
2553
2554
2555 /* Local variables */
2556 Standard_Integer iter, ipnt;
2557 Standard_Real ecart, emold, erold, tpara;
2558 Standard_Real vocri[4], j1cibl, vtest, vseuil;
2559 Standard_Integer i, numint, flag;
2560 TColStd_Array1OfReal tbpoid(myFirstPoint, myLastPoint);
2561 Standard_Boolean loptim, lrejet;
2562 Handle(AppDef_SmoothCriterion) JNew;
2563 Handle(FEmTool_Curve) CNew;
2564 Standard_Real E1, E2, E3;
2565
2566
2567 /* (0.b) Initialisations */
2568
2569 loptim = Standard_True;
2570 iter = 0;
2571 tbpoid.Init(1.);
2572
2573
2574 /* ============ boucle sur le moteur de lissage ============== */
2575
2576 vtest = WQuality * .9;
2577 j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
2578
2579 while(loptim) {
2580
2581 ++iter;
2582
2583 /* (1) Sauvegarde de l'etat precedents */
2584
2585 vocri[0] = myCriterium[0];
2586 vocri[1] = myCriterium[1];
2587 vocri[2] = myCriterium[2];
2588 vocri[3] = myCriterium[3];
2589 erold = myMaxError;
2590 emold = myAverageError;
2591
2592 /* (2) Augmentation du poids des moindre carre */
2593
2594 if (j1cibl > vtest) {
2595 WQuadratic = j1cibl / vtest * WQuadratic;
2596 }
2597
2598 /* (3) Augmentation du poid associe aux points a problemes */
2599
2600 vseuil = WQuality * .88;
2601
2602 for (ipnt = myFirstPoint; ipnt <= myLastPoint; ++ipnt) {
2603 if (Ecarts(ipnt) > vtest) {
2604 ecart = (Ecarts(ipnt) - vseuil) / WQuality;
2605 tbpoid(ipnt) = (ecart * 3 + 1.) * tbpoid(ipnt);
2606 }
2607 }
2608
2609 /* (4) Decoupe force */
2610
2611 if (TheCurve->NbElements() < myMaxSegment && myWithCutting) {
2612
2613 numint = NearIndex(myParameters->Value(myMaxErrorIndex), TheCurve->Knots(), 0, flag);
2614
2615 tpara = (TheCurve->Knots()(numint) + TheCurve->Knots()(numint + 1) +
2616 myParameters->Value(myMaxErrorIndex) * 2) / 4;
2617
2618 CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements() + 1,
2619 TheCurve->Base(), CurvTol);
2620
2621 for(i = 1; i <= numint; i++) CNew->Knots()(i) = TheCurve->Knots()(i);
2622 for(i = numint + 1; i <= TheCurve->Knots().Length(); i++)
2623 CNew->Knots()(i + 1) = TheCurve->Knots()(i);
2624
2625 CNew->Knots()(numint + 1) = tpara;
2626
2627 } else {
2628
2629 CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements(), TheCurve->Base(), CurvTol);
2630
2631 CNew->Knots() = TheCurve->Knots();
2632 }
2633
2634
2635 JNew = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint);
2636
2637 JNew->EstLength() = J->EstLength();
2638
2639 J->GetEstimation(E1, E2, E3);
2640
2641 JNew->SetEstimation(E1, E2, E3);
2642
2643 JNew->SetParameters(myParameters);
2644
2645 JNew->SetWeight(WQuadratic, WQuality, myPercent[0], myPercent[1], myPercent[2]);
2646
2647 JNew->SetWeight(tbpoid);
2648
2649 JNew->SetCurve(CNew);
2650
2651 /* (5) Relissage */
2652
2653 TheMotor(JNew, WQuadratic, WQuality, CNew, Ecarts);
2654
2655 /* (6) Tests de rejet */
2656
2657 j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
2658 vseuil = Sqrt(vocri[1]) + (erold - myMaxError) * 4;
2659
2660 lrejet = ((myMaxError > WQuality) && (myMaxError > erold * 1.01)) || (Sqrt(myCriterium[1]) > vseuil * 1.05);
2661
2662 if (lrejet) {
2663 myCriterium[0] = vocri[0];
2664 myCriterium[1] = vocri[1];
2665 myCriterium[2] = vocri[2];
2666 myCriterium[3] = vocri[3];
2667 myMaxError = erold;
2668 myAverageError = emold;
2669
2670 loptim = Standard_False;
2671 }
2672 else {
2673 J = JNew;
2674 TheCurve = CNew;
2675 J->SetCurve(TheCurve);
2676 }
2677
2678 /* (7) Test de convergence */
2679
2680 if (((iter >= mxiter) && (myMaxSegment == CNew->NbElements())) || myMaxError < WQuality) {
2681 loptim = Standard_False;
2682 }
2683 }
2684 }
2685
NotParallel(gp_Vec & T,gp_Vec & V)2686 static Standard_Boolean NotParallel(gp_Vec& T, gp_Vec& V)
2687 {
2688 V = T;
2689 V.SetX(V.X() + 1.);
2690 if (V.CrossMagnitude(T) > 1.e-12)
2691 return Standard_True;
2692 V.SetY(V.Y() + 1.);
2693 if (V.CrossMagnitude(T) > 1.e-12)
2694 return Standard_True;
2695 V.SetZ(V.Z() + 1.);
2696 if (V.CrossMagnitude(T) > 1.e-12)
2697 return Standard_True;
2698 return Standard_False;
2699 }
2700
AssemblingConstraints(const Handle (FEmTool_Curve)& Curve,const TColStd_Array1OfReal & Parameters,const Standard_Real CBLONG,FEmTool_Assembly & A) const2701 void AppDef_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve,
2702 const TColStd_Array1OfReal& Parameters,
2703 const Standard_Real CBLONG,
2704 FEmTool_Assembly& A ) const
2705 {
2706
2707 Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
2708 NbElm = Curve->NbElements(),
2709 NbDim = Curve->Dimension();
2710
2711 TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg);
2712 math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg),
2713 V1((Standard_Real*)&G1(0), 0, MxDeg),
2714 V2((Standard_Real*)&G2(0), 0, MxDeg);
2715
2716 Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1,
2717 NTang3d, NTang2d,
2718 Point, TypOfConstr,
2719 p0 = Parameters.Lower() - myFirstPoint,
2720 curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim,
2721 jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d;
2722 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2723
2724 // Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints;
2725 // Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints;
2726 Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints;
2727 Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints;
2728 NBeg2d = Ng3d * myNbP3d;
2729 // NgPC1 = NbConstr + myNbCurvPoints;
2730 NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints;
2731 NPass = 0;
2732 NTang3d = 3 * NgPC1;
2733 NTang2d = 2 * NgPC1;
2734
2735 TColStd_Array1OfReal& Intervals = Curve->Knots();
2736
2737 Standard_Real t, R1, R2;
2738
2739 Handle(PLib_Base) myBase = Curve->Base();
2740 Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase);
2741 Standard_Integer Order = myHermitJacobi->NivConstr() + 1;
2742
2743 Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1;
2744
2745 A.NullifyConstraint();
2746
2747 ipnt = -1;
2748 ityp = 0;
2749 for(i = 1; i <= NbConstr; i++) {
2750
2751 ipnt += 2; ityp += 2;
2752
2753 Point = myTypConstraints->Value(ipnt);
2754 TypOfConstr = myTypConstraints->Value(ityp);
2755
2756 t = Parameters(p0 + Point);
2757
2758 for(el = curel; el <= NbElm; ) {
2759 if( t <= Intervals(++el) ) {
2760 curel = el - 1;
2761 break;
2762 }
2763 }
2764
2765
2766 UFirst = Intervals(curel); ULast = Intervals(curel + 1);
2767 coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.;
2768
2769 t = (t - c0) / coeff;
2770
2771 if(TypOfConstr == 0) {
2772 myBase->D0(t, G0);
2773 for(k = 1; k < Order; k++) {
2774 mfact = Pow(coeff, k);
2775 G0(k) *= mfact;
2776 G0(k + Order) *= mfact;
2777 }
2778 }
2779 else if(TypOfConstr == 1) {
2780 myBase->D1(t, G0, G1);
2781 for(k = 1; k < Order; k++) {
2782 mfact = Pow(coeff, k);
2783 G0(k) *= mfact;
2784 G0(k + Order) *= mfact;
2785 G1(k) *= mfact;
2786 G1(k + Order) *= mfact;
2787 }
2788 mfact = 1./coeff;
2789 for(k = 0; k <= MxDeg; k++) {
2790 G1(k) *= mfact;
2791 }
2792 }
2793 else {
2794 myBase->D2(t, G0, G1, G2);
2795 for(k = 1; k < Order; k++) {
2796 mfact = Pow(coeff, k);
2797 G0(k) *= mfact;
2798 G0(k + Order) *= mfact;
2799 G1(k) *= mfact;
2800 G1(k + Order) *= mfact;
2801 G2(k) *= mfact;
2802 G2(k + Order) *= mfact;
2803 }
2804 mfact = 1. / coeff;
2805 mfact1 = mfact / coeff;
2806 for(k = 0; k <= MxDeg; k++) {
2807 G1(k) *= mfact;
2808 G2(k) *= mfact1;
2809 }
2810 }
2811
2812 NPass++;
2813
2814 j = NbDim * (Point - myFirstPoint);
2815 Standard_Integer n0 = NPass;
2816 curdim = 0;
2817 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2818 IndexOfConstraint = n0;
2819 for(k = 1; k <= 3; k++) {
2820 curdim++;
2821 A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
2822 IndexOfConstraint += NgPC1;
2823 }
2824 j +=3;
2825 n0 += Ng3d;
2826 }
2827
2828 n0 = NPass + NBeg2d;
2829 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2830 IndexOfConstraint = n0;
2831 for(k = 1; k <= 2; k++) {
2832 curdim++;
2833 A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
2834 IndexOfConstraint += NgPC1;
2835 }
2836 j +=2;
2837 n0 += Ng2d;
2838 }
2839
2840 /* if(TypOfConstr == 1) {
2841
2842 IndexOfConstraint = NTang3d + 1;
2843 jt = Ntheta * (i - 1);
2844 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2845 for(k = 1; k <= 3; k++) {
2846 A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
2847 A.AddConstraint(IndexOfConstraint + 1, curel, k, myTtheta->Value(jt + 3 + k) * V1, 0.);
2848 }
2849 IndexOfConstraint += Ng3d;
2850 jt += 6;
2851 }
2852
2853 IndexOfConstraint = NBeg2d + NTang2d + 1;
2854 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2855 for(k = 1; k <= 2; k++) {
2856 A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
2857 }
2858 IndexOfConstraint += Ng2d;
2859 jt += 2;
2860 }
2861
2862 NTang3d += 2;
2863 NTang2d += 1;
2864 } */
2865 if(TypOfConstr == 1) {
2866
2867 NPass++;
2868 n0 = NPass;
2869 j = 2 * NbDim * (i - 1);
2870 curdim = 0;
2871 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2872 IndexOfConstraint = n0;
2873 for(k = 1; k <= 3; k++) {
2874 curdim++;
2875 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2876 IndexOfConstraint += NgPC1;
2877 }
2878 n0 += Ng3d;
2879 j += 6;
2880 }
2881
2882 n0 = NPass + NBeg2d;
2883 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2884 IndexOfConstraint = n0;
2885 for(k = 1; k <= 2; k++) {
2886 curdim++;
2887 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2888 IndexOfConstraint += NgPC1;
2889 }
2890 n0 += Ng2d;
2891 j += 4;
2892 }
2893 }
2894 if(TypOfConstr == 2) {
2895
2896 NPass++;
2897 n0 = NPass;
2898 j = 2 * NbDim * (i - 1);
2899 curdim = 0;
2900 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2901 IndexOfConstraint = n0;
2902 for(k = 1; k <= 3; k++) {
2903 curdim++;
2904 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2905 IndexOfConstraint += NgPC1;
2906 }
2907 n0 += Ng3d;
2908 j += 6;
2909 }
2910
2911 n0 = NPass + NBeg2d;
2912 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2913 IndexOfConstraint = n0;
2914 for(k = 1; k <= 2; k++) {
2915 curdim++;
2916 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2917 IndexOfConstraint += NgPC1;
2918 }
2919 n0 += Ng2d;
2920 j += 4;
2921 }
2922
2923 j = 2 * NbDim * (i - 1) + 3;
2924 jt = Ntheta * (i - 1);
2925 IndexOfConstraint = NTang3d + 1;
2926 curdim = 0;
2927 for(pnt = 1; pnt <= myNbP3d; pnt++) {
2928 R1 = 0.; R2 = 0.;
2929 for(k = 1; k <= 3; k++) {
2930 R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
2931 R2 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + 3 + k);
2932 }
2933 R1 *= CBLONG * CBLONG;
2934 R2 *= CBLONG * CBLONG;
2935 for(k = 1; k <= 3; k++) {
2936 curdim++;
2937 if(k > 1) R1 = R2 = 0.;
2938 A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
2939 A.AddConstraint(IndexOfConstraint + 1, curel, curdim, myTfthet->Value(jt + 3 + k) * V2, R2);
2940 }
2941 IndexOfConstraint += Ng3d;
2942 j += 6;
2943 jt += 6;
2944 }
2945
2946 j--;
2947 IndexOfConstraint = NBeg2d + NTang2d + 1;
2948 for(pnt = 1; pnt <= myNbP2d; pnt++) {
2949 R1 = 0.;
2950 for(k = 1; k <= 2; k++) {
2951 R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
2952 }
2953 R1 *= CBLONG * CBLONG;
2954 for(k = 1; k <= 2; k++) {
2955 curdim++;
2956 if(k > 1) R1 = 0.;
2957 A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
2958 }
2959 IndexOfConstraint += Ng2d;
2960 j += 4;
2961 jt += 2;
2962 }
2963
2964 NTang3d += 2;
2965 NTang2d += 1;
2966 }
2967
2968 }
2969
2970 }
2971
InitTthetaF(const Standard_Integer ndimen,const AppParCurves_Constraint typcon,const Standard_Integer begin,const Standard_Integer jndex)2972 Standard_Boolean AppDef_Variational::InitTthetaF(const Standard_Integer ndimen,
2973 const AppParCurves_Constraint typcon,
2974 const Standard_Integer begin,
2975 const Standard_Integer jndex)
2976 {
2977 if ((ndimen < 2)||(ndimen >3))
2978 return Standard_False;
2979 gp_Vec T, V;
2980 gp_Vec theta1, theta2;
2981 gp_Vec F;
2982 Standard_Real XX, XY, YY, XZ, YZ, ZZ;
2983
2984 if ((typcon == AppParCurves_TangencyPoint)||(typcon == AppParCurves_CurvaturePoint))
2985 {
2986 T.SetX(myTabConstraints->Value(jndex));
2987 T.SetY(myTabConstraints->Value(jndex + 1));
2988 if (ndimen == 3)
2989 T.SetZ(myTabConstraints->Value(jndex + 2));
2990 else
2991 T.SetZ(0.);
2992 if (ndimen == 2)
2993 {
2994 V.SetX(0.);
2995 V.SetY(0.);
2996 V.SetZ(1.);
2997 }
2998 if (ndimen == 3)
2999 if (!NotParallel(T, V))
3000 return Standard_False;
3001 theta1 = V ^ T;
3002 theta1.Normalize();
3003 myTtheta->SetValue(begin, theta1.X());
3004 myTtheta->SetValue(begin + 1, theta1.Y());
3005 if (ndimen == 3)
3006 {
3007 theta2 = T ^ theta1;
3008 theta2.Normalize();
3009 myTtheta->SetValue(begin + 2, theta1.Z());
3010 myTtheta->SetValue(begin + 3, theta2.X());
3011 myTtheta->SetValue(begin + 4, theta2.Y());
3012 myTtheta->SetValue(begin + 5, theta2.Z());
3013 }
3014
3015 // Calculation of myTfthet
3016 if (typcon == AppParCurves_CurvaturePoint)
3017 {
3018 XX = Pow(T.X(), 2);
3019 XY = T.X() * T.Y();
3020 YY = Pow(T.Y(), 2);
3021 if (ndimen == 2)
3022 {
3023 F.SetX(YY * theta1.X() - XY * theta1.Y());
3024 F.SetY(XX * theta1.Y() - XY * theta1.X());
3025 myTfthet->SetValue(begin, F.X());
3026 myTfthet->SetValue(begin + 1, F.Y());
3027 }
3028 if (ndimen == 3)
3029 {
3030 XZ = T.X() * T.Z();
3031 YZ = T.Y() * T.Z();
3032 ZZ = Pow(T.Z(), 2);
3033
3034 F.SetX((ZZ + YY) * theta1.X() - XY * theta1.Y() - XZ * theta1.Z());
3035 F.SetY((XX + ZZ) * theta1.Y() - XY * theta1.X() - YZ * theta1.Z());
3036 F.SetZ((XX + YY) * theta1.Z() - XZ * theta1.X() - YZ * theta1.Y());
3037 myTfthet->SetValue(begin, F.X());
3038 myTfthet->SetValue(begin + 1, F.Y());
3039 myTfthet->SetValue(begin + 2, F.Z());
3040 F.SetX((ZZ + YY) * theta2.X() - XY * theta2.Y() - XZ * theta2.Z());
3041 F.SetY((XX + ZZ) * theta2.Y() - XY * theta2.X() - YZ * theta2.Z());
3042 F.SetZ((XX + YY) * theta2.Z() - XZ * theta2.X() - YZ * theta2.Y());
3043 myTfthet->SetValue(begin + 3, F.X());
3044 myTfthet->SetValue(begin + 4, F.Y());
3045 myTfthet->SetValue(begin + 5, F.Z());
3046 }
3047 }
3048 }
3049 return Standard_True;
3050 }
3051
3052