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