projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Integration of OCCT 6.5.0 from SVN
[occt.git] / src / AppParCurves / AppParCurves_Variational_6.gxx
CommitLineData
7fd59977 1// File: AppParCurves_Variational_6.gxx
2// Created: Mon Dec 21 10:31:30 1998
3// Author: Igor FEOKTISTOV
4// <ifv@paradox.nnov.matra-dtv.fr>
5
6
7// Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559
8
9#include <PLib_Base.hxx>
10#include <PLib_JacobiPolynomial.hxx>
11#include <PLib_HermitJacobi.hxx>
12#include <FEmTool_Curve.hxx>
13#include <math_Vector.hxx>
14#include <TColStd_Array1OfReal.hxx>
15
16//
17//=======================================================================
18//function : InitSmoothCriterion
19//purpose : Initializes the SmoothCriterion
20//
21//=======================================================================
22//
23void AppParCurves_Variational::InitSmoothCriterion()
24{
25
26 const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9;
27// const Standard_Real J1 = .01, J2 = .001, J3 = .001;
28
29
30
31 Standard_Real Length;
32
33 InitParameters(Length);
34
35 mySmoothCriterion->SetParameters(myParameters);
36
37 Standard_Real E1, E2, E3;
38
39 InitCriterionEstimations(Length, E1, E2, E3);
40/*
41 J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6;
42
43 if(E1 < J1) E1 = J1;
44 if(E2 < J2) E2 = J2;
45 if(E3 < J3) E3 = J3;
46*/
47 mySmoothCriterion->EstLength() = Length;
48 mySmoothCriterion->SetEstimation(E1, E2, E3);
49
50 Standard_Real WQuadratic, WQuality;
51
52 if(!myWithMinMax && myTolerance != 0.)
53 WQuality = myTolerance;
54 else if (myTolerance == 0.)
55 WQuality = 1.;
56 else
57 WQuality = Max(myTolerance, Eps2 * Length);
58
59 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
60 WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality;
61 if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic;
62
63 if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.);
64
65
66 mySmoothCriterion->SetWeight(WQuadratic, WQuality,
67 myPercent[0], myPercent[1], myPercent[2]);
68
69
70 Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity);
71 Handle(FEmTool_Curve) TheCurve;
72 Standard_Integer NbElem;
73 Standard_Real CurvTol = Eps2 * Length / myNbPoints;
74
75 // Decoupe de l'intervalle en fonction des contraintes
76 if(myWithCutting == Standard_True && NbConstr != 0) {
77
78 InitCutting(TheBase, CurvTol, TheCurve);
79
80 }
81 else {
82
83 NbElem = 1;
84 TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol);
85 TheCurve->Knots().SetValue(TheCurve->Knots().Lower(),
86 myParameters->Value(myFirstPoint));
87 TheCurve->Knots().SetValue(TheCurve->Knots().Upper(),
88 myParameters->Value(myLastPoint));
89
90 }
91
92 mySmoothCriterion->SetCurve(TheCurve);
93
94 return;
95
96}
97
98//
99//=======================================================================
100//function : InitParameters
101//purpose : Calculation of initial estimation of the multicurve's length
102// and parameters for multipoints.
103//=======================================================================
104//
105void AppParCurves_Variational::InitParameters(Standard_Real& Length)
106{
107
108 const Standard_Real Eps1 = Precision::Confusion() * .01;
109
110 Standard_Real aux, dist;
111 Standard_Integer i, i0, i1 = 0, ipoint;
112
113
114 Length = 0.;
115 myParameters->SetValue(myFirstPoint, Length);
116
117 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) {
118 i0 = i1; i1 += myDimension;
119 dist = 0;
120 for(i = 1; i <= myDimension; i++) {
121 aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i);
122 dist += aux * aux;
123 }
124 Length += Sqrt(dist);
125 myParameters->SetValue(ipoint, Length);
126 }
127
128
129 if(Length <= Eps1)
130 Standard_ConstructionError::Raise("AppParCurves_Variational::InitParameters");
131
132
133 for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++)
134 myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length);
135
136 myParameters->SetValue(myLastPoint, 1.);
137
138 // Avec peu de point il y a sans doute sous estimation ...
139 if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1));
140}
141
142//
143//=======================================================================
144//function : InitCriterionEstimations
145//purpose : Calculation of initial estimation of the linear criteria
146//
147//=======================================================================
148//
149void AppParCurves_Variational::InitCriterionEstimations(const Standard_Real Length,
150 Standard_Real& E1,
151 Standard_Real& E2,
152 Standard_Real& E3) const
153{
154 E1 = Length * Length;
155
156 const Standard_Real Eps1 = Precision::Confusion() * .01;
157
158 math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension),
159 VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension);
160
161 // ========== Treatment of first point =================
162
163 Standard_Integer ipnt = myFirstPoint;
164
165 EstTangent(ipnt, VTang1);
166 ipnt++;
167 EstTangent(ipnt, VTang2);
168 ipnt++;
169 EstTangent(ipnt, VTang3);
170
171 ipnt = myFirstPoint;
172 EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1);
173 ipnt++;
174 EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2);
175
176// Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 Begin
177// Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt));
178 Standard_Integer anInd = ipnt;
179 Standard_Real Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt));
180// Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 End
181
182 if(Delta <= Eps1) Delta = 1.;
183
184 E2 = VScnd1.Norm2() * Delta;
185
186 E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
187 // ========== Treatment of internal points =================
188
189 Standard_Integer CurrPoint = 2;
190
191 for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) {
192
193 Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1));
194
195 if(CurrPoint == 1) {
196 if(ipnt + 1 != myLastPoint) {
197 EstTangent(ipnt + 2, VTang3);
198 EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2);
199 }
200 else
201 EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2);
202
203 E2 += VScnd1.Norm2() * Delta;
204 E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.;
205
206 }
207 else if(CurrPoint == 2) {
208 if(ipnt + 1 != myLastPoint) {
209 EstTangent(ipnt + 2, VTang1);
210 EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3);
211 }
212 else
213 EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3);
214
215 E2 += VScnd2.Norm2() * Delta;
216 E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
217
218 }
219 else {
220 if(ipnt + 1 != myLastPoint) {
221 EstTangent(ipnt + 2, VTang2);
222 EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1);
223 }
224 else
225 EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1);
226
227 E2 += VScnd3.Norm2() * Delta;
228 E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.;
229
230 }
231
232 CurrPoint++; if(CurrPoint == 4) CurrPoint = 1;
233 }
234
235 // ========== Treatment of last point =================
236
237 Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1));
238 if(Delta <= Eps1) Delta = 1.;
239
240 Standard_Real aux;
241
242 if(CurrPoint == 1) {
243
244 E2 += VScnd1.Norm2() * Delta;
245 aux = VScnd1.Subtracted(VScnd3).Norm2();
246 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
247
248 }
249 else if(CurrPoint == 2) {
250
251 E2 += VScnd2.Norm2() * Delta;
252 aux = VScnd2.Subtracted(VScnd1).Norm2();
253 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
254
255 }
256 else {
257
258 E2 += VScnd3.Norm2() * Delta;
259 aux = VScnd3.Subtracted(VScnd2).Norm2();
260 E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
261
262 }
263
264 aux = Length * Length;
265
266 E2 *= aux; E3 *= aux;
267
268
269}
270
271//
272//=======================================================================
273//function : EstTangent
274//purpose : Calculation of estimation of the Tangent
275// (see fortran routine MMLIPRI)
276//=======================================================================
277//
278
279void AppParCurves_Variational::EstTangent(const Standard_Integer ipnt,
280 math_Vector& VTang) const
281
282{
283 Standard_Integer i ;
284 const Standard_Real Eps1 = Precision::Confusion() * .01;
285 const Standard_Real EpsNorm = 1.e-9;
286
287 Standard_Real Wpnt = 1.;
288
289
290 if(ipnt == myFirstPoint) {
291 // Estimation at first point
292 if(myNbPoints < 3)
293 Wpnt = 0.;
294 else {
295
296 Standard_Integer adr1 = 1, adr2 = adr1 + myDimension,
297 adr3 = adr2 + myDimension;
298
299 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
300 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
301 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
302
303 // Parabolic interpolation
304 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
305 // first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d))
306 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d
307 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
308 Standard_Real V2 = 0.;
309 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
310 if(V2 > Eps1) {
311 Standard_Real d = V1 / (V1 + V2), d1;
312 d1 = 1. / (d * (1 - d)); d *= d;
313 VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1;
314 }
315 else {
316 // Simple 2-point estimation
317
318 VTang = Pnt2 - Pnt1;
319 }
320 }
321 }
322 else if (ipnt == myLastPoint) {
323 // Estimation at last point
324 if(myNbPoints < 3)
325 Wpnt = 0.;
326 else {
327
328 Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1,
329 adr2 = adr1 + myDimension,
330 adr3 = adr2 + myDimension;
331
332 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
333 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
334 math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
335
336 // Parabolic interpolation
337 // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
338 // first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d))
339 // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2)
340 Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
341 Standard_Real V2 = 0.;
342 if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
343 if(V2 > Eps1) {
344 Standard_Real d = V1 / (V1 + V2), d1;
345 d1 = 1. / (d * (1 - d)); d *= d - 2;
346 VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1;
347 }
348 else {
349 // Simple 2-point estimation
350
351 VTang = Pnt3 - Pnt2;
352 }
353 }
354 }
355 else {
356
357 Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1,
358 adr2 = adr1 + 2 * myDimension;
359
360 math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
361 math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
362
363 VTang = Pnt2 - Pnt1;
364
365 }
366
367 Standard_Real Vnorm = VTang.Norm();
368
369 if(Vnorm <= EpsNorm)
370 VTang.Init(0.);
371 else
372 VTang /= Vnorm;
373
374 // Estimation with constraints
375
376 Standard_Real Wcnt = 0.;
377 Standard_Integer IdCnt = 1;
378
379// Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
380
381 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
382 math_Vector VCnt(1, myDimension, 0.);
383
384 if(NbConstr > 0) {
385
386 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
387 IdCnt <= NbConstr) IdCnt++;
388 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
389 (myTypConstraints->Value(2 * IdCnt) >= 1)) {
390 Wcnt = 1.;
391 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0;
392 for( i = 1; i <= myNbP3d; i++) {
393 for(Standard_Integer j = 1; j <= 3; j++)
394 VCnt(++k) = myTabConstraints->Value(++i0);
395 i0 += 3;
396 }
397 for(i = 1; i <= myNbP2d; i++) {
398 for(Standard_Integer j = 1; j <= 2; j++)
399 VCnt(++k) = myTabConstraints->Value(++i0);
400 i0 += 2;
401 }
402 }
403 }
404
405 // Averaging of estimation
406
407 Standard_Real Denom = Wpnt + Wcnt;
408 if(Denom == 0.) Denom = 1.;
409 else Denom = 1. / Denom;
410
411 VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom;
412
413 Vnorm = VTang.Norm();
414
415 if(Vnorm <= EpsNorm)
416 VTang.Init(0.);
417 else
418 VTang /= Vnorm;
419
420
421}
422
423//
424//=======================================================================
425//function : EstSecnd
426//purpose : Calculation of estimation of the second derivative
427// (see fortran routine MLIMSCN)
428//=======================================================================
429//
430void AppParCurves_Variational::EstSecnd(const Standard_Integer ipnt,
431 const math_Vector& VTang1,
432 const math_Vector& VTang2,
433 const Standard_Real Length,
434 math_Vector& VScnd) const
435{
436 Standard_Integer i ;
437
438 const Standard_Real Eps = 1.e-9;
439
440 Standard_Real Wpnt = 1.;
441
442 Standard_Real aux;
443
444 if(ipnt == myFirstPoint)
445 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt);
446 else if(ipnt == myLastPoint)
447 aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1);
448 else
449 aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1);
450
451 if(aux <= Eps)
452 aux = 1.;
453 else
454 aux = 1. / aux;
455
456 VScnd = (VTang2 - VTang1) * aux;
457
458 // Estimation with constraints
459
460 Standard_Real Wcnt = 0.;
461 Standard_Integer IdCnt = 1;
462
463// Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
464
465 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
466 math_Vector VCnt(1, myDimension, 0.);
467
468 if(NbConstr > 0) {
469
470 while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
471 IdCnt <= NbConstr) IdCnt++;
472
473 if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
474 (myTypConstraints->Value(2 * IdCnt) >= 2)) {
475 Wcnt = 1.;
476 Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0;
477 for( i = 1; i <= myNbP3d; i++) {
478 for(Standard_Integer j = 1; j <= 3; j++)
479 VCnt(++k) = myTabConstraints->Value(++i0);
480 i0 += 3;
481 }
482 i0--;
483 for(i = 1; i <= myNbP2d; i++) {
484 for(Standard_Integer j = 1; j <= 2; j++)
485 VCnt(++k) = myTabConstraints->Value(++i0);
486 i0 += 2;
487 }
488 }
489 }
490
491 // Averaging of estimation
492
493 Standard_Real Denom = Wpnt + Wcnt;
494 if(Denom == 0.) Denom = 1.;
495 else Denom = 1. / Denom;
496
497 VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom;
498
499}
500
501
502//
503//=======================================================================
504//function : InitCutting
505//purpose : Realisation of curve's cutting a priory accordingly to
506// constraints (see fortran routine MLICUT)
507//=======================================================================
508//
509void AppParCurves_Variational::InitCutting(const Handle(PLib_Base)& aBase,
510 const Standard_Real CurvTol,
511 Handle(FEmTool_Curve)& aCurve) const
512{
513
514 // Definition of number of elements
515 Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem;
516 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
517
518 for(i = 1; i <= NbConstr; i++) {
519 kk = Abs(myTypConstraints->Value(2 * i)) + 1;
520 ORCMx = Max(ORCMx, kk);
521 NCont += kk;
522 }
523
524 if(ORCMx > myMaxDegree - myNivCont)
525 Standard_ConstructionError::Raise("AppParCurves_Variational::InitCutting");
526
527 Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4,
528 myNivCont + 1);
529
530 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
531
532 while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) {
533
534 NLibre++;
535 NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
536
537 }
538
539
540 if(NbElem > myMaxSegment)
541 Standard_ConstructionError::Raise("AppParCurves_Variational::InitCutting");
542
543
544 aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol);
545
546 Standard_Integer NCnt = (NCont - 1) / NbElem + 1;
547 Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont);
548
549 TColStd_Array1OfReal& Knot = aCurve->Knots();
550
551 Standard_Integer IDeb = 0, IFin = NbConstr + 1,
552 NDeb = 0, NFin = 0,
553 IndEl = Knot.Lower(), IUpper = Knot.Upper(),
554 NbEl = 0;
555
556
557 Knot(IndEl) = myParameters->Value(myFirstPoint);
558 Knot(IUpper) = myParameters->Value(myLastPoint);
559
560 while(NbElem - NbEl > 1) {
561
562 IndEl++; NbEl++;
563 if(NPlus == 0) NCnt--;
564
565 while(NDeb < NCnt && IDeb < IFin) {
566 IDeb++;
567 NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1;
568 }
569
570 if(NDeb == NCnt) {
571 NDeb = 0;
572 if(NPlus == 1 &&
573 myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1))
574
575 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
576 else
577 Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) +
578 myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2;
579
580 }
581 else {
582 NDeb -= NCnt;
583 Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
584 }
585
586 NPlus--;
587 if(NPlus == 0) NCnt--;
588
589 if(NbElem - NbEl == 1) break;
590
591 NbEl++;
592
593 while(NFin < NCnt && IDeb < IFin) {
594 IFin--;
595 NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1;
596 }
597
598 if(NFin == NCnt) {
599 NFin = 0;
600 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
601 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
602 }
603 else {
604 NFin -= NCnt;
605 if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1))
606 Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1));
607 else
608 Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
609 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
610 }
611
612 }
613
614
615}
Open CASCADE Technology repository