1 // Created on: 1998-07-02
2 // Created by: Stephanie HUMEAU
3 // Copyright (c) 1998-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
23 #include <GeomFill_GuideTrihedronPlan.ixx>
26 #include <gp_Pnt2d.hxx>
27 //#include <gp_Trsf2d.hxx>
28 //#include <Bnd_Box2d.hxx>
31 #include <Adaptor3d_Curve.hxx>
32 #include <GeomAdaptor_HCurve.hxx>
33 #include <GeomAdaptor_HSurface.hxx>
35 #include <Geom_Plane.hxx>
37 #include <IntCurveSurface_IntersectionPoint.hxx>
38 #include <IntCurveSurface_HInter.hxx>
40 #include <GeomFill_Frenet.hxx>
41 #include <GeomFill_PlanFunc.hxx>
43 #include <math_Vector.hxx>
45 #include <math_FunctionRoot.hxx>
46 #include <math_Matrix.hxx>
48 #include <Precision.hxx>
51 #include <DrawTrSurf.hxx>
55 static void TracePlan(const Handle(Geom_Surface)& /*Plan*/)
57 cout << "Pas d'intersection Guide/Plan" << endl;
59 char* Temp = "ThePlan" ;
60 DrawTrSurf::Set(Temp, Plan);
61 // DrawTrSurf::Set("ThePlan", Plan);
66 //==================================================================
67 //Function: InGoodPeriod
68 //Purpose : Recadre un paramtere
69 //==================================================================
70 static void InGoodPeriod(const Standard_Real Prec,
71 const Standard_Real Period,
72 Standard_Real& Current)
74 Standard_Real Diff=Current-Prec;
75 Standard_Integer nb = (Standard_Integer ) IntegerPart(Diff/Period);
78 if (Diff > Period/2) Current -= Period;
79 else if (Diff < -Period/2) Current += Period;
82 //=======================================================================
83 //function : GuideTrihedronPlan
84 //purpose : Constructor
85 //=======================================================================
86 GeomFill_GuideTrihedronPlan::GeomFill_GuideTrihedronPlan (const Handle(Adaptor3d_HCurve)& theGuide) :
90 myStatus(GeomFill_PipeOk)
93 myGuide = theGuide; // guide
95 myNbPts = 20; // nb points pour calculs
96 Pole = new (TColgp_HArray2OfPnt2d)(1,1,1,myNbPts);//tab pr stocker Pprime (pt sur guide)
97 frenet = new (GeomFill_Frenet)();
99 XTol(1) = myGuide->Resolution(1.e-6);
102 //=======================================================================
104 //purpose : calcule myNbPts points sur la courbe guide (<=> normale)
105 //=======================================================================
106 void GeomFill_GuideTrihedronPlan::Init()
108 myStatus = GeomFill_PipeOk;
111 // Box.Update(-0.1, -0.1, 0.1, 0.1); // Taille minimal
112 gp_Vec Tangent,Normal,BiNormal;
114 Standard_Real t, DeltaG, w = 0.;
115 Standard_Real f = myCurve->FirstParameter();
116 Standard_Real l = myCurve->LastParameter();
120 Handle(Geom_Plane) Plan;
121 Handle(GeomAdaptor_HSurface) Pl;
122 IntCurveSurface_IntersectionPoint PInt;
123 IntCurveSurface_HInter Int;
124 frenet->SetCurve(myCurve);
125 DeltaG = (myGuide->LastParameter() - myGuide->FirstParameter())/2;
127 Inf(1) = myGuide->FirstParameter() - DeltaG;
128 Sup(1) = myGuide->LastParameter() + DeltaG;
130 if (!myGuide->IsPeriodic()) {
131 myTrimG = myGuide->Trim(myGuide->FirstParameter()- DeltaG/100,
132 myGuide->LastParameter() + DeltaG/100,
138 // Standard_Real Step = DeltaG/100;
140 for (ii=1; ii<=myNbPts; ii++)
142 t = Standard_Real(myNbPts - ii)*f + Standard_Real(ii - 1)*l;
145 frenet->D0(t, Tangent, Normal, BiNormal);
146 Plan = new (Geom_Plane) (P, Tangent);
147 Pl = new(GeomAdaptor_HSurface) (Plan);
149 Int.Perform(myTrimG, Pl); // intersection plan / guide
150 if (Int.NbPoints() == 0) {
154 w = (fabs(myGuide->LastParameter() -w) > fabs(myGuide->FirstParameter()-w) ? myGuide->FirstParameter() : myGuide->LastParameter());
156 myStatus = GeomFill_PlaneNotIntersectGuide;
164 Standard_Real Dmin = P.Distance(Pmin);
165 for (Standard_Integer jj=2;jj<=Int.NbPoints();jj++)
167 Pmin = Int.Point(jj).Pnt();
168 if (P.Distance(Pmin) < Dmin)
170 PInt = Int.Point(jj);
171 Dmin = P.Distance(Pmin);
178 Standard_Real Diff = w - Pole->Value(1, ii-1).Y();
179 if (Abs(Diff) > DeltaG) {
180 if (myGuide->IsPeriodic()) {
181 InGoodPeriod (Pole->Value(1, ii-1).Y(),
182 myGuide->Period(), w);
184 Diff = w - Pole->Value(1, ii-1).Y();
189 if (Abs(Diff) > DeltaG) {
190 cout << "Trihedron Plan Diff on Guide : " <<
196 gp_Pnt2d p1(t, w); // on stocke les parametres
197 Pole->SetValue(1, ii, p1);
202 //=======================================================================
203 //function : SetCurve
204 //purpose : calculation of trihedron
205 //=======================================================================
206 void GeomFill_GuideTrihedronPlan::SetCurve(const Handle(Adaptor3d_HCurve)& C)
209 if (!myCurve.IsNull()) Init();
212 //=======================================================================
214 //purpose : calculation of trihedron
215 //=======================================================================
217 Handle(Adaptor3d_HCurve) GeomFill_GuideTrihedronPlan::Guide()const
222 //=======================================================================
224 //purpose : calculation of trihedron
225 //=======================================================================
226 Standard_Boolean GeomFill_GuideTrihedronPlan::D0(const Standard_Real Param,
234 myCurve->D0(Param, P);
236 frenet->D0(Param,Tangent,Normal,BiNormal);
238 //initialisation de la recherche
241 Standard_Integer Iter = 50;
243 // fonction dont il faut trouver la racine : G(W)-Pl(U,V)=0
244 GeomFill_PlanFunc E(P, Tangent, myGuide);
247 math_FunctionRoot Result(E, X(1), XTol(1),
248 Inf(1), Sup(1), Iter);
252 Standard_Real Res = Result.Root();
253 // R = Result.Root(); // solution
255 Pprime = myTrimG->Value(Res); // pt sur courbe guide
256 gp_Vec n (P, Pprime); // vecteur definissant la normale du triedre
258 Normal = n.Normalized();
259 BiNormal = Tangent.Crossed(Normal);
260 BiNormal.Normalized();
265 // plan ortho a la trajectoire pour determiner Pprime
266 Handle(Geom_Plane) Plan = new (Geom_Plane)(P, Tangent);
269 myStatus = GeomFill_PlaneNotIntersectGuide;
270 return Standard_False;
273 return Standard_True;
276 //=======================================================================
278 //purpose : calculation of trihedron and first derivative
279 //=======================================================================
280 Standard_Boolean GeomFill_GuideTrihedronPlan::D1(const Standard_Real Param,
288 // return Standard_False;
294 // triedre de frenet sur la trajectoire
295 myCurve->D1(Param, P, To);
296 frenet->D1(Param,Tangent,DTangent,Normal,DNormal,BiNormal,DBiNormal);
300 Standard_Integer Iter = 50;
302 // fonction dont il faut trouver la racine : G(W)-Pl(U,V)=0
304 GeomFill_PlanFunc E(P, Tangent, myGuide);
307 math_FunctionRoot Result(E, X(1), XTol(1),
308 Inf(1), Sup(1), Iter);
312 Standard_Real Res = Result.Root();
313 // R = Result.Root(); // solution
314 myTrimG->D1(Res, PG, TG);
315 gp_Vec n (P, PG), dn; // vecteur definissant la normale du triedre
316 Standard_Real Norm = n.Magnitude();
324 BiNormal = Tangent.Crossed(Normal);
326 // derivee premiere du triedre
327 Standard_Real dedx, dedt, dtg_dt;
328 E.Derivative(Res, dedx);
329 E.DEDT(Res, To, DTangent, dedt);
333 /* Standard_Real h=1.e-7, e, etg, etc;
336 if ( Abs( (etg-e)/h - dedx) > 1.e-4) {
337 cout << "err :" << (etg-e)/h - dedx << endl;
341 myCurve->D0(Param+h, pdbg);
342 frenet->D0(Param+h,td, nb, bnb);
344 GeomFill_PlanFunc Edeb(pdbg, td, myGuide);
345 Edeb.Value(Res, etc);
346 if ( Abs( (etc-e)/h - dedt) > 1.e-4) {
347 cout << "err :" << (etc-e)/h - dedt << endl;
350 dn.SetLinearForm(dtg_dt, TG, -1, To);
352 DNormal.SetLinearForm(-(n*dn), n, dn);
354 DBiNormal.SetLinearForm(Tangent.Crossed(DNormal),
355 DTangent.Crossed(Normal));
360 // plan ortho a la trajectoire
361 Handle(Geom_Plane) Plan = new (Geom_Plane)(P, Tangent);
364 myStatus = GeomFill_PlaneNotIntersectGuide;
365 return Standard_False;
368 return Standard_True;
372 //=======================================================================
374 //purpose : calculation of trihedron and derivatives
375 //=======================================================================
376 Standard_Boolean GeomFill_GuideTrihedronPlan::D2(const Standard_Real Param,
389 // gp_Vec To,DTo,TG,DTG;
392 myCurve->D2(Param, P, To, DTo);
394 // triedre de Frenet sur la trajectoire
395 frenet->D2(Param,Tangent,DTangent,D2Tangent,
396 Normal,DNormal,D2Normal,
397 BiNormal,DBiNormal,D2BiNormal);
400 // plan ortho a Tangent pour trouver la pt Pprime sur le guide
401 Handle(Geom_Plane) Plan = new (Geom_Plane)(P, Tangent);
402 Handle(GeomAdaptor_HSurface) Pl= new(GeomAdaptor_HSurface)(Plan);
405 Standard_Integer Iter = 50;
406 // fonction dont il faut trouver la racine : G(W) - Pl(U,V)=0
407 GeomFill_FunctionPipe E(Pl , myGuide);
411 math_FunctionSetRoot Result(E, X, XTol,
416 R = Result.Root(); // solution
417 myTrimG->D2(R(1), PG, TG, DTG);
419 gp_Vec n (P, PG); // vecteur definissant la normale du triedre
420 Standard_Real Norm = n.Magnitude();
422 Normal = n.Normalized();
423 BiNormal = Tangent.Crossed(Normal);
427 // derivee premiere du triedre
428 Standard_Real dtp_dt;
429 dtp_dt = (To*Tangent - Norm*(n*DTangent))/(Tangent*TG);
431 dn.SetLinearForm(dtp_dt, TG, -1, To);
433 DNormal.SetLinearForm(-(n*dn), n, dn);
435 DBiNormal = Tangent.Crossed(DNormal) + DTangent.Crossed(Normal);
437 // derivee seconde du triedre
438 Standard_Real d2tp_dt2;
439 d2tp_dt2 = (DTo*Tangent+To*DTangent - dn*DTangent-Norm*n*D2Tangent)/(TG*Tangent)
440 - (To*Tangent-Norm*n*DTangent) * (DTG*dtp_dt*Tangent+TG*DTangent)
441 / ((TG*Tangent)*(TG*Tangent));
444 d2n.SetLinearForm(dtp_dt*dtp_dt, DTG, d2tp_dt2, TG, -DTo);
448 D2Normal.SetLinearForm(3*Pow(n*dn,2)- (dn.SquareMagnitude() + n*d2n), n,
452 D2BiNormal.SetLinearForm(1, D2Tangent.Crossed(Normal),
453 2, DTangent.Crossed(DNormal),
454 Tangent.Crossed(D2Normal));
461 myStatus = GeomFill_PlaneNotIntersectGuide;
462 return Standard_False;
465 // return Standard_True;
466 return Standard_False;
470 //=======================================================================
473 //=======================================================================
474 Handle(GeomFill_TrihedronLaw) GeomFill_GuideTrihedronPlan::Copy() const
476 Handle(GeomFill_GuideTrihedronPlan) copy =
477 new (GeomFill_GuideTrihedronPlan) (myGuide);
478 copy->SetCurve(myCurve);
482 //=======================================================================
483 //function : ErrorStatus
485 //=======================================================================
486 GeomFill_PipeError GeomFill_GuideTrihedronPlan::ErrorStatus() const
492 //=======================================================================
493 //function : NbIntervals
494 //purpose : Version provisoire : Il faut tenir compte du guide
495 //=======================================================================
496 Standard_Integer GeomFill_GuideTrihedronPlan::NbIntervals(const GeomAbs_Shape S)const
501 case GeomAbs_C0: tmpS = GeomAbs_C1; break;
502 case GeomAbs_C1: tmpS = GeomAbs_C2; break;
503 case GeomAbs_C2: tmpS = GeomAbs_C3; break;
504 default: tmpS = GeomAbs_CN;
507 Nb = myCurve->NbIntervals(tmpS);
510 //======================================================================
511 //function :Intervals
513 //=======================================================================
514 void GeomFill_GuideTrihedronPlan::Intervals(TColStd_Array1OfReal& TT,
515 const GeomAbs_Shape S) const
519 case GeomAbs_C0: tmpS = GeomAbs_C1; break;
520 case GeomAbs_C1: tmpS = GeomAbs_C2; break;
521 case GeomAbs_C2: tmpS = GeomAbs_C3; break;
522 default: tmpS = GeomAbs_CN;
524 myCurve->Intervals(TT, tmpS);
527 //======================================================================
528 //function :SetInterval
530 //=======================================================================
531 void GeomFill_GuideTrihedronPlan::SetInterval(const Standard_Real First,
532 const Standard_Real Last)
534 myTrimmed = myCurve->Trim(First, Last, Precision::Confusion());
538 //=======================================================================
539 //function : GetAverageLaw
541 //=======================================================================
542 void GeomFill_GuideTrihedronPlan::GetAverageLaw(gp_Vec& ATangent,
547 Standard_Real t, Delta = (myCurve->LastParameter() -
548 myCurve->FirstParameter())/20.001;
550 ATangent.SetCoord(0.,0.,0.);
551 ANormal.SetCoord(0.,0.,0.);
552 ABiNormal.SetCoord(0.,0.,0.);
555 for (ii=1, T; ii<=20; ii++) {
556 t = myCurve->FirstParameter() +(ii-1)*Delta;
567 //=======================================================================
568 //function : IsConstant
570 //=======================================================================
571 Standard_Boolean GeomFill_GuideTrihedronPlan::IsConstant() const
573 if ((myCurve->GetType() == GeomAbs_Line) &&
574 (myGuide->GetType() == GeomAbs_Line)) {
576 Angle = myCurve->Line().Angle(myGuide->Line());
577 if ((Angle<1.e-12) || ((2*M_PI-Angle)<1.e-12) )
578 return Standard_True;
581 return Standard_False;
584 //=======================================================================
585 //function : IsOnlyBy3dCurve
587 //=======================================================================
588 Standard_Boolean GeomFill_GuideTrihedronPlan::IsOnlyBy3dCurve() const
590 return Standard_False;
593 //=======================================================================
595 //purpose : Nothing!!
596 //=======================================================================
597 void GeomFill_GuideTrihedronPlan::Origine(const Standard_Real ,
598 const Standard_Real )
602 //==================================================================
604 //Purpose : recherche par interpolation d'une valeur initiale
605 //==================================================================
606 void GeomFill_GuideTrihedronPlan::InitX(const Standard_Real Param)
609 Standard_Integer Ideb = 1, Ifin = Pole->RowLength(), Idemi;
610 Standard_Real Valeur, t1, t2;
613 Valeur = Pole->Value(1, Ideb).X();
614 if (Param == Valeur) {
618 Valeur = Pole->Value(1, Ifin).X();
619 if (Param == Valeur) {
623 while ( Ideb+1 != Ifin) {
624 Idemi = (Ideb+Ifin)/2;
625 Valeur = Pole->Value(1, Idemi).X();
626 if (Valeur < Param) {
630 if ( Valeur > Param) { Ifin = Idemi;}
638 t1 = Pole->Value(1,Ideb).X();
639 t2 = Pole->Value(1,Ifin).X();
640 Standard_Real diff = t2-t1;
642 Standard_Real b = (Param-t1) / diff,
643 a = (t2-Param) / diff;
645 X(1) = Pole->Value(1,Ideb).Coord(2) * a
646 + Pole->Value(1,Ifin).Coord(2) * b; //param guide
649 X(1) = (Pole->Value(1, Ideb).Coord(2) +
650 Pole->Value(1, Ifin).Coord(2)) / 2;
652 if (myGuide->IsPeriodic()) {
653 X(1) = ElCLib::InPeriod(X(1), myGuide->FirstParameter(),
654 myGuide->LastParameter());