1 // Created by: Julia GERASIMOVA
2 // Copyright (c) 2015 OPEN CASCADE SAS
3 //
4 // This file is part of Open CASCADE Technology software library.
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
11 //
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
14 
15 
16 #include <Adaptor2d_Curve2d.hxx>
17 #include <Adaptor3d_Curve.hxx>
18 #include <Adaptor3d_Surface.hxx>
19 #include <BlendFunc.hxx>
20 #include <BlendFunc_ConstThroatInv.hxx>
21 #include <math_Matrix.hxx>
22 #include <Precision.hxx>
23 
24 //=======================================================================
25 //function : BlendFunc_ConstThroatInv
26 //purpose  :
27 //=======================================================================
28 
BlendFunc_ConstThroatInv(const Handle (Adaptor3d_Surface)& S1,const Handle (Adaptor3d_Surface)& S2,const Handle (Adaptor3d_Curve)& C)29 BlendFunc_ConstThroatInv::BlendFunc_ConstThroatInv(const Handle(Adaptor3d_Surface)& S1,
30                                                    const Handle(Adaptor3d_Surface)& S2,
31                                                    const Handle(Adaptor3d_Curve)&   C)
32 : BlendFunc_GenChamfInv(S1,S2,C),
33   Throat(0.0),
34   param(0.0),
35   sign1(0.0),
36   sign2(0.0),
37   normtg(0.0),
38   theD(0.0)
39 {
40 }
41 
42 
43 //=======================================================================
44 //function : Set
45 //purpose  :
46 //=======================================================================
47 
Set(const Standard_Real theThroat,const Standard_Real,const Standard_Integer Choix)48 void BlendFunc_ConstThroatInv::Set(const Standard_Real    theThroat,
49                                    const Standard_Real,
50                                    const Standard_Integer Choix)
51 {
52   //Standard_Real dis1,dis2;
53 
54   Throat = theThroat;
55 
56   choix = Choix;
57   switch (choix) {
58   case 1:
59   case 2:
60     {
61       sign1 = -1;
62       sign2 = -1;
63     }
64     break;
65   case 3:
66   case 4:
67     {
68       sign1 = 1;
69       sign2 = -1;
70     }
71     break;
72   case 5:
73   case 6:
74     {
75       sign1 = 1;
76       sign2 = 1;
77     }
78     break;
79   case 7:
80   case 8:
81     {
82       sign1 = -1;
83       sign2 = 1;
84     }
85     break;
86   default:
87     sign1 = -1;
88     sign2 = -1;
89   }
90 }
91 
92 //=======================================================================
93 //function : IsSolution
94 //purpose  :
95 //=======================================================================
96 
IsSolution(const math_Vector & Sol,const Standard_Real Tol)97 Standard_Boolean BlendFunc_ConstThroatInv::IsSolution(const math_Vector& Sol, const Standard_Real Tol)
98 {
99   math_Vector valsol(1,4);
100   Value(Sol, valsol);
101 
102   if (Abs(valsol(1)) <= Tol &&
103       Abs(valsol(2)) <= Tol &&
104       Abs(valsol(3)) <= Tol*Tol &&
105       Abs(valsol(4)) <= Tol*Tol)
106     return Standard_True;
107 
108   return Standard_False;
109 }
110 
111 //=======================================================================
112 //function : Value
113 //purpose  :
114 //=======================================================================
115 
Value(const math_Vector & X,math_Vector & F)116 Standard_Boolean BlendFunc_ConstThroatInv::Value(const math_Vector& X, math_Vector& F)
117 {
118   gp_Pnt2d p2d;
119   gp_Vec2d v2d;
120   csurf->D1(X(1),p2d,v2d);
121   param = X(2);
122   curv->D2(param,ptgui,d1gui,d2gui);
123   normtg = d1gui.Magnitude();
124   nplan  = d1gui.Normalized();
125   theD = - (nplan.XYZ().Dot(ptgui.XYZ()));
126 
127   math_Vector XX(1,4);
128 
129   if(first){
130     XX(1) = p2d.X(); XX(2) = p2d.Y();
131     XX(3) = X(3); XX(4) = X(4);
132   }
133 
134   else{
135     XX(1) = X(3); XX(2) = X(4);
136     XX(3) = p2d.X(); XX(4) = p2d.Y();
137   }
138 
139   surf1->D0( XX(1), XX(2), pts1 );
140   surf2->D0( XX(3), XX(4), pts2 );
141 
142   F(1) = nplan.XYZ().Dot(pts1.XYZ()) + theD;
143   F(2) = nplan.XYZ().Dot(pts2.XYZ()) + theD;
144 
145   const gp_Pnt ptmid((pts1.XYZ() + pts2.XYZ())/2);
146   const gp_Vec vmid(ptgui, ptmid);
147 
148   F(3) = vmid.SquareMagnitude() - Throat*Throat;
149 
150   const gp_Vec vref1(ptgui, pts1);
151   const gp_Vec vref2(ptgui, pts2);
152 
153   F(4) = vref1.SquareMagnitude() - vref2.SquareMagnitude();
154 
155   return Standard_True;
156 }
157 
158 //=======================================================================
159 //function : Derivatives
160 //purpose  :
161 //=======================================================================
162 
Derivatives(const math_Vector & X,math_Matrix & D)163 Standard_Boolean BlendFunc_ConstThroatInv::Derivatives(const math_Vector& X, math_Matrix& D)
164 {
165   //Standard_Integer i, j;
166   gp_Pnt2d p2d;
167   gp_Vec2d v2d; //, df1, df2;
168   //gp_Pnt pts, ptgui;
169   gp_Vec dnplan, temp, temp1, temp2, tempmid; //, d1u, d1v, nplan;
170   math_Vector XX(1,4); //x1(1,2), x2(1,2);
171   //math_Matrix d1(1,2,1,2), d2(1,2,1,2);
172 
173   csurf->D1(X(1), p2d, v2d);
174   //corde1.SetParam(X(2));
175   //corde2.SetParam(X(2));
176   param = X(2);
177   curv->D2(param,ptgui,d1gui,d2gui);
178   normtg = d1gui.Magnitude();
179   nplan  = d1gui.Normalized();
180   theD = - (nplan.XYZ().Dot(ptgui.XYZ()));
181 
182   dnplan.SetLinearForm(1./normtg,d2gui,
183                        -1./normtg*(nplan.Dot(d2gui)),nplan);
184 
185   temp1.SetXYZ(pts1.XYZ() - ptgui.XYZ());
186   temp2.SetXYZ(pts2.XYZ() - ptgui.XYZ());
187   tempmid.SetXYZ((pts1.XYZ() + pts2.XYZ())/2 - ptgui.XYZ());
188 
189   //x1(1) = p2d.X(); x1(2) = p2d.Y();
190   //x2(1) = X(3); x2(2) = X(4);
191   if (first)
192   {
193     XX(1) = p2d.X(); XX(2) = p2d.Y();
194     XX(3) = X(3); XX(4) = X(4);
195   }
196   else
197   {
198     XX(1) = X(3); XX(2) = X(4);
199     XX(3) = p2d.X(); XX(4) = p2d.Y();
200   }
201 
202   surf1->D1(XX(1), XX(2), pts1, d1u1, d1v1);
203   surf2->D1(XX(3), XX(4), pts2, d1u2, d1v2);
204 
205   if( first ){
206     // p2d = pts est sur surf1
207     //ptgui = corde1.PointOnGuide();
208     //nplan = corde1.NPlan();
209     temp.SetLinearForm(v2d.X(),d1u1, v2d.Y(),d1v1);
210 
211     D(1,1) = nplan.Dot(temp);
212     D(2,1) = 0.;
213     D(3,1) = gp_Vec(ptgui,pts1).Dot(temp);
214     D(4,1) = 2*(gp_Vec(ptgui,pts1).Dot(temp));
215 
216     D(1,3) = 0.;
217     D(1,4) = 0.;
218     D(2,3) = nplan.Dot(d1u2);
219     D(2,4) = nplan.Dot(d1v2);
220     D(3,3) = gp_Vec((pts1.XYZ() + pts2.XYZ())/2 - ptgui.XYZ()).Dot(d1u2);
221     D(3,4) = gp_Vec((pts1.XYZ() + pts2.XYZ())/2 - ptgui.XYZ()).Dot(d1v2);
222     D(4,3) = -2.*gp_Vec(ptgui,pts2).Dot(d1u2);
223     D(4,4) = -2.*gp_Vec(ptgui,pts2).Dot(d1v2);
224 
225     //surf1->D1(x1(1),x1(2),pts,d1u,d1v);
226   }
227   else{
228     //  p2d = pts est sur surf2
229     //ptgui = corde2.PointOnGuide();
230     //nplan = corde2.NPlan();
231     temp.SetLinearForm(v2d.X(),d1u2, v2d.Y(),d1v2);
232 
233     D(1,1) = 0.;
234     D(2,1) = nplan.Dot(temp);
235     D(3,1) = gp_Vec(ptgui,pts2).Dot(temp);
236     D(4,1) = -2*(gp_Vec(ptgui,pts2).Dot(temp));
237 
238     D(1,3) = nplan.Dot(d1u1);
239     D(1,4) = nplan.Dot(d1v1);
240     D(2,3) = 0.;
241     D(2,4) = 0.;
242     D(3,3) = gp_Vec((pts1.XYZ() + pts2.XYZ())/2 - ptgui.XYZ()).Dot(d1u1);
243     D(3,4) = gp_Vec((pts1.XYZ() + pts2.XYZ())/2 - ptgui.XYZ()).Dot(d1v1);
244     D(4,3) = 2.*gp_Vec(ptgui,pts1).Dot(d1u1);
245     D(4,4) = 2.*gp_Vec(ptgui,pts1).Dot(d1v1);
246 
247     //surf2->D1(x1(1),x1(2),pts,d1u,d1v);
248   }
249 
250   D(1,2) = dnplan.Dot(temp1) - nplan.Dot(d1gui);
251   D(2,2) = dnplan.Dot(temp2) - nplan.Dot(d1gui);
252   D(3,2) = -2.*d1gui.Dot(tempmid);
253   D(4,2) = 2.*d1gui.Dot(temp1) - 2.*d1gui.Dot(temp2);
254 
255   return Standard_True;
256 }
257