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