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