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