0030597: BRepOffsetAPI_MakePipeShell raises Standard_RangeError
[occt.git] / src / GeomFill / GeomFill_SectionPlacement.cxx
CommitLineData
b311480e 1// Created on: 1997-12-15
2// Created by: Philippe MANGIN
3// Copyright (c) 1997-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//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 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
7fd59977 17
42cf5bc1 18#include <Adaptor3d_HCurve.hxx>
7fd59977 19#include <Bnd_Box.hxx>
20#include <BndLib_Add3dCurve.hxx>
7fd59977 21#include <Extrema_ExtCC.hxx>
22#include <Extrema_POnCurv.hxx>
7fd59977 23#include <Geom_BezierCurve.hxx>
42cf5bc1 24#include <Geom_BSplineCurve.hxx>
7fd59977 25#include <Geom_CartesianPoint.hxx>
42cf5bc1 26#include <Geom_Conic.hxx>
ec357c5c 27#include <Geom_Curve.hxx>
42cf5bc1 28#include <Geom_Geometry.hxx>
29#include <Geom_Line.hxx>
30#include <Geom_Plane.hxx>
31#include <Geom_TrimmedCurve.hxx>
32#include <GeomAbs_CurveType.hxx>
33#include <GeomAdaptor_HSurface.hxx>
34#include <GeomFill_LocationLaw.hxx>
35#include <GeomFill_SectionPlacement.hxx>
36#include <GeomLib.hxx>
37#include <GeomLProp_CLProps.hxx>
38#include <gp.hxx>
39#include <gp_Ax2.hxx>
40#include <gp_Ax3.hxx>
41#include <gp_Dir.hxx>
42#include <gp_Mat.hxx>
43#include <gp_Pnt.hxx>
44#include <gp_Trsf.hxx>
45#include <gp_Vec.hxx>
46#include <IntCurveSurface_HInter.hxx>
47#include <IntCurveSurface_IntersectionPoint.hxx>
48#include <Precision.hxx>
49#include <StdFail_NotDone.hxx>
50#include <TColgp_HArray1OfPnt.hxx>
7fd59977 51
52//===============================================================
53// Function :
54// Purpose :
55//===============================================================
56static void Tangente(const Adaptor3d_Curve& Path,
57 const Standard_Real Param,
58 gp_Pnt& P,
59 gp_Vec& Tang)
60{
61 Path.D1(Param, P, Tang);
62 Standard_Real Norm = Tang.Magnitude();
63
64 for (Standard_Integer ii=2; (ii<12) && (Norm < Precision::Confusion());
65 ii++) {
66 Tang = Path.DN(Param, ii);
67 Norm = Tang.Magnitude();
68 }
69
70 if (Norm > 100*gp::Resolution()) Tang /= Norm;
71}
72
73
74static Standard_Real Penalite(const Standard_Real angle,
75 const Standard_Real dist)
76{
77 Standard_Real penal;
78
79 if (dist < 1)
80 penal = Sqrt(dist);
81 else if (dist <2)
82 penal = Pow(dist, 2);
83 else
84 penal = dist + 2;
85
86 if (angle > 1.e-3) {
c6541a0c 87 penal += 1./angle -2./M_PI;
7fd59977 88 }
89 else {
90 penal += 1.e3;
91 }
92
93 return penal;
94}
95
96static Standard_Real EvalAngle(const gp_Vec& V1,
97 const gp_Vec& V2)
98{
99 Standard_Real angle;
100 angle = V1.Angle(V2);
c6541a0c 101 if (angle > M_PI/2) angle = M_PI - angle;
7fd59977 102 return angle;
103}
104
105static Standard_Integer NbSamples(const Handle(Geom_Curve)& aCurve)
106{
107 Standard_Real nbs = 100.; //on default
108
109 Handle(Geom_Curve) theCurve = aCurve;
110 if (aCurve->IsInstance(STANDARD_TYPE(Geom_TrimmedCurve)))
111 theCurve = (Handle(Geom_TrimmedCurve)::DownCast(aCurve))->BasisCurve();
112
113 if (theCurve->IsInstance(STANDARD_TYPE(Geom_Line)))
114 nbs = 1;
115 else if (theCurve->IsKind(STANDARD_TYPE(Geom_Conic)))
116 nbs = 4;
117 else if (theCurve->IsInstance(STANDARD_TYPE(Geom_BezierCurve)))
118 {
119 Handle(Geom_BezierCurve) BC = Handle(Geom_BezierCurve)::DownCast(theCurve);
120 nbs = 3 + BC->NbPoles();
121 }
122 else if (theCurve->IsInstance(STANDARD_TYPE(Geom_BSplineCurve)))
123 {
124 Handle(Geom_BSplineCurve) BC = Handle(Geom_BSplineCurve)::DownCast(theCurve);
125 nbs = BC->NbKnots();
126 nbs *= BC->Degree();
127 Standard_Real ratio =
128 (aCurve->LastParameter() - aCurve->FirstParameter())/(BC->LastParameter() - BC->FirstParameter());
129 nbs *= ratio;
130 if(nbs < 4.0)
131 nbs = 4;
132 }
133
134 if (nbs > 300.)
135 nbs = 300;
136
137 return ((Standard_Integer)nbs);
138}
139
140//===============================================================
141// Function :DistMini
142// Purpose : Examine un extrema pour updater <Dist> & <Param>
143//===============================================================
144static void DistMini(const Extrema_ExtPC& Ext,
145 const Adaptor3d_Curve& C,
146 Standard_Real& Dist,
147 Standard_Real& Param)
148{
149 Standard_Real dist1, dist2;
150 Standard_Integer ii;
151 gp_Pnt P1, P2;
152 Standard_Real Dist2 = RealLast();
153
154 Ext.TrimmedSquareDistances(dist1, dist2, P1, P2);
155 if ( (dist1<Dist2) || (dist2<Dist2) ) {
156 if (dist1 < dist2) {
157 Dist2 = dist1;
158 Param = C.FirstParameter();
159 }
160 else {
161 Dist2 = dist2;
162 Param = C.LastParameter();
163 }
164 }
165
166 if (Ext.IsDone())
167 {
168 for (ii=1; ii<= Ext.NbExt(); ii++) {
169 if (Ext.SquareDistance(ii) < Dist2) {
170 Dist2 = Ext.SquareDistance(ii);
171 Param = Ext.Point(ii).Parameter();
172 }
173 }
174 }
175 Dist = sqrt (Dist2);
176}
177
178
179//===============================================================
180// Function : Constructeur
181// Purpose :
182//===============================================================
183GeomFill_SectionPlacement::
184GeomFill_SectionPlacement(const Handle(GeomFill_LocationLaw)& L,
185 const Handle(Geom_Geometry)& Section) :
186 myLaw(L), /* myAdpSection(Section), mySection(Section), */
187 Dist(RealLast()), AngleMax(0.)
188{
189
190 done = Standard_False;
191 isplan = Standard_False;
192 myIsPoint = Standard_False;
193
194 if (Section->IsInstance(STANDARD_TYPE(Geom_CartesianPoint)))
195 {
196 myIsPoint = Standard_True;
197 myPoint = Handle(Geom_CartesianPoint)::DownCast(Section)->Pnt();
198 isplan = Standard_True;
199 }
200 else
201 {
202 Handle(Geom_Curve) CurveSection = Handle(Geom_Curve)::DownCast(Section);
203 myAdpSection.Load(CurveSection);
204 mySection = CurveSection;
205 }
206
207 Standard_Integer i, j, NbPoles=0;
208 Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
209
210 // Boite d'encombrement de la section pour en deduire le gabarit
211 Bnd_Box box;
212 if (myIsPoint)
213 box.Add(myPoint);
214 else
215 BndLib_Add3dCurve::Add(myAdpSection, 1.e-4, box);
216 box.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
217
218 Standard_Real DX = aXmax - aXmin ;
219 Standard_Real DY = aYmax - aYmin ;
220 Standard_Real DZ = aZmax - aZmin ;
221 Gabarit = Sqrt( DX*DX+ DY*DY + DZ*DZ )/2. ;
222
223 Gabarit += Precision::Confusion(); // Cas des toute petite
224
225
226 // Initialisation de TheAxe pour les cas singulier
227 if (!myIsPoint)
228 {
229 gp_Pnt P;
230 gp_Vec V;
231 Tangente(myAdpSection,
232 (myAdpSection.FirstParameter()+myAdpSection.LastParameter())/2,
233 P, V);
234 TheAxe.SetLocation(P);
235 TheAxe.SetDirection(V);
236
237 // y a t'il un Plan moyen ?
238 GeomAbs_CurveType TheType = myAdpSection.GetType();
239 switch (TheType) {
240 case GeomAbs_Circle:
241 {
242 isplan = Standard_True;
243 TheAxe = myAdpSection.Circle().Axis();
244 break;
245 }
246 case GeomAbs_Ellipse:
247 {
248 isplan = Standard_True;
249 TheAxe = myAdpSection.Ellipse().Axis();
250 break;
251 }
252 case GeomAbs_Hyperbola:
253 {
254 isplan = Standard_True;
255 TheAxe = myAdpSection.Hyperbola().Axis();
256 break;
257 }
258 case GeomAbs_Parabola:
259 {
260 isplan = Standard_True;
261 TheAxe = myAdpSection.Parabola().Axis();
262 break;
263 }
264 case GeomAbs_Line:
265 {
266 NbPoles = 0; // Pas de Plan !!
267 break;
268 }
269 case GeomAbs_BezierCurve:
270 case GeomAbs_BSplineCurve:
271 {
272 NbPoles = myAdpSection.NbPoles();
273 break;
274 }
275 default:
276 NbPoles = 21;
277 }
278
279
280 if (!isplan && NbPoles>2)
281 {
282 // Calcul d'un plan moyen.
283 Handle(TColgp_HArray1OfPnt) Pnts;
284 Standard_Real first = myAdpSection.FirstParameter();
285 Standard_Real last = myAdpSection.LastParameter();
a5278fc1 286 if (myAdpSection.IsPeriodic())
287 {
288 //Correct boundaries to avoid mistake of LocateU
289 Handle(Geom_Curve) aCurve = myAdpSection.Curve();
290 if (aCurve->IsInstance(STANDARD_TYPE(Geom_TrimmedCurve)))
291 aCurve = (Handle(Geom_TrimmedCurve)::DownCast(aCurve))->BasisCurve();
292 Standard_Real Ufirst = aCurve->FirstParameter();
293 Standard_Real aPeriod = aCurve->Period();
294 Standard_Real U1 = Ufirst + Floor((first - Ufirst)/aPeriod) * aPeriod;
295 Standard_Real U2 = U1 + aPeriod;
296 if (Abs(first - U1) <= Precision::PConfusion())
297 first = U1;
298 if (Abs(last - U2) <= Precision::PConfusion())
299 last = U2;
300 }
7fd59977 301 Standard_Real t, delta;
302 if (myAdpSection.GetType() == GeomAbs_BSplineCurve)
303 {
304 Handle(Geom_BSplineCurve) BC =
305 Handle(Geom_BSplineCurve)::DownCast(myAdpSection.Curve());
306 Standard_Integer I1, I2, I3, I4;
307 BC->LocateU( first, Precision::Confusion(), I1, I2 );
308 BC->LocateU( last, Precision::Confusion(), I3, I4 );
309 Standard_Integer NbKnots = I3 - I2 + 1;
310
311 Standard_Integer NbLocalPnts = 10;
312 Standard_Integer NbPnts = (NbKnots-1) * NbLocalPnts;
313 if (I1 != I2)
314 NbPnts += NbLocalPnts;
315 if (I3 != I4 && first < BC->Knot(I3))
316 NbPnts += NbLocalPnts;
317 if (!myAdpSection.IsClosed())
318 NbPnts++;
319 Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
320 Standard_Integer nb = 1;
321 if (I1 != I2)
322 {
323 Standard_Real locallast = (BC->Knot(I2) < last)? BC->Knot(I2) : last;
324 delta = (locallast - first) / NbLocalPnts;
325 for (j = 0; j < NbLocalPnts; j++)
326 {
327 t = first + j*delta;
328 Pnts->SetValue( nb++, myAdpSection.Value(t) );
329 }
330 }
331 for (i = I2; i < I3; i++)
332 {
333 t = BC->Knot(i);
334 delta = (BC->Knot(i+1) - t) / NbLocalPnts;
335 for (j = 0; j < NbLocalPnts; j++)
336 {
7fd59977 337 Pnts->SetValue( nb++, myAdpSection.Value(t) );
c8ea5b8e 338 t += delta;
7fd59977 339 }
340 }
341 if (I3 != I4 && first < BC->Knot(I3))
342 {
343 t = BC->Knot(I3);
344 delta = (last - t) / NbLocalPnts;
345 for (j = 0; j < NbLocalPnts; j++)
346 {
7fd59977 347 Pnts->SetValue( nb++, myAdpSection.Value(t) );
c8ea5b8e 348 t += delta;
7fd59977 349 }
350 }
351 if (!myAdpSection.IsClosed())
352 Pnts->SetValue( nb, myAdpSection.Value(last) );
353 }
354 else // other type
355 {
356 Standard_Integer NbPnts = NbPoles-1;
357 if (!myAdpSection.IsClosed())
358 NbPnts++;
359 Pnts = new TColgp_HArray1OfPnt(1, NbPnts);
360 delta = (last - first) / (NbPoles-1);
361 for (i = 0; i < NbPoles-1; i++)
362 {
363 t = first + i*delta;
364 Pnts->SetValue( i+1, myAdpSection.Value(t) );
365 }
366 if (!myAdpSection.IsClosed())
367 Pnts->SetValue( NbPnts, myAdpSection.Value(last) );
368 }
369
370 Standard_Boolean issing;
371 gp_Ax2 axe;
372 GeomLib::AxeOfInertia(Pnts->Array1(), axe, issing, Precision::Confusion());
373 if (!issing) {
374 isplan = Standard_True;
375 TheAxe.SetLocation(axe.Location());
376 TheAxe.SetDirection(axe.Direction());
377 }
378 }
379
380
381 myExt.Initialize(myAdpSection,
382 myAdpSection.FirstParameter(),
383 myAdpSection.LastParameter(),
384 Precision::Confusion());
385 }
386}
387
388//===============================================================
389// Function :SetLocation
390// Purpose :
391//===============================================================
392void GeomFill_SectionPlacement::
393SetLocation(const Handle(GeomFill_LocationLaw)& L)
394{
395 myLaw = L;
396}
397
398//===============================================================
399// Function : Perform
400// Purpose : Le plus simple
401//===============================================================
402void GeomFill_SectionPlacement::Perform(const Standard_Real Tol)
403{
404 Handle(Adaptor3d_HCurve) Path;
405 Path = myLaw->GetCurve();
406 Perform(Path, Tol);
407}
408
409//===============================================================
410// Function :Perform
411// Purpose : Recherche automatique
412//===============================================================
413void GeomFill_SectionPlacement::Perform(const Handle(Adaptor3d_HCurve)& Path,
414 const Standard_Real Tol)
415{
416 Standard_Real IntTol = 1.e-5;
417 Standard_Real DistCenter = Precision::Infinite();
418
419 if (myIsPoint)
420 {
421 Extrema_ExtPC Projector(myPoint, Path->Curve(), Precision::Confusion());
422 DistMini( Projector, Path->Curve(), Dist, PathParam );
c6541a0c 423 AngleMax = M_PI/2;
7fd59977 424 }
425 else
426 {
427 PathParam = Path->FirstParameter();
428 SecParam = myAdpSection.FirstParameter();
429
576e3066 430 Standard_Real distaux, taux = 0.0, alpha;
7fd59977 431 gp_Pnt PonPath, PonSec, P;
432 gp_Vec VRef, dp1;
433 VRef.SetXYZ(TheAxe.Direction().XYZ());
434
435 Tangente( Path->Curve(), PathParam, PonPath, dp1);
436 PonSec = myAdpSection.Value(SecParam);
437 Dist = PonPath.Distance(PonSec);
438 if (Dist > Tol) { // On Cherche un meilleur point sur la section
439 myExt.Perform(PonPath);
440 if ( myExt.IsDone() ) {
441 DistMini(myExt, myAdpSection, Dist, SecParam);
442 PonSec = myAdpSection.Value(SecParam);
443 }
444 }
445 AngleMax = EvalAngle(VRef, dp1);
c6541a0c 446 if (isplan) AngleMax = M_PI/2 - AngleMax;
7fd59977 447
448 Standard_Boolean Trouve = Standard_False;
449 Standard_Integer ii;
450
451 if (isplan) {
452 // (1.1) Distances Point-Plan
453 Standard_Real DistPlan;
454 gp_Vec V1 (PonPath, TheAxe.Location());
455 DistPlan = Abs(V1.Dot(VRef));
456 if (DistPlan <= IntTol)
457 DistCenter = V1.Magnitude();
458
459 gp_Pnt Plast = Path->Value( Path->LastParameter() );
460 V1.SetXYZ( TheAxe.Location().XYZ() - Plast.XYZ() );
461 DistPlan = Abs(V1.Dot(VRef));
462 if (DistPlan <= IntTol)
463 {
464 Standard_Real aDist = V1.Magnitude();
465 if (aDist < DistCenter)
466 {
467 DistCenter = aDist;
468 PonPath = Plast;
469 PathParam = Path->LastParameter();
470 }
471 }
472
473 // (1.2) Intersection Plan-courbe
474 gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
475 Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
476 Handle(GeomAdaptor_HSurface) adplan =
477 new (GeomAdaptor_HSurface)(plan);
478 IntCurveSurface_HInter Intersector;
479 Intersector.Perform(Path, adplan);
480 if (Intersector.IsDone())
481 {
482 Standard_Real w;
7fd59977 483 Standard_Real aDist;
484 for (ii=1; ii<=Intersector.NbPoints(); ii++)
485 {
486 w = Intersector.Point(ii).W();
51740958 487 P = Path->Value( w );
7fd59977 488 aDist = P.Distance( TheAxe.Location() );
489 if (aDist < DistCenter)
490 {
491 DistCenter = aDist;
492 PonPath = P;
493 PathParam = w;
494 }
495 }
496 }
497 if (!Intersector.IsDone() || Intersector.NbPoints() == 0)
498 {
499 Standard_Integer NbPnts = NbSamples( mySection );
500 TColgp_Array1OfPnt Pnts( 1, NbPnts+1 );
501 Standard_Real delta = (mySection->LastParameter()-mySection->FirstParameter())/NbPnts;
502 for (ii = 0; ii <= NbPnts; ii++)
503 Pnts(ii+1) = mySection->Value( mySection->FirstParameter() + ii*delta );
504
505 gp_Pnt BaryCenter;
506 gp_Dir Xdir, Ydir;
507 Standard_Real Xgap, Ygap, Zgap;
508 GeomLib::Inertia( Pnts, BaryCenter, Xdir, Ydir, Xgap, Ygap, Zgap );
509
510 gp_Pnt Pfirst = Path->Value( Path->FirstParameter() );
511 if (Pfirst.Distance(BaryCenter) < Plast.Distance(BaryCenter))
512 PathParam = Path->FirstParameter();
513 else
514 {
515 PathParam = Path->LastParameter();
516 Tangente( Path->Curve(), PathParam, PonPath, dp1);
517 PonSec = myAdpSection.Value(SecParam);
518 Dist = PonPath.Distance(PonSec);
519 if (Dist > Tol) { // On Cherche un meilleur point sur la section
520 myExt.Perform(PonPath);
521 if ( myExt.IsDone() ) {
522 DistMini(myExt, myAdpSection, Dist, SecParam);
523 PonSec = myAdpSection.Value(SecParam);
524 }
525 }
526 AngleMax = EvalAngle(VRef, dp1);
c6541a0c 527 AngleMax = M_PI/2 - AngleMax;
7fd59977 528 }
529 }
530
531 /*
532 // (1.1) Distances Point-Plan
533 Standard_Real DistPlan;
534 gp_Vec V1 (PonPath, TheAxe.Location());
535 DistPlan = Abs(V1.Dot(VRef));
536
537 // On examine l'autre extremite
538 gp_Pnt P;
539 Tangente(Path->Curve(), Path->LastParameter(), P, dp1);
540 V1.SetXYZ(TheAxe.Location().XYZ()-P.XYZ());
541 if (Abs(V1.Dot(VRef)) <= DistPlan ) { // On prend l'autre extremite
c6541a0c 542 alpha = M_PI/2 - EvalAngle(VRef, dp1);
7fd59977 543 distaux = PonPath.Distance(PonSec);
544 if (distaux > Tol) {
545 myExt.Perform(P);
546 if ( myExt.IsDone() )
547 DistMini(myExt, myAdpSection, distaux, taux);
548 }
549 else
550 taux = SecParam;
551
552 if (Choix(distaux, alpha)) {
553 Dist = distaux;
554 AngleMax = alpha;
555 PonPath = P;
556 PathParam = Path->LastParameter();
557 }
558 }
559
560 // (1.2) Intersection Plan-courbe
561 gp_Ax3 axe (TheAxe.Location(), TheAxe.Direction());
562 Handle(Geom_Plane) plan = new (Geom_Plane)(axe);
563 Handle(GeomAdaptor_HSurface) adplan =
564 new (GeomAdaptor_HSurface)(plan);
565 IntCurveSurface_HInter Intersector;
566 Intersector.Perform(Path, adplan);
567 if (Intersector.IsDone()) {
568 Standard_Real w;
569 gp_Vec V;
570 for (ii=1; ii<=Intersector.NbPoints(); ii++) {
571 w = Intersector.Point(ii).W();
572 //(1.3) test d'angle
573 Tangente( Path->Curve(), w, P, V);
c6541a0c 574 alpha = M_PI/2 - EvalAngle(V, VRef);
7fd59977 575 //(1.4) Test de distance Point-Courbe
576 myExt.Perform(P);
577 if ( myExt.IsDone() ) {
578 DistMini(myExt, myAdpSection, distaux, taux);
579 if (Choix(distaux, alpha)) {
580 Dist = distaux;
581 SecParam = taux;
582 AngleMax = alpha;
583 PonPath = P;
584 PathParam = w;
585 PonSec = myAdpSection.Value(SecParam);
586 }
587 }
588 else {
589 distaux = P.Distance(PonSec);
590 if (Choix(distaux, alpha)) {
591 Dist = distaux;
592 AngleMax = alpha;
593 PonPath = P;
594 PathParam = w;
595 }
596 }
597 }
598 }
599 */
0797d9d3 600#ifdef OCCT_DEBUG
7fd59977 601 if (Intersector.NbPoints() == 0) {
602 Intersector.Dump();
603 }
604#endif
605 }
606
607 // Cas General
608 if (! isplan) {
609 // (2.1) Distance avec les extremites ...
610 myExt.Perform(PonPath);
611 if ( myExt.IsDone() ) {
612 DistMini(myExt, myAdpSection, distaux, taux);
613 if (distaux < Dist) {
614 Dist = distaux;
615 SecParam = taux;
616 }
617 }
618 Trouve = (Dist <= Tol);
619 if (!Trouve) {
620 Tangente( Path->Curve(), Path->LastParameter(), P, dp1);
621 alpha = EvalAngle(VRef, dp1);
622 myExt.Perform(P);
623 if ( myExt.IsDone() ) {
624 if ( myExt.IsDone() ) {
625 DistMini(myExt, myAdpSection, distaux, taux);
626 if (Choix(distaux, alpha)) {
627 Dist = distaux;
628 SecParam = taux;
629 AngleMax = alpha;
630 PonPath = P;
631 PathParam = Path->LastParameter();
632 }
633 }
634 }
635 Trouve = (Dist <= Tol);
636 }
637
638 // (2.2) Distance courbe-courbe
639 if (!Trouve) {
640 Extrema_ExtCC Ext(Path->Curve(), myAdpSection,
641 Path->FirstParameter(), Path->LastParameter(),
642 myAdpSection.FirstParameter(),
643 myAdpSection.LastParameter(),
644 Path->Resolution(Tol/100),
645 myAdpSection.Resolution(Tol/100));
646 if (Ext.IsDone()) {
647 Extrema_POnCurv P1, P2;
648 for (ii=1; ii<=Ext.NbExt(); ii++) {
649 distaux = sqrt (Ext.SquareDistance(ii));
650 Ext.Points(ii, P1, P2);
651 Tangente(Path->Curve(), P1.Parameter(), P, dp1);
652 alpha = EvalAngle(VRef, dp1);
653 if (Choix(distaux, alpha)) {
654 Trouve = Standard_True;
655 Dist = distaux;
656 PathParam = P1.Parameter();
657 SecParam = P2.Parameter();
658 PonSec = P2.Value();
659 PonPath = P;
660 AngleMax = alpha;
661 }
662 }
663 }
664 if (!Trouve){
665 // Si l'on a toujours rien, on essai une distance point/path
666 // c'est la derniere chance.
667 Extrema_ExtPC PExt;
668 PExt.Initialize(Path->Curve(),
669 Path->FirstParameter(),
670 Path->LastParameter(),
671 Precision::Confusion());
672 PExt.Perform(PonSec);
673 if ( PExt.IsDone() ) {
674 // modified for OCC13595 Tue Oct 17 15:00:08 2006.BEGIN
675 // DistMini(PExt, myAdpSection, distaux, taux);
676 DistMini(PExt, Path->Curve(), distaux, taux);
677 // modified for OCC13595 Tue Oct 17 15:00:11 2006.END
678 Tangente(Path->Curve(), taux, P, dp1);
679 alpha = EvalAngle(VRef, dp1);
680 if (Choix(distaux, alpha)) {
681 Dist = distaux;
682 PonPath = P;
683 AngleMax = alpha;
684 PathParam = taux;
685 }
686 }
687 }
688 }
689 }
690 }
691
692 done = Standard_True;
693}
694
695
696//===============================================================
697// Function : Perform
698// Purpose : Calcul le placement pour un parametre donne.
699//===============================================================
700void GeomFill_SectionPlacement::Perform(const Standard_Real Param,
701 const Standard_Real Tol)
702{
703 done = Standard_True;
704 Handle(Adaptor3d_HCurve) Path;
705 Path = myLaw->GetCurve();
706
707 PathParam = Param;
708 if (myIsPoint)
709 {
710 gp_Pnt PonPath = Path->Value( PathParam );
711 Dist = PonPath.Distance(myPoint);
c6541a0c 712 AngleMax = M_PI/2;
7fd59977 713 }
714 else
715 {
716 SecParam = myAdpSection.FirstParameter();
717
718 // Standard_Real distaux, taux, alpha;
719 // gp_Pnt PonPath, PonSec, P;
720 gp_Pnt PonPath, PonSec;
721 gp_Vec VRef, dp1;
722 VRef.SetXYZ(TheAxe.Direction().XYZ());
723
724 Tangente( Path->Curve(), PathParam, PonPath, dp1);
725 PonSec = myAdpSection.Value(SecParam);
726 Dist = PonPath.Distance(PonSec);
727 if (Dist > Tol) { // On Cherche un meilleur point sur la section
728 myExt.Perform(PonPath);
729 if ( myExt.IsDone() ) {
730 DistMini(myExt, myAdpSection, Dist, SecParam);
731 PonSec = myAdpSection.Value(SecParam);
732 }
733 }
734 AngleMax = EvalAngle(VRef, dp1);
c6541a0c 735 if (isplan) AngleMax = M_PI/2 - AngleMax;
7fd59977 736 }
737
738 done = Standard_True;
739}
740
741//===============================================================
742// Function :IsDone
743// Purpose :
744//===============================================================
745 Standard_Boolean GeomFill_SectionPlacement::IsDone() const
746{
747 return done;
748}
749
750//===============================================================
751// Function : ParameterOnPath
752// Purpose :
753//===============================================================
754 Standard_Real GeomFill_SectionPlacement::ParameterOnPath() const
755{
756 return PathParam;
757}
758
759//===============================================================
760// Function : ParameterOnSection
761// Purpose :
762//===============================================================
763 Standard_Real GeomFill_SectionPlacement::ParameterOnSection() const
764{
765 return SecParam;
766}
767
768//===============================================================
769// Function : Distance
770// Purpose :
771//===============================================================
772 Standard_Real GeomFill_SectionPlacement::Distance() const
773{
774 return Dist;
775}
776
777//===============================================================
778// Function : Angle
779// Purpose :
780//===============================================================
781 Standard_Real GeomFill_SectionPlacement::Angle() const
782{
783 return AngleMax;
784}
785
786//===============================================================
787// Function : Transformation
788// Purpose :
789//===============================================================
790 gp_Trsf GeomFill_SectionPlacement::Transformation(
791 const Standard_Boolean WithTranslation,
792 const Standard_Boolean WithCorrection) const
793{
794 gp_Vec V;
795 gp_Mat M;
796 gp_Dir DT, DN, D;
797// modified by NIZHNY-MKK Fri Oct 17 15:27:07 2003
798 gp_Pnt P(0., 0., 0.), PSection(0., 0., 0.);
799
800 // Calcul des reperes
801 myLaw->D0(PathParam, M, V);
802
803 P.SetXYZ(V.XYZ());
804 D.SetXYZ (M.Column(3));
805 DN.SetXYZ(M.Column(1));
806 gp_Ax3 Paxe(P, D, DN);
807
808
809 if (WithTranslation || WithCorrection) {
810 if (myIsPoint)
811 PSection = myPoint;
812 else
813 PSection = mySection->Value(SecParam);
814 }
815 // modified by NIZHNY-MKK Wed Oct 8 15:03:19 2003.BEGIN
816 gp_Trsf Rot;
817
818 if (WithCorrection && !myIsPoint) {
819 Standard_Real angle;
820
821 if (!isplan)
9775fa61 822 throw Standard_Failure("Illegal usage: can't rotate non-planar profile");
7fd59977 823
824 gp_Dir ProfileNormal = TheAxe.Direction();
825 gp_Dir SpineStartDir = Paxe.Direction();
826
827 if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() )) {
828 gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
829 angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
830 gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
831 Rot.SetRotation( AxeOfRotation, angle );
832 }
833 PSection.Transform(Rot);
834 }
835 // modified by NIZHNY-MKK Wed Oct 8 15:03:21 2003.END
836
837 if (WithTranslation) {
838 P.ChangeCoord().SetLinearForm(-1,PSection.XYZ(),
839 V.XYZ());
840 }
841 else {
842 P.SetCoord(0., 0., 0.);
843 }
844
845 gp_Ax3 Saxe(P, gp::DZ(), gp::DX());
846
847 // Transfo
848 gp_Trsf Tf;
849 Tf.SetTransformation(Saxe, Paxe);
850
851 if (WithCorrection) {
852 // modified by NIZHNY-MKK Fri Oct 17 15:26:36 2003.BEGIN
853 // Standard_Real angle;
854 // gp_Trsf Rot;
855
856 // if (!isplan)
9775fa61 857 // throw Standard_Failure("Illegal usage: can't rotate non-planar profile");
7fd59977 858
859 // gp_Dir ProfileNormal = TheAxe.Direction();
860 // gp_Dir SpineStartDir = Paxe.Direction();
861 // if (! ProfileNormal.IsParallel( SpineStartDir, Precision::Angular() ))
862 // {
863 // gp_Dir DirAxeOfRotation = ProfileNormal ^ SpineStartDir;
864 // angle = ProfileNormal.AngleWithRef( SpineStartDir, DirAxeOfRotation );
865 // gp_Ax1 AxeOfRotation( TheAxe.Location(), DirAxeOfRotation );
866 // Rot.SetRotation( AxeOfRotation, angle );
867 // }
868 // modified by NIZHNY-MKK Fri Oct 17 15:26:42 2003.END
869
870 Tf *= Rot;
871 }
872
873 return Tf;
874}
875
876//===============================================================
877// Function : Section
878// Purpose :
879//===============================================================
880 Handle(Geom_Curve) GeomFill_SectionPlacement::
881 Section(const Standard_Boolean WithTranslation) const
882{
883 Handle(Geom_Curve) TheSection =
884 Handle(Geom_Curve)::DownCast(mySection->Copy());
885 TheSection->Transform(Transformation(WithTranslation, Standard_False));
886 return TheSection;
887}
888
889
890//===============================================================
891// Function :
892// Purpose :
893//===============================================================
894 Handle(Geom_Curve) GeomFill_SectionPlacement::
895 ModifiedSection(const Standard_Boolean WithTranslation) const
896{
897 Handle(Geom_Curve) TheSection =
898 Handle(Geom_Curve)::DownCast(mySection->Copy());
899 TheSection->Transform(Transformation(WithTranslation, Standard_True));
900 return TheSection;
901}
902
903//===============================================================
904// Function :SectionAxis
905// Purpose :
906//===============================================================
907 void GeomFill_SectionPlacement::SectionAxis(const gp_Mat& M,
908 gp_Vec& T,
909 gp_Vec& N,
910 gp_Vec& BN) const
911{
912 Standard_Real Eps = 1.e-10;
913 gp_Dir D;
914 gp_Vec PathNormal;
915 GeomLProp_CLProps CP(mySection, SecParam, 2, Eps);
916 if (CP.IsTangentDefined()) {
917 CP.Tangent(D);
918 T.SetXYZ(D.XYZ());
919 T.Normalize();
920 if (CP.Curvature() > Eps) {
921 CP.Normal(D);
922 N.SetXYZ(D.XYZ());
923 }
924 else {
925 // Cas ambigu, on essai de recuperer la normal a la trajectoire
926 // Reste le probleme des points d'inflexions, qui n'est pas
927 // bien traiter par LProp (pas de recuperation de la normal) :
928 // A voir...
929 PathNormal.SetXYZ(M.Column(1));
930 PathNormal.Normalize();
931 BN = T ^ PathNormal;
932 if (BN.Magnitude() > Eps) {
933 BN.Normalize();
934 //N = BN ^ T;
935 }
936 N = BN ^ T;
937 }
938 }
939 else { // Cas indefinie, on prend le triedre
940 // complet sur la trajectoire
941 T.SetXYZ(M.Column(3));
942 N.SetXYZ(M.Column(2));
943 }
944 BN = T ^ N;
945}
946
947
948//===============================================================
949// Function :Choix
950// Purpose : Decide si le couple (dist, angle) est "meilleur" que
951// couple courrant.
952//===============================================================
953 Standard_Boolean GeomFill_SectionPlacement::Choix(const Standard_Real dist,
954 const Standard_Real angle) const
955{
956 Standard_Real evoldist, evolangle;
957 evoldist = dist - Dist;
958 evolangle = angle - AngleMax;
959 // (1) Si la gain en distance est > que le gabarit, on prend
960 if (evoldist < - Gabarit)
961 return Standard_True;
962
963 // (2) si l'ecart en distance est de l'ordre du gabarit
964 if (Abs(evoldist) < Gabarit) {
965// (2.1) si le gain en angle est important on garde
966 if (evolangle > 0.5)
967 return Standard_True;
968 // (2.2) si la variation d'angle est moderee on evalue
969 // une fonction de penalite
970 if (Penalite(angle, dist/Gabarit) < Penalite(AngleMax, Dist/Gabarit) )
971 return Standard_True;
972 }
973
974 return Standard_False;
975}
976
977
978
979
980
981
982
983