projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0026360: Useless global functions IsSimilar() in Standard
[occt.git] / src / GeomLib / GeomLib.cxx
CommitLineData
b311480e 1// Created on: 1993-07-07
2// Created by: Jean Claude VAUTHIER
3// Copyright (c) 1993-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 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
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
7fd59977 17// Version:
b311480e 18//pmn 24/09/96 Ajout du prolongement de courbe.
7fd59977 19// jct 15/04/97 Ajout du prolongement de surface.
20// jct 24/04/97 simplification ou suppression de calculs
21// inutiles dans ExtendSurfByLength
22// correction de Tbord et Continuity=0 accepte
23// correction du calcul de lambda et appel a
24// TangExtendToConstraint avec lambmin au lieu de 1.
25// correction du passage Sr rat --> BSp nD
26// xab 26/06/97 treatement partiel anulation des derivees
27// partiels du denonimateur des Surfaces BSplines Rationnelles
28// dans le cas de valeurs proportionnelles des denominateurs
29// en umin umax et/ou vmin vmax.
30// pmn 4/07/97 Gestion de la continuite dans BuildCurve3d (PRO9097)
7fd59977 31// xab 10/07/97 on revient en arriere sur l'ajout du 26/06/97
32// pmn 26/09/97 Ajout des parametres d'approx dans BuildCurve3d
33// xab 29/09/97 on reintegre l'ajout du 26/06/97
34// pmn 31/10/97 Ajoute AdjustExtremity
35// jct 26/11/98 blindage dans ExtendSurf qd NTgte = 0 (CTS21288)
36// jct 19/01/99 traitement de la periodicite dans ExtendSurf
37// Design:
38// Warning: None
39// References: None
40// Language: C++2.0
41// Purpose:
7fd59977 42// Declarations:
43
7fd59977 44#include <Adaptor2d_HCurve2d.hxx>
42cf5bc1 45#include <Adaptor3d_Curve.hxx>
46#include <Adaptor3d_CurveOnSurface.hxx>
7fd59977 47#include <Adaptor3d_HCurve.hxx>
48#include <Adaptor3d_HSurface.hxx>
42cf5bc1 49#include <AdvApprox_ApproxAFunction.hxx>
50#include <AdvApprox_PrefAndRec.hxx>
51#include <BSplCLib.hxx>
52#include <BSplSLib.hxx>
53#include <CSLib.hxx>
54#include <CSLib_NormalStatus.hxx>
55#include <ElCLib.hxx>
7fd59977 56#include <Geom2d_BezierCurve.hxx>
42cf5bc1 57#include <Geom2d_BSplineCurve.hxx>
7fd59977 58#include <Geom2d_Circle.hxx>
42cf5bc1 59#include <Geom2d_Curve.hxx>
7fd59977 60#include <Geom2d_Ellipse.hxx>
7fd59977 61#include <Geom2d_Hyperbola.hxx>
42cf5bc1 62#include <Geom2d_Line.hxx>
7fd59977 63#include <Geom2d_OffsetCurve.hxx>
42cf5bc1 64#include <Geom2d_Parabola.hxx>
65#include <Geom2d_TrimmedCurve.hxx>
66#include <Geom2dAdaptor_Curve.hxx>
67#include <Geom2dAdaptor_GHCurve.hxx>
68#include <Geom2dAdaptor_HCurve.hxx>
69#include <Geom2dConvert.hxx>
70#include <Geom_BezierCurve.hxx>
7fd59977 71#include <Geom_BezierSurface.hxx>
42cf5bc1 72#include <Geom_BoundedCurve.hxx>
73#include <Geom_BoundedSurface.hxx>
74#include <Geom_BSplineCurve.hxx>
7fd59977 75#include <Geom_BSplineSurface.hxx>
42cf5bc1 76#include <Geom_Circle.hxx>
77#include <Geom_Curve.hxx>
78#include <Geom_Ellipse.hxx>
79#include <Geom_Hyperbola.hxx>
80#include <Geom_Line.hxx>
81#include <Geom_OffsetCurve.hxx>
82#include <Geom_Parabola.hxx>
83#include <Geom_Plane.hxx>
84#include <Geom_RectangularTrimmedSurface.hxx>
85#include <Geom_Surface.hxx>
86#include <Geom_TrimmedCurve.hxx>
87#include <GeomAdaptor_HSurface.hxx>
88#include <GeomAdaptor_Surface.hxx>
89#include <GeomConvert.hxx>
90#include <GeomConvert_ApproxSurface.hxx>
91#include <GeomConvert_CompCurveToBSplineCurve.hxx>
92#include <GeomLib.hxx>
93#include <GeomLib_DenominatorMultiplier.hxx>
94#include <GeomLib_DenominatorMultiplierPtr.hxx>
95#include <GeomLib_LogSample.hxx>
96#include <GeomLib_MakeCurvefromApprox.hxx>
97#include <GeomLib_PolyFunc.hxx>
98#include <gp_Ax2.hxx>
7fd59977 99#include <gp_Circ.hxx>
100#include <gp_Circ2d.hxx>
42cf5bc1 101#include <gp_Dir.hxx>
7fd59977 102#include <gp_Elips.hxx>
103#include <gp_Elips2d.hxx>
42cf5bc1 104#include <gp_GTrsf2d.hxx>
7fd59977 105#include <gp_Hypr.hxx>
106#include <gp_Hypr2d.hxx>
42cf5bc1 107#include <gp_Lin.hxx>
108#include <gp_Lin2d.hxx>
7fd59977 109#include <gp_Parab.hxx>
110#include <gp_Parab2d.hxx>
42cf5bc1 111#include <gp_Pnt.hxx>
112#include <gp_Pnt2d.hxx>
7fd59977 113#include <gp_Trsf2d.hxx>
42cf5bc1 114#include <gp_TrsfForm.hxx>
115#include <gp_Vec.hxx>
116#include <Hermit.hxx>
117#include <math.hxx>
118#include <math_FunctionAllRoots.hxx>
119#include <math_FunctionSample.hxx>
120#include <math_Jacobi.hxx>
121#include <math_Matrix.hxx>
122#include <math_Vector.hxx>
123#include <PLib.hxx>
124#include <Precision.hxx>
7fd59977 125#include <Standard_ConstructionError.hxx>
42cf5bc1 126#include <Standard_NotImplemented.hxx>
127#include <TColgp_Array1OfPnt.hxx>
128#include <TColgp_Array1OfPnt2d.hxx>
129#include <TColgp_Array1OfVec.hxx>
130#include <TColgp_Array1OfXYZ.hxx>
131#include <TColgp_Array2OfPnt.hxx>
132#include <TColgp_HArray2OfPnt.hxx>
133#include <TColStd_Array1OfInteger.hxx>
134#include <TColStd_Array1OfReal.hxx>
135#include <TColStd_Array2OfReal.hxx>
136#include <TColStd_HArray1OfReal.hxx>
137#include <TColStd_HArray2OfReal.hxx>
7fd59977 138
139//=======================================================================
140//function : ComputeLambda
141//purpose : Calcul le facteur lambda qui minimise la variation de vittesse
142// sur une interpolation d'hermite d'ordre (i,0)
143//=======================================================================
144static void ComputeLambda(const math_Matrix& Constraint,
145 const math_Matrix& Hermit,
146 const Standard_Real Length,
147 Standard_Real& Lambda )
148{
149 Standard_Integer size = Hermit.RowNumber();
150 Standard_Integer Continuity = size-2;
151 Standard_Integer ii, jj, ip, pp;
152
153 //Minimization
154 math_Matrix HDer(1, size-1, 1, size);
155 for (jj=1; jj<=size; jj++) {
156 for (ii=1; ii<size;ii++) {
157 HDer(ii, jj) = ii*Hermit(jj, ii+1);
158 }
159 }
160
161 math_Vector V(1, size);
162 math_Vector Vec1(1, Constraint.RowNumber());
163 math_Vector Vec2(1, Constraint.RowNumber());
164 math_Vector Vec3(1, Constraint.RowNumber());
165 math_Vector Vec4(1, Constraint.RowNumber());
166
167 Standard_Real * polynome = &HDer(1,1);
168 Standard_Real * valhder = &V(1);
169 Vec2 = Constraint.Col(2);
170 Vec2 /= Length;
171 Standard_Real t, squared1 = Vec2.Norm2(), GW;
172// math_Matrix Vec(1, Constraint.RowNumber(), 1, size-1);
173// gp_Vec Vfirst(p0.XYZ()), Vlast(Point.XYZ());
174// TColgp_Array1OfVec Der(2, 4);
175// Der(2) = d1; Der(3) = d2; Der(4) = d3;
176
177 Standard_Integer GOrdre = 4 + 4*Continuity,
178 DDim=Continuity*(Continuity+2);
179 math_Vector GaussP(1, GOrdre), GaussW(1, GOrdre),
180 pol2(1, 2*Continuity+1),
181 pol4(1, 4*Continuity+1);
182 math::GaussPoints(GOrdre, GaussP);
183 math::GaussWeights (GOrdre, GaussW);
184 pol4.Init(0.);
185
186 for (ip=1; ip<=GOrdre; ip++) {
187 t = (GaussP(ip)+1.)/2;
188 GW = GaussW(ip);
189 PLib::NoDerivativeEvalPolynomial(t , Continuity, Continuity+2, DDim,
190 polynome[0], valhder[0]);
191 V /= Length; //Normalisation
192
193 // i
194 // C'(t) = SUM Vi*Lambda
195 Vec1 = Constraint.Col(1);
196 Vec1 *= V(1);
197 Vec1 += V(size)*Constraint.Col(size);
198 Vec2 = Constraint.Col(2);
199 Vec2 *= V(2);
200 if (Continuity > 1) {
201 Vec3 = Constraint.Col(3);
202 Vec3 *= V(3);
203 if (Continuity > 2) {
204 Vec4 = Constraint.Col(4);
205 Vec4 *= V(4);
206 }
207 }
208
209
210 // 2 2
211 // C'(t) - C'(0)
212
213 pol2(1) = Vec1.Norm2();
214 pol2(2) = 2*(Vec1.Multiplied(Vec2));
215 pol2(3) = Vec2.Norm2() - squared1;
216 if (Continuity>1) {
217 pol2(3) += 2*(Vec1.Multiplied(Vec3));
218 pol2(4) = 2*(Vec2.Multiplied(Vec3));
219 pol2(5) = Vec3.Norm2();
220 if (Continuity>2) {
221 pol2(4)+= 2*(Vec1.Multiplied(Vec4));
222 pol2(5)+= 2*(Vec2.Multiplied(Vec4));
223 pol2(6) = 2*(Vec3.Multiplied(Vec4));
224 pol2(7) = Vec4.Norm2();
225 }
226 }
227
228 // 2 2 2
229 // Integrale de ( C'(t) - C'(0) )
230 for (ii=1; ii<=pol2.Length(); ii++) {
231 pp = ii;
232 for(jj=1; jj<ii; jj++, pp++) {
233 pol4(pp) += 2*GW*pol2(ii)*pol2(jj);
234 }
235 pol4(2*ii-1) += GW*Pow(pol2(ii), 2);
236 }
237 }
238
239 Standard_Real EMin, E;
240 PLib::NoDerivativeEvalPolynomial(Lambda , pol4.Length()-1, 1,
241 pol4.Length()-1,
242 pol4(1), EMin);
243
244 if (EMin > Precision::Confusion()) {
245 // Recheche des extrema de la fonction
246 GeomLib_PolyFunc FF(pol4);
247 GeomLib_LogSample S(Lambda/1000, 50*Lambda, 100);
248 math_FunctionAllRoots Solve(FF, S, Precision::Confusion(),
249 Precision::Confusion()*(Length+1),
250 1.e-15);
251 if (Solve.IsDone()) {
252 for (ii=1; ii<=Solve.NbPoints(); ii++) {
253 t = Solve.GetPoint(ii);
254 PLib::NoDerivativeEvalPolynomial(t , pol4.Length()-1, 1,
255 pol4.Length()-1,
256 pol4(1), E);
257 if (E < EMin) {
258 Lambda = t;
259 EMin = E;
260 }
261 }
262 }
263 }
264}
265
266#include <Extrema_LocateExtPC.hxx>
ec357c5c 267#include <Geom2d_Curve.hxx>
7fd59977 268//=======================================================================
269//function : RemovePointsFromArray
270//purpose :
271//=======================================================================
272
273void GeomLib::RemovePointsFromArray(const Standard_Integer NumPoints,
274 const TColStd_Array1OfReal& InParameters,
275 Handle(TColStd_HArray1OfReal)& OutParameters)
276{
277 Standard_Integer ii,
278 jj,
279 add_one_point,
280 loc_num_points,
281 num_points,
282 index ;
283 Standard_Real delta,
284 current_parameter ;
285
286 loc_num_points = Max(0,NumPoints-2) ;
287 delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ;
288 delta /= (Standard_Real) (loc_num_points + 1) ;
289 num_points = 1 ;
290 current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ;
291 ii = InParameters.Lower() + 1 ;
292 for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) {
293 add_one_point = 0 ;
294 while ( ii < InParameters.Upper() && InParameters(ii) < current_parameter) {
295 ii += 1 ;
296 add_one_point = 1 ;
297 }
298 num_points += add_one_point ;
299 current_parameter += delta ;
300 }
301 if (NumPoints <= 2) {
302 num_points = 2 ;
303 }
304 index = 2 ;
305 current_parameter = InParameters(InParameters.Lower()) + delta * 0.5e0 ;
306 OutParameters =
307 new TColStd_HArray1OfReal(1,num_points) ;
308 OutParameters->ChangeArray1()(1) = InParameters(InParameters.Lower()) ;
309 ii = InParameters.Lower() + 1 ;
310 for (jj = 0 ; ii < InParameters.Upper() && jj < NumPoints ; jj++) {
311 add_one_point = 0 ;
312 while (ii < InParameters.Upper() && InParameters(ii) < current_parameter) {
313 ii += 1 ;
314 add_one_point = 1 ;
315 }
316 if (add_one_point && index <= num_points) {
317 OutParameters->ChangeArray1()(index) = InParameters(ii-1) ;
318 index += 1 ;
319 }
320 current_parameter += delta ;
321 }
322 OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ;
323}
324//=======================================================================
325//function : DensifyArray1OfReal
326//purpose :
327//=======================================================================
328
329void GeomLib::DensifyArray1OfReal(const Standard_Integer MinNumPoints,
330 const TColStd_Array1OfReal& InParameters,
331 Handle(TColStd_HArray1OfReal)& OutParameters)
332{
333 Standard_Integer ii,
334 in_order,
335 num_points,
336 num_parameters_to_add,
337 index ;
338 Standard_Real delta,
339 current_parameter ;
340
341 in_order = 1 ;
342 if (MinNumPoints > InParameters.Length()) {
343
344 //
345 // checks the paramaters are in increasing order
346 //
347 for (ii = InParameters.Lower() ; ii < InParameters.Upper() ; ii++) {
348 if (InParameters(ii) > InParameters(ii+1)) {
349 in_order = 0 ;
350 break ;
351 }
352 }
353 if (in_order) {
354 num_parameters_to_add = MinNumPoints - InParameters.Length() ;
355 delta = InParameters(InParameters.Upper()) - InParameters(InParameters.Lower()) ;
356 delta /= (Standard_Real) (num_parameters_to_add + 1) ;
357 num_points = MinNumPoints ;
358 OutParameters =
359 new TColStd_HArray1OfReal(1,num_points) ;
360 index = 1 ;
361 current_parameter = InParameters(InParameters.Lower()) ;
362 OutParameters->ChangeArray1()(index) = current_parameter ;
363 index += 1 ;
364 current_parameter += delta ;
365 for (ii = InParameters.Lower() + 1 ; index <= num_points && ii <= InParameters.Upper() ; ii++) {
366 while (current_parameter < InParameters(ii) && index <= num_points) {
367 OutParameters->ChangeArray1()(index) = current_parameter ;
368 index += 1 ;
369 current_parameter += delta ;
370 }
371 if (index <= num_points) {
372 OutParameters->ChangeArray1()(index) = InParameters(ii) ;
373 }
374 index += 1 ;
375 }
376 //
377 // beware of roundoff !
378 //
379 OutParameters->ChangeArray1()(num_points) = InParameters(InParameters.Upper()) ;
380 }
381 else {
382 index = 1 ;
383 num_points = InParameters.Length() ;
384 OutParameters =
385 new TColStd_HArray1OfReal(1,num_points) ;
386 for (ii = InParameters.Lower() ; ii <= InParameters.Upper() ; ii++) {
387 OutParameters->ChangeArray1()(index) = InParameters(ii) ;
388 index += 1 ;
389 }
390 }
391 }
392 else {
393 index = 1 ;
394 num_points = InParameters.Length() ;
395 OutParameters =
396 new TColStd_HArray1OfReal(1,num_points) ;
397 for (ii = InParameters.Lower() ; ii <= InParameters.Upper() ; ii++) {
398 OutParameters->ChangeArray1()(index) = InParameters(ii) ;
399 index += 1 ;
400 }
401 }
402}
403
404//=======================================================================
405//function : FuseIntervals
406//purpose :
407//=======================================================================
408void GeomLib::FuseIntervals(const TColStd_Array1OfReal& I1,
409 const TColStd_Array1OfReal& I2,
410 TColStd_SequenceOfReal& Seq,
411 const Standard_Real Epspar)
412{
413 Standard_Integer ind1=1, ind2=1;
414 Standard_Real v1, v2;
415// Initialisations : les IND1 et IND2 pointent sur le 1er element
416// de chacune des 2 tables a traiter.INDS pointe sur le dernier
417// element cree de TABSOR
418
419
420//--- On remplit TABSOR en parcourant TABLE1 et TABLE2 simultanement ---
421//------------------ en eliminant les occurrences multiples ------------
422
423 while ((ind1<=I1.Upper()) && (ind2<=I2.Upper())) {
424 v1 = I1(ind1);
425 v2 = I2(ind2);
426 if (Abs(v1-v2)<= Epspar) {
427// Ici les elements de I1 et I2 conviennent .
428 Seq.Append((v1+v2)/2);
429 ind1++;
430 ind2++;
431 }
432 else if (v1 < v2) {
433 // Ici l' element de I1 convient.
434 Seq.Append(v1);
435 ind1++;
436 }
437 else {
438// Ici l' element de TABLE2 convient.
439 Seq.Append(v2);
440 ind2++;
441 }
442 }
443
444 if (ind1>I1.Upper()) {
445//----- Ici I1 est epuise, on complete avec la fin de TABLE2 -------
446
447 for (; ind2<=I2.Upper(); ind2++) {
448 Seq.Append(I2(ind2));
449 }
450 }
451
452 if (ind2>I2.Upper()) {
453//----- Ici I2 est epuise, on complete avec la fin de I1 -------
454 for (; ind1<=I1.Upper(); ind1++) {
455 Seq.Append(I1(ind1));
456 }
457 }
458}
459
460
461//=======================================================================
462//function : EvalMaxParametricDistance
463//purpose :
464//=======================================================================
465
466void GeomLib::EvalMaxParametricDistance(const Adaptor3d_Curve& ACurve,
467 const Adaptor3d_Curve& AReferenceCurve,
468// const Standard_Real Tolerance,
469 const Standard_Real ,
470 const TColStd_Array1OfReal& Parameters,
471 Standard_Real& MaxDistance)
472{
473 Standard_Integer ii ;
474
475 Standard_Real max_squared = 0.0e0,
476// tolerance_squared,
477 local_distance_squared ;
478
479// tolerance_squared = Tolerance * Tolerance ;
480 gp_Pnt Point1 ;
481 gp_Pnt Point2 ;
482 for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
483 ACurve.D0(Parameters(ii),
484 Point1) ;
485 AReferenceCurve.D0(Parameters(ii),
486 Point2) ;
487 local_distance_squared =
488 Point1.SquareDistance (Point2) ;
489 max_squared = Max(max_squared,local_distance_squared) ;
490 }
491 if (max_squared > 0.0e0) {
492 MaxDistance = sqrt(max_squared) ;
493 }
494 else {
495 MaxDistance = 0.0e0 ;
496 }
497
498}
499//=======================================================================
500//function : EvalMaxDistanceAlongParameter
501//purpose :
502//=======================================================================
503
504void GeomLib::EvalMaxDistanceAlongParameter(const Adaptor3d_Curve& ACurve,
505 const Adaptor3d_Curve& AReferenceCurve,
506 const Standard_Real Tolerance,
507 const TColStd_Array1OfReal& Parameters,
508 Standard_Real& MaxDistance)
509{
510 Standard_Integer ii ;
511 Standard_Real max_squared = 0.0e0,
512 tolerance_squared = Tolerance * Tolerance,
513 other_parameter,
514 para_tolerance,
515 local_distance_squared ;
516 gp_Pnt Point1 ;
517 gp_Pnt Point2 ;
518
519
520
521 para_tolerance =
522 AReferenceCurve.Resolution(Tolerance) ;
523 other_parameter = Parameters(Parameters.Lower()) ;
524 ACurve.D0(other_parameter,
525 Point1) ;
526 Extrema_LocateExtPC a_projector(Point1,
527 AReferenceCurve,
528 other_parameter,
529 para_tolerance) ;
530 for (ii = Parameters.Lower() ; ii <= Parameters.Upper() ; ii++) {
531 ACurve.D0(Parameters(ii),
532 Point1) ;
533 AReferenceCurve.D0(Parameters(ii),
534 Point2) ;
535 local_distance_squared =
536 Point1.SquareDistance (Point2) ;
537
538 local_distance_squared =
539 Point1.SquareDistance (Point2) ;
540
541
542 if (local_distance_squared > tolerance_squared) {
543
544
545 a_projector.Perform(Point1,
546 other_parameter) ;
547 if (a_projector.IsDone()) {
548 other_parameter =
549 a_projector.Point().Parameter() ;
550 AReferenceCurve.D0(other_parameter,
551 Point2) ;
552 local_distance_squared =
553 Point1.SquareDistance (Point2) ;
554 }
555 else {
556 local_distance_squared = 0.0e0 ;
557 other_parameter = Parameters(ii) ;
558 }
559 }
560 else {
561 other_parameter = Parameters(ii) ;
562 }
563
564
565 max_squared = Max(max_squared,local_distance_squared) ;
566 }
567 if (max_squared > tolerance_squared) {
568 MaxDistance = sqrt(max_squared) ;
569 }
570 else {
571 MaxDistance = Tolerance ;
572 }
573}
574
575
576
577// Aliases:
578
579// Global data definitions:
580
581// Methods :
582
583
584//=======================================================================
585//function : To3d
586//purpose :
587//=======================================================================
588
589Handle(Geom_Curve) GeomLib::To3d (const gp_Ax2& Position,
590 const Handle(Geom2d_Curve)& Curve2d ) {
591 Handle(Geom_Curve) Curve3d;
592 Handle(Standard_Type) KindOfCurve = Curve2d->DynamicType();
593
594 if (KindOfCurve == STANDARD_TYPE (Geom2d_TrimmedCurve)) {
595 Handle(Geom2d_TrimmedCurve) Ct =
596 Handle(Geom2d_TrimmedCurve)::DownCast(Curve2d);
597 Standard_Real U1 = Ct->FirstParameter ();
598 Standard_Real U2 = Ct->LastParameter ();
599 Handle(Geom2d_Curve) CBasis2d = Ct->BasisCurve();
600 Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d);
601 Curve3d = new Geom_TrimmedCurve (CC, U1, U2);
602 }
603 else if (KindOfCurve == STANDARD_TYPE (Geom2d_OffsetCurve)) {
604 Handle(Geom2d_OffsetCurve) Co =
605 Handle(Geom2d_OffsetCurve)::DownCast(Curve2d);
606 Standard_Real Offset = Co->Offset();
607 Handle(Geom2d_Curve) CBasis2d = Co->BasisCurve();
608 Handle(Geom_Curve) CC = GeomLib::To3d(Position, CBasis2d);
609 Curve3d = new Geom_OffsetCurve (CC, Offset, Position.Direction());
610 }
611 else if (KindOfCurve == STANDARD_TYPE (Geom2d_BezierCurve)) {
612 Handle(Geom2d_BezierCurve) CBez2d =
613 Handle(Geom2d_BezierCurve)::DownCast (Curve2d);
614 Standard_Integer Nbpoles = CBez2d->NbPoles ();
615 TColgp_Array1OfPnt2d Poles2d (1, Nbpoles);
616 CBez2d->Poles (Poles2d);
617 TColgp_Array1OfPnt Poles3d (1, Nbpoles);
618 for (Standard_Integer i = 1; i <= Nbpoles; i++) {
619 Poles3d (i) = ElCLib::To3d (Position, Poles2d (i));
620 }
621 Handle(Geom_BezierCurve) CBez3d;
622 if (CBez2d->IsRational()) {
623 TColStd_Array1OfReal TheWeights (1, Nbpoles);
624 CBez2d->Weights (TheWeights);
625 CBez3d = new Geom_BezierCurve (Poles3d, TheWeights);
626 }
627 else {
628 CBez3d = new Geom_BezierCurve (Poles3d);
629 }
630 Curve3d = CBez3d;
631 }
632 else if (KindOfCurve == STANDARD_TYPE (Geom2d_BSplineCurve)) {
633 Handle(Geom2d_BSplineCurve) CBSpl2d =
634 Handle(Geom2d_BSplineCurve)::DownCast (Curve2d);
635 Standard_Integer Nbpoles = CBSpl2d->NbPoles ();
636 Standard_Integer Nbknots = CBSpl2d->NbKnots ();
637 Standard_Integer TheDegree = CBSpl2d->Degree ();
638 Standard_Boolean IsPeriodic = CBSpl2d->IsPeriodic();
639 TColgp_Array1OfPnt2d Poles2d (1, Nbpoles);
640 CBSpl2d->Poles (Poles2d);
641 TColgp_Array1OfPnt Poles3d (1, Nbpoles);
642 for (Standard_Integer i = 1; i <= Nbpoles; i++) {
643 Poles3d (i) = ElCLib::To3d (Position, Poles2d (i));
644 }
645 TColStd_Array1OfReal TheKnots (1, Nbknots);
646 TColStd_Array1OfInteger TheMults (1, Nbknots);
647 CBSpl2d->Knots (TheKnots);
648 CBSpl2d->Multiplicities (TheMults);
649 Handle(Geom_BSplineCurve) CBSpl3d;
650 if (CBSpl2d->IsRational()) {
651 TColStd_Array1OfReal TheWeights (1, Nbpoles);
652 CBSpl2d->Weights (TheWeights);
653 CBSpl3d = new Geom_BSplineCurve (Poles3d, TheWeights, TheKnots, TheMults, TheDegree, IsPeriodic);
654 }
655 else {
656 CBSpl3d = new Geom_BSplineCurve (Poles3d, TheKnots, TheMults, TheDegree, IsPeriodic);
657 }
658 Curve3d = CBSpl3d;
659 }
660 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Line)) {
661 Handle(Geom2d_Line) Line2d = Handle(Geom2d_Line)::DownCast (Curve2d);
662 gp_Lin2d L2d = Line2d->Lin2d();
663 gp_Lin L3d = ElCLib::To3d (Position, L2d);
664 Handle(Geom_Line) GeomL3d = new Geom_Line (L3d);
665 Curve3d = GeomL3d;
666 }
667 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Circle)) {
668 Handle(Geom2d_Circle) Circle2d =
669 Handle(Geom2d_Circle)::DownCast (Curve2d);
670 gp_Circ2d C2d = Circle2d->Circ2d();
671 gp_Circ C3d = ElCLib::To3d (Position, C2d);
672 Handle(Geom_Circle) GeomC3d = new Geom_Circle (C3d);
673 Curve3d = GeomC3d;
674 }
675 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Ellipse)) {
676 Handle(Geom2d_Ellipse) Ellipse2d =
677 Handle(Geom2d_Ellipse)::DownCast (Curve2d);
678 gp_Elips2d E2d = Ellipse2d->Elips2d ();
679 gp_Elips E3d = ElCLib::To3d (Position, E2d);
680 Handle(Geom_Ellipse) GeomE3d = new Geom_Ellipse (E3d);
681 Curve3d = GeomE3d;
682 }
683 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Parabola)) {
684 Handle(Geom2d_Parabola) Parabola2d =
685 Handle(Geom2d_Parabola)::DownCast (Curve2d);
686 gp_Parab2d Prb2d = Parabola2d->Parab2d ();
687 gp_Parab Prb3d = ElCLib::To3d (Position, Prb2d);
688 Handle(Geom_Parabola) GeomPrb3d = new Geom_Parabola (Prb3d);
689 Curve3d = GeomPrb3d;
690 }
691 else if (KindOfCurve == STANDARD_TYPE (Geom2d_Hyperbola)) {
692 Handle(Geom2d_Hyperbola) Hyperbola2d =
693 Handle(Geom2d_Hyperbola)::DownCast (Curve2d);
694 gp_Hypr2d H2d = Hyperbola2d->Hypr2d ();
695 gp_Hypr H3d = ElCLib::To3d (Position, H2d);
696 Handle(Geom_Hyperbola) GeomH3d = new Geom_Hyperbola (H3d);
697 Curve3d = GeomH3d;
698 }
699 else {
700 Standard_NotImplemented::Raise();
701 }
702
703 return Curve3d;
704}
705
706
707
708//=======================================================================
709//function : GTransform
710//purpose :
711//=======================================================================
712
713Handle(Geom2d_Curve) GeomLib::GTransform(const Handle(Geom2d_Curve)& Curve,
714 const gp_GTrsf2d& GTrsf)
715{
716 gp_TrsfForm Form = GTrsf.Form();
717
718 if ( Form != gp_Other) {
719
720 // Alors, la GTrsf est en fait une Trsf.
721 // La geometrie des courbes sera alors inchangee.
722
723 Handle(Geom2d_Curve) C =
724 Handle(Geom2d_Curve)::DownCast(Curve->Transformed(GTrsf.Trsf2d()));
725 return C;
726 }
727 else {
728
729 // Alors, la GTrsf est une other Transformation.
730 // La geometrie des courbes est alors changee, et les conics devront
731 // etre converties en BSplines.
732
733 Handle(Standard_Type) TheType = Curve->DynamicType();
734
735 if ( TheType == STANDARD_TYPE(Geom2d_TrimmedCurve)) {
736
737 // On va recurer sur la BasisCurve
738
739 Handle(Geom2d_TrimmedCurve) C =
740 Handle(Geom2d_TrimmedCurve)::DownCast(Curve->Copy());
741
742 Handle(Standard_Type) TheBasisType = (C->BasisCurve())->DynamicType();
743
744 if (TheBasisType == STANDARD_TYPE(Geom2d_BSplineCurve) ||
745 TheBasisType == STANDARD_TYPE(Geom2d_BezierCurve) ) {
746
747 // Dans ces cas le parametrage est conserve sur la courbe transformee
748 // on peut donc la trimmer avec les parametres de la courbe de base.
749
750 Standard_Real U1 = C->FirstParameter();
751 Standard_Real U2 = C->LastParameter();
752
753 Handle(Geom2d_TrimmedCurve) result =
754 new Geom2d_TrimmedCurve(GTransform(C->BasisCurve(), GTrsf), U1,U2);
755 return result;
756 }
757 else if ( TheBasisType == STANDARD_TYPE(Geom2d_Line)) {
758
759 // Dans ce cas, le parametrage n`est plus conserve.
760 // Il faut recalculer les parametres de Trimming sur la courbe
761 // resultante. ( Calcul par projection ( ElCLib) des points debut
762 // et fin transformes)
763
764 Handle(Geom2d_Line) L =
765 Handle(Geom2d_Line)::DownCast(GTransform(C->BasisCurve(), GTrsf));
766 gp_Lin2d Lin = L->Lin2d();
767
768 gp_Pnt2d P1 = C->StartPoint();
769 gp_Pnt2d P2 = C->EndPoint();
770 P1.SetXY(GTrsf.Transformed(P1.XY()));
771 P2.SetXY(GTrsf.Transformed(P2.XY()));
772 Standard_Real U1 = ElCLib::Parameter(Lin,P1);
773 Standard_Real U2 = ElCLib::Parameter(Lin,P2);
774
775 Handle(Geom2d_TrimmedCurve) result =
776 new Geom2d_TrimmedCurve(L,U1,U2);
777 return result;
778 }
779 else if (TheBasisType == STANDARD_TYPE(Geom2d_Circle) ||
780 TheBasisType == STANDARD_TYPE(Geom2d_Ellipse) ||
781 TheBasisType == STANDARD_TYPE(Geom2d_Parabola) ||
782 TheBasisType == STANDARD_TYPE(Geom2d_Hyperbola) ) {
783
784 // Dans ces cas, la geometrie de la courbe n`est pas conservee
785 // on la convertir en BSpline avant de lui appliquer la Trsf.
786
787 Handle(Geom2d_BSplineCurve) BS =
788 Geom2dConvert::CurveToBSplineCurve(C);
789 return GTransform(BS,GTrsf);
790 }
791 else {
792
793 // La transformee d`une OffsetCurve vaut ????? Sais pas faire !!
794
795 Handle(Geom2d_Curve) dummy;
796 return dummy;
797 }
798 }
799 else if ( TheType == STANDARD_TYPE(Geom2d_Line)) {
800
801 Handle(Geom2d_Line) L =
802 Handle(Geom2d_Line)::DownCast(Curve->Copy());
803 gp_Lin2d Lin = L->Lin2d();
804 gp_Pnt2d P = Lin.Location();
805 gp_Pnt2d PP = L->Value(10.); // pourquoi pas !!
806 P.SetXY(GTrsf.Transformed(P.XY()));
807 PP.SetXY(GTrsf.Transformed(PP.XY()));
808 L->SetLocation(P);
809 gp_Vec2d V(P,PP);
810 L->SetDirection(gp_Dir2d(V));
811 return L;
812 }
813 else if ( TheType == STANDARD_TYPE(Geom2d_BezierCurve)) {
814
815 // Les GTrsf etant des operation lineaires, la transformee d`une courbe
816 // a poles est la courbe dont les poles sont la transformee des poles
817 // de la courbe de base.
818
819 Handle(Geom2d_BezierCurve) C =
820 Handle(Geom2d_BezierCurve)::DownCast(Curve->Copy());
821 Standard_Integer NbPoles = C->NbPoles();
822 TColgp_Array1OfPnt2d Poles(1,NbPoles);
823 C->Poles(Poles);
824 for ( Standard_Integer i = 1; i <= NbPoles; i++) {
825 Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY()));
826 C->SetPole(i,Poles(i));
827 }
828 return C;
829 }
830 else if ( TheType == STANDARD_TYPE(Geom2d_BSplineCurve)) {
831
832 // Voir commentaire pour les Bezier.
833
834 Handle(Geom2d_BSplineCurve) C =
835 Handle(Geom2d_BSplineCurve)::DownCast(Curve->Copy());
836 Standard_Integer NbPoles = C->NbPoles();
837 TColgp_Array1OfPnt2d Poles(1,NbPoles);
838 C->Poles(Poles);
839 for ( Standard_Integer i = 1; i <= NbPoles; i++) {
840 Poles(i).SetXY(GTrsf.Transformed(Poles(i).XY()));
841 C->SetPole(i,Poles(i));
842 }
843 return C;
844 }
845 else if ( TheType == STANDARD_TYPE(Geom2d_Circle) ||
846 TheType == STANDARD_TYPE(Geom2d_Ellipse) ) {
847
848 // Dans ces cas, la geometrie de la courbe n`est pas conservee
849 // on la convertir en BSpline avant de lui appliquer la Trsf.
850
851 Handle(Geom2d_BSplineCurve) C =
852 Geom2dConvert::CurveToBSplineCurve(Curve);
853 return GTransform(C, GTrsf);
854 }
855 else if ( TheType == STANDARD_TYPE(Geom2d_Parabola) ||
856 TheType == STANDARD_TYPE(Geom2d_Hyperbola) ||
857 TheType == STANDARD_TYPE(Geom2d_OffsetCurve) ) {
858
859 // On ne sait pas faire : return a null Handle;
860
861 Handle(Geom2d_Curve) dummy;
862 return dummy;
863 }
864 }
865
866 Handle(Geom2d_Curve) WNT__; // portage Windows.
867 return WNT__;
868}
869
870
871//=======================================================================
872//function : SameRange
873//purpose :
874//=======================================================================
875void GeomLib::SameRange(const Standard_Real Tolerance,
876 const Handle(Geom2d_Curve)& CurvePtr,
877 const Standard_Real FirstOnCurve,
878 const Standard_Real LastOnCurve,
879 const Standard_Real RequestedFirst,
880 const Standard_Real RequestedLast,
881 Handle(Geom2d_Curve)& NewCurvePtr)
882{
883 if(CurvePtr.IsNull()) Standard_Failure::Raise();
884 if (Abs(LastOnCurve - RequestedLast) <= Tolerance &&
885 Abs(FirstOnCurve - RequestedFirst) <= Tolerance) {
886 NewCurvePtr = CurvePtr;
887 return;
888 }
889
890 // the parametrisation lentgh must at least be the same.
891 if (Abs(LastOnCurve - FirstOnCurve - RequestedLast + RequestedFirst)
892 <= Tolerance) {
893 if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Line))) {
894 Handle(Geom2d_Line) Line =
895 Handle(Geom2d_Line)::DownCast(CurvePtr->Copy());
896 Standard_Real dU = FirstOnCurve - RequestedFirst;
897 gp_Dir2d D = Line->Direction() ;
898 Line->Translate(dU * gp_Vec2d(D));
899 NewCurvePtr = Line;
900 }
901 else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_Circle))) {
902 gp_Trsf2d Trsf;
903 NewCurvePtr = Handle(Geom2d_Curve)::DownCast(CurvePtr->Copy());
904 Handle(Geom2d_Circle) Circ =
905 Handle(Geom2d_Circle)::DownCast(NewCurvePtr);
906 gp_Pnt2d P = Circ->Location();
907 Standard_Real dU;
908 if (Circ->Circ2d().IsDirect()) {
909 dU = FirstOnCurve - RequestedFirst;
910 }
911 else {
912 dU = RequestedFirst - FirstOnCurve;
913 }
914 Trsf.SetRotation(P,dU);
915 NewCurvePtr->Transform(Trsf) ;
916 }
917 else if (CurvePtr->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve))) {
918 Handle(Geom2d_TrimmedCurve) TC =
919 Handle(Geom2d_TrimmedCurve)::DownCast(CurvePtr);
920 GeomLib::SameRange(Tolerance,
921 TC->BasisCurve(),
922 FirstOnCurve , LastOnCurve,
923 RequestedFirst, RequestedLast,
924 NewCurvePtr);
925 NewCurvePtr = new Geom2d_TrimmedCurve( NewCurvePtr, RequestedFirst, RequestedLast );
926 }
927//
928// attention a des problemes de limitation : utiliser le MEME test que dans
929// Geom2d_TrimmedCurve::SetTrim car sinon comme on risque de relimite sur
930// RequestedFirst et RequestedLast on aura un probleme
931//
932//
933 else if (Abs(LastOnCurve - FirstOnCurve) > Precision::PConfusion() ||
934 Abs(RequestedLast + RequestedFirst) > Precision::PConfusion()) {
935
936 Handle(Geom2d_TrimmedCurve) TC =
937 new Geom2d_TrimmedCurve(CurvePtr,FirstOnCurve,LastOnCurve);
938
939 Handle(Geom2d_BSplineCurve) BS =
940 Geom2dConvert::CurveToBSplineCurve(TC);
941 TColStd_Array1OfReal Knots(1,BS->NbKnots());
942 BS->Knots(Knots);
943
944 BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots);
945
946 BS->SetKnots(Knots);
947 NewCurvePtr = BS;
948 }
949
950 }
951 else { // On segmente le resultat
952 Handle(Geom2d_TrimmedCurve) TC =
953 new Geom2d_TrimmedCurve( CurvePtr, FirstOnCurve, LastOnCurve );
905522ee 954 //
7fd59977 955 Handle(Geom2d_BSplineCurve) BS =
956 Geom2dConvert::CurveToBSplineCurve(TC);
7fd59977 957 TColStd_Array1OfReal Knots(1,BS->NbKnots());
958 BS->Knots(Knots);
959
960 BSplCLib::Reparametrize(RequestedFirst,RequestedLast,Knots);
961
962 BS->SetKnots(Knots);
963 NewCurvePtr = BS;
964 }
965}
966
967//=======================================================================
968//class : GeomLib_CurveOnSurfaceEvaluator
969//purpose: The evaluator for the Curve 3D building
970//=======================================================================
971
972class GeomLib_CurveOnSurfaceEvaluator : public AdvApprox_EvaluatorFunction
973{
974 public:
975 GeomLib_CurveOnSurfaceEvaluator (Adaptor3d_CurveOnSurface& theCurveOnSurface,
976 Standard_Real theFirst, Standard_Real theLast)
977 : CurveOnSurface(theCurveOnSurface), FirstParam(theFirst), LastParam(theLast) {}
978
979 virtual void Evaluate (Standard_Integer *Dimension,
980 Standard_Real StartEnd[2],
981 Standard_Real *Parameter,
982 Standard_Integer *DerivativeRequest,
983 Standard_Real *Result, // [Dimension]
984 Standard_Integer *ErrorCode);
985
986 private:
987 Adaptor3d_CurveOnSurface& CurveOnSurface;
988 Standard_Real FirstParam;
989 Standard_Real LastParam;
990
991 Handle(Adaptor3d_HCurve) TrimCurve;
992};
993
994void GeomLib_CurveOnSurfaceEvaluator::Evaluate (Standard_Integer *,/*Dimension*/
995 Standard_Real DebutFin[2],
996 Standard_Real *Parameter,
997 Standard_Integer *DerivativeRequest,
998 Standard_Real *Result,// [Dimension]
999 Standard_Integer *ReturnCode)
a7f510bf 1000{
1001 gp_Pnt Point;
7fd59977 1002
1003 //Gestion des positionnements gauche / droite
1004 if ((DebutFin[0] != FirstParam) || (DebutFin[1] != LastParam))
1005 {
1006 TrimCurve = CurveOnSurface.Trim(DebutFin[0], DebutFin[1], Precision::PConfusion());
1007 FirstParam = DebutFin[0];
1008 LastParam = DebutFin[1];
1009 }
1010
1011 //Positionemment
1012 if (*DerivativeRequest == 0)
1013 {
1014 TrimCurve->D0((*Parameter), Point) ;
1015
a7f510bf 1016 for (Standard_Integer ii = 0 ; ii < 3 ; ii++)
7fd59977 1017 Result[ii] = Point.Coord(ii + 1);
1018 }
1019 if (*DerivativeRequest == 1)
1020 {
1021 gp_Vec Vector;
1022 TrimCurve->D1((*Parameter), Point, Vector);
a7f510bf 1023 for (Standard_Integer ii = 0 ; ii < 3 ; ii++)
7fd59977 1024 Result[ii] = Vector.Coord(ii + 1) ;
1025 }
1026 if (*DerivativeRequest == 2)
1027 {
1028 gp_Vec Vector, VecBis;
1029 TrimCurve->D2((*Parameter), Point, VecBis, Vector);
a7f510bf 1030 for (Standard_Integer ii = 0 ; ii < 3 ; ii++)
7fd59977 1031 Result[ii] = Vector.Coord(ii + 1) ;
1032 }
1033 ReturnCode[0] = 0;
1034}
1035
1036//=======================================================================
1037//function : BuildCurve3d
1038//purpose :
1039//=======================================================================
1040
1041void GeomLib::BuildCurve3d(const Standard_Real Tolerance,
1042 Adaptor3d_CurveOnSurface& Curve,
1043 const Standard_Real FirstParameter,
1044 const Standard_Real LastParameter,
857ffd5e 1045 Handle(Geom_Curve)& NewCurvePtr,
7fd59977 1046 Standard_Real& MaxDeviation,
1047 Standard_Real& AverageDeviation,
1048 const GeomAbs_Shape Continuity,
1049 const Standard_Integer MaxDegree,
1050 const Standard_Integer MaxSegment)
1051
1052{
1053
1054
1055 Standard_Integer curve_not_computed = 1 ;
1056 MaxDeviation = 0.0e0 ;
1057 AverageDeviation = 0.0e0 ;
c5f3a425 1058 Handle(GeomAdaptor_HSurface) geom_adaptor_surface_ptr (Handle(GeomAdaptor_HSurface)::DownCast(Curve.GetSurface()) );
1059 Handle(Geom2dAdaptor_HCurve) geom_adaptor_curve_ptr (Handle(Geom2dAdaptor_HCurve)::DownCast(Curve.GetCurve()) );
7fd59977 1060
1061 if (! geom_adaptor_curve_ptr.IsNull() &&
1062 ! geom_adaptor_surface_ptr.IsNull()) {
1063 Handle(Geom_Plane) P ;
1064 const GeomAdaptor_Surface & geom_surface =
1065 * (GeomAdaptor_Surface *) &geom_adaptor_surface_ptr->Surface() ;
1066
1067 Handle(Geom_RectangularTrimmedSurface) RT =
1068 Handle(Geom_RectangularTrimmedSurface)::
1069 DownCast(geom_surface.Surface());
1070 if ( RT.IsNull()) {
1071 P = Handle(Geom_Plane)::DownCast(geom_surface.Surface());
1072 }
1073 else {
1074 P = Handle(Geom_Plane)::DownCast(RT->BasisSurface());
1075 }
1076
1077
1078 if (! P.IsNull()) {
1079 // compute the 3d curve
1080 gp_Ax2 axes = P->Position().Ax2();
1081 const Geom2dAdaptor_Curve & geom2d_curve =
1082 * (Geom2dAdaptor_Curve *) & geom_adaptor_curve_ptr->Curve2d() ;
1083 NewCurvePtr =
1084 GeomLib::To3d(axes,
1085 geom2d_curve.Curve());
1086 curve_not_computed = 0 ;
1087
1088 }
1089 }
1090 if (curve_not_computed) {
1091
1092 //
1093 // Entree
1094 //
1095 Handle(TColStd_HArray1OfReal) Tolerance1DPtr,Tolerance2DPtr;
1096 Handle(TColStd_HArray1OfReal) Tolerance3DPtr =
1097 new TColStd_HArray1OfReal(1,1) ;
1098 Tolerance3DPtr->SetValue(1,Tolerance);
1099
1100 // Recherche des discontinuitees
1101 Standard_Integer NbIntervalC2 = Curve.NbIntervals(GeomAbs_C2);
1102 TColStd_Array1OfReal Param_de_decoupeC2 (1, NbIntervalC2+1);
1103 Curve.Intervals(Param_de_decoupeC2, GeomAbs_C2);
1104
1105 Standard_Integer NbIntervalC3 = Curve.NbIntervals(GeomAbs_C3);
1106 TColStd_Array1OfReal Param_de_decoupeC3 (1, NbIntervalC3+1);
1107 Curve.Intervals(Param_de_decoupeC3, GeomAbs_C3);
1108
1109 // Note extension of the parameteric range
1110 // Pour forcer le Trim au premier appel de l'evaluateur
1111 GeomLib_CurveOnSurfaceEvaluator ev (Curve, FirstParameter - 1., LastParameter + 1.);
1112
1113 // Approximation avec decoupe preferentiel
1114 AdvApprox_PrefAndRec Preferentiel(Param_de_decoupeC2,
1115 Param_de_decoupeC3);
1116 AdvApprox_ApproxAFunction anApproximator(0,
1117 0,
1118 1,
1119 Tolerance1DPtr,
1120 Tolerance2DPtr,
1121 Tolerance3DPtr,
1122 FirstParameter,
1123 LastParameter,
1124 Continuity,
1125 MaxDegree,
1126 MaxSegment,
1127 ev,
1128// CurveOnSurfaceEvaluator,
1129 Preferentiel) ;
1130
1131 if (anApproximator.HasResult()) {
1132 GeomLib_MakeCurvefromApprox
1133 aCurveBuilder(anApproximator) ;
1134
1135 Handle(Geom_BSplineCurve) aCurvePtr =
1136 aCurveBuilder.Curve(1) ;
1137 // On rend les resultats de l'approx
1138 MaxDeviation = anApproximator.MaxError(3,1) ;
1139 AverageDeviation = anApproximator.AverageError(3,1) ;
1140 NewCurvePtr = aCurvePtr ;
1141 }
1142 }
1143 }
1144
1145//=======================================================================
1146//function : AdjustExtremity
1147//purpose :
1148//=======================================================================
1149
1150void GeomLib::AdjustExtremity(Handle(Geom_BoundedCurve)& Curve,
1151 const gp_Pnt& P1,
1152 const gp_Pnt& P2,
1153 const gp_Vec& T1,
1154 const gp_Vec& T2)
1155{
1156// il faut Convertir l'entree (en preservant si possible le parametrage)
1157 Handle(Geom_BSplineCurve) aIn, aDef;
1158 aIn = GeomConvert::CurveToBSplineCurve(Curve, Convert_QuasiAngular);
1159
1160 Standard_Integer ii, jj;
1161 gp_Pnt P;
1162 gp_Vec V, Vtan, DV;
1163 TColgp_Array1OfPnt PolesDef(1,4), Coeffs(1,4);
1164 TColStd_Array1OfReal FK(1, 8);
1165 TColStd_Array1OfReal Ti(1, 4);
1166 TColStd_Array1OfInteger Contact(1, 4);
1167
1168 Ti(1) = Ti(2) = aIn->FirstParameter();
1169 Ti(3) = Ti(4) = aIn->LastParameter();
1170 Contact(1) = Contact(3) = 0;
1171 Contact(2) = Contact(4) = 1;
1172 for (ii=1; ii<=4; ii++) {
1173 FK(ii) = aIn->FirstParameter();
1174 FK(ii) = aIn->LastParameter();
1175 }
1176
1177 // Calculs des contraintes de deformations
1178 aIn->D1(Ti(1), P, V);
1179 PolesDef(1).ChangeCoord() = P1.XYZ()-P.XYZ();
1180 Vtan = T1;
1181 Vtan.Normalize();
1182 DV = Vtan * (Vtan * V) - V;
1183 PolesDef(2).ChangeCoord() = (Ti(4)-Ti(1))*DV.XYZ();
1184
1185 aIn->D1(Ti(4), P, V);
1186 PolesDef(3).ChangeCoord() = P2.XYZ()-P.XYZ();
1187 Vtan = T2;
1188 Vtan.Normalize();
1189 DV = Vtan * (Vtan * V) - V;
1190 PolesDef(4).ChangeCoord() = (Ti(4)-Ti(1))* DV.XYZ();
1191
1192 // Interpolation des contraintes
1193 math_Matrix Mat(1, 4, 1, 4);
1194 if (!PLib::HermiteCoefficients(0., 1., 1, 1, Mat))
1195 Standard_ConstructionError::Raise();
1196
1197 for (jj=1; jj<=4; jj++) {
1198 gp_XYZ aux(0.,0.,0.);
1199 for (ii=1; ii<=4; ii++) {
1200 aux.SetLinearForm(Mat(ii,jj), PolesDef(ii).XYZ(), aux);
1201 }
1202 Coeffs(jj).SetXYZ(aux);
1203 }
1204
1205 PLib::CoefficientsPoles(Coeffs, PLib::NoWeights(),
1206 PolesDef, PLib::NoWeights());
1207
1208 // Ajout de la deformation
1209 TColStd_Array1OfReal K(1, 2);
1210 TColStd_Array1OfInteger M(1, 2);
1211 K(1) = Ti(1);
1212 K(2) = Ti(4);
1213 M.Init(4);
1214
1215 aDef = new (Geom_BSplineCurve) (PolesDef, K, M, 3);
1216 if (aIn->Degree() < 3) aIn->IncreaseDegree(3);
1217 else aDef->IncreaseDegree(aIn->Degree());
1218
1219 for (ii=2; ii<aIn->NbKnots(); ii++) {
1220 aDef->InsertKnot(aIn->Knot(ii), aIn->Multiplicity(ii));
1221 }
1222
1223 if (aDef->NbPoles() != aIn->NbPoles())
1224 Standard_ConstructionError::Raise("Inconsistent poles's number");
1225
1226 for (ii=1; ii<=aDef->NbPoles(); ii++) {
1227 P = aIn->Pole(ii);
1228 P.ChangeCoord() += aDef->Pole(ii).XYZ();
1229 aIn->SetPole(ii, P);
1230 }
1231 Curve = aIn;
1232}
1233//=======================================================================
1234//function : ExtendCurveToPoint
1235//purpose :
1236//=======================================================================
1237
1238void GeomLib::ExtendCurveToPoint(Handle(Geom_BoundedCurve)& Curve,
1239 const gp_Pnt& Point,
1240 const Standard_Integer Continuity,
1241 const Standard_Boolean After)
1242{
1243 if(Continuity < 1 || Continuity > 3) return;
1244 Standard_Integer size = Continuity + 2;
1245 Standard_Real Ubord, Tol=1.e-6;
1246 math_Matrix MatCoefs(1,size, 1,size);
1247 Standard_Real Lambda, L1;
1248 Standard_Integer ii, jj;
1249 gp_Vec d1, d2, d3;
1250 gp_Pnt p0;
1251// il faut Convertir l'entree (en preservant si possible le parametrage)
1252 GeomConvert_CompCurveToBSplineCurve Concat(Curve, Convert_QuasiAngular);
1253
1254// Les contraintes de constructions
1255 TColgp_Array1OfXYZ Cont(1,size);
1256 if (After) {
1257 Ubord = Curve->LastParameter();
1258
1259 }
1260 else {
1261 Ubord = Curve->FirstParameter();
1262 }
1263 PLib::HermiteCoefficients(0, 1, // Les Bornes
1264 Continuity, 0, // Les Ordres de contraintes
1265 MatCoefs);
1266
1267 Curve->D3(Ubord, p0, d1, d2, d3);
1268 if (!After) { // Inversion du parametrage
1269 d1 *= -1;
1270 d3 *= -1;
1271 }
1272
1273 L1 = p0.Distance(Point);
1274 if (L1 > Tol) {
1275 // Lambda est le ratio qu'il faut appliquer a la derive de la courbe
1276 // pour obtenir la derive du prolongement (fixe arbitrairement a la
1277 // longueur du segment bout de la courbe - point cible.
1278 // On essai d'avoir sur le prolongement la vitesse moyenne que l'on
1279 // a sur la courbe.
1280 gp_Vec daux;
1281 gp_Pnt pp;
1282 Standard_Real f= Curve->FirstParameter(), t, dt, norm;
1283 dt = (Curve->LastParameter()-f)/9;
1284 norm = d1.Magnitude();
1285 for (ii=1, t=f+dt; ii<=8; ii++, t+=dt) {
1286 Curve->D1(t, pp, daux);
1287 norm += daux.Magnitude();
1288 }
1289 norm /= 9;
1290 dt = d1.Magnitude() / norm;
1291 if ((dt<1.5) && (dt>0.75)) { // Le bord est dans la moyenne on le garde
1292 Lambda = ((Standard_Real)1) / Max (d1.Magnitude() / L1, Tol);
1293 }
1294 else {
1295 Lambda = ((Standard_Real)1) / Max (norm / L1, Tol);
1296 }
1297 }
1298 else {
1299 return; // Pas d'extension
1300 }
1301
1302 // Optimisation du Lambda
1303 math_Matrix Cons(1, 3, 1, size);
1304 Cons(1,1) = p0.X(); Cons(2,1) = p0.Y(); Cons(3,1) = p0.Z();
1305 Cons(1,2) = d1.X(); Cons(2,2) = d1.Y(); Cons(3,2) = d1.Z();
1306 Cons(1,size) = Point.X(); Cons(2,size) = Point.Y(); Cons(3,size) = Point.Z();
1307 if (Continuity >= 2) {
1308 Cons(1,3) = d2.X(); Cons(2,3) = d2.Y(); Cons(3,3) = d2.Z();
1309 }
1310 if (Continuity >= 3) {
1311 Cons(1,4) = d3.X(); Cons(2,4) = d3.Y(); Cons(3,4) = d3.Z();
1312 }
1313 ComputeLambda(Cons, MatCoefs, L1, Lambda);
1314
1315 // Construction dans la Base Polynomiale
1316 Cont(1) = p0.XYZ();
1317 Cont(2) = d1.XYZ() * Lambda;
1318 if(Continuity >= 2) Cont(3) = d2.XYZ() * Pow(Lambda,2);
1319 if(Continuity >= 3) Cont(4) = d3.XYZ() * Pow(Lambda,3);
1320 Cont(size) = Point.XYZ();
1321
1322
1323 TColgp_Array1OfPnt ExtrapPoles(1, size);
1324 TColgp_Array1OfPnt ExtraCoeffs(1, size);
1325
1326 gp_Pnt PNull(0.,0.,0.);
1327 ExtraCoeffs.Init(PNull);
1328 for (ii=1; ii<=size; ii++) {
1329 for (jj=1; jj<=size; jj++) {
1330 ExtraCoeffs(jj).ChangeCoord() += MatCoefs(ii,jj)*Cont(ii);
1331 }
1332 }
1333
1334 // Convertion Dans la Base de Bernstein
1335 PLib::CoefficientsPoles(ExtraCoeffs, PLib::NoWeights(),
1336 ExtrapPoles, PLib::NoWeights());
1337
1338 Handle(Geom_BezierCurve) Bezier = new (Geom_BezierCurve) (ExtrapPoles);
1339
1340 Standard_Real dist = ExtrapPoles(1).Distance(p0);
1341 Standard_Boolean Ok;
1342 Tol += dist;
1343
1344 // Concatenation
1345 Ok = Concat.Add(Bezier, Tol, After);
1346 if (!Ok) Standard_ConstructionError::Raise("ExtendCurveToPoint");
1347
1348 Curve = Concat.BSplineCurve();
1349}
1350
1351
1352//=======================================================================
1353//function : ExtendKPart
1354//purpose : Extension par longueur des surfaces cannonique
1355//=======================================================================
1356static Standard_Boolean
1357ExtendKPart(Handle(Geom_RectangularTrimmedSurface)& Surface,
1358 const Standard_Real Length,
1359 const Standard_Boolean InU,
1360 const Standard_Boolean After)
1361{
1362
1363 if (Surface.IsNull()) return Standard_False;
1364
1365 Standard_Boolean Ok=Standard_True;
1366 Standard_Real Uf, Ul, Vf, Vl;
1367 Handle(Geom_Surface) Support = Surface->BasisSurface();
1368 GeomAbs_SurfaceType Type;
1369
1370 Surface->Bounds(Uf, Ul, Vf, Vl);
1371 GeomAdaptor_Surface AS(Surface);
1372 Type = AS.GetType();
1373
1374 if (InU) {
1375 switch(Type) {
1376 case GeomAbs_Plane :
1377 {
1378 if (After) Ul+=Length;
1379 else Uf-=Length;
1380 Surface = new (Geom_RectangularTrimmedSurface)
1381 (Support, Uf, Ul, Vf, Vl);
1382 break;
1383 }
1384
1385 default:
1386 Ok = Standard_False;
1387 }
1388 }
1389 else {
1390 switch(Type) {
1391 case GeomAbs_Plane :
1392 case GeomAbs_Cylinder :
1393 case GeomAbs_SurfaceOfExtrusion :
1394 {
1395 if (After) Vl+=Length;
1396 else Vf-=Length;
1397 Surface = new (Geom_RectangularTrimmedSurface)
1398 (Support, Uf, Ul, Vf, Vl);
1399 break;
1400 }
1401 default:
1402 Ok = Standard_False;
1403 }
1404 }
1405
1406 return Ok;
1407}
1408
1409//=======================================================================
1410//function : ExtendSurfByLength
1411//purpose :
1412//=======================================================================
1413void GeomLib::ExtendSurfByLength(Handle(Geom_BoundedSurface)& Surface,
1414 const Standard_Real Length,
1415 const Standard_Integer Continuity,
1416 const Standard_Boolean InU,
1417 const Standard_Boolean After)
1418{
1419 if(Continuity < 0 || Continuity > 3) return;
1420 Standard_Integer Cont = Continuity;
1421
1422 // Kpart ?
1423 Handle(Geom_RectangularTrimmedSurface) TS =
1424 Handle(Geom_RectangularTrimmedSurface)::DownCast (Surface);
1425 if (ExtendKPart(TS,Length, InU, After) ) {
1426 Surface = TS;
1427 return;
1428 }
1429
1430// format BSplineSurface avec un degre suffisant pour la continuite voulue
1431 Handle(Geom_BSplineSurface) BS =
1432 Handle(Geom_BSplineSurface)::DownCast (Surface);
1433 if (BS.IsNull()) {
1434 //BS = GeomConvert::SurfaceToBSplineSurface(Surface);
1435 Standard_Real Tol = Precision::Confusion(); //1.e-4;
1436 GeomAbs_Shape UCont = GeomAbs_C1, VCont = GeomAbs_C1;
1437 Standard_Integer degU = 14, degV = 14;
1438 Standard_Integer nmax = 16;
543a9964 1439 Standard_Integer thePrec = 1;
1440 const Handle(Geom_Surface)& aSurf = Surface; // to resolve ambiguity
1441 GeomConvert_ApproxSurface theApprox(aSurf,Tol,UCont,VCont,degU,degV,nmax,thePrec);
7fd59977 1442 if (theApprox.HasResult())
1443 BS = theApprox.Surface();
1444 else
1445 BS = GeomConvert::SurfaceToBSplineSurface(Surface);
1446 }
1447 if (InU&&(BS->UDegree()<Continuity+1))
1448 BS->IncreaseDegree(Continuity+1,BS->VDegree());
1449 if (!InU&&(BS->VDegree()<Continuity+1))
1450 BS->IncreaseDegree(BS->UDegree(),Continuity+1);
1451
1452 // si BS etait periodique dans le sens de l'extension, elle ne le sera plus
1453 if ( (InU&&(BS->IsUPeriodic())) || (!InU&&(BS->IsVPeriodic())) ) {
1454 Standard_Real U0,U1,V0,V1;
1455 BS->Bounds(U0,U1,V0,V1);
1456 BS->Segment(U0,U1,V0,V1);
1457 }
1458
1459
47c580a7
A		 1460// IFV Fix OCC bug 0022694 - wrong result extrapolating rational surfaces
1461// Standard_Boolean rational = ( InU && BS->IsURational() )
1462// || ( !InU && BS->IsVRational() ) ;
1463 Standard_Boolean rational = (BS->IsURational() || BS->IsVRational());
7fd59977 1464 Standard_Boolean NullWeight;
1465 Standard_Real EpsW = 10*Precision::PConfusion();
1466 Standard_Integer gap = 3;
1467 if ( rational ) gap++;
1468
1469
1470
1d47d8d0 1471 Standard_Integer Cdeg = 0, Cdim = 0, NbP = 0, Ksize = 0, Psize = 1;
7fd59977 1472 Standard_Integer ii, jj, ipole, Kount;
1473 Standard_Real Tbord, lambmin=Length;
1d47d8d0 1474 Standard_Real * Padr = NULL;
7fd59977 1475 Standard_Boolean Ok;
1476 Handle(TColStd_HArray1OfReal) FKnots, Point, lambda, Tgte, Poles;
1477
1478
1479
1480
1481 for (Kount=0, Ok=Standard_False; Kount<=2 && !Ok; Kount++) {
1482 // transformation de la surface en une BSpline non rationnelle a une variable
1483 // de degre UDegree ou VDegree et de dimension 3 ou 4 x NbVpoles ou NbUpoles
1484 // le nombre de poles egal a NbUpoles ou NbVpoles
1485 // ATTENTION : dans le cas rationnel, un point de coordonnees (x,y,z)
1486 // et de poids w devient un point de coordonnees (wx, wy, wz, w )
1487
1488
1489 if (InU) {
1490 Cdeg = BS->UDegree();
1491 NbP = BS->NbUPoles();
1492 Cdim = BS->NbVPoles() * gap;
1493 }
1494 else {
1495 Cdeg = BS->VDegree();
1496 NbP = BS->NbVPoles();
1497 Cdim = BS->NbUPoles() * gap;
1498 }
1499
1500 // les noeuds plats
1501 Ksize = NbP + Cdeg + 1;
1502 FKnots = new (TColStd_HArray1OfReal) (1,Ksize);
1503 if (InU)
1504 BS->UKnotSequence(FKnots->ChangeArray1());
1505 else
1506 BS->VKnotSequence(FKnots->ChangeArray1());
1507
1508 // le parametre du noeud de raccord
1509 if (After)
1510 Tbord = FKnots->Value(FKnots->Upper()-Cdeg);
1511 else
1512 Tbord = FKnots->Value(FKnots->Lower()+Cdeg);
1513
1514 // les poles
1515 Psize = Cdim * NbP;
1516 Poles = new (TColStd_HArray1OfReal) (1,Psize);
1517
1518 if (InU) {
1519 for (ii=1,ipole=1; ii<=NbP; ii++) {
1520 for (jj=1;jj<=BS->NbVPoles();jj++) {
1521 Poles->SetValue(ipole, BS->Pole(ii,jj).X());
1522 Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y());
1523 Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z());
1524 if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj));
1525 ipole+=gap;
1526 }
1527 }
1528 }
1529 else {
1530 for (jj=1,ipole=1; jj<=NbP; jj++) {
1531 for (ii=1;ii<=BS->NbUPoles();ii++) {
1532 Poles->SetValue(ipole, BS->Pole(ii,jj).X());
1533 Poles->SetValue(ipole+1, BS->Pole(ii,jj).Y());
1534 Poles->SetValue(ipole+2, BS->Pole(ii,jj).Z());
1535 if (rational) Poles->SetValue(ipole+3, BS->Weight(ii,jj));
1536 ipole+=gap;
1537 }
1538 }
1539 }
1540 Padr = (Standard_Real *) &Poles->ChangeValue(1);
1541
1542 // calcul du point de raccord et de la tangente
1543 Point = new (TColStd_HArray1OfReal)(1,Cdim);
1544 Tgte = new (TColStd_HArray1OfReal)(1,Cdim);
1545 lambda = new (TColStd_HArray1OfReal)(1,Cdim);
1546
1547 Standard_Boolean periodic_flag = Standard_False ;
1548 Standard_Integer extrap_mode[2], derivative_request = Max(Continuity,1);
1549 extrap_mode[0] = extrap_mode[1] = Cdeg;
1550 TColStd_Array1OfReal Result(1, Cdim * (derivative_request+1)) ;
1551
1552 TColStd_Array1OfReal& tgte = Tgte->ChangeArray1();
1553 TColStd_Array1OfReal& point = Point->ChangeArray1();
1554 TColStd_Array1OfReal& lamb = lambda->ChangeArray1();
1555
1556 Standard_Real * Radr = (Standard_Real *) &Result(1) ;
1557
1558 BSplCLib::Eval(Tbord,periodic_flag,derivative_request,extrap_mode[0],
1559 Cdeg,FKnots->Array1(),Cdim,*Padr,*Radr);
1560 Ok = Standard_True;
1561 for (ii=1;ii<=Cdim;ii++) {
1562 point(ii) = Result(ii);
1563 tgte(ii) = Result(ii+Cdim);
1564 }
1565
1566 // calcul de la contrainte a atteindre
1567
1568 gp_Vec CurT, OldT;
1569
1570 Standard_Real NTgte, val, Tgtol = 1.e-12, OldN = 0.0;
1571 if (rational) {
1572 for (ii=gap;ii<=Cdim;ii+=gap) {
1573 tgte(ii) = 0.;
1574 }
1575 for (ii=gap;ii<=Cdim;ii+=gap) {
1576 CurT.SetCoord(tgte(ii-3),tgte(ii-2), tgte(ii-1));
1577 NTgte=CurT.Magnitude();
1578 if (NTgte>Tgtol) {
1579 val = Length/NTgte;
1580 // Attentions aux Cas ou le segment donne par les poles
1581 // est oppose au sens de la derive
1582 // Exemple: Certaine portions de tore.
1583 if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) {
1584 Ok = Standard_False;
1585 }
1586
1587 lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = val;
1588 lamb(ii) = 0.;
1589 lambmin = Min(lambmin, val);
1590 }
1591 else {
1592 lamb(ii-1) = lamb(ii-2) = lamb(ii-3) = 0.;
1593 lamb(ii) = 0.;
1594 }
1595 OldT = CurT;
1596 OldN = NTgte;
1597 }
1598 }
1599 else {
1600 for (ii=gap;ii<=Cdim;ii+=gap) {
1601 CurT.SetCoord(tgte(ii-2),tgte(ii-1), tgte(ii));
1602 NTgte=CurT.Magnitude();
1603 if (NTgte>Tgtol) {
1604 val = Length/NTgte;
1605 // Attentions aux Cas ou le segment donne par les poles
1606 // est oppose au sens de la derive
1607 // Exemple: Certaine portion de tore.
1608 if ( (OldN > Tgtol) && (CurT.Angle(OldT) > 2)) {
1609 Ok = Standard_False;
1610 }
1611 lamb(ii) = lamb(ii-1) = lamb(ii-2) = val;
1612 lambmin = Min(lambmin, val);
1613 }
1614 else {
1615 lamb(ii) =lamb(ii-1) = lamb(ii-2) = 0.;
1616 }
1617 OldT = CurT;
1618 OldN = NTgte;
1619 }
1620 }
1621 if (!Ok && Kount<2) {
1622 // On augmente le degre de l'iso bord afin de rapprocher les poles de la surface
1623 // Et on ressaye
1624 if (InU) BS->IncreaseDegree(BS->UDegree(), BS->VDegree()+2);
1625 else BS->IncreaseDegree(BS->UDegree()+2, BS->VDegree());
1626 }
1627 }
1628
1629
1630 TColStd_Array1OfReal ConstraintPoint(1,Cdim);
1631 if (After) {
1632 for (ii=1;ii<=Cdim;ii++) {
1633 ConstraintPoint(ii) = Point->Value(ii) + lambda->Value(ii)*Tgte->Value(ii);
1634 }
1635 }
1636 else {
1637 for (ii=1;ii<=Cdim;ii++) {
1638 ConstraintPoint(ii) = Point->Value(ii) - lambda->Value(ii)*Tgte->Value(ii);
1639 }
1640 }
1641
1642// cas particulier du rationnel
1643 if (rational) {
1644 for (ipole=1;ipole<=Psize;ipole+=gap) {
1645 Poles->ChangeValue(ipole) *= Poles->Value(ipole+3);
1646 Poles->ChangeValue(ipole+1) *= Poles->Value(ipole+3);
1647 Poles->ChangeValue(ipole+2) *= Poles->Value(ipole+3);
1648 }
1649 for (ii=1;ii<=Cdim;ii+=gap) {
1650 ConstraintPoint(ii) *= ConstraintPoint(ii+3);
1651 ConstraintPoint(ii+1) *= ConstraintPoint(ii+3);
1652 ConstraintPoint(ii+2) *= ConstraintPoint(ii+3);
1653 }
1654 }
1655
1656// tableaux necessaires pour l'extension
1d47d8d0 1657 Standard_Integer Ksize2 = Ksize+Cdeg, NbPoles, NbKnots = 0;
7fd59977 1658 TColStd_Array1OfReal FK(1, Ksize2) ;
1659 Standard_Real * FKRadr = &FK(1);
1660
1661 Standard_Integer Psize2 = Psize+Cdeg*Cdim;
1662 TColStd_Array1OfReal PRes(1, Psize2) ;
1663 Standard_Real * PRadr = &PRes(1);
1664 Standard_Real ww;
1665 Standard_Boolean ExtOk = Standard_False;
1666 Handle(TColgp_HArray2OfPnt) NewPoles;
1667 Handle(TColStd_HArray2OfReal) NewWeights;
1668
1669
1670 for (Kount=1; Kount<=5 && !ExtOk; Kount++) {
1671 // extension
1672 BSplCLib::TangExtendToConstraint(FKnots->Array1(),
1673 lambmin,NbP,*Padr,
1674 Cdim,Cdeg,
1675 ConstraintPoint, Cont, After,
1676 NbPoles, NbKnots,*FKRadr, *PRadr);
1677
1678 // recopie des poles du resultat sous forme de points 3D et de poids
1679 Standard_Integer NU, NV, indice ;
1680 if (InU) {
1681 NU = NbPoles;
1682 NV = BS->NbVPoles();
1683 }
1684 else {
1685 NU = BS->NbUPoles();
1686 NV = NbPoles;
1687 }
1688
1689 NewPoles = new (TColgp_HArray2OfPnt)(1,NU,1,NV);
1690 TColgp_Array2OfPnt& NewP = NewPoles->ChangeArray2();
1691 NewWeights = new (TColStd_HArray2OfReal) (1,NU,1,NV);
1692 TColStd_Array2OfReal& NewW = NewWeights->ChangeArray2();
1693
1694 if (!rational) NewW.Init(1.);
1695 NullWeight= Standard_False;
1696
1697 if (InU) {
1698 for (ii=1; ii<=NU && !NullWeight; ii++) {
1699 for (jj=1; jj<=NV && !NullWeight; jj++) {
1700 indice = 1+(ii-1)*Cdim+(jj-1)*gap;
1701 NewP(ii,jj).SetCoord(1,PRes(indice));
1702 NewP(ii,jj).SetCoord(2,PRes(indice+1));
1703 NewP(ii,jj).SetCoord(3,PRes(indice+2));
1704 if (rational) {
1705 ww = PRes(indice+3);
94f71cad 1706 if (Abs(ww - 1.0) < EpsW)
1707 ww = 1.0;
7fd59977 1708 if (ww < EpsW) {
1709 NullWeight = Standard_True;
1710 }
1711 else {
1712 NewW(ii,jj) = ww;
1713 NewP(ii,jj).ChangeCoord() /= ww;
1714 }
1715 }
1716 }
1717 }
1718 }
1719 else {
1720 for (jj=1; jj<=NV && !NullWeight; jj++) {
1721 for (ii=1; ii<=NU && !NullWeight; ii++) {
1722 indice = 1+(ii-1)*gap+(jj-1)*Cdim;
1723 NewP(ii,jj).SetCoord(1,PRes(indice));
1724 NewP(ii,jj).SetCoord(2,PRes(indice+1));
1725 NewP(ii,jj).SetCoord(3,PRes(indice+2));
1726 if (rational) {
1727 ww = PRes(indice+3);
94f71cad 1728 if (Abs(ww - 1.0) < EpsW)
1729 ww = 1.0;
7fd59977 1730 if (ww < EpsW) {
1731 NullWeight = Standard_True;
1732 }
1733 else {
1734 NewW(ii,jj) = ww;
1735 NewP(ii,jj).ChangeCoord() /= ww;
1736 }
1737 }
1738 }
1739 }
1740 }
1741
1742 if (NullWeight) {
0797d9d3 1743#ifdef OCCT_DEBUG
7fd59977 1744 cout << "Echec de l'Extension rationnelle" << endl;
1745#endif
1746 lambmin /= 3.;
1747 NullWeight = Standard_False;
1748 }
1749 else {
1750 ExtOk = Standard_True;
1751 }
1752 }
1753
1754
1755// recopie des noeuds plats sous forme de noeuds avec leurs multiplicites
1756// calcul des degres du resultat
1757 Standard_Integer Usize = BS->NbUKnots(), Vsize = BS->NbVKnots(), UDeg, VDeg;
1758 if (InU)
1759 Usize++;
1760 else
1761 Vsize++;
1762 TColStd_Array1OfReal UKnots(1,Usize);
1763 TColStd_Array1OfReal VKnots(1,Vsize);
1764 TColStd_Array1OfInteger UMults(1,Usize);
1765 TColStd_Array1OfInteger VMults(1,Vsize);
1766 TColStd_Array1OfReal FKRes(1, NbKnots);
1767
1768 for (ii=1; ii<=NbKnots; ii++)
1769 FKRes(ii) = FK(ii);
1770
1771 if (InU) {
1772 BSplCLib::Knots(FKRes, UKnots, UMults);
1773 UDeg = Cdeg;
1774 UMults(Usize) = UDeg+1; // Petite verrue utile quand la continuite
1775 // n'est pas ok.
1776 BS->VKnots(VKnots);
1777 BS->VMultiplicities(VMults);
1778 VDeg = BS->VDegree();
1779 }
1780 else {
1781 BSplCLib::Knots(FKRes, VKnots, VMults);
1782 VDeg = Cdeg;
1783 VMults(Vsize) = VDeg+1;
1784 BS->UKnots(UKnots);
1785 BS->UMultiplicities(UMults);
1786 UDeg = BS->UDegree();
1787 }
1788
1789// construction de la surface BSpline resultat
1790 Handle(Geom_BSplineSurface) Res =
1791 new (Geom_BSplineSurface) (NewPoles->Array2(),
1792 NewWeights->Array2(),
1793 UKnots,VKnots,
1794 UMults,VMults,
1795 UDeg,VDeg,
1796 BS->IsUPeriodic(),
1797 BS->IsVPeriodic());
1798 Surface = Res;
1799}
1800
1801//=======================================================================
1802//function : Inertia
1803//purpose :
1804//=======================================================================
1805void GeomLib::Inertia(const TColgp_Array1OfPnt& Points,
1806 gp_Pnt& Bary,
1807 gp_Dir& XDir,
1808 gp_Dir& YDir,
1809 Standard_Real& Xgap,
1810 Standard_Real& Ygap,
1811 Standard_Real& Zgap)
1812{
1813 gp_XYZ GB(0., 0., 0.), Diff;
1814// gp_Vec A,B,C,D;
1815
1816 Standard_Integer i,nb=Points.Length();
1817 GB.SetCoord(0.,0.,0.);
1818 for (i=1; i<=nb; i++)
1819 GB += Points(i).XYZ();
1820
1821 GB /= nb;
1822
1823 math_Matrix M (1, 3, 1, 3);
1824 M.Init(0.);
1825 for (i=1; i<=nb; i++) {
1826 Diff.SetLinearForm(-1, Points(i).XYZ(), GB);
1827 M(1,1) += Diff.X() * Diff.X();
1828 M(2,2) += Diff.Y() * Diff.Y();
1829 M(3,3) += Diff.Z() * Diff.Z();
1830 M(1,2) += Diff.X() * Diff.Y();
1831 M(1,3) += Diff.X() * Diff.Z();
1832 M(2,3) += Diff.Y() * Diff.Z();
1833 }
1834
1835 M(2,1)=M(1,2) ;
1836 M(3,1)=M(1,3) ;
1837 M(3,2)=M(2,3) ;
1838
1839 M /= nb;
1840
1841 math_Jacobi J(M);
1842 if (!J.IsDone()) {
0797d9d3 1843#ifdef OCCT_DEBUG
7fd59977 1844 cout << "Erreur dans Jacobbi" << endl;
1845 M.Dump(cout);
1846#endif
1847 }
1848
1849 Standard_Real n1,n2,n3;
1850
1851 n1=J.Value(1);
1852 n2=J.Value(2);
1853 n3=J.Value(3);
1854
1855 Standard_Real r1 = Min(Min(n1,n2),n3), r2;
1856 Standard_Integer m1, m2, m3;
1857 if (r1==n1) {
1858 m1 = 1;
1859 r2 = Min(n2,n3);
1860 if (r2==n2) {
1861 m2 = 2;
1862 m3 = 3;
1863 }
1864 else {
1865 m2 = 3;
1866 m3 = 2;
1867 }
1868 }
1869 else {
1870 if (r1==n2) {
1871 m1 = 2 ;
1872 r2 = Min(n1,n3);
1873 if (r2==n1) {
1874 m2 = 1;
1875 m3 = 3;
1876 }
1877 else {
1878 m2 = 3;
1879 m3 = 1;
1880 }
1881 }
1882 else {
1883 m1 = 3 ;
1884 r2 = Min(n1,n2);
1885 if (r2==n1) {
1886 m2 = 1;
1887 m3 = 2;
1888 }
1889 else {
1890 m2 = 2;
1891 m3 = 1;
1892 }
1893 }
1894 }
1895
1896 math_Vector V2(1,3),V3(1,3);
1897 J.Vector(m2,V2);
1898 J.Vector(m3,V3);
1899
1900 Bary.SetXYZ(GB);
1901 XDir.SetCoord(V3(1),V3(2),V3(3));
1902 YDir.SetCoord(V2(1),V2(2),V2(3));
1903
1904 Zgap = sqrt(Abs(J.Value(m1)));
1905 Ygap = sqrt(Abs(J.Value(m2)));
1906 Xgap = sqrt(Abs(J.Value(m3)));
1907}
1908//=======================================================================
1909//function : AxeOfInertia
1910//purpose :
1911//=======================================================================
1912void GeomLib::AxeOfInertia(const TColgp_Array1OfPnt& Points,
1913 gp_Ax2& Axe,
1914 Standard_Boolean& IsSingular,
1915 const Standard_Real Tol)
1916{
1917 gp_Pnt Bary;
1918 gp_Dir OX,OY,OZ;
1919 Standard_Real gx, gy, gz;
1920
1921 GeomLib::Inertia(Points, Bary, OX, OY, gx, gy, gz);
1922
1923 if (gy*Points.Length()<=Tol) {
1924 gp_Ax2 axe (Bary, OX);
1925 OY = axe.XDirection();
1926 IsSingular = Standard_True;
1927 }
1928 else {
1929 IsSingular = Standard_False;
1930 }
1931
1932 OZ = OX^OY;
1933 gp_Ax2 TheAxe(Bary, OZ, OX);
1934 Axe = TheAxe;
1935}
1936
1937//=======================================================================
1938//function : CanBeTreated
1939//purpose : indicates if the surface can be treated(if the conditions are
1940// filled) and need to be treated(if the surface hasn't been yet
1941// treated or if the surface is rationnal and non periodic)
1942//=======================================================================
1943
1944static Standard_Boolean CanBeTreated(Handle(Geom_BSplineSurface)& BSurf)
1945
1946{Standard_Integer i;
1947 Standard_Real lambda; //proportionnality coefficient
1948 Standard_Boolean AlreadyTreated=Standard_True;
1949
1950 if (!BSurf->IsURational()||(BSurf->IsUPeriodic()))
1951 return Standard_False;
1952 else {
1953 lambda=(BSurf->Weight(1,1)/BSurf->Weight(BSurf->NbUPoles(),1));
1954 for (i=1;i<=BSurf->NbVPoles();i++) //test of the proportionnality of the denominator on the boundaries
1955 if ((BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))<(1-Precision::Confusion()))||
1956 (BSurf->Weight(1,i)/(lambda*BSurf->Weight(BSurf->NbUPoles(),i))>(1+Precision::Confusion())))
1957 return Standard_False;
1958 i=1;
1959 while ((AlreadyTreated) && (i<=BSurf->NbVPoles())){ //tests if the surface has already been treated
1960 if (((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))<(1-Precision::Confusion()))||
1961 ((BSurf->Weight(1,i)/(BSurf->Weight(2,i)))>(1+Precision::Confusion()))||
1962 ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))<(1-Precision::Confusion()))||
1963 ((BSurf->Weight(BSurf->NbUPoles()-1,i)/(BSurf->Weight(BSurf->NbUPoles(),i)))>(1+Precision::Confusion())))
1964 AlreadyTreated=Standard_False;
1965 i++;
1966 }
1967 if (AlreadyTreated)
1968 return Standard_False;
1969 }
1970 return Standard_True;
1971}
1972
1973//=======================================================================
41194117
K		 1974//class : law_evaluator
1975//purpose : usefull to estimate the value of a function of 2 variables
7fd59977 1976//=======================================================================
1977
41194117
K		 1978class law_evaluator : public BSplSLib_EvaluatorFunction
1979{
7fd59977 1980
41194117 1981public:
7fd59977 1982
41194117
K		 1983 law_evaluator (const GeomLib_DenominatorMultiplierPtr theDenominatorPtr)
1984 : myDenominator (theDenominatorPtr) {}
1985
1986 virtual void Evaluate (const Standard_Integer theDerivativeRequest,
1987 const Standard_Real theUParameter,
1988 const Standard_Real theVParameter,
1989 Standard_Real& theResult,
1990 Standard_Integer& theErrorCode) const
1991 {
1992 if ((myDenominator != NULL) && (theDerivativeRequest == 0))
1993 {
1994 theResult = myDenominator->Value (theUParameter, theVParameter);
1995 theErrorCode = 0;
1996 }
1997 else
1998 {
1999 theErrorCode = 1;
2000 }
7fd59977 2001 }
41194117
K		 2002
2003private:
2004
2005 GeomLib_DenominatorMultiplierPtr myDenominator;
2006
2007};
2008
7fd59977 2009//=======================================================================
2010//function : CheckIfKnotExists
2011//purpose : true if the knot already exists in the knot sequence
2012//=======================================================================
2013
2014static Standard_Boolean CheckIfKnotExists(const TColStd_Array1OfReal& surface_knots,
2015 const Standard_Real knot)
2016
2017{Standard_Integer i;
2018 for (i=1;i<=surface_knots.Length();i++)
2019 if ((surface_knots(i)-Precision::Confusion()<=knot)&&(surface_knots(i)+Precision::Confusion()>=knot))
2020 return Standard_True;
2021 return Standard_False;
2022}
2023
2024//=======================================================================
2025//function : AddAKnot
2026//purpose : add a knot and its multiplicity to the knot sequence. This knot
2027// will be C2 and the degree is increased of deltasurface_degree
2028//=======================================================================
2029
2030static void AddAKnot(const TColStd_Array1OfReal& knots,
2031 const TColStd_Array1OfInteger& mults,
2032 const Standard_Real knotinserted,
2033 const Standard_Integer deltasurface_degree,
2034 const Standard_Integer finalsurfacedegree,
2035 Handle(TColStd_HArray1OfReal) & newknots,
2036 Handle(TColStd_HArray1OfInteger) & newmults)
2037
2038{Standard_Integer i;
2039
2040 newknots=new TColStd_HArray1OfReal(1,knots.Length()+1);
2041 newmults=new TColStd_HArray1OfInteger(1,knots.Length()+1);
2042 i=1;
2043 while (knots(i)<knotinserted){
2044 newknots->SetValue(i,knots(i));
2045 newmults->SetValue(i,mults(i)+deltasurface_degree);
2046 i++;
2047 }
2048 newknots->SetValue(i,knotinserted); //insertion of the new knot
2049 newmults->SetValue(i,finalsurfacedegree-2);
2050 i++;
2051 while (i<=newknots->Length()){
2052 newknots->SetValue(i,knots(i-1));
2053 newmults->SetValue(i,mults(i-1)+deltasurface_degree);
2054 i++;
2055 }
2056}
2057
2058//=======================================================================
2059//function : Sort
2060//purpose : give the new flat knots(u or v) of the surface
2061//=======================================================================
2062
2063static void BuildFlatKnot(const TColStd_Array1OfReal& surface_knots,
2064 const TColStd_Array1OfInteger& surface_mults,
2065 const Standard_Integer deltasurface_degree,
2066 const Standard_Integer finalsurface_degree,
2067 const Standard_Real knotmin,
2068 const Standard_Real knotmax,
2069 Handle(TColStd_HArray1OfReal)& ResultKnots,
2070 Handle(TColStd_HArray1OfInteger)& ResultMults)
2071
2072{
2073 Standard_Integer i;
2074
2075 if (CheckIfKnotExists(surface_knots,knotmin) &&
2076 CheckIfKnotExists(surface_knots,knotmax)){
2077 ResultKnots=new TColStd_HArray1OfReal(1,surface_knots.Length());
2078 ResultMults=new TColStd_HArray1OfInteger(1,surface_knots.Length());
2079 for (i=1;i<=surface_knots.Length();i++){
2080 ResultKnots->SetValue(i,surface_knots(i));
2081 ResultMults->SetValue(i,surface_mults(i)+deltasurface_degree);
2082 }
2083 }
2084 else{
2085 if ((CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax)))
2086 AddAKnot(surface_knots,surface_mults,knotmax,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2087 else{
2088 if ((!CheckIfKnotExists(surface_knots,knotmin))&&(CheckIfKnotExists(surface_knots,knotmax)))
2089 AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2090 else{
2091 if ((!CheckIfKnotExists(surface_knots,knotmin))&&(!CheckIfKnotExists(surface_knots,knotmax))&&
2092 (knotmin==knotmax)){
2093 AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,ResultKnots,ResultMults);
2094 }
2095 else{
2096 Handle(TColStd_HArray1OfReal) IntermedKnots;
2097 Handle(TColStd_HArray1OfInteger) IntermedMults;
2098 AddAKnot(surface_knots,surface_mults,knotmin,deltasurface_degree,finalsurface_degree,IntermedKnots,IntermedMults);
2099 AddAKnot(IntermedKnots->ChangeArray1(),IntermedMults->ChangeArray1(),knotmax,0,finalsurface_degree,ResultKnots,ResultMults);
2100 }
2101 }
2102 }
2103 }
2104}
2105
2106//=======================================================================
2107//function : FunctionMultiply
2108//purpose : multiply the surface BSurf by a(u,v) (law_evaluator) on its
2109// numerator and denominator
2110//=======================================================================
2111
2112static void FunctionMultiply(Handle(Geom_BSplineSurface)& BSurf,
2113 const Standard_Real knotmin,
2114 const Standard_Real knotmax)
2115
2116{TColStd_Array1OfReal surface_u_knots(1,BSurf->NbUKnots()) ;
2117 TColStd_Array1OfInteger surface_u_mults(1,BSurf->NbUKnots()) ;
2118 TColStd_Array1OfReal surface_v_knots(1,BSurf->NbVKnots()) ;
2119 TColStd_Array1OfInteger surface_v_mults(1,BSurf->NbVKnots()) ;
2120 TColgp_Array2OfPnt surface_poles(1,BSurf->NbUPoles(),
2121 1,BSurf->NbVPoles()) ;
2122 TColStd_Array2OfReal surface_weights(1,BSurf->NbUPoles(),
2123 1,BSurf->NbVPoles()) ;
2124 Standard_Integer i,j,k,status,new_num_u_poles,new_num_v_poles,length=0;
2125 Handle(TColStd_HArray1OfReal) newuknots,newvknots;
2126 Handle(TColStd_HArray1OfInteger) newumults,newvmults;
2127
2128 BSurf->UKnots(surface_u_knots) ;
2129 BSurf->UMultiplicities(surface_u_mults) ;
2130 BSurf->VKnots(surface_v_knots) ;
2131 BSurf->VMultiplicities(surface_v_mults) ;
2132 BSurf->Poles(surface_poles) ;
2133 BSurf->Weights(surface_weights) ;
2134
2135 TColStd_Array1OfReal Knots(1,2);
2136 TColStd_Array1OfInteger Mults(1,2);
2137 Handle(TColStd_HArray1OfReal) NewKnots;
2138 Handle(TColStd_HArray1OfInteger) NewMults;
2139
2140 Knots(1)=0;
2141 Knots(2)=1;
2142 Mults(1)=4;
2143 Mults(2)=4;
2144 BuildFlatKnot(Knots,Mults,0,3,knotmin,knotmax,NewKnots,NewMults);
2145
2146 for (i=1;i<=NewMults->Length();i++)
2147 length+=NewMults->Value(i);
2148 TColStd_Array1OfReal FlatKnots(1,length);
2149 BSplCLib::KnotSequence(NewKnots->ChangeArray1(),NewMults->ChangeArray1(),FlatKnots);
2150
41194117 2151 GeomLib_DenominatorMultiplier aDenominator (BSurf, FlatKnots);
7fd59977 2152
2153 BuildFlatKnot(surface_u_knots,
2154 surface_u_mults,
2155 3,
2156 BSurf->UDegree()+3,
2157 knotmin,
2158 knotmax,
2159 newuknots,
2160 newumults);
2161 BuildFlatKnot(surface_v_knots,
2162 surface_v_mults,
2163 BSurf->VDegree(),
2164 2*(BSurf->VDegree()),
2165 1.0,
2166 0.0,
2167 newvknots,
2168 newvmults);
2169 length=0;
2170 for (i=1;i<=newumults->Length();i++)
2171 length+=newumults->Value(i);
2172 new_num_u_poles=(length-BSurf->UDegree()-3-1);
2173 TColStd_Array1OfReal newuflatknots(1,length);
2174 length=0;
2175 for (i=1;i<=newvmults->Length();i++)
2176 length+=newvmults->Value(i);
2177 new_num_v_poles=(length-2*BSurf->VDegree()-1);
2178 TColStd_Array1OfReal newvflatknots(1,length);
2179
2180 TColgp_Array2OfPnt NewNumerator(1,new_num_u_poles,1,new_num_v_poles);
2181 TColStd_Array2OfReal NewDenominator(1,new_num_u_poles,1,new_num_v_poles);
2182
2183 BSplCLib::KnotSequence(newuknots->ChangeArray1(),newumults->ChangeArray1(),newuflatknots);
2184 BSplCLib::KnotSequence(newvknots->ChangeArray1(),newvmults->ChangeArray1(),newvflatknots);
2185//POP pour WNT
41194117 2186 law_evaluator ev (&aDenominator);
7fd59977 2187// BSplSLib::FunctionMultiply(law_evaluator, //multiplication
2188 BSplSLib::FunctionMultiply(ev, //multiplication
2189 BSurf->UDegree(),
2190 BSurf->VDegree(),
2191 surface_u_knots,
2192 surface_v_knots,
2193 surface_u_mults,
2194 surface_v_mults,
2195 surface_poles,
2196 surface_weights,
2197 newuflatknots,
2198 newvflatknots,
2199 BSurf->UDegree()+3,
2200 2*(BSurf->VDegree()),
2201 NewNumerator,
2202 NewDenominator,
2203 status);
2204 if (status!=0)
2205 Standard_ConstructionError::Raise("GeomLib Multiplication Error") ;
2206 for (i = 1 ; i <= new_num_u_poles ; i++) {
2207 for (j = 1 ; j <= new_num_v_poles ; j++) {
2208 for (k = 1 ; k <= 3 ; k++) {
2209 NewNumerator(i,j).SetCoord(k,NewNumerator(i,j).Coord(k)/NewDenominator(i,j)) ;
2210 }
2211 }
2212 }
2213 BSurf= new Geom_BSplineSurface(NewNumerator,
2214 NewDenominator,
2215 newuknots->ChangeArray1(),
2216 newvknots->ChangeArray1(),
2217 newumults->ChangeArray1(),
2218 newvmults->ChangeArray1(),
2219 BSurf->UDegree()+3,
2220 2*(BSurf->VDegree()) );
2221}
2222
2223//=======================================================================
2224//function : CancelDenominatorDerivative1D
2225//purpose : cancel the denominator derivative in one direction
2226//=======================================================================
2227
2228static void CancelDenominatorDerivative1D(Handle(Geom_BSplineSurface) & BSurf)
2229
2230{Standard_Integer i,j;
2231 Standard_Real uknotmin=1.0,uknotmax=0.0,
2232 x,y,
2233 startu_value,
2234 endu_value;
2235 TColStd_Array1OfReal BSurf_u_knots(1,BSurf->NbUKnots()) ;
2236
2237 startu_value=BSurf->UKnot(1);
2238 endu_value=BSurf->UKnot(BSurf->NbUKnots());
2239 BSurf->UKnots(BSurf_u_knots) ;
2240 BSplCLib::Reparametrize(0.0,1.0,BSurf_u_knots);
2241 BSurf->SetUKnots(BSurf_u_knots); //reparametrisation of the surface
2242 Handle(Geom_BSplineCurve) BCurve;
2243 TColStd_Array1OfReal BCurveWeights(1,BSurf->NbUPoles());
2244 TColgp_Array1OfPnt BCurvePoles(1,BSurf->NbUPoles());
2245 TColStd_Array1OfReal BCurveKnots(1,BSurf->NbUKnots());
2246 TColStd_Array1OfInteger BCurveMults(1,BSurf->NbUKnots());
2247
2248 if (CanBeTreated(BSurf)){
2249 for (i=1;i<=BSurf->NbVPoles();i++){ //loop on each pole function
2250 x=1.0;y=0.0;
2251 for (j=1;j<=BSurf->NbUPoles();j++){
2252 BCurveWeights(j)=BSurf->Weight(j,i);
2253 BCurvePoles(j)=BSurf->Pole(j,i);
2254 }
2255 BSurf->UKnots(BCurveKnots);
2256 BSurf->UMultiplicities(BCurveMults);
2257 BCurve = new Geom_BSplineCurve(BCurvePoles, //building of a pole function
2258 BCurveWeights,
2259 BCurveKnots,
2260 BCurveMults,
2261 BSurf->UDegree());
2262 Hermit::Solutionbis(BCurve,x,y,Precision::Confusion(),Precision::Confusion());
2263 if (x<uknotmin)
2264 uknotmin=x; //uknotmin,uknotmax:extremal knots
2265 if ((x!=1.0)&&(x>uknotmax))
2266 uknotmax=x;
2267 if ((y!=0.0)&&(y<uknotmin))
2268 uknotmin=y;
2269 if (y>uknotmax)
2270 uknotmax=y;
2271 }
2272
2273 FunctionMultiply(BSurf,uknotmin,uknotmax); //multiplication
2274
2275 BSurf->UKnots(BSurf_u_knots) ;
2276 BSplCLib::Reparametrize(startu_value,endu_value,BSurf_u_knots);
2277 BSurf->SetUKnots(BSurf_u_knots);
2278 }
2279}
2280
2281//=======================================================================
2282//function : CancelDenominatorDerivative
2283//purpose :
2284//=======================================================================
2285
2286void GeomLib::CancelDenominatorDerivative(Handle(Geom_BSplineSurface) & BSurf,
2287 const Standard_Boolean udirection,
2288 const Standard_Boolean vdirection)
2289
2290{if (udirection && !vdirection)
2291 CancelDenominatorDerivative1D(BSurf);
2292 else{
2293 if (!udirection && vdirection) {
2294 BSurf->ExchangeUV();
2295 CancelDenominatorDerivative1D(BSurf);
2296 BSurf->ExchangeUV();
2297 }
2298 else{
2299 if (udirection && vdirection){ //optimize the treatment
2300 if (BSurf->UDegree()<=BSurf->VDegree()){
2301 CancelDenominatorDerivative1D(BSurf);
2302 BSurf->ExchangeUV();
2303 CancelDenominatorDerivative1D(BSurf);
2304 BSurf->ExchangeUV();
2305 }
2306 else{
2307 BSurf->ExchangeUV();
2308 CancelDenominatorDerivative1D(BSurf);
2309 BSurf->ExchangeUV();
2310 CancelDenominatorDerivative1D(BSurf);
2311 }
2312 }
2313 }
2314 }
2315}
2316
2317//=======================================================================
2318//function : NormEstim
2319//purpose :
2320//=======================================================================
2321
2322Standard_Integer GeomLib::NormEstim(const Handle(Geom_Surface)& S,
2323 const gp_Pnt2d& UV,
2324 const Standard_Real Tol, gp_Dir& N)
2325{
2326 gp_Vec DU, DV;
2327 gp_Pnt DummyPnt;
2328 Standard_Real aTol2 = Square(Tol);
2329
2330 S->D1(UV.X(), UV.Y(), DummyPnt, DU, DV);
2331
2332 Standard_Real MDU = DU.SquareMagnitude(), MDV = DV.SquareMagnitude();
2333
7fd59977 2334 if(MDU >= aTol2 && MDV >= aTol2) {
2335 gp_Vec Norm = DU^DV;
2336 Standard_Real Magn = Norm.SquareMagnitude();
2337 if(Magn < aTol2) return 3;
2338
2339 //Magn = sqrt(Magn);
2340 N.SetXYZ(Norm.XYZ());
2341
2342 return 0;
2343 }
7fd59977 2344 else {
2b21c641 2345 gp_Vec D2U, D2V, D2UV;
2346 Standard_Boolean isDone;
2347 CSLib_NormalStatus aStatus;
2348 gp_Dir aNormal;
2349
2350 S->D2(UV.X(), UV.Y(), DummyPnt, DU, DV, D2U, D2V, D2UV);
2351 CSLib::Normal(DU, DV, D2U, D2V, D2UV, Tol, isDone, aStatus, aNormal);
2352
2353 if (isDone) {
2354 Standard_Real Umin, Umax, Vmin, Vmax;
2355 Standard_Real step = 1.0e-5;
2356 Standard_Real eps = 1.0e-16;
23b894f7 2357 Standard_Real sign = -1.0;
2b21c641 2358
2359 S->Bounds(Umin, Umax, Vmin, Vmax);
23b894f7 2360
2361 // check for cone apex singularity point
2362 if ((UV.Y() > Vmin + step) && (UV.Y() < Vmax - step))
2363 {
2364 gp_Dir aNormal1, aNormal2;
2365 Standard_Real aConeSingularityAngleEps = 1.0e-4;
2366 S->D1(UV.X(), UV.Y() - sign * step, DummyPnt, DU, DV);
2367 if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) {
2368 aNormal1 = DU^DV;
2369 S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV);
2370 if ((DU.XYZ().SquareModulus() > eps) && (DV.XYZ().SquareModulus() > eps)) {
2371 aNormal2 = DU^DV;
2372 if (aNormal1.IsOpposite(aNormal2, aConeSingularityAngleEps))
2373 return 2;
2374 }
2375 }
2376 }
2377
2b21c641 2378 // Along V
2379 if(MDU < aTol2 && MDV >= aTol2) {
23b894f7 2380 if ((Vmax - UV.Y()) > (UV.Y() - Vmin))
2381 sign = 1.0;
2b21c641 2382 S->D1(UV.X(), UV.Y() + sign * step, DummyPnt, DU, DV);
2383 gp_Vec Norm = DU^DV;
23b894f7 2384 if (Norm.SquareMagnitude() < eps) {
2385 Standard_Real sign1 = -1.0;
2386 if ((Umax - UV.X()) > (UV.X() - Umin))
2387 sign1 = 1.0;
2388 S->D1(UV.X() + sign1 * step, UV.Y() + sign * step, DummyPnt, DU, DV);
2389 Norm = DU^DV;
2390 }
2b21c641 2391 if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0))
23b894f7 2392 aNormal.Reverse();
2b21c641 2393 }
23b894f7 2394
2b21c641 2395 // Along U
2396 if(MDV < aTol2 && MDU >= aTol2) {
23b894f7 2397 if ((Umax - UV.X()) > (UV.X() - Umin))
2398 sign = 1.0;
2b21c641 2399 S->D1(UV.X() + sign * step, UV.Y(), DummyPnt, DU, DV);
2400 gp_Vec Norm = DU^DV;
23b894f7 2401 if (Norm.SquareMagnitude() < eps) {
2402 Standard_Real sign1 = -1.0;
2403 if ((Vmax - UV.Y()) > (UV.Y() - Vmin))
2404 sign1 = 1.0;
2405 S->D1(UV.X() + sign * step, UV.Y() + sign1 * step, DummyPnt, DU, DV);
2406 Norm = DU^DV;
2407 }
2b21c641 2408 if ((Norm.SquareMagnitude() >= eps) && (Norm.Dot(aNormal) < 0.0))
2409 aNormal.Reverse();
2410 }
7fd59977 2411
2b21c641 2412 // quasysingular
2413 if ((aStatus == CSLib_D1NuIsNull) || (aStatus == CSLib_D1NvIsNull) ||
2414 (aStatus == CSLib_D1NuIsParallelD1Nv)) {
2415 N.SetXYZ(aNormal.XYZ());
2416 return 1;
2417 }
2418 // conical
2419 if (aStatus == CSLib_InfinityOfSolutions)
2420 return 2;
7fd59977 2421 }
2b21c641 2422 // computation is impossible
7fd59977 2423 else {
2b21c641 2424 // conical
2425 if (aStatus == CSLib_D1NIsNull) {
2426 return 2;
2427 }
2428 return 3;
7fd59977 2429 }
7fd59977 2430 }
2b21c641 2431 return 3;
7fd59977 2432}
2433
2434
Open CASCADE Technology repository