1 // Created on: 1992-01-20
2 // Created by: Remi GILET
3 // Copyright (c) 1992-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.
19 #include <Geom2dAdaptor_Curve.hxx>
20 #include <Geom2dGcc_CurveTool.hxx>
21 #include <Geom2dGcc_FunctionTanCuCuCu.hxx>
23 #include <gp_Circ2d.hxx>
24 #include <gp_Lin2d.hxx>
25 #include <gp_Pnt2d.hxx>
26 #include <gp_Vec2d.hxx>
27 #include <math_Matrix.hxx>
28 #include <Standard_ConstructionError.hxx>
30 void Geom2dGcc_FunctionTanCuCuCu::
31 InitDerivative(const math_Vector& X,
42 case Geom2dGcc_CiCuCu:
44 ElCLib::D2(X(1),Circ1,Point1,Tan1,D21);
45 Geom2dGcc_CurveTool::D2(Curv2,X(2),Point2,Tan2,D22);
46 Geom2dGcc_CurveTool::D2(Curv3,X(3),Point3,Tan3,D23);
49 case Geom2dGcc_CiCiCu:
51 ElCLib::D2(X(1),Circ1,Point1,Tan1,D21);
52 ElCLib::D2(X(2),Circ2,Point2,Tan2,D22);
53 Geom2dGcc_CurveTool::D2(Curv3,X(3),Point3,Tan3,D23);
56 case Geom2dGcc_CiLiCu:
58 ElCLib::D2(X(1),Circ1,Point1,Tan1,D21);
59 ElCLib::D1(X(2),Lin2,Point2,Tan2);
60 D22 = gp_Vec2d(0.,0.);
61 Geom2dGcc_CurveTool::D2(Curv3,X(3),Point3,Tan3,D23);
64 case Geom2dGcc_LiCuCu:
66 ElCLib::D1(X(1),Lin1,Point1,Tan1);
67 D21 = gp_Vec2d(0.,0.);
68 Geom2dGcc_CurveTool::D2(Curv2,X(2),Point2,Tan2,D22);
69 Geom2dGcc_CurveTool::D2(Curv3,X(3),Point3,Tan3,D23);
72 case Geom2dGcc_LiLiCu:
74 ElCLib::D1(X(1),Lin1,Point1,Tan1);
75 D21 = gp_Vec2d(0.,0.);
76 ElCLib::D1(X(2),Lin2,Point2,Tan2);
77 D22 = gp_Vec2d(0.,0.);
78 Geom2dGcc_CurveTool::D2(Curv3,X(3),Point3,Tan3,D23);
81 case Geom2dGcc_CuCuCu:
83 Geom2dGcc_CurveTool::D2(Curv1,X(1),Point1,Tan1,D21);
84 Geom2dGcc_CurveTool::D2(Curv2,X(2),Point2,Tan2,D22);
85 Geom2dGcc_CurveTool::D2(Curv3,X(3),Point3,Tan3,D23);
90 throw Standard_ConstructionError();
95 Geom2dGcc_FunctionTanCuCuCu::
96 Geom2dGcc_FunctionTanCuCuCu(const Geom2dAdaptor_Curve& C1 ,
97 const Geom2dAdaptor_Curve& C2 ,
98 const Geom2dAdaptor_Curve& C3 ) {
102 TheType = Geom2dGcc_CuCuCu;
105 Geom2dGcc_FunctionTanCuCuCu::
106 Geom2dGcc_FunctionTanCuCuCu(const gp_Circ2d& C1 ,
107 const Geom2dAdaptor_Curve& C2 ,
108 const Geom2dAdaptor_Curve& C3 ) {
112 TheType = Geom2dGcc_CiCuCu;
115 Geom2dGcc_FunctionTanCuCuCu::
116 Geom2dGcc_FunctionTanCuCuCu(const gp_Circ2d& C1 ,
117 const gp_Circ2d& C2 ,
118 const Geom2dAdaptor_Curve& C3 ) {
122 TheType = Geom2dGcc_CiCiCu;
125 Geom2dGcc_FunctionTanCuCuCu::
126 Geom2dGcc_FunctionTanCuCuCu(const gp_Circ2d& C1 ,
128 const Geom2dAdaptor_Curve& C3 ) {
132 TheType = Geom2dGcc_CiLiCu;
135 Geom2dGcc_FunctionTanCuCuCu::
136 Geom2dGcc_FunctionTanCuCuCu(const gp_Lin2d& L1 ,
138 const Geom2dAdaptor_Curve& C3 ) {
142 TheType = Geom2dGcc_LiLiCu;
145 Geom2dGcc_FunctionTanCuCuCu::
146 Geom2dGcc_FunctionTanCuCuCu(const gp_Lin2d& L1 ,
147 const Geom2dAdaptor_Curve& C2 ,
148 const Geom2dAdaptor_Curve& C3 ) {
152 TheType = Geom2dGcc_LiCuCu;
155 //==========================================================================
157 Standard_Integer Geom2dGcc_FunctionTanCuCuCu::
158 NbVariables () const { return 3; }
160 Standard_Integer Geom2dGcc_FunctionTanCuCuCu::
161 NbEquations () const { return 3; }
163 Standard_Boolean Geom2dGcc_FunctionTanCuCuCu::
164 Value (const math_Vector& X ,
165 math_Vector& Fval ) {
175 InitDerivative(X,Point1,Point2,Point3,Tan1,Tan2,Tan3,D21,D22,D23);
176 //pipj (normes) et PiPj (non Normes).
177 gp_XY P1P2(gp_Vec2d(Point1,Point2).XY());
178 gp_XY P2P3(gp_Vec2d(Point2,Point3).XY());
179 gp_XY P3P1(gp_Vec2d(Point3,Point1).XY());
180 Standard_Real NorP1P2 = P1P2.Modulus();
181 Standard_Real NorP2P3 = P2P3.Modulus();
182 Standard_Real NorP3P1 = P3P1.Modulus();
183 gp_XY p1p2,p2p3,p3p1;
184 if (NorP1P2 >= gp::Resolution()) { p1p2 = P1P2/NorP1P2; }
185 else { p1p2 = gp_XY(0.,0.); }
186 if (NorP2P3 >= gp::Resolution()) { p2p3 = P2P3/NorP2P3; }
187 else { p2p3 = gp_XY(0.,0.); }
188 if (NorP3P1 >= gp::Resolution()) { p3p1 = P3P1/NorP3P1; }
189 else { p3p1 = gp_XY(0.,0.); }
190 //derivees premieres non normees Deriv1ui.
191 gp_XY Deriv1u1(Tan1.XY());
192 gp_XY Deriv1u2(Tan2.XY());
193 gp_XY Deriv1u3(Tan3.XY());
194 //normales aux courbes.
195 Standard_Real nnor1 = Deriv1u1.Modulus();
196 Standard_Real nnor2 = Deriv1u2.Modulus();
197 Standard_Real nnor3 = Deriv1u3.Modulus();
198 gp_XY Nor1(-Deriv1u1.Y(),Deriv1u1.X());
199 gp_XY Nor2(-Deriv1u2.Y(),Deriv1u2.X());
200 gp_XY Nor3(-Deriv1u3.Y(),Deriv1u3.X());
201 gp_XY nor1,nor2,nor3;
202 if (nnor1 >= gp::Resolution()) { nor1 = Nor1/nnor1; }
203 else { nor1 = gp_XY(0.,0.); }
204 if (nnor2 >= gp::Resolution()) { nor2 = Nor2/nnor2; }
205 else { nor2 = gp_XY(0.,0.); }
206 if (nnor3 >= gp::Resolution()) { nor3 = Nor3/nnor3; }
207 else { nor3 = gp_XY(0.,0.); }
208 //determination des signes pour les produits scalaires.
209 Standard_Real signe1 = 1.;
210 Standard_Real signe2 = 1.;
211 Standard_Real signe3 = 1.;
212 gp_Pnt2d Pcenter(gp_XY(Point1.XY()+Point2.XY()+Point3.XY())/3.);
213 gp_XY fic1(Pcenter.XY()-Point1.XY());
214 gp_XY fic2(Pcenter.XY()-Point2.XY());
215 gp_XY fic3(Pcenter.XY()-Point3.XY());
216 Standard_Real pscal11 = nor1.Dot(fic1);
217 Standard_Real pscal22 = nor2.Dot(fic2);
218 Standard_Real pscal33 = nor3.Dot(fic3);
219 if (pscal11 <= 0.) { signe1 = -1; }
220 if (pscal22 <= 0.) { signe2 = -1; }
221 if (pscal33 <= 0.) { signe3 = -1; }
224 Fval(1) = signe1*nor1.Dot(p1p2)+signe2*nor2.Dot(p1p2);
225 Fval(2) = signe2*nor2.Dot(p2p3)+signe3*nor3.Dot(p2p3);
226 Fval(3) = signe3*nor3.Dot(p3p1)+signe1*nor1.Dot(p3p1);
227 return Standard_True;
232 Standard_Boolean Geom2dGcc_FunctionTanCuCuCu::Derivatives (const math_Vector&,
245 InitDerivative(X,Point1,Point2,Point3,Tan1,Tan2,Tan3,D21,D22,D23);
246 //derivees premieres non normees Deriv1ui.
247 gp_XY Deriv1u1(Tan1.XY());
248 gp_XY Deriv1u2(Tan2.XY());
249 gp_XY Deriv1u3(Tan3.XY());
250 //pipj (normes) et PiPj (non Normes).
251 gp_XY P1P2(gp_Vec2d(Point1,Point2).XY());
252 gp_XY P2P3(gp_Vec2d(Point2,Point3).XY());
253 gp_XY P3P1(gp_Vec2d(Point3,Point1).XY());
254 Standard_Real NorP1P2 = P1P2.Modulus();
255 Standard_Real NorP2P3 = P2P3.Modulus();
256 Standard_Real NorP3P1 = P3P1.Modulus();
257 gp_XY p1p2,p2p3,p3p1;
258 if (NorP1P2 >= gp::Resolution()) { p1p2 = P1P2/NorP1P2; }
259 else { p1p2 = gp_XY(0.,0.); }
260 if (NorP2P3 >= gp::Resolution()) { p2p3 = P2P3/NorP2P3; }
261 else { p2p3 = gp_XY(0.,0.); }
262 if (NorP3P1 >= gp::Resolution()) { p3p1 = P3P1/NorP3P1; }
263 else { p3p1 = gp_XY(0.,0.); }
264 //normales au courbes normees Nori et non nromees nori et norme des nori.
265 Standard_Real nnor1 = Deriv1u1.Modulus();
266 Standard_Real nnor2 = Deriv1u2.Modulus();
267 Standard_Real nnor3 = Deriv1u3.Modulus();
268 gp_XY Nor1(-Deriv1u1.Y(),Deriv1u1.X());
269 gp_XY Nor2(-Deriv1u2.Y(),Deriv1u2.X());
270 gp_XY Nor3(-Deriv1u3.Y(),Deriv1u3.X());
271 gp_XY nor1,nor2,nor3;
272 if (nnor1 >= gp::Resolution()) { nor1 = Nor1/nnor1; }
273 else { nor1 = gp_XY(0.,0.); }
274 if (nnor2 >= gp::Resolution()) { nor2 = Nor2/nnor2; }
275 else { nor2 = gp_XY(0.,0.); }
276 if (nnor3 >= gp::Resolution()) { nor3 = Nor3/nnor3; }
277 else { nor3 = gp_XY(0.,0.); }
278 //derivees des normales.
279 gp_XY NorD21,NorD22,NorD23;
280 NorD21 = gp_XY(-D21.Y(),D21.X());
281 NorD22 = gp_XY(-D22.Y(),D22.X());
282 NorD23 = gp_XY(-D23.Y(),D23.X());
283 //determination des signes pour les produits scalaires.
284 Standard_Real signe1 = 1.;
285 Standard_Real signe2 = 1.;
286 Standard_Real signe3 = 1.;
287 gp_XY P = Point1.XY();
292 gp_XY fic1 = Pcenter.XY();
294 gp_XY fic2 = Pcenter.XY();
296 gp_XY fic3 = Pcenter.XY();
298 Standard_Real pscal11 = nor1.Dot(fic1);
299 Standard_Real pscal22 = nor2.Dot(fic2);
300 Standard_Real pscal33 = nor3.Dot(fic3);
301 if (pscal11 <= 0.) { signe1 = -1; }
302 if (pscal22 <= 0.) { signe2 = -1; }
303 if (pscal33 <= 0.) { signe3 = -1; }
305 // Derivees dFui/uj 1 <= ui <= 3 , 1 <= uj <= 3
306 // =============================================
307 Standard_Real partie1,partie2;
308 if (nnor1 <= gp::Resolution()) { partie1 = 0.; }
309 else { partie1 = (signe1*NorD21/nnor1-(Nor1.Dot(NorD21)/(nnor1*nnor1*nnor1))
311 if (NorP1P2 <= gp::Resolution()) { partie2 = 0.; }
312 else {partie2=((Deriv1u1.Dot(p1p2)/(NorP1P2*NorP1P2))*P1P2-Deriv1u1/NorP1P2)
313 .Dot(signe1*nor1+signe2*nor2); }
314 Deriv(1,1) = partie1 + partie2;
315 if (nnor2 <= gp::Resolution()) { partie1 = 0.; }
316 else { partie1=(signe2*NorD22/(nnor2)-(Nor2.Dot(NorD22)/(nnor2*nnor2*nnor2))
318 if (NorP1P2 <= gp::Resolution()) { partie2 = 0.; }
319 else{partie2=((-Deriv1u2.Dot(p1p2)/(NorP1P2*NorP1P2))*P1P2+Deriv1u2/NorP1P2)
320 .Dot(signe1*nor1+signe2*nor2); }
321 Deriv(1,2) = partie1 + partie2;
324 if (nnor2 <= gp::Resolution()) { partie1 = 0.; }
325 else { partie1=(signe2*NorD22/(nnor2)-(Nor2.Dot(NorD22)/(nnor2*nnor2*nnor2))
327 if (NorP2P3 <= gp::Resolution()) { partie2 = 0.; }
328 else { partie2=((Deriv1u2.Dot(p2p3)/(NorP2P3*NorP2P3))*P2P3-Deriv1u2/NorP2P3)
329 .Dot(signe2*nor2+signe3*nor3); }
330 Deriv(2,2) = partie1 +partie2;
331 if (nnor3 <= gp::Resolution()) { partie1 = 0.; }
332 else { partie1=(signe3*NorD23/(nnor3)-(Nor3.Dot(NorD23)/(nnor3*nnor3*nnor3))
334 if (NorP2P3 <= gp::Resolution()) { partie2 = 0.; }
335 else {partie2=((-Deriv1u3.Dot(p2p3)/(NorP2P3*NorP2P3))*P2P3+Deriv1u3/NorP2P3)
336 .Dot(signe2*nor2+signe3*nor3); }
337 Deriv(2,3) = partie1 + partie2;
338 if (nnor1 <= gp::Resolution()) { partie1 = 0.; }
339 else { partie1 =(signe1*NorD21/(nnor1)-(Nor1.Dot(NorD21)/(nnor1*nnor1*nnor1))
341 if (NorP3P1 <= gp::Resolution()) { partie2 = 0.; }
342 else {partie2=((-Deriv1u1.Dot(p3p1)/(NorP3P1*NorP3P1))*P3P1+Deriv1u1/NorP3P1)
343 .Dot(signe1*nor1+signe3*nor3); }
344 Deriv(3,1) = partie1 + partie2;
346 if (nnor3 <= gp::Resolution()) { partie1 = 0.; }
347 else { partie1=(signe3*NorD23/(nnor3)-(Nor3.Dot(NorD23)/(nnor3*nnor3*nnor3))
349 if (NorP3P1 <= gp::Resolution()) { partie2 = 0.; }
350 else {partie2=((Deriv1u3.Dot(p3p1)/(NorP3P1*NorP3P1))*P3P1-Deriv1u3/NorP3P1)
351 .Dot(signe1*nor1+signe3*nor3); }
352 Deriv(3,3) = partie1+partie2;
354 return Standard_True;
358 Standard_Boolean Geom2dGcc_FunctionTanCuCuCu:: Values (const math_Vector&,
372 InitDerivative(X,Point1,Point2,Point3,Tan1,Tan2,Tan3,D21,D22,D23);
373 //derivees premieres non normees Deriv1ui.
374 gp_XY Deriv1u1(Tan1.XY());
375 gp_XY Deriv1u2(Tan2.XY());
376 gp_XY Deriv1u3(Tan3.XY());
377 //pipj (normes) et PiPj (non Normes).
378 gp_XY P1P2(gp_Vec2d(Point1,Point2).XY());
379 gp_XY P2P3(gp_Vec2d(Point2,Point3).XY());
380 gp_XY P3P1(gp_Vec2d(Point3,Point1).XY());
381 Standard_Real NorP1P2 = P1P2.Modulus();
382 Standard_Real NorP2P3 = P2P3.Modulus();
383 Standard_Real NorP3P1 = P3P1.Modulus();
384 gp_XY p1p2,p2p3,p3p1;
385 if (NorP1P2 >= gp::Resolution()) { p1p2 = P1P2/NorP1P2; }
386 else { p1p2 = gp_XY(0.,0.); }
387 if (NorP2P3 >= gp::Resolution()) { p2p3 = P2P3/NorP2P3; }
388 else { p2p3 = gp_XY(0.,0.); }
389 if (NorP3P1 >= gp::Resolution()) { p3p1 = P3P1/NorP3P1; }
390 else { p3p1 = gp_XY(0.,0.); }
391 //normales au courbes normees Nori et non nromees nori et norme des nori.
392 Standard_Real nnor1 = Deriv1u1.Modulus();
393 Standard_Real nnor2 = Deriv1u2.Modulus();
394 Standard_Real nnor3 = Deriv1u3.Modulus();
395 gp_XY Nor1(-Deriv1u1.Y(),Deriv1u1.X());
396 gp_XY Nor2(-Deriv1u2.Y(),Deriv1u2.X());
397 gp_XY Nor3(-Deriv1u3.Y(),Deriv1u3.X());
398 gp_XY nor1,nor2,nor3;
399 if (nnor1 >= gp::Resolution()) { nor1 = Nor1/nnor1; }
400 else { nor1 = gp_XY(0.,0.); }
401 if (nnor2 >= gp::Resolution()) { nor2 = Nor2/nnor2; }
402 else { nor2 = gp_XY(0.,0.); }
403 if (nnor3 >= gp::Resolution()) { nor3 = Nor3/nnor3; }
404 else { nor3 = gp_XY(0.,0.); }
405 //derivees des normales.
406 gp_XY NorD21,NorD22,NorD23;
407 NorD21 = gp_XY(-D21.Y(),D21.X());
408 NorD22 = gp_XY(-D22.Y(),D22.X());
409 NorD23 = gp_XY(-D23.Y(),D23.X());
410 //determination des signes pour les produits scalaires.
411 Standard_Real signe1 = 1.;
412 Standard_Real signe2 = 1.;
413 Standard_Real signe3 = 1.;
414 gp_XY P = Point1.XY();
419 gp_XY fic1 = Pcenter.XY();
421 gp_XY fic2 = Pcenter.XY();
423 gp_XY fic3 = Pcenter.XY();
425 Standard_Real pscal11 = nor1.Dot(fic1);
426 Standard_Real pscal22 = nor2.Dot(fic2);
427 Standard_Real pscal33 = nor3.Dot(fic3);
428 if (pscal11 <= 0.) { signe1 = -1; }
429 if (pscal22 <= 0.) { signe2 = -1; }
430 if (pscal33 <= 0.) { signe3 = -1; }
434 Fval(1) = signe1*nor1.Dot(p1p2)+signe2*nor2.Dot(p1p2);
435 Fval(2) = signe2*nor2.Dot(p2p3)+signe3*nor3.Dot(p2p3);
436 Fval(3) = signe3*nor3.Dot(p3p1)+signe1*nor1.Dot(p3p1);
437 // Derivees dFui/uj 1 <= ui <= 3 , 1 <= uj <= 3
438 // =============================================
439 Standard_Real partie1,partie2;
440 if (nnor1 <= gp::Resolution()) { partie1 = 0.; }
441 else { partie1 = signe1*(NorD21/nnor1-(Nor1.Dot(NorD21)/(nnor1*nnor1*nnor1))
443 if (NorP1P2 <= gp::Resolution()) { partie2 = 0.; }
444 else {partie2=((Deriv1u1.Dot(p1p2)/(NorP1P2*NorP1P2))*P1P2-Deriv1u1/NorP1P2)
445 .Dot(signe1*nor1+signe2*nor2); }
446 Deriv(1,1) = partie1 + partie2;
447 if (nnor2 <= gp::Resolution()) { partie1 = 0.; }
448 else { partie1=signe2*(NorD22/(nnor2)-(Nor2.Dot(NorD22)/(nnor2*nnor2*nnor2))
450 if (NorP1P2 <= gp::Resolution()) { partie2 = 0.; }
451 else{partie2=((-Deriv1u2.Dot(p1p2)/(NorP1P2*NorP1P2))*P1P2+Deriv1u2/NorP1P2)
452 .Dot(signe1*nor1+signe2*nor2); }
453 Deriv(1,2) = partie1 + partie2;
456 if (nnor2 <= gp::Resolution()) { partie1 = 0.; }
457 else { partie1=signe2*(NorD22/(nnor2)-(Nor2.Dot(NorD22)/(nnor2*nnor2*nnor2))
459 if (NorP2P3 <= gp::Resolution()) { partie2 = 0.; }
460 else { partie2=((Deriv1u2.Dot(p2p3)/(NorP2P3*NorP2P3))*P2P3-Deriv1u2/NorP2P3)
461 .Dot(signe2*nor2+signe3*nor3); }
462 Deriv(2,2) = partie1 +partie2;
463 if (nnor3 <= gp::Resolution()) { partie1 = 0.; }
464 else { partie1=signe3*(NorD23/(nnor3)-(Nor3.Dot(NorD23)/(nnor3*nnor3*nnor3))
466 if (NorP2P3 <= gp::Resolution()) { partie2 = 0.; }
467 else {partie2=((-Deriv1u3.Dot(p2p3)/(NorP2P3*NorP2P3))*P2P3+Deriv1u3/NorP2P3)
468 .Dot(signe2*nor2+signe3*nor3); }
469 Deriv(2,3) = partie1 + partie2;
470 if (nnor1 <= gp::Resolution()) { partie1 = 0.; }
471 else { partie1 =signe1*(NorD21/(nnor1)-(Nor1.Dot(NorD21)/(nnor1*nnor1*nnor1))
473 if (NorP3P1 <= gp::Resolution()) { partie2 = 0.; }
474 else {partie2=((-Deriv1u1.Dot(p3p1)/(NorP3P1*NorP3P1))*P3P1+Deriv1u1/NorP3P1)
475 .Dot(signe1*nor1+signe3*nor3); }
476 Deriv(3,1) = partie1 + partie2;
478 if (nnor3 <= gp::Resolution()) { partie1 = 0.; }
479 else { partie1=signe3*(NorD23/(nnor3)-(Nor3.Dot(NorD23)/(nnor3*nnor3*nnor3))
481 if (NorP3P1 <= gp::Resolution()) { partie2 = 0.; }
482 else {partie2=((Deriv1u3.Dot(p3p1)/(NorP3P1*NorP3P1))*P3P1-Deriv1u3/NorP3P1)
483 .Dot(signe1*nor1+signe3*nor3); }
484 Deriv(3,3) = partie1+partie2;
486 return Standard_True;