1 // Created on: 1998-04-21
2 // Created by: Stephanie HUMEAU
3 // Copyright (c) 1998-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
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
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
18 #include <Adaptor3d_HCurve.hxx>
19 #include <Adaptor3d_HSurface.hxx>
20 #include <Geom_Line.hxx>
21 #include <Geom_Surface.hxx>
22 #include <GeomAdaptor_HCurve.hxx>
23 #include <GeomAdaptor_HSurface.hxx>
24 #include <GeomFill_DraftTrihedron.hxx>
25 #include <GeomFill_FunctionDraft.hxx>
26 #include <GeomFill_LocationDraft.hxx>
27 #include <GeomFill_LocationLaw.hxx>
28 #include <GeomFill_Tensor.hxx>
29 #include <GeomFill_TrihedronLaw.hxx>
34 #include <IntCurveSurface_HInter.hxx>
35 #include <IntCurveSurface_Intersection.hxx>
36 #include <IntCurveSurface_IntersectionPoint.hxx>
37 #include <math_FunctionSetWithDerivatives.hxx>
38 #include <math_Gauss.hxx>
39 #include <math_Matrix.hxx>
40 #include <math_NewtonFunctionSetRoot.hxx>
41 #include <math_Vector.hxx>
42 #include <Standard_NotImplemented.hxx>
43 #include <Standard_OutOfRange.hxx>
44 #include <Standard_Type.hxx>
46 //==================================================================
47 //Function: GeomFill_LocationDraft
48 //Purpose : constructor
49 //==================================================================
50 GeomFill_LocationDraft::GeomFill_LocationDraft
51 (const gp_Dir& Direction,
52 const Standard_Real Angle)
54 myDir = Direction; // direction de depouille
56 myAngle = Angle; // angle de depouille (teta prime)
59 myLaw = new (GeomFill_DraftTrihedron)(myDir, Angle); // triedre
60 myNbPts = 41; // nb de points utilises pour les calculs
61 myPoles2d = new (TColgp_HArray1OfPnt2d)(1, 2*myNbPts);
62 Intersec = Standard_False; //intersection avec surface d'arret ?
63 WithTrans = Standard_False;
68 //==================================================================
71 //==================================================================
72 Handle(GeomFill_LocationLaw) GeomFill_LocationDraft::Copy() const
74 Handle(GeomFill_TrihedronLaw) law;
76 Handle(GeomFill_LocationDraft) copy =
77 new (GeomFill_LocationDraft) (myDir,myAngle);
78 copy->SetCurve(myCurve);
79 copy->SetStopSurf(mySurf);
80 if (WithTrans) copy->SetTrsf(Trans);
85 //==================================================================
88 //==================================================================
89 void GeomFill_LocationDraft::SetTrsf(const gp_Mat& Transfo)
95 WithTrans = Standard_False; // Au cas ou Trans = I
96 for (Standard_Integer ii=1; ii<=3 && !WithTrans ; ii++)
97 for (Standard_Integer jj=1; jj<=3 && !WithTrans; jj++)
98 if (Abs(Aux.Value(ii, jj)) > 1.e-14) WithTrans = Standard_True;
101 //==================================================================
103 //Purpose : Calcul des poles sur la surfaces d'arret (intersection
104 // entre la generatrice et la surface en myNbPts points de la section)
105 //==================================================================
106 void GeomFill_LocationDraft::SetCurve(const Handle(Adaptor3d_HCurve)& C)
115 //==================================================================
116 //Function: SetStopSurf
118 //==================================================================
119 void GeomFill_LocationDraft::SetStopSurf(const Handle(Adaptor3d_HSurface)& Surf)
125 //==================================================================
128 //==================================================================
129 void GeomFill_LocationDraft::SetAngle(const Standard_Real Angle)
132 myLaw->SetAngle(myAngle);
136 //==================================================================
138 //Purpose : Poses les jalon de l'intersection : depouille / Surface
139 //==================================================================
140 void GeomFill_LocationDraft::Prepare()
142 if (mySurf.IsNull()) {
143 Intersec = Standard_False;
147 Intersec = Standard_True;
149 Standard_Integer ii,jj;
150 Standard_Real f, l, t;
154 IntCurveSurface_IntersectionPoint P1,P2;
155 f = myCurve->FirstParameter();
156 l = myCurve->LastParameter();
158 for (ii=1; ii<=myNbPts; ii++)
160 t = Standard_Real(myNbPts - ii)*f + Standard_Real(ii - 1)*l;
167 D = Cos(myAngle)*B + Sin(myAngle)*N;
169 L = new (Geom_Line) (P, D);
171 IntCurveSurface_HInter Int; // intersection surface / generatrice
172 Handle(GeomAdaptor_HCurve) AC = new (GeomAdaptor_HCurve) (L);
173 Int.Perform(AC, mySurf); // calcul de l'intersection
175 if (Int.NbPoints() > 0) // il y a au moins 1 intersection
177 P1 = Int.Point(1); // 1er point d'intersection
179 for (jj=2 ; jj<=Int.NbPoints() ; jj++)
182 if(P1.W() > P2.W()) P1 = P2; // point le plus proche
185 gp_Pnt2d p (P1.W(), t); // point de la courbe
186 gp_Pnt2d q (P1.U(),P1.V()); // point sur la surface
187 myPoles2d->SetValue(2*ii-1,p); // point de la courbe (indice impair)
188 myPoles2d->SetValue(2*ii,q); // point sur la surface (indice pair)
191 {// au moins un point ou il n'y a pas intersection
192 Intersec = Standard_False;
198 //==================================================================
200 //Purpose : return the path
201 //==================================================================
202 const Handle(Adaptor3d_HCurve)& GeomFill_LocationDraft::GetCurve() const
208 //==================================================================
211 //==================================================================
212 Standard_Boolean GeomFill_LocationDraft::D0(const Standard_Real Param,
220 myTrimmed->D0(Param, P);
223 Ok = myLaw->D0(Param, T, N, B);
225 M.SetCols(N.XYZ(), B.XYZ(), T.XYZ());
231 return Standard_True;
234 //==================================================================
236 //Purpose : calcul de l'intersection (C0) sur la surface
237 //==================================================================
238 Standard_Boolean GeomFill_LocationDraft::D0(const Standard_Real Param,
241 TColgp_Array1OfPnt2d& Poles2d)
244 // gp_Vec D,T,N,B,DT,DN,DB;
248 myCurve->D0(Param, P);
250 Ok = myLaw->D0(Param, T, N, B);
252 M.SetCols(N.XYZ(), B.XYZ(), T.XYZ());
258 if (Intersec == Standard_True) {
259 // la generatrice intersecte la surface d'arret
261 D = Cos(myAngle)*B + Sin(myAngle)*N;
263 Handle(Geom_Line) L = new (Geom_Line) (P, D);
264 Handle(GeomAdaptor_HCurve) G = new (GeomAdaptor_HCurve) (L);
266 Standard_Real t1,t2,Paramt1,t2Param;
267 Standard_Real U0=0,V0=0,W0=0;
269 Standard_Integer ii = 1;
271 // on recherche l'intervalle auquel appartient Param
272 while (ii<2*myNbPts && myPoles2d->Value(ii).Coord(2) < Param) ii=ii+2;
274 if (ii<2*myNbPts && !IsEqual(myPoles2d->Value(ii).Coord(2),Param))
277 // interpolation lineaire pour initialiser le germe de la recherche
278 t1 = myPoles2d->Value(ii).Coord(2);
279 t2 = myPoles2d->Value(ii-2).Coord(2);
281 Paramt1 = (Param-t1) / (t2-t1);
282 t2Param = (t2-Param) / (t2-t1);
284 W0 = myPoles2d->Value(ii-2).Coord(1)*Paramt1
285 + myPoles2d->Value(ii).Coord(1)*t2Param;
286 U0 = myPoles2d->Value(ii-1).Coord(1)*Paramt1
287 + myPoles2d->Value(ii+1).Coord(1)*t2Param;
288 V0 = myPoles2d->Value(ii-1).Coord(2)*Paramt1
289 + myPoles2d->Value(ii+1).Coord(2)*t2Param;
291 // on est sur un param ou les points ont deja ete calcules
292 else if (ii<2*myNbPts && IsEqual(myPoles2d->Value(ii).Coord(2) ,Param))
294 W0 = myPoles2d->Value(ii).Coord(1);
295 U0 = myPoles2d->Value(ii+1).Coord(1);
296 V0 = myPoles2d->Value(ii+1).Coord(2);
299 // recherche de la solution (pt d'intersection generatrice / surface)
307 math_Vector XTol(1,3);
311 Standard_Real FTol = 0.0000001;
312 Standard_Integer Iter = 100;
314 // fonction dont il faut trouver la racine : G(W)-S(U,V)=0
315 GeomFill_FunctionDraft E(mySurf , G);
318 math_NewtonFunctionSetRoot Result(E, XTol, FTol, Iter);
319 Result.Perform(E, X);
324 Result.Root(R); // solution
326 gp_Pnt2d p (R(2), R(3)); // point sur la surface
327 gp_Pnt2d q (R(1), Param); // point de la courbe
328 Poles2d.SetValue(1,p);
329 Poles2d.SetValue(2,q);
332 return Standard_False;
336 // la generatrice n'intersecte pas la surface d'arret
337 return Standard_True;
340 //==================================================================
342 //Purpose : calcul de l'intersection (C1) sur la surface
343 //==================================================================
344 Standard_Boolean GeomFill_LocationDraft::D1(const Standard_Real Param,
349 TColgp_Array1OfPnt2d& Poles2d,
350 TColgp_Array1OfVec2d& DPoles2d)
353 gp_Vec D,T,N,B,DT,DN,DB;
356 myCurve->D1(Param, P, DV);
359 Ok = myLaw->D1(Param, T, DT, N, DN, B, DB);
360 if (!Ok) return Standard_False;
362 M.SetCols(N.XYZ(), B.XYZ(), T.XYZ());
363 DM.SetCols(DN.XYZ(), DB.XYZ(), DT.XYZ());
370 if (Intersec == Standard_True)
371 { // la generatrice intersecte la surface d'arret
373 D = Cos(myAngle)*B + Sin(myAngle)*N;
375 Handle(Geom_Line) L = new (Geom_Line) (P, D);
376 Handle(GeomAdaptor_HCurve) G = new (GeomAdaptor_HCurve) (L);
378 Standard_Real t1,t2,Paramt1,t2Param;
379 Standard_Real U0=0,V0=0,W0=0;
381 Standard_Integer ii = 1;
383 // recherche de la solution (pt d'intersection generatrice / surface)
385 // on recherche l'intervalle auquel appartient Param
386 while (ii < 2*myNbPts && myPoles2d->Value(ii).Coord(2) < Param) ii=ii+2;
389 if (ii < 2*myNbPts && !IsEqual(myPoles2d->Value(ii).Coord(2) ,Param))
391 // interpolation lineaire pour initialiser le germe de la recherche
392 t1 = myPoles2d->Value(ii).Coord(2);
393 t2 = myPoles2d->Value(ii-2).Coord(2);
395 Paramt1 = (Param-t1) / (t2-t1);
396 t2Param = (t2-Param) / (t2-t1);
398 W0 = myPoles2d->Value(ii-2).Coord(1)*Paramt1
399 + myPoles2d->Value(ii).Coord(1)*t2Param;
400 U0 = myPoles2d->Value(ii-1).Coord(1)*Paramt1
401 + myPoles2d->Value(ii+1).Coord(1)*t2Param;
402 V0 = myPoles2d->Value(ii-1).Coord(2)*Paramt1
403 + myPoles2d->Value(ii+1).Coord(2)*t2Param;
405 else if (ii<2*myNbPts && IsEqual(myPoles2d->Value(ii).Coord(2) ,Param))
407 W0 = myPoles2d->Value(ii).Coord(1);
408 U0 = myPoles2d->Value(ii+1).Coord(1);
409 V0 = myPoles2d->Value(ii+1).Coord(2);
419 math_Vector XTol(1,3);
423 Standard_Real FTol = 0.000001;
424 Standard_Integer Iter = 100;
427 // fonction dont il faut trouver la racine : G(W)-S(U,V)=0
428 GeomFill_FunctionDraft E(mySurf,G);
432 math_NewtonFunctionSetRoot Result(E, XTol, FTol, Iter);
433 Result.Perform(E, X);
438 Result.Root(R); // solution
440 gp_Pnt2d p (R(2), R(3)); // point sur la surface
441 gp_Pnt2d q (R(1), Param); // point de la courbe
442 Poles2d.SetValue(1,p);
443 Poles2d.SetValue(2,q);
445 // derivee de la fonction par rapport a Param
446 math_Vector DEDT(1,3,0);
447 E.DerivT(myTrimmed, Param, R(1), DN, myAngle, DEDT); // dE/dt => DEDT
449 math_Vector DSDT (1,3,0);
450 math_Matrix DEDX (1,3,1,3,0);
451 E.Derivatives(R, DEDX); // dE/dx au point R => DEDX
453 // resolution du syst. lin. : DEDX*DSDT = -DEDT
457 Ga.Solve (DEDT.Opposite(), DSDT); // resolution du syst. lin.
458 gp_Vec2d dp (DSDT(2), DSDT(3)); // surface
459 gp_Vec2d dq (DSDT(1), 1); // courbe
460 DPoles2d.SetValue(1, dp);
461 DPoles2d.SetValue(2, dq);
466 else {// la generatrice n'intersecte pas la surface d'arret
467 return Standard_False;
470 return Standard_True;
473 //==================================================================
475 //Purpose : calcul de l'intersection (C2) sur la surface
476 //==================================================================
477 Standard_Boolean GeomFill_LocationDraft::D2(const Standard_Real Param,
484 TColgp_Array1OfPnt2d& Poles2d,
485 TColgp_Array1OfVec2d& DPoles2d,
486 TColgp_Array1OfVec2d& D2Poles2d)
489 gp_Vec D,T,N,B,DT,DN,DB,D2T,D2N,D2B;
492 myCurve->D2(Param, P, DV, D2V);
495 Ok = myLaw->D2(Param, T, DT, D2T, N, DN, D2N, B, DB, D2B);
498 M.SetCols(N.XYZ(), B.XYZ(), T.XYZ());
499 DM.SetCols(DN.XYZ(), DB.XYZ(), DT.XYZ());
500 D2M.SetCols(D2N.XYZ(), D2B.XYZ(), D2T.XYZ());
507 if (Intersec == Standard_True)
508 {// la generatrice intersecte la surface d'arret
511 D = Cos(myAngle) * B + Sin(myAngle) * N;
513 Handle(Geom_Line) L = new (Geom_Line) (P, D);
514 Handle(GeomAdaptor_HCurve) G = new (GeomAdaptor_HCurve) (L);
516 Standard_Real t1,t2,Paramt1,t2Param;
517 Standard_Real U0=0,V0=0,W0=0;
519 Standard_Integer ii = 1;
521 // on recherche l'intervalle auquel appartient Param
522 while (ii<2*myNbPts && myPoles2d->Value(ii).Coord(2) < Param) ii=ii+2;
524 if (ii<2*myNbPts && !IsEqual(myPoles2d->Value(ii).Coord(2) ,Param))
527 // interpolation lineaire pour initialiser le germe de la recherche
528 t1 = myPoles2d->Value(ii).Coord(2);
529 t2 = myPoles2d->Value(ii-2).Coord(2);
531 Paramt1 = (Param-t1) / (t2-t1);
532 t2Param = (t2-Param) / (t2-t1);
534 W0 = myPoles2d->Value(ii-2).Coord(1)*Paramt1 +
535 myPoles2d->Value(ii).Coord(1)*t2Param;
536 U0 = myPoles2d->Value(ii-1).Coord(1)*Paramt1 +
537 myPoles2d->Value(ii+1).Coord(1)*t2Param;
538 V0 = myPoles2d->Value(ii-1).Coord(2)*Paramt1 +
539 myPoles2d->Value(ii+1).Coord(2)*t2Param;
541 else if (ii<2*myNbPts && IsEqual(myPoles2d->Value(ii).Coord(2) ,Param))
543 W0 = myPoles2d->Value(ii).Coord(1);
544 U0 = myPoles2d->Value(ii+1).Coord(1);
545 V0 = myPoles2d->Value(ii+1).Coord(2);
548 // recherche de la solution (pt d'intersection generatrice / surface)
556 math_Vector XTol(1,3);
560 Standard_Real FTol = 0.000001;
561 Standard_Integer Iter = 150;
564 // fonction dont il faut trouver la racine : G(W)-S(U,V)=0
565 GeomFill_FunctionDraft E(mySurf,G);
569 math_NewtonFunctionSetRoot Result (E, XTol, FTol, Iter);
570 Result.Perform(E, X);
575 Result.Root(R); // solution
578 gp_Pnt2d p (R(2), R(3)); // point sur la surface
579 gp_Pnt2d q (R(1), Param); // point de la courbe
580 Poles2d.SetValue(1,p);
581 Poles2d.SetValue(2,q);
584 // premiere derivee de la fonction
585 math_Vector DEDT(1,3,0);
586 E.DerivT(myTrimmed, Param, R(1), DN, myAngle, DEDT); // dE/dt => DEDT
588 math_Vector DSDT (1,3,0);
589 math_Matrix DEDX (1,3,1,3,0);
590 E.Derivatives(R, DEDX); // dE/dx => DEDX
592 // resolution du syst. lin.
593 math_Gauss Ga (DEDX);
596 Ga.Solve (DEDT.Opposite(), DSDT);
597 gp_Vec2d dp (DSDT(2), DSDT(3)); // surface
598 gp_Vec2d dq (DSDT(1), 1); // courbe
599 DPoles2d.SetValue(1, dp);
600 DPoles2d.SetValue(2, dq);
605 GeomFill_Tensor D2EDX2(3,3,3);
606 E.Deriv2X(R, D2EDX2); // d2E/dx2
608 math_Vector D2EDT2(1,3,0);
609 E.Deriv2T(myTrimmed, Param, R(1), D2N, myAngle ,D2EDT2); // d2E/dt2
611 math_Matrix D2EDTDX(1,3,1,3,0);
612 E.DerivTX(DN, myAngle, D2EDTDX); // d2E/dtdx
614 math_Vector D2SDT2(1,3,0); // d2s/dt2
615 math_Matrix T(1,3,1,3,0);
616 D2EDX2.Multiply(DSDT,T);
618 // resolution du syst. lin.
619 math_Gauss Ga1 (DEDX);
622 Ga1.Solve ( -T*DSDT - 2*D2EDTDX*DSDT - D2EDT2 , D2SDT2);
623 gp_Vec2d d2p (D2SDT2(2), D2SDT2(3)); // surface
624 gp_Vec2d d2q (D2SDT2(1), 0); // courbe
625 D2Poles2d.SetValue(1, d2p);
626 D2Poles2d.SetValue(2, d2q);
628 else {// la generatrice n'intersecte pas la surface d'arret
629 return Standard_False;
634 return Standard_True;
637 //==================================================================
638 //Function : HasFirstRestriction
640 //==================================================================
641 Standard_Boolean GeomFill_LocationDraft::HasFirstRestriction() const
643 return Standard_False;
646 //==================================================================
647 //Function : HasLastRestriction
649 //==================================================================
650 Standard_Boolean GeomFill_LocationDraft::HasLastRestriction() const
653 if (Intersec == Standard_True) return Standard_True;
654 else return Standard_False;
657 //==================================================================
658 //Function : TraceNumber
660 //==================================================================
661 Standard_Integer GeomFill_LocationDraft::TraceNumber() const
663 if (Intersec == Standard_True) return 1;
667 //==================================================================
668 //Function:NbIntervals
670 //==================================================================
671 Standard_Integer GeomFill_LocationDraft::NbIntervals
672 (const GeomAbs_Shape S) const
674 return myLaw->NbIntervals(S);
677 //==================================================================
680 //==================================================================
681 void GeomFill_LocationDraft::Intervals(TColStd_Array1OfReal& T,
682 const GeomAbs_Shape S) const
684 myLaw->Intervals(T, S);
687 //==================================================================
688 //Function:SetInterval
690 //==================================================================
691 void GeomFill_LocationDraft::SetInterval(const Standard_Real First,
692 const Standard_Real Last)
694 myLaw->SetInterval( First, Last);
695 myTrimmed = myCurve->Trim( First, Last, 0);
697 //==================================================================
698 //Function: GetInterval
700 //==================================================================
701 void GeomFill_LocationDraft::GetInterval(Standard_Real& First,
702 Standard_Real& Last) const
704 First = myTrimmed->FirstParameter();
705 Last = myTrimmed->LastParameter();
708 //==================================================================
709 //Function: GetDomain
711 //==================================================================
712 void GeomFill_LocationDraft::GetDomain(Standard_Real& First,
713 Standard_Real& Last) const
715 First = myCurve->FirstParameter();
716 Last = myCurve->LastParameter();
719 //==================================================================
720 //function : Resolution
722 //==================================================================
723 void GeomFill_LocationDraft::Resolution (const Standard_Integer Index,
724 const Standard_Real Tol,
726 Standard_Real& TolV) const
730 TolU = mySurf->UResolution(Tol);
731 TolV = mySurf->VResolution(Tol);
739 //==================================================================
740 //Function:GetMaximalNorm
741 //Purpose : On suppose les triedres normes => return 1
742 //==================================================================
743 Standard_Real GeomFill_LocationDraft::GetMaximalNorm()
748 //==================================================================
749 //Function:GetAverageLaw
751 //==================================================================
752 void GeomFill_LocationDraft::GetAverageLaw(gp_Mat& AM,
756 Standard_Real U, delta;
759 myLaw->GetAverageLaw(V1, V2, V3);
760 AM.SetCols(V1.XYZ(), V2.XYZ(), V3.XYZ());
762 AV.SetCoord(0., 0., 0.);
763 delta = (myTrimmed->LastParameter() - myTrimmed->FirstParameter())/10;
764 U= myTrimmed->FirstParameter();
765 for (ii=0; ii<=10; ii++, U+=delta) {
766 V.SetXYZ( myTrimmed->Value(U).XYZ() );
772 //==================================================================
773 //Function : IsTranslation
775 //==================================================================
776 // Standard_Boolean GeomFill_LocationDraft::IsTranslation(Standard_Real& Error) const
777 Standard_Boolean GeomFill_LocationDraft::IsTranslation(Standard_Real& ) const
779 return myLaw->IsConstant();
782 //==================================================================
783 //Function : IsRotation
785 //==================================================================
786 Standard_Boolean GeomFill_LocationDraft::IsRotation(Standard_Real& Error) const
788 GeomAbs_CurveType Type;
790 Type = myCurve->GetType();
791 if (Type == GeomAbs_Circle) {
792 return myLaw->IsOnlyBy3dCurve();
794 return Standard_False;
797 //==================================================================
798 //Function : Rotation
800 //==================================================================
801 void GeomFill_LocationDraft::Rotation(gp_Pnt& Centre) const
803 Centre = myCurve->Circle().Location();
807 //==================================================================
808 //Function : IsIntersec
810 //==================================================================
811 Standard_Boolean GeomFill_LocationDraft::IsIntersec() const
817 //==================================================================
818 //Function : Direction
820 //==================================================================
821 gp_Dir GeomFill_LocationDraft::Direction() const