1 // Copyright (c) 1995-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
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.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
19 #include <gp_Cone.hxx>
20 #include <gp_Cylinder.hxx>
23 #include <gp_Sphere.hxx>
24 #include <gp_Torus.hxx>
26 #include <IntSurf_Quadric.hxx>
27 #include <StdFail_NotDone.hxx>
29 // ============================================================
30 IntSurf_Quadric::IntSurf_Quadric ():typ(GeomAbs_OtherSurface),
31 prm1(0.), prm2(0.), prm3(0.), prm4(0.)
33 // ============================================================
34 IntSurf_Quadric::IntSurf_Quadric (const gp_Pln& P):
35 ax3(P.Position()),typ(GeomAbs_Plane)
37 ax3direc = ax3.Direct();
38 P.Coefficients(prm1,prm2,prm3,prm4);
40 // ============================================================
41 IntSurf_Quadric::IntSurf_Quadric (const gp_Cylinder& C):
43 ax3(C.Position()),lin(ax3.Axis()),typ(GeomAbs_Cylinder)
46 ax3direc=ax3.Direct();
49 // ============================================================
50 IntSurf_Quadric::IntSurf_Quadric (const gp_Sphere& S):
52 ax3(S.Position()),lin(ax3.Axis()),typ(GeomAbs_Sphere)
55 ax3direc = ax3.Direct();
58 // ============================================================
59 IntSurf_Quadric::IntSurf_Quadric (const gp_Cone& C):
61 ax3(C.Position()),typ(GeomAbs_Cone)
63 ax3direc = ax3.Direct();
64 lin.SetPosition(ax3.Axis());
70 // ============================================================
71 IntSurf_Quadric::IntSurf_Quadric (const gp_Torus& T):
73 ax3(T.Position()),typ(GeomAbs_Torus)
75 ax3direc = ax3.Direct();
76 lin.SetPosition(ax3.Axis());
77 prm1 = T.MajorRadius();
78 prm2 = T.MinorRadius();
82 // ============================================================
83 void IntSurf_Quadric::SetValue (const gp_Pln& P)
87 ax3direc = ax3.Direct();
88 P.Coefficients(prm1,prm2,prm3,prm4);
90 // ============================================================
91 void IntSurf_Quadric::SetValue (const gp_Cylinder& C)
93 typ = GeomAbs_Cylinder;
95 ax3direc = ax3.Direct();
96 lin.SetPosition(ax3.Axis());
100 // ============================================================
101 void IntSurf_Quadric::SetValue (const gp_Sphere& S)
103 typ = GeomAbs_Sphere;
105 ax3direc = ax3.Direct();
106 lin.SetPosition(ax3.Axis());
110 // ============================================================
111 void IntSurf_Quadric::SetValue (const gp_Cone& C)
115 ax3direc = ax3.Direct();
116 lin.SetPosition(ax3.Axis());
117 prm1 = C.RefRadius();
118 prm2 = C.SemiAngle();
122 // ============================================================
123 void IntSurf_Quadric::SetValue (const gp_Torus& T)
127 ax3direc = ax3.Direct();
128 lin.SetPosition(ax3.Axis());
129 prm1 = T.MajorRadius();
130 prm2 = T.MinorRadius();
134 // ============================================================
135 Standard_Real IntSurf_Quadric::Distance (const gp_Pnt& P) const {
137 case GeomAbs_Plane: // plan
138 return prm1*P.X() + prm2*P.Y() + prm3*P.Z() + prm4;
139 case GeomAbs_Cylinder: // cylindre
140 return (lin.Distance(P) - prm1);
141 case GeomAbs_Sphere: // sphere
142 return (lin.Location().Distance(P) - prm1);
143 case GeomAbs_Cone: // cone
145 Standard_Real dist = lin.Distance(P);
147 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
148 gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
149 Standard_Real distp = lin.Distance(Pp);
150 dist = (dist-distp)/prm3;
153 case GeomAbs_Torus: // torus
158 gp_Vec OZ (ax3.Direction());
159 Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
161 gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
162 ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
163 PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
165 Standard_Real dist = P.Distance(PT) - prm2;
175 // ============================================================
176 gp_Vec IntSurf_Quadric::Gradient (const gp_Pnt& P) const {
179 case GeomAbs_Plane: // plan
180 grad.SetCoord(prm1,prm2,prm3);
182 case GeomAbs_Cylinder: // cylindre
184 gp_XYZ PP(lin.Location().XYZ());
185 PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
186 grad.SetXYZ(P.XYZ()-PP);
187 Standard_Real N = grad.Magnitude();
188 if(N>1e-14) { grad.Divide(N); }
189 else { grad.SetCoord(0.0,0.0,0.0); }
192 case GeomAbs_Sphere: // sphere
195 grad.SetXYZ((PP-lin.Location().XYZ()));
196 Standard_Real N = grad.Magnitude();
197 if(N>1e-14) { grad.Divide(N); }
198 else { grad.SetCoord(0.0,0.0,0.0); }
201 case GeomAbs_Cone: // cone
204 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
205 gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
207 ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
208 grad=D1u.Crossed(D1v);
209 if(ax3direc==Standard_False) {
215 case GeomAbs_Torus: // torus
220 gp_Vec OZ (ax3.Direction());
221 Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
223 gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
224 ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
225 PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
227 grad.SetXYZ(P.XYZ() - PT.XYZ());
228 Standard_Real N = grad.Magnitude();
229 if(N>1e-14) { grad.Divide(N); }
230 else { grad.SetCoord(0., 0., 0.); }
239 // ============================================================
240 void IntSurf_Quadric::ValAndGrad (const gp_Pnt& P,
248 Dist = prm1*P.X() + prm2*P.Y() + prm3*P.Z() + prm4;
249 Grad.SetCoord(prm1,prm2,prm3);
252 case GeomAbs_Cylinder:
254 Dist = lin.Distance(P) - prm1;
255 gp_XYZ PP(lin.Location().XYZ());
256 PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
257 Grad.SetXYZ((P.XYZ()-PP));
258 Standard_Real N = Grad.Magnitude();
259 if(N>1e-14) { Grad.Divide(N); }
260 else { Grad.SetCoord(0.0,0.0,0.0); }
265 Dist = lin.Location().Distance(P) - prm1;
267 Grad.SetXYZ((PP-lin.Location().XYZ()));
268 Standard_Real N = Grad.Magnitude();
269 if(N>1e-14) { Grad.Divide(N); }
270 else { Grad.SetCoord(0.0,0.0,0.0); }
275 Standard_Real dist = lin.Distance(P);
279 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
280 ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
281 Standard_Real distp = lin.Distance(Pp);
282 dist = (dist-distp)/prm3;
284 Grad=D1u.Crossed(D1v);
285 if(ax3direc==Standard_False) {
288 //-- lbr le 7 mars 96
289 //-- Si le gardient est nul, on est sur l axe
290 //-- et dans ce cas dist vaut 0
291 //-- On peut donc renvoyer une valeur quelconque.
294 || Grad.Z() > 1e-13) {
304 gp_Vec OZ (ax3.Direction());
305 Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
307 gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
308 ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
309 PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
311 Dist = P.Distance(PT) - prm2;
313 Grad.SetXYZ(P.XYZ()-PT.XYZ());
314 Standard_Real N = Grad.Magnitude();
315 if(N>1e-14) { Grad.Divide(N); }
316 else { Grad.SetCoord(0., 0., 0.); }
324 // ============================================================
325 gp_Pnt IntSurf_Quadric::Value(const Standard_Real U,
326 const Standard_Real V) const
331 return ElSLib::PlaneValue(U,V,ax3);
332 case GeomAbs_Cylinder:
333 return ElSLib::CylinderValue(U,V,ax3,prm1);
335 return ElSLib::SphereValue(U,V,ax3,prm1);
337 return ElSLib::ConeValue(U,V,ax3,prm1,prm2);
339 return ElSLib::TorusValue(U,V,ax3,prm1,prm2);
348 // return gp_Pnt(0,0,0);
350 // ============================================================
351 void IntSurf_Quadric::D1(const Standard_Real U,
352 const Standard_Real V,
359 ElSLib::PlaneD1(U,V,ax3,P,D1U,D1V);
361 case GeomAbs_Cylinder:
362 ElSLib::CylinderD1(U,V,ax3,prm1,P,D1U,D1V);
365 ElSLib::SphereD1(U,V,ax3,prm1,P,D1U,D1V);
368 ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
371 ElSLib::TorusD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
379 // ============================================================
380 gp_Vec IntSurf_Quadric::DN(const Standard_Real U,
381 const Standard_Real V,
382 const Standard_Integer Nu,
383 const Standard_Integer Nv) const
387 return ElSLib::PlaneDN(U,V,ax3,Nu,Nv);
388 case GeomAbs_Cylinder:
389 return ElSLib::CylinderDN(U,V,ax3,prm1,Nu,Nv);
391 return ElSLib::SphereDN(U,V,ax3,prm1,Nu,Nv);
393 return ElSLib::ConeDN(U,V,ax3,prm1,prm2,Nu,Nv);
395 return ElSLib::TorusDN(U,V,ax3,prm1,prm2,Nu,Nv);
404 // return gp_Vec(0,0,0);
406 // ============================================================
407 gp_Vec IntSurf_Quadric::Normale(const Standard_Real U,
408 const Standard_Real V) const
413 return ax3.Direction();
415 return ax3.Direction().Reversed();
416 case GeomAbs_Cylinder:
417 return Normale(Value(U,V));
419 return Normale(Value(U,V));
424 ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1u,D1v);
425 if(D1u.Magnitude()<0.0000001) {
426 gp_Vec Vn(0.0,0.0,0.0);
429 return(D1u.Crossed(D1v));
432 return Normale(Value(U,V));
441 // return gp_Vec(0,0,0);
443 // ============================================================
444 gp_Vec IntSurf_Quadric::Normale (const gp_Pnt& P) const
449 return ax3.Direction();
451 return ax3.Direction().Reversed();
452 case GeomAbs_Cylinder:
455 return lin.Normal(P).Direction();
458 gp_Dir D(lin.Normal(P).Direction());
466 gp_Vec ax3P(ax3.Location(),P);
470 gp_Vec Pax3(P,ax3.Location());
477 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
485 gp_Vec OZ (ax3.Direction());
486 Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
488 gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
489 ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
490 PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
491 if (PT.SquareDistance(P) < 1e-14) {
494 gp_Dir aD(ax3direc ? gp_Vec(PT, P) : gp_Vec(P, PT));
504 // ============================================================
505 void IntSurf_Quadric::Parameters (const gp_Pnt& P,
507 Standard_Real& V) const
511 ElSLib::PlaneParameters(ax3,P,U,V);
513 case GeomAbs_Cylinder:
514 ElSLib::CylinderParameters(ax3,prm1,P,U,V);
517 ElSLib::SphereParameters(ax3,prm1,P,U,V);
520 ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
523 ElSLib::TorusParameters(ax3,prm1,prm2,P,U,V);
529 // ============================================================