0022922: Clean up warnings on uninitialized / unused variables
[occt.git] / src / GeomPlate / GeomPlate_BuildPlateSurface.cxx
CommitLineData
7fd59977 1// File: GeomPlate_BuildPlateSurface.cxx
2// Created: Mon May 5 17:07:03 1997
3// Author: Jerome LEMONIER
4// <jlr@sgi64>
5// Modification de l'interface
6// Amelioration de l'aglo de remplissage
7// 29-09-97 ; jct; correction du calcul des points doublons (tol 2d et non 3d)
8// 16-10-97 ; jct; on remet un Raise au lieu du cout sous debug (sinon plantage dans k1fab)
9// 03-11-97 ; jct; plus de reference a BRepAdaptor et BRepTools
10
11#include <stdio.h>
12
13#include <GeomPlate_BuildPlateSurface.ixx>
14
15#include <TColStd_HArray1OfReal.hxx>
16#include <TColgp_HArray1OfPnt.hxx>
17
18#include <gp_Pnt.hxx>
19#include <gp_Pnt2d.hxx>
20#include <gp_Vec.hxx>
21#include <gp_Vec2d.hxx>
22
23#include <Geom_RectangularTrimmedSurface.hxx>
24
25#include <GCPnts_AbscissaPoint.hxx>
26#include <Law_Interpol.hxx>
27
28#include <Plate_PinpointConstraint.hxx>
29#include <Plate_Plate.hxx>
30#include <Plate_GtoCConstraint.hxx>
31#include <Plate_D1.hxx>
32#include <Plate_D2.hxx>
33#include <Plate_FreeGtoCConstraint.hxx>
34
35#include <Precision.hxx>
36
37#include <GeomPlate_BuildAveragePlane.hxx>
38#include <GeomPlate_HArray1OfSequenceOfReal.hxx>
39// pour la verif G2
40#include <LocalAnalysis_SurfaceContinuity.hxx>
41
42// pour les intersection entre les courbes
43#include <Geom2dInt_GInter.hxx>
44#include <IntRes2d_IntersectionPoint.hxx>
45#include <Geom2dAdaptor_Curve.hxx>
46#include <GeomAdaptor_Surface.hxx>
47
48// projection
49#include <Extrema_ExtPS.hxx>
50#include <Extrema_POnSurf.hxx>
51#include <ProjLib_CompProjectedCurve.hxx>
52#include <ProjLib_HCompProjectedCurve.hxx>
53#include <Approx_CurveOnSurface.hxx>
54#include <GeomAdaptor.hxx>
55#include <GeomAdaptor_HSurface.hxx>
56#include <GeomLProp_SLProps.hxx>
57#include <Geom2d_BezierCurve.hxx>
58#include <TColgp_Array1OfPnt2d.hxx>
59
60// pour mes tests
61#ifdef DEB
62#include <OSD_Chronometer.hxx>
63
64static Standard_Integer Affich=0;
65// 0 : Pas de display
66// 1 : Display des Geometries et controle intermediaire
67// 2 : Display du nombre de contrainte par courbe + Intersection
68// 3 : Dump des contraintes dans Plate
69static Standard_Integer NbPlan = 0;
70//static Standard_Integer NbCurv2d = 0;
71static Standard_Integer NbMark = 0;
72static Standard_Integer NbProj = 0;
73#endif
74
75#ifdef DRAW
76#include <DrawTrSurf.hxx>
77#include <Draw_Marker3D.hxx>
78#include <Draw_Marker2D.hxx>
79#include <Draw.hxx>
80#endif
81
82#include <TColStd_SequenceOfInteger.hxx>
83#include <TColgp_SequenceOfVec.hxx>
84#include <TColgp_HArray2OfPnt.hxx>
85#include <Geom_BSplineSurface.hxx>
86#include <GeomPlate_SequenceOfAij.hxx>
87#include <GeomPlate_MakeApprox.hxx>
88
89
90//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
91// =========================================================
92// C O N S T R U C T E U R S
93// =========================================================
94//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
95
96//---------------------------------------------------------
97// Constructeur compatible avec l'ancienne version
98//---------------------------------------------------------
99GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface (
100 const Handle(TColStd_HArray1OfInteger)& NPoints,
101 const Handle(GeomPlate_HArray1OfHCurveOnSurface)& TabCurve,
102 const Handle(TColStd_HArray1OfInteger)& Tang,
103 const Standard_Integer Degree,
104 const Standard_Integer NbIter,
105 const Standard_Real Tol2d,
106 const Standard_Real Tol3d,
107 const Standard_Real TolAng,
108 const Standard_Real TolCurv,
109 const Standard_Boolean Anisotropie
110) :
111myAnisotropie(Anisotropie),
112myDegree(Degree),
113myNbIter(NbIter),
114myProj(),
115myTol2d(Tol2d),
116myTol3d(Tol3d),
117myTolAng(TolAng),
118myTolCurv(TolCurv),
119myNbBounds(0)
120{ Standard_Integer NTCurve=TabCurve->Length();// Nombre de contraintes lineaires
121 myNbPtsOnCur = 0; // Debrayage du calcul du nombre de points
122 // en fonction de la longueur
123 myLinCont = new GeomPlate_HSequenceOfCurveConstraint;
124 myPntCont = new GeomPlate_HSequenceOfPointConstraint;
125 if (myNbIter<1)
126 Standard_ConstructionError::Raise("GeomPlate : Number of iteration must be >= 1");
127
128 if (NTCurve==0)
129 Standard_ConstructionError::Raise("GeomPlate : the bounds Array is null");
130 if (Tang->Length()==0)
131 Standard_ConstructionError::Raise("GeomPlate : the constraints Array is null");
132 Standard_Integer nbp = 0;
133 Standard_Integer i ;
134 for ( i=1;i<=NTCurve;i++)
135 { nbp+=NPoints->Value(i);
136 }
137 if (nbp==0)
138 Standard_ConstructionError::Raise("GeomPlate : the resolution is impossible if the number of constraints points is 0");
139 if (myDegree<2)
140 Standard_ConstructionError::Raise("GeomPlate ; the degree resolution must be upper of 2");
141 // Remplissage des champs passage de l'ancien constructeur au nouveau
142 for(i=1;i<=NTCurve;i++)
143 { Handle(GeomPlate_CurveConstraint) Cont = new GeomPlate_CurveConstraint
144 ( TabCurve->Value(i),
145 Tang->Value(i),
146 NPoints->Value(i));
147 myLinCont->Append(Cont);
148 }
149 mySurfInitIsGive=Standard_False;
150
151 myIsLinear = Standard_True;
152 myFree = Standard_False;
153}
154
155//------------------------------------------------------------------
156// Constructeur avec surface d'init et degre de resolution de plate
157//------------------------------------------------------------------
158GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface (
159 const Handle(Geom_Surface)& Surf,
160 const Standard_Integer Degree,
161 const Standard_Integer NbPtsOnCur,
162 const Standard_Integer NbIter,
163 const Standard_Real Tol2d,
164 const Standard_Real Tol3d,
165 const Standard_Real TolAng,
166 const Standard_Real TolCurv,
167 const Standard_Boolean Anisotropie ) :
168mySurfInit(Surf),
169myAnisotropie(Anisotropie),
170myDegree(Degree),
171myNbPtsOnCur(NbPtsOnCur),
172myNbIter(NbIter),
173myProj(),
174myTol2d(Tol2d),
175myTol3d(Tol3d),
176myTolAng(TolAng),
177myTolCurv(TolCurv),
178myNbBounds(0)
179{ if (myNbIter<1)
180 Standard_ConstructionError::Raise("GeomPlate : Number of iteration must be >= 1");
181 if (myDegree<2)
182 Standard_ConstructionError::Raise("GeomPlate : the degree resolution must be upper of 2");
183 myLinCont = new GeomPlate_HSequenceOfCurveConstraint;
184 myPntCont = new GeomPlate_HSequenceOfPointConstraint;
185 mySurfInitIsGive=Standard_True;
186
187 myIsLinear = Standard_True;
188 myFree = Standard_False;
189}
190
191
192//---------------------------------------------------------
193// Constructeur avec degre de resolution de plate
194//---------------------------------------------------------
195GeomPlate_BuildPlateSurface::GeomPlate_BuildPlateSurface (
196 const Standard_Integer Degree,
197 const Standard_Integer NbPtsOnCur,
198 const Standard_Integer NbIter,
199 const Standard_Real Tol2d,
200 const Standard_Real Tol3d,
201 const Standard_Real TolAng,
202 const Standard_Real TolCurv,
203 const Standard_Boolean Anisotropie ) :
204myAnisotropie(Anisotropie),
205myDegree(Degree),
206myNbPtsOnCur(NbPtsOnCur),
207myNbIter(NbIter),
208myProj(),
209myTol2d(Tol2d),
210myTol3d(Tol3d),
211myTolAng(TolAng),
212myTolCurv(TolCurv),
213myNbBounds(0)
214{ if (myNbIter<1)
215 Standard_ConstructionError::Raise("GeomPlate : Number of iteration must be >= 1");
216 if (myDegree<2)
217 Standard_ConstructionError::Raise("GeomPlate : the degree resolution must be upper of 2");
218 myLinCont = new GeomPlate_HSequenceOfCurveConstraint;
219 myPntCont = new GeomPlate_HSequenceOfPointConstraint;
220 mySurfInitIsGive=Standard_False;
221
222 myIsLinear = Standard_True;
223 myFree = Standard_False;
224}
225
226//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
227// =========================================================
228// M E T H O D E S P U B L I Q U E S
229// =========================================================
230//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
231
232
233
234//-------------------------------------------------------------------------
235// Fonction : TrierTab, Fonction interne, ne fait partie de la classe
236//-------------------------------------------------------------------------
237// Reordonne le tableau des permutations
238// Apres l'appel a CourbeJointive l'ordre des courbes est modifie
239// Ex : ordre init des courbe ==> A B C D E F
240// Dans TabInit on note ==> 1 2 3 4 5 6
241// apres CourbeJointive ==> A E C B D F
242// TabInit ==> 1 5 3 2 4 6
243// apres TrierTab le Tableau contient ==> 1 4 3 5 2 6
244// On peut ainsi acceder a la 2eme courbe en prenant TabInit[2]
245// c'est a dire la 4eme du tableau des courbes classees
246//-------------------------------------------------------------------------
247void TrierTab(Handle(TColStd_HArray1OfInteger)& Tab)
248{ // Trie le tableau des permutations pour retrouver l'ordre initial
249 Standard_Integer Nb=Tab->Length();
250 TColStd_Array1OfInteger TabTri(1,Nb);
251 for (Standard_Integer i=1;i<=Nb;i++)
252 TabTri.SetValue(Tab->Value(i),i);
253 Tab->ChangeArray1()=TabTri;
254}
255
256//---------------------------------------------------------
257// Fonction : ProjectCurve
258//---------------------------------------------------------
259Handle(Geom2d_Curve) GeomPlate_BuildPlateSurface::ProjectCurve(const Handle(Adaptor3d_HCurve)& Curv)
260{ // Projection d'une courbe sur un plan
261 Handle(Geom2d_Curve) Curve2d ;
262 Handle(GeomAdaptor_HSurface) hsur = new GeomAdaptor_HSurface(mySurfInit);
263 gp_Pnt2d P2d;
264
265 ProjLib_CompProjectedCurve Projector(hsur, Curv, myTol3d/10, myTol3d/10);
266
267 Standard_Real Udeb, Ufin, ProjUdeb, ProjUfin;
268 Udeb = Curv->FirstParameter();
269 Ufin = Curv->LastParameter();
270 Projector.Bounds( 1, ProjUdeb, ProjUfin );
271
272 if (Projector.NbCurves() != 1 ||
273 Abs( Udeb-ProjUdeb ) > Precision::PConfusion() ||
274 Abs( Ufin-ProjUfin ) > Precision::PConfusion())
275 {
276 if (Projector.IsSinglePnt(1, P2d))
277 {
278 //solution ponctuelles
279 TColgp_Array1OfPnt2d poles(1, 2);
280 poles.Init(P2d);
281 Curve2d = new (Geom2d_BezierCurve) (poles);
282 }
283 else
284 {
285 Curve2d.Nullify(); // Pas de solution continue
9230430c 286#if PLATE_DEB
7fd59977 287 cout << "BuildPlateSurace :: Pas de projection continue" << endl;
288#endif
289 }
290 }
291 else {
292 GeomAbs_Shape Continuity = GeomAbs_C1;
293 Standard_Integer MaxDegree = 10, MaxSeg;
294 Standard_Real Udeb, Ufin;
295 Handle(ProjLib_HCompProjectedCurve) HProjector =
296 new ProjLib_HCompProjectedCurve();
297 HProjector->Set(Projector);
298 Projector.Bounds(1, Udeb, Ufin);
299
300 MaxSeg = 20 + HProjector->NbIntervals(GeomAbs_C3);
301 Approx_CurveOnSurface appr(HProjector, hsur, Udeb, Ufin, myTol3d,
302 Continuity, MaxDegree, MaxSeg,
303 Standard_False, Standard_True);
304
305 Curve2d = appr.Curve2d();
306 }
307#if DRAW
308 if (Affich) {
309 char* name = new char[100];
310 sprintf(name,"proj_%d",++NbProj);
311 DrawTrSurf::Set(name, Curve2d);
312 }
313#endif
314 return Curve2d;
315}
316//---------------------------------------------------------
317// Fonction : ProjectedCurve
318//---------------------------------------------------------
319Handle(Adaptor2d_HCurve2d) GeomPlate_BuildPlateSurface::ProjectedCurve( Handle(Adaptor3d_HCurve)& Curv)
320{ // Projection d'une courbe sur la surface d'init
321
322 Handle(GeomAdaptor_HSurface) hsur = new GeomAdaptor_HSurface(mySurfInit);
323
324 ProjLib_CompProjectedCurve Projector(hsur, Curv, myTolU/10, myTolV/10);
325 Handle(ProjLib_HCompProjectedCurve) HProjector =
326 new ProjLib_HCompProjectedCurve();
327
328 if (Projector.NbCurves() != 1) {
329
330 HProjector.Nullify(); // Pas de solution continue
9230430c 331#if PLATE_DEB
7fd59977 332 cout << "BuildPlateSurace :: Pas de projection continue" << endl;
333#endif
334 }
335 else
336 {
337 Standard_Real First1,Last1,First2,Last2;
338 First1=Curv->FirstParameter();
339 Last1 =Curv->LastParameter();
340 Projector.Bounds(1,First2,Last2);
341
342 if (Abs(First1-First2) <= Max(myTolU,myTolV) &&
343 Abs(Last1-Last2) <= Max(myTolU,myTolV))
344 {
345
346 HProjector->Set(Projector);
347 HProjector = Handle(ProjLib_HCompProjectedCurve)::
348 DownCast(HProjector->Trim(First2,Last2,Precision::PConfusion()));
349 }
350 else
351 {
352 HProjector.Nullify(); // Pas de solution continue
9230430c 353#if PLATE_DEB
7fd59977 354 cout << "BuildPlateSurace :: Pas de projection complete" << endl;
355#endif
356 }
357 }
358 return HProjector;
359}
360
361//---------------------------------------------------------
362// Fonction : ProjectPoint
363//---------------------------------------------------------
364// Projete une point sur la surface d'init
365//---------------------------------------------------------
366gp_Pnt2d GeomPlate_BuildPlateSurface::ProjectPoint(const gp_Pnt &p3d)
367{ Extrema_POnSurf P;
368 myProj.Perform(p3d);
369 Standard_Integer nearest = 1;
370 if( myProj.NbExt() > 1)
371 {
372 Standard_Real dist2mini = myProj.SquareDistance(1);
373 for (Standard_Integer i=2; i<=myProj.NbExt();i++)
374 {
375 if (myProj.SquareDistance(i) < dist2mini)
376 {
377 dist2mini = myProj.SquareDistance(i);
378 nearest = i;
379 }
380 }
381 }
382 P=myProj.Point(nearest);
383 Standard_Real u,v;
384 P.Parameter(u,v);
385 gp_Pnt2d p2d;
386 p2d.SetCoord(u,v);
387
388 return p2d;
389}
390//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
391// =========================================================
392// M E T H O D E S P U B L I Q U E S
393// =========================================================
394//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
395//---------------------------------------------------------
396// Fonction : Init
397//---------------------------------------------------------
398// Initialise les contraintes lineaires et ponctuelles
399//---------------------------------------------------------
400void GeomPlate_BuildPlateSurface::Init()
401{ myLinCont->Clear();
402 myPntCont->Clear();
403 myPntCont = new GeomPlate_HSequenceOfPointConstraint;
404 myLinCont = new GeomPlate_HSequenceOfCurveConstraint;
405}
406
407//---------------------------------------------------------
408// Fonction : LoadInitSurface
409//---------------------------------------------------------
410// Charge la surface initiale
411//---------------------------------------------------------
412void GeomPlate_BuildPlateSurface::LoadInitSurface(const Handle(Geom_Surface)& Surf)
413{ mySurfInit = Surf;
414 mySurfInitIsGive=Standard_True;
415}
416
417//---------------------------------------------------------
418//fonction : Add
419//---------------------------------------------------------
420//---------------------------------------------------------
421// Ajout d'une contrainte lineaire
422//---------------------------------------------------------
423void GeomPlate_BuildPlateSurface::
424 Add(const Handle(GeomPlate_CurveConstraint)& Cont)
425{ myLinCont->Append(Cont);
426}
427
428void GeomPlate_BuildPlateSurface::SetNbBounds( const Standard_Integer NbBounds )
429{
430 myNbBounds = NbBounds;
431}
432
433//---------------------------------------------------------
434//fonction : Add
435//---------------------------------------------------------
436//---------------------------------------------------------
437// Ajout d'une contrainte ponctuelle
438//---------------------------------------------------------
439void GeomPlate_BuildPlateSurface::
440 Add(const Handle(GeomPlate_PointConstraint)& Cont)
441{
442 myPntCont->Append(Cont);
443}
444
445//---------------------------------------------------------
446//fonction : Perform
447// Calcul la surface de remplissage avec les contraintes chargees
448//---------------------------------------------------------
449void GeomPlate_BuildPlateSurface::Perform()
450{
451#ifdef DEB
452 // Chronmetrage
453 OSD_Chronometer Chrono;
454 Chrono.Reset();
455 Chrono.Start();
456#endif
457
458 if (myNbBounds == 0)
459 myNbBounds = myLinCont->Length();
460
461 myPlate.Init();
462 //=====================================================================
463 //Declaration des variables.
464 //=====================================================================
465 Standard_Integer NTLinCont = myLinCont->Length(),
466 NTPntCont = myPntCont->Length(), NbBoucle=0;
467 // La variable NTPoint peut etre enlevee
468 Standard_Boolean Fini=Standard_True;
469 if ((NTLinCont+NTPntCont)==0)
470 Standard_RangeError::Raise("GeomPlate : The number of constraints is null.");
471
472 //======================================================================
473 // Surface Initiale
474 //======================================================================
475 if (!mySurfInitIsGive)
476 ComputeSurfInit();
477
478 else {
479 if (NTLinCont>=2)
480 { // Tableau des permutations pour conserver l'ordre initial voir TrierTab
481 myInitOrder = new TColStd_HArray1OfInteger(1,NTLinCont);
482 for (Standard_Integer l=1;l<=NTLinCont;l++)
483 myInitOrder->SetValue(l,l);
484 if (!CourbeJointive(myTol3d))
485 {// Standard_Failure::Raise("Curves are not joined");
9230430c 486#ifdef PLATE_DEB
7fd59977 487 cout<<"WARNING : Courbes non jointives a "<<myTol3d<<" pres"<<endl;
488#endif
489 }
490 TrierTab(myInitOrder); // Reordonne le tableau des permutations
491 }
492 else if(NTLinCont > 0)//Patch
493 {
494 mySense = new TColStd_HArray1OfInteger( 1, NTLinCont, 0 );
495 myInitOrder = new TColStd_HArray1OfInteger( 1, NTLinCont, 1 );
496 }
497 }
498
499 Standard_Real u1,v1,u2,v2;
500 mySurfInit->Bounds(u1,v1,u2,v2);
501 GeomAdaptor_Surface Surf(mySurfInit);
502 myTolU = Surf.UResolution(myTol3d);
503 myTolV = Surf.VResolution(myTol3d);
504 myProj.Initialize(Surf,u1,v1,u2,v2,
505 myTolU,myTolV);
506
507 //======================================================================
508 // Projection des courbes
509 //======================================================================
510 Standard_Integer i;
511 Standard_Boolean Ok = Standard_True;
512 for (i = 1; i <= NTLinCont; i++)
513 if(myLinCont->Value(i)->Curve2dOnSurf().IsNull())
514 {
515 Handle( Geom2d_Curve ) Curve2d = ProjectCurve( myLinCont->Value(i)->Curve3d() );
516 if (Curve2d.IsNull())
517 {
518 Ok = Standard_False;
519 break;
520 }
521 myLinCont->ChangeValue(i)->SetCurve2dOnSurf( Curve2d );
522 }
523 if (!Ok)
524 {
525 GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d,
526 1, 3,
527 15*myTol3d,
528 -1, GeomAbs_C0,
529 1.3);
530 mySurfInit = App.Surface();
531
532 mySurfInit->Bounds(u1,v1,u2,v2);
533 GeomAdaptor_Surface Surf(mySurfInit);
534 myTolU = Surf.UResolution(myTol3d);
535 myTolV = Surf.VResolution(myTol3d);
536 myProj.Initialize(Surf,u1,v1,u2,v2,
537 myTolU,myTolV);
538
539 Ok = Standard_True;
540 for (i = 1; i <= NTLinCont; i++)
541 {
542 Handle( Geom2d_Curve ) Curve2d = ProjectCurve( myLinCont->Value(i)->Curve3d() );
543 if (Curve2d.IsNull())
544 {
545 Ok = Standard_False;
546 break;
547 }
548 myLinCont->ChangeValue(i)->SetCurve2dOnSurf( Curve2d );
549 }
550 if (!Ok)
551 {
552 mySurfInit = myPlanarSurfInit;
553
554 mySurfInit->Bounds(u1,v1,u2,v2);
555 GeomAdaptor_Surface Surf(mySurfInit);
556 myTolU = Surf.UResolution(myTol3d);
557 myTolV = Surf.VResolution(myTol3d);
558 myProj.Initialize(Surf,u1,v1,u2,v2,
559 myTolU,myTolV);
560
561 for (i = 1; i <= NTLinCont; i++)
562 myLinCont->ChangeValue(i)->
563 SetCurve2dOnSurf(ProjectCurve( myLinCont->Value(i)->Curve3d() ) );
564 }
565 else { // Project the points
566 for ( i=1;i<=NTPntCont;i++) {
567 gp_Pnt P;
568 myPntCont->Value(i)->D0(P);
569 myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
570 }
571 }
572 }
573
574 //======================================================================
575 // Projection des points
576 //======================================================================
577 for ( i=1;i<=NTPntCont;i++) {
578 if (! myPntCont->Value(i)->HasPnt2dOnSurf()) {
579 gp_Pnt P;
580 myPntCont->Value(i)->D0(P);
581 myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
582 }
583 }
584
585 //======================================================================
586 // Nbre de point par courbe
587 //======================================================================
588 if ((NTLinCont !=0)&&(myNbPtsOnCur!=0))
589 CalculNbPtsInit();
590
591 //======================================================================
592 // Gestion des incompatibilites entre courbes
593 //======================================================================
594 Handle(GeomPlate_HArray1OfSequenceOfReal) PntInter;
595 Handle(GeomPlate_HArray1OfSequenceOfReal) PntG1G1;
596 if (NTLinCont != 0) {
597 PntInter = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
598 PntG1G1 = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
599 Intersect(PntInter, PntG1G1);
600 }
601
602 //======================================================================
603 // Boucle pour obtenir une meilleur surface
604 //======================================================================
605
606 myFree = !myIsLinear;
607
608 do
609 {
9230430c 610#if PLATE_DEB
7fd59977 611 if (Affich && NbBoucle) {
612 cout<<"Resultats boucle"<< NbBoucle << endl;
613 cout<<"DistMax="<<myG0Error<<endl;
614 if (myG1Error!=0)
615 cout<<"AngleMax="<<myG1Error<<endl;
616 if (myG2Error!=0)
617 cout<<"CourbMax="<<myG2Error<<endl;
618 }
619#endif
620 NbBoucle++;
621 if (NTLinCont!=0)
622 { //====================================================================
623 // Calcul le nombre de point total et le maximum de points par courbe
624 //====================================================================
625 Standard_Integer NPointMax=0;
626 for (Standard_Integer i=1;i<=NTLinCont;i++)
627 if ((myLinCont->Value(i)->NbPoints())>NPointMax)
628 NPointMax=(Standard_Integer )( myLinCont->Value(i)->NbPoints());
629 //====================================================================
630 // Discretisation des courbes
631 //====================================================================
632 Discretise(PntInter, PntG1G1);
633 //====================================================================
634 //Preparation des points de contrainte pour plate
635 //====================================================================
636 LoadCurve( NbBoucle );
637 if ( myPntCont->Length() != 0)
638 LoadPoint( NbBoucle );
639 //====================================================================
640 //Resolution de la surface
641 //====================================================================
642 myPlate.SolveTI(myDegree, ComputeAnisotropie());
643 if (!myPlate.IsDone())
644 Standard_Failure::Raise("GeomPlate : abort calcul of Plate.");
645 myGeomPlateSurface = new GeomPlate_Surface(mySurfInit,myPlate);
646 Standard_Real Umin,Umax,Vmin,Vmax;
647 myPlate.UVBox(Umin,Umax,Vmin,Vmax);
648 myGeomPlateSurface->SetBounds(Umin,Umax,Vmin,Vmax);
649
650 Fini = VerifSurface(NbBoucle);
651 if ((NbBoucle >= myNbIter)&&(!Fini))
652 {
9230430c 653#ifdef PLATE_DEB
7fd59977 654 cout<<"Warning objectif non atteint"<<endl;
655#endif
656 Fini = Standard_True;
657 }
658
659 if ((NTPntCont != 0)&&(Fini))
660 { Standard_Real di,an,cu;
661 VerifPoints(di,an,cu);
662 }
663 }
664 else
665 { LoadPoint( NbBoucle );
666 //====================================================================
667 //Resolution de la surface
668 //====================================================================
669 myPlate.SolveTI(myDegree, ComputeAnisotropie());
670 if (!myPlate.IsDone())
671 Standard_Failure::Raise("GeomPlate : abort calcul of Plate.");
672 myGeomPlateSurface = new GeomPlate_Surface(mySurfInit,myPlate);
673 Standard_Real Umin,Umax,Vmin,Vmax;
674 myPlate.UVBox(Umin,Umax,Vmin,Vmax);
675 myGeomPlateSurface->SetBounds(Umin,Umax,Vmin,Vmax);
676 Fini = Standard_True;
677 Standard_Real di,an,cu;
678 VerifPoints(di,an,cu);
679 }
680 } while (!Fini); // Fin boucle pour meilleur surface
9230430c 681#ifdef PLATE_DEB
7fd59977 682 if (NTLinCont != 0)
683 { cout<<"======== Resultats globaux ==========="<<endl;
684 cout<<"DistMax="<<myG0Error<<endl;
685 if (myG1Error!=0)
686 cout<<"AngleMax="<<myG1Error<<endl;
687 if (myG2Error!=0)
688 cout<<"CourbMax="<<myG2Error<<endl;
689 }
690 Chrono.Stop();
691 Standard_Real Tps;
692 Chrono.Show(Tps);
693 cout<<"*** FIN DE GEOMPLATE ***"<<endl;
694 cout<<"Temps de calcul : "<<Tps<<endl;
695 cout<<"Nombre de boucle : "<<NbBoucle<<endl;
696#endif
697}
698
699//---------------------------------------------------------
700// fonction : EcartContraintesMIL
701//---------------------------------------------------------
702void GeomPlate_BuildPlateSurface::
703EcartContraintesMil ( const Standard_Integer c,
704 Handle(TColStd_HArray1OfReal)& d,
705 Handle(TColStd_HArray1OfReal)& an,
706 Handle(TColStd_HArray1OfReal)& courb)
707{
708 Standard_Integer NbPt=myParCont->Value(c).Length();
709 Standard_Real U;
710 if (NbPt<3)
711 NbPt=4;
712 else
713 NbPt=myParCont->Value(c).Length();
714 gp_Vec v1i, v1f, v2i, v2f, v3i, v3f;
715 gp_Pnt Pi, Pf;
716 gp_Pnt2d P2d;Standard_Integer i;
717 Handle(GeomPlate_CurveConstraint) LinCont = myLinCont->Value(c);
718 switch (LinCont->Order())
719 { case 0 :
720 for (i=1; i<NbPt; i++)
721 { U = ( myParCont->Value(c).Value(i) +
722 myParCont->Value(c).Value(i+1) )/2;
723 LinCont->D0(U, Pi);
724 if (!LinCont->ProjectedCurve().IsNull())
725 P2d = LinCont->ProjectedCurve()->Value(U);
726 else
727 { if (!LinCont->Curve2dOnSurf().IsNull())
728 P2d = LinCont->Curve2dOnSurf()->Value(U);
729 else
730 P2d = ProjectPoint(Pi);
731 }
732 myGeomPlateSurface->D0( P2d.Coord(1), P2d.Coord(2), Pf);
733 an->Init(0);
734 courb->Init(0);
735 d->ChangeValue(i) = Pf.Distance(Pi);
736 }
737 break;
738 case 1 :
739 for (i=1; i<NbPt; i++)
740 { U = ( myParCont->Value(c).Value(i) +
741 myParCont->Value(c).Value(i+1) )/2;
742 LinCont->D1(U, Pi, v1i, v2i);
743 if (!LinCont->ProjectedCurve().IsNull())
744 P2d = LinCont->ProjectedCurve()->Value(U);
745 else
746 { if (!LinCont->Curve2dOnSurf().IsNull())
747 P2d = LinCont->Curve2dOnSurf()->Value(U);
748 else
749 P2d = ProjectPoint(Pi);
750 }
751 myGeomPlateSurface->D1( P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f);
752 d->ChangeValue(i) = Pf.Distance(Pi);
753 v3i = v1i^v2i; v3f=v1f^v2f;
754 Standard_Real angle=v3f.Angle(v3i);
c6541a0c
D
755 if (angle>(M_PI/2))
756 an->ChangeValue(i) = M_PI -angle;
7fd59977 757 else
758 an->ChangeValue(i) = angle;
759 courb->Init(0);
760 }
761 break;
762 case 2 :
763 Handle(Geom_Surface) Splate;
764 Splate = Handle(Geom_Surface)::DownCast(myGeomPlateSurface);
765 LocalAnalysis_SurfaceContinuity CG2;
766 for (i=1; i<NbPt; i++)
767 { U = ( myParCont->Value(c).Value(i) +
768 myParCont->Value(c).Value(i+1) )/2;
769 LinCont->D0(U, Pi);
770 if (!LinCont->ProjectedCurve().IsNull())
771 P2d = LinCont->ProjectedCurve()->Value(U);
772 else
773 { if (!LinCont->Curve2dOnSurf().IsNull())
774 P2d = LinCont->Curve2dOnSurf()->Value(U);
775 else
776 P2d = ProjectPoint(Pi);
777 }
778 GeomLProp_SLProps Prop (Splate, P2d.Coord(1), P2d.Coord(2),
779 2, 0.001);
780 CG2.ComputeAnalysis(Prop, myLinCont->Value(c)->LPropSurf(U),
781 GeomAbs_G2);
782 d->ChangeValue(i)=CG2.C0Value();
783 an->ChangeValue(i)=CG2.G1Angle();
784 courb->ChangeValue(i)=CG2.G2CurvatureGap();
785 }
786 break;
787 }
788}
789
790
791
792//---------------------------------------------------------
793// fonction : Disc2dContour
794//---------------------------------------------------------
795void GeomPlate_BuildPlateSurface::
796 Disc2dContour ( const Standard_Integer nbp,
797 TColgp_SequenceOfXY& Seq2d)
798{
9230430c
M
799#ifdef PLATE_DEB
800 if (nbp!=4)
801 cout<<"nbp doit etre egal a 4 pour Disc2dContour"<<endl;
802#endif
7fd59977 803 // initialisation
804 Seq2d.Clear();
805
806 // echantillonnage en "cosinus" + 3 points sur chaque intervalle
807 Standard_Integer NTCurve = myLinCont->Length();
808 Standard_Integer NTPntCont = myPntCont->Length();
809 gp_Pnt2d P2d;
810 gp_XY UV;
811 gp_Pnt PP;
812 Standard_Real u1,v1,u2,v2;
813 Standard_Integer i ;
814
815 mySurfInit->Bounds(u1,v1,u2,v2);
816 GeomAdaptor_Surface Surf(mySurfInit);
817 myProj.Initialize(Surf,u1,v1,u2,v2,
818 myTolU,myTolV);
819
820 for( i=1; i<=NTPntCont; i++)
821 if (myPntCont->Value(i)->Order()!=-1)
822 { P2d = myPntCont->Value(i)->Pnt2dOnSurf();
823 UV.SetX(P2d.Coord(1));
824 UV.SetY(P2d.Coord(2));
825 Seq2d.Append(UV);
826 }
827 for(i=1; i<=NTCurve; i++)
828 {
829 Handle(GeomPlate_CurveConstraint) LinCont = myLinCont->Value(i);
830 if (LinCont->Order()!=-1)
831 { Standard_Integer NbPt=myParCont->Value(i).Length();
832 // premier point de contrainte (j=0)
833 if (!LinCont->ProjectedCurve().IsNull())
834 P2d = LinCont->ProjectedCurve()->Value(myParCont->Value(i).Value(1));
835
836 else
837 { if (!LinCont->Curve2dOnSurf().IsNull())
838 P2d = LinCont->Curve2dOnSurf()->Value(myParCont->Value(i).Value(1));
839
840 else
841 {
842 LinCont->D0(myParCont->Value(i).Value(1),PP);
843 P2d = ProjectPoint(PP);
844 }
845 }
846
847 UV.SetX(P2d.Coord(1));
848 UV.SetY(P2d.Coord(2));
849 Seq2d.Append(UV);
850 for (Standard_Integer j=2; j<NbPt; j++)
851 { Standard_Real Uj=myParCont->Value(i).Value(j),
852 Ujp1=myParCont->Value(i).Value(j+1);
853 if (!LinCont->ProjectedCurve().IsNull())
854 P2d = LinCont->ProjectedCurve()->Value((Ujp1+3*Uj)/4);
855
856 else
857 { if (!LinCont->Curve2dOnSurf().IsNull())
858 P2d = LinCont->Curve2dOnSurf()->Value((Ujp1+3*Uj)/4);
859
860 else
861 {
862 LinCont->D0((Ujp1+3*Uj)/4,PP);
863 P2d = ProjectPoint(PP);
864 }
865 }
866 UV.SetX(P2d.Coord(1));
867 UV.SetY(P2d.Coord(2));
868 Seq2d.Append(UV);
869 // point milieu precedent
870 if (!LinCont->ProjectedCurve().IsNull())
871 P2d = LinCont->ProjectedCurve()->Value((Ujp1+Uj)/2);
872
873 else
874 { if (!LinCont->Curve2dOnSurf().IsNull())
875 P2d = LinCont->Curve2dOnSurf()->Value((Ujp1+Uj)/2);
876
877 else
878 {
879 LinCont->D0((Ujp1+Uj)/2,PP);
880 P2d = ProjectPoint(PP);
881 }
882 }
883
884 UV.SetX(P2d.Coord(1));
885 UV.SetY(P2d.Coord(2));
886 Seq2d.Append(UV);
887 // point 3/4 precedent
888 if (!LinCont->ProjectedCurve().IsNull())
889 P2d = LinCont->ProjectedCurve()->Value((3*Ujp1+Uj)/4);
890
891 else
892 { if (!LinCont->Curve2dOnSurf().IsNull())
893 P2d = LinCont->Curve2dOnSurf()->Value((3*Ujp1+Uj)/4);
894
895 else
896 {
897 LinCont->D0((3*Ujp1+Uj)/4,PP);
898 P2d = ProjectPoint(PP);
899 }
900 }
901
902 UV.SetX(P2d.Coord(1));
903 UV.SetY(P2d.Coord(2));
904 Seq2d.Append(UV);
905 // point de contrainte courant
906 if (!LinCont->ProjectedCurve().IsNull())
907 P2d = LinCont->ProjectedCurve()->Value(Ujp1);
908
909 else
910 { if (!LinCont->Curve2dOnSurf().IsNull())
911 P2d = LinCont->Curve2dOnSurf()->Value(Ujp1);
912
913 else
914 {
915 LinCont->D0(Ujp1,PP);
916 P2d = ProjectPoint(PP);
917 }
918 }
919
920 UV.SetX(P2d.Coord(1));
921 UV.SetY(P2d.Coord(2));
922 Seq2d.Append(UV);
923 }
924 }
925 }
926}
927
928//---------------------------------------------------------
929// fonction : Disc3dContour
930//---------------------------------------------------------
931void GeomPlate_BuildPlateSurface::
932Disc3dContour ( const Standard_Integer nbp,
933 const Standard_Integer iordre,
934 TColgp_SequenceOfXYZ& Seq3d)
935{
9230430c
M
936#ifdef PLATE_DEB
937 if (nbp!=4)
938 cout<<"nbp doit etre egal a 4 pour Disc3dContour"<<endl;
939 if (iordre!=0&&iordre!=1)
940 cout<<"iordre incorrect pour Disc3dContour"<<endl;
941#endif
7fd59977 942 // initialisation
943 Seq3d.Clear();
944 // echantillonnage en "cosinus" + 3 points sur chaque intervalle
945 Standard_Real u1,v1,u2,v2;
946 mySurfInit->Bounds(u1,v1,u2,v2);
947 GeomAdaptor_Surface Surf(mySurfInit);
948 myProj.Initialize(Surf,u1,v1,u2,v2,
949 Surf.UResolution(myTol3d),
950 Surf.VResolution(myTol3d));
951 Standard_Integer NTCurve = myLinCont->Length();
952 Standard_Integer NTPntCont = myPntCont->Length();
953// gp_Pnt2d P2d;
954 gp_Pnt P3d;
955 gp_Vec v1h,v2h,v3h;
956 gp_XYZ Pos;
957 Standard_Integer i ;
958
959 for( i=1; i<=NTPntCont; i++)
960 if (myPntCont->Value(i)->Order()!=-1)
961 { if (iordre==0)
962 { myPntCont->Value(i)->D0(P3d);
963 Pos.SetX(P3d.X());
964 Pos.SetY(P3d.Y());
965 Pos.SetZ(P3d.Z());
966 Seq3d.Append(Pos);
967 }
968 else {
969 myPntCont->Value(i)->D1(P3d,v1h,v2h);
970 v3h=v1h^v2h;
971 Pos.SetX(v3h.X());
972 Pos.SetY(v3h.Y());
973 Pos.SetZ(v3h.Z());
974 Seq3d.Append(Pos);
975 }
976 }
977
978 for(i=1; i<=NTCurve; i++)
979 if (myLinCont->Value(i)->Order()!=-1)
980
981 { Standard_Integer NbPt=myParCont->Value(i).Length();;
982 // premier point de contrainte (j=0)
983 // Standard_Integer NbPt=myParCont->Length();
984 if (iordre==0) {
985 myLinCont->Value(i)->D0(myParCont->Value(i).Value(1),P3d);
986 Pos.SetX(P3d.X());
987 Pos.SetY(P3d.Y());
988 Pos.SetZ(P3d.Z());
989 Seq3d.Append(Pos);
990 }
991 else {
992 myLinCont->Value(i)->D1(myParCont->Value(i).Value(1),P3d,v1h,v2h);
993 v3h=v1h^v2h;
994 Pos.SetX(v3h.X());
995 Pos.SetY(v3h.Y());
996 Pos.SetZ(v3h.Z());
997 Seq3d.Append(Pos);
998 }
999
1000 for (Standard_Integer j=2; j<NbPt; j++)
1001 { Standard_Real Uj=myParCont->Value(i).Value(j),
1002 Ujp1=myParCont->Value(i).Value(j+1);
1003 if (iordre==0) {
1004 // point 1/4 precedent
1005 myLinCont->Value(i)->D0((Ujp1+3*Uj)/4,P3d);
1006 Pos.SetX(P3d.X());
1007 Pos.SetY(P3d.Y());
1008 Pos.SetZ(P3d.Z());
1009 Seq3d.Append(Pos);
1010 // point milieu precedent
1011 myLinCont->Value(i)->D0((Ujp1+Uj)/2,P3d);
1012 Pos.SetX(P3d.X());
1013 Pos.SetY(P3d.Y());
1014 Pos.SetZ(P3d.Z());
1015 Seq3d.Append(Pos);
1016 // point 3/4 precedent
1017 myLinCont->Value(i)->D0((3*Ujp1+Uj)/4,P3d);
1018 Pos.SetX(P3d.X());
1019 Pos.SetY(P3d.Y());
1020 Pos.SetZ(P3d.Z());
1021 Seq3d.Append(Pos);
1022 // point de contrainte courant
1023 myLinCont->Value(i)->D0(Ujp1,P3d);
1024 Pos.SetX(P3d.X());
1025 Pos.SetY(P3d.Y());
1026 Pos.SetZ(P3d.Z());
1027 Seq3d.Append(Pos);
1028 }
1029 else {
1030 // point 1/4 precedent
1031 myLinCont->Value(i)->D1((Ujp1+3*Uj)/4,P3d,v1h,v2h);
1032 v3h=v1h^v2h;
1033 Pos.SetX(v3h.X());
1034 Pos.SetY(v3h.Y());
1035 Pos.SetZ(v3h.Z());
1036 Seq3d.Append(Pos);
1037 // point milieu precedent
1038 myLinCont->Value(i)->D1((Ujp1+Uj)/2,P3d,v1h,v2h);
1039 v3h=v1h^v2h;
1040 Pos.SetX(v3h.X());
1041 Pos.SetY(v3h.Y());
1042 Pos.SetZ(v3h.Z());
1043 Seq3d.Append(Pos);
1044 // point 3/4 precedent
1045 myLinCont->Value(i)->D1((3*Ujp1+Uj)/4,P3d,v1h,v2h);
1046 v3h=v1h^v2h;
1047 Pos.SetX(v3h.X());
1048 Pos.SetY(v3h.Y());
1049 Pos.SetZ(v3h.Z());
1050 Seq3d.Append(Pos);
1051 // point de contrainte courant
1052 myLinCont->Value(i)->D1(Ujp1,P3d,v1h,v2h);
1053 v3h=v1h^v2h;
1054 Pos.SetX(v3h.X());
1055 Pos.SetY(v3h.Y());
1056 Pos.SetZ(v3h.Z());
1057 Seq3d.Append(Pos);
1058 }
1059 }
1060 }
1061
1062}
1063
1064
1065//---------------------------------------------------------
1066// fonction : IsDone
1067//---------------------------------------------------------
1068Standard_Boolean GeomPlate_BuildPlateSurface::IsDone() const
1069{ return myPlate.IsDone();
1070}
1071
1072
1073
1074//---------------------------------------------------------
1075// fonction : Surface
1076//---------------------------------------------------------
1077Handle_GeomPlate_Surface GeomPlate_BuildPlateSurface::Surface() const
1078{ return myGeomPlateSurface ;
1079}
1080
1081//---------------------------------------------------------
1082// fonction : SurfInit
1083//---------------------------------------------------------
1084Handle_Geom_Surface GeomPlate_BuildPlateSurface::SurfInit() const
1085{ return mySurfInit ;
1086}
1087
1088
1089//---------------------------------------------------------
1090// fonction : Sense
1091//---------------------------------------------------------
1092Handle_TColStd_HArray1OfInteger GeomPlate_BuildPlateSurface::Sense() const
1093{ Standard_Integer NTCurve = myLinCont->Length();
1094 Handle(TColStd_HArray1OfInteger) Sens = new TColStd_HArray1OfInteger(1,
1095 NTCurve);
1096 for (Standard_Integer i=1; i<=NTCurve; i++)
1097 Sens->SetValue(i,mySense->Value(myInitOrder->Value(i)));
1098 return Sens;
1099}
1100
1101
1102
1103//---------------------------------------------------------
1104// fonction : Curve2d
1105//---------------------------------------------------------
1106Handle_TColGeom2d_HArray1OfCurve GeomPlate_BuildPlateSurface::Curves2d() const
1107{ Standard_Integer NTCurve = myLinCont->Length();
1108 Handle(TColGeom2d_HArray1OfCurve) C2dfin =
1109 new TColGeom2d_HArray1OfCurve(1,NTCurve);
1110
1111 for (Standard_Integer i=1; i<=NTCurve; i++)
1112 C2dfin->SetValue(i,myLinCont->Value(myInitOrder->Value(i))->Curve2dOnSurf());
1113 return C2dfin;
1114
1115}
1116
1117
1118//---------------------------------------------------------
1119//fonction : Order
1120//---------------------------------------------------------
1121Handle_TColStd_HArray1OfInteger GeomPlate_BuildPlateSurface::Order() const
1122{ Handle_TColStd_HArray1OfInteger result=
1123 new TColStd_HArray1OfInteger(1,myLinCont->Length());
1124 for (Standard_Integer i=1;i<=myLinCont->Length();i++)
1125 result->SetValue(myInitOrder->Value(i),i);
1126 return result;
1127}
1128
1129
1130//---------------------------------------------------------
1131// fonction : G0Error
1132//---------------------------------------------------------
1133Standard_Real GeomPlate_BuildPlateSurface::G0Error() const
1134{ return myG0Error;
1135}
1136
1137//---------------------------------------------------------
1138// fonction : G1Error
1139//---------------------------------------------------------
1140Standard_Real GeomPlate_BuildPlateSurface::G1Error() const
1141{ return myG1Error;
1142}
1143
1144//---------------------------------------------------------
1145// fonction : G2Error
1146//---------------------------------------------------------
1147Standard_Real GeomPlate_BuildPlateSurface::G2Error() const
1148{ return myG2Error;
1149}
1150
1151//=======================================================================
1152//function : G0Error
1153//purpose :
1154//=======================================================================
1155
1156Standard_Real GeomPlate_BuildPlateSurface::G0Error( const Standard_Integer Index )
1157{
1158 Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1159 Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1160 Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1161 EcartContraintesMil( Index, tdistance, tangle, tcurvature );
1162 Standard_Real MaxDistance = 0.;
1163 for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++)
1164 if (tdistance->Value(i) > MaxDistance)
1165 MaxDistance = tdistance->Value(i);
1166 return MaxDistance;
1167}
1168
1169//=======================================================================
1170//function : G1Error
1171//purpose :
1172//=======================================================================
1173
1174Standard_Real GeomPlate_BuildPlateSurface::G1Error( const Standard_Integer Index )
1175{
1176 Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1177 Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1178 Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1179 EcartContraintesMil( Index, tdistance, tangle, tcurvature );
1180 Standard_Real MaxAngle = 0.;
1181 for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++)
1182 if (tangle->Value(i) > MaxAngle)
1183 MaxAngle = tangle->Value(i);
1184 return MaxAngle;
1185}
1186
1187//=======================================================================
1188//function : G2Error
1189//purpose :
1190//=======================================================================
1191
1192Standard_Real GeomPlate_BuildPlateSurface::G2Error( const Standard_Integer Index )
1193{
1194 Handle( TColStd_HArray1OfReal ) tdistance = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1195 Handle( TColStd_HArray1OfReal ) tangle = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1196 Handle( TColStd_HArray1OfReal ) tcurvature = new TColStd_HArray1OfReal( 1, myNbPtsOnCur );
1197 EcartContraintesMil( Index, tdistance, tangle, tcurvature );
1198 Standard_Real MaxCurvature = 0.;
1199 for (Standard_Integer i = 1; i <= myNbPtsOnCur; i++)
1200 if (tcurvature->Value(i) > MaxCurvature)
1201 MaxCurvature = tcurvature->Value(i);
1202 return MaxCurvature;
1203}
1204
1205Handle(GeomPlate_CurveConstraint) GeomPlate_BuildPlateSurface::CurveConstraint( const Standard_Integer order) const
1206{
1207 return myLinCont->Value(order);
1208}
1209Handle(GeomPlate_PointConstraint) GeomPlate_BuildPlateSurface::PointConstraint( const Standard_Integer order) const
1210{
1211 return myPntCont->Value(order);
1212}
1213//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
1214// =========================================================
1215// M E T H O D E S P R I V E S
1216// =========================================================
1217//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//\\//
1218
1219//=======================================================================
1220// function : CourbeJointive
1221// purpose : Effectue un chainage des courbes pour pouvoir calculer
1222// la surface d'init avec la methode du flux maxi.
1223// Retourne vrai si il s'agit d'un contour ferme.
1224//=======================================================================
1225
1226Standard_Boolean GeomPlate_BuildPlateSurface::
1227 CourbeJointive(const Standard_Real tolerance)
1228{ Standard_Integer nbf = myLinCont->Length();
1229 Standard_Real Ufinal1,Uinit1,Ufinal2,Uinit2;
1230 mySense = new TColStd_HArray1OfInteger(1,nbf,0);
1231 Standard_Boolean result=Standard_True;
1232 Standard_Integer j=1,i;
1233 gp_Pnt P1,P2;
1234
1235 while (j <= (myNbBounds-1))
1236 {
1237 Standard_Integer a=0;
1238 i=j+1;
1239 if (i > myNbBounds)
1240 { result = Standard_False;
1241 a=2;
1242 }
1243 while (a<1)
1244 { if (i > myNbBounds)
1245 { result = Standard_False;
1246 a=2;
1247 }
1248 else
1249 {
1250 Uinit1=myLinCont->Value(j)->FirstParameter();
1251 Ufinal1=myLinCont->Value(j)->LastParameter();
1252 Uinit2=myLinCont->Value(i)->FirstParameter();
1253 Ufinal2=myLinCont->Value(i)->LastParameter();
1254 if (mySense->Value(j)==1)
1255 Ufinal1=Uinit1;
1256 myLinCont->Value(j)->D0(Ufinal1,P1);
1257 myLinCont->Value(i)->D0(Uinit2,P2);
1258 if (P1.Distance(P2)<tolerance)
1259 { if (i!=j+1) {
1260 Handle(GeomPlate_CurveConstraint) tampon = myLinCont->Value(j+1);
1261 myLinCont->SetValue(j+1,myLinCont->Value(i));
1262 myLinCont->SetValue(i,tampon);
1263 Standard_Integer Tmp=myInitOrder->Value(j+1);
1264 //Voir fonction TrierTab pour le fonctionnement de myInitOrder
1265 myInitOrder->SetValue(j+1,myInitOrder->Value(i));
1266 myInitOrder->SetValue(i,Tmp);
1267
1268
1269 };
1270 a=2;
1271 mySense->SetValue(j+1,0);
1272
1273 }
1274 else
1275 { myLinCont->Value(i)->D0(Ufinal2,P2);
1276
1277 if (P1.Distance(P2)<tolerance)
1278 { if (i!=j+1) {
1279 Handle(GeomPlate_CurveConstraint) tampon =
1280 myLinCont->Value(j+1);
1281 myLinCont->SetValue(j+1,myLinCont->Value(i));
1282 myLinCont->SetValue(i,tampon);
1283 Standard_Integer Tmp=myInitOrder->Value(j+1);
1284 //Voir fonction TrierTab pour le fonctionnement de myInitOrder
1285 myInitOrder->SetValue(j+1,myInitOrder->Value(i));
1286 myInitOrder->SetValue(i,Tmp);
1287 };
1288 a=2;
1289 mySense->SetValue(j+1,1);
1290 }
1291 }
1292 }
1293 i++;
1294 }
1295 j++;
1296 }
1297 Uinit1=myLinCont->Value( myNbBounds )->FirstParameter();
1298 Ufinal1=myLinCont->Value( myNbBounds )->LastParameter();
1299 Uinit2=myLinCont->Value(1)->FirstParameter();
1300 Ufinal2=myLinCont->Value(1)->LastParameter();
1301 myLinCont->Value( myNbBounds )->D0(Ufinal1,P1);
1302 myLinCont->Value(1)->D0(Uinit2,P2);
1303 if ((mySense->Value( myNbBounds )==0)
1304 &&(P1.Distance(P2)<tolerance))
1305 {
1306 return ((Standard_True)&&(result));
1307 }
1308 myLinCont->Value( myNbBounds )->D0(Uinit1,P1);
1309 if ((mySense->Value( myNbBounds )==1)
1310 &&(P1.Distance(P2)<tolerance))
1311 {
1312 return ((Standard_True)&&(result));
1313 }
1314 else return Standard_False;
1315}
1316
1317
1318//-------------------------------------------------------------------------
1319// fonction : ComputeSurfInit
1320//-------------------------------------------------------------------------
1321// Calcul la surface d'init soit par la methode du flux maxi soit par
1322// la methode du plan d'inertie si le contour n'est pas ferme ou si
1323// il y a des contraintes ponctuelles
1324//-------------------------------------------------------------------------
1325
1326void GeomPlate_BuildPlateSurface::ComputeSurfInit()
1327{
1328 Standard_Integer nopt=2, popt=2, Np=1;
1329 Standard_Boolean isHalfSpace = Standard_True;
1330 Standard_Real LinTol = 0.001, AngTol = 0.001; //AngTol = 0.0001; //LinTol = 0.0001
1331
1332 // Option pour le calcul du plan initial
1333 Standard_Integer NTLinCont = myLinCont->Length(), NTPntCont = myPntCont->Length();
1334 // Tableau des permutation pour conserver l'ordre initial voir TrierTab
1335 if (NTLinCont != 0) {
1336 myInitOrder = new TColStd_HArray1OfInteger(1,NTLinCont);
1337 for (Standard_Integer i = 1; i <= NTLinCont; i++)
1338 myInitOrder->SetValue( i, i );
1339 }
1340
1341 Standard_Boolean CourbeJoint=Standard_False;
1342
1343 if (NTLinCont != 0 && (CourbeJoint = CourbeJointive( myTol3d )) && IsOrderG1())
1344 {
1345 nopt = 3;
1346 // Tableau contenant le nuage de point pour le calcul du plan
1347 Standard_Integer NbPoint = 20, Discr = NbPoint/4, pnum = 0;
1348 Handle( TColgp_HArray1OfPnt ) Pts = new TColgp_HArray1OfPnt( 1, (NbPoint+1)*NTLinCont+NTPntCont );
1349 TColgp_SequenceOfVec Vecs, NewVecs;
1350 GeomPlate_SequenceOfAij Aset;
1351 Standard_Real Uinit, Ufinal, Uif;
1352 gp_Vec LastVec;
1353 Standard_Integer i ;
1354 for ( i = 1; i <= NTLinCont; i++)
1355 {
1356 Standard_Integer Order = myLinCont->Value(i)->Order();
1357
1358 NewVecs.Clear();
1359
1360 Uinit = myLinCont->Value(i)->FirstParameter();
1361 Ufinal = myLinCont->Value(i)->LastParameter();
1362 Uif = Ufinal - Uinit;
1363 if (mySense->Value(i) == 1)
1364 {
1365 Uinit = Ufinal;
1366 Uif = -Uif;
1367 }
1368
1369 gp_Vec Vec1, Vec2, Normal;
1370 Standard_Boolean ToReverse = Standard_False;
1371 if (i > 1 && Order >= GeomAbs_G1)
1372 {
1373 gp_Pnt P;
1374 myLinCont->Value(i)->D1( Uinit, P, Vec1, Vec2 );
1375 Normal = Vec1 ^ Vec2;
1376 if (LastVec.IsOpposite( Normal, AngTol ))
1377 ToReverse = Standard_True;
1378 }
1379
1380 for (Standard_Integer j = 0; j <= NbPoint; j++)
1381 { // Nombre de point par courbe = 20
1382 // Repartition lineaire
1383 Standard_Real Inter = j*Uif/(NbPoint);
1384 if (Order < GeomAbs_G1 || j % Discr != 0)
1385 myLinCont->Value(i)->D0( Uinit+Inter, Pts->ChangeValue(++pnum) );
1386 else
1387 {
1388 myLinCont->Value(i)->D1( Uinit+Inter, Pts->ChangeValue(++pnum), Vec1, Vec2 );
1389 Normal = Vec1 ^ Vec2;
1390 Normal.Normalize();
1391 if (ToReverse)
1392 Normal.Reverse();
1393 Standard_Boolean isNew = Standard_True;
1394 Standard_Integer k ;
1395 for ( k = 1; k <= Vecs.Length(); k++)
1396 if (Vecs(k).IsEqual( Normal, LinTol, AngTol ))
1397 {
1398 isNew = Standard_False;
1399 break;
1400 }
1401 if (isNew)
1402 for (k = 1; k <= NewVecs.Length(); k++)
1403 if (NewVecs(k).IsEqual( Normal, LinTol, AngTol ))
1404 {
1405 isNew = Standard_False;
1406 break;
1407 }
1408 if (isNew)
1409 NewVecs.Append( Normal );
1410 }
1411 }
1412 if (Order >= GeomAbs_G1)
1413 {
1414 isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol );
1415 if (! isHalfSpace)
1416 break;
1417 LastVec = Normal;
1418 }
1419 } //for (i = 1; i <= NTLinCont; i++)
1420
1421 if (isHalfSpace)
1422 {
1423 for (i = 1; i <= NTPntCont; i++)
1424 {
1425 Standard_Integer Order = myPntCont->Value(i)->Order();
1426
1427 NewVecs.Clear();
1428 gp_Vec Vec1, Vec2, Normal;
1429 if (Order < GeomAbs_G1)
1430 myPntCont->Value(i)->D0( Pts->ChangeValue(++pnum) );
1431 else
1432 {
1433 myPntCont->Value(i)->D1( Pts->ChangeValue(++pnum), Vec1, Vec2 );
1434 Normal = Vec1 ^ Vec2;
1435 Normal.Normalize();
1436 Standard_Boolean isNew = Standard_True;
1437 for (Standard_Integer k = 1; k <= Vecs.Length(); k++)
1438 if (Vecs(k).IsEqual( Normal, LinTol, AngTol ))
1439 {
1440 isNew = Standard_False;
1441 break;
1442 }
1443 if (isNew)
1444 {
1445 NewVecs.Append( Normal );
1446 isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol );
1447 if (! isHalfSpace)
1448 {
1449 NewVecs(1).Reverse();
1450 isHalfSpace = GeomPlate_BuildAveragePlane::HalfSpace( NewVecs, Vecs, Aset, LinTol, AngTol );
1451 }
1452 if (! isHalfSpace)
1453 break;
1454 }
1455 }
1456 } //for (i = 1; i <= NTPntCont; i++)
1457
1458 if (isHalfSpace)
1459 {
1460 Standard_Boolean NullExist = Standard_True;
1461 while (NullExist)
1462 {
1463 NullExist = Standard_False;
1464 for (i = 1; i <= Vecs.Length(); i++)
1465 if (Vecs(i).SquareMagnitude() == 0.)
1466 {
1467 NullExist = Standard_True;
1468 Vecs.Remove(i);
1469 break;
1470 }
1471 }
1472 GeomPlate_BuildAveragePlane BAP( Vecs, Pts );
1473 Standard_Real u1,u2,v1,v2;
1474 BAP.MinMaxBox(u1,u2,v1,v2);
1475 // On agrandit le bazar pour les projections
1476 Standard_Real du = u2 - u1;
1477 Standard_Real dv = v2 - v1;
1478 u1 -= du; u2 += du; v1 -= dv; v2 += dv;
1479 mySurfInit = new Geom_RectangularTrimmedSurface( BAP.Plane(), u1,u2,v1,v2 );
1480 }
1481 } //if (isHalfSpace)
1482 if (!isHalfSpace)
1483 {
9230430c 1484#ifdef PLATE_DEB
7fd59977 1485 cout<<endl<<"Normals are not in half space"<<endl<<endl;
9230430c 1486#endif
7fd59977 1487 myIsLinear = Standard_False;
1488 nopt = 2;
1489 }
1490 } //if (NTLinCont != 0 && (CourbeJoint = CourbeJointive( myTol3d )) && IsOrderG1())
1491
1492 if (NTLinCont != 0)
1493 TrierTab( myInitOrder ); // Reordonne le tableau des permutations
1494
1495 if (nopt != 3)
1496 {
1497 if ( NTPntCont != 0)
1498 nopt = 1; //Calcul par la methode du plan d'inertie
1499 else if (!CourbeJoint || NTLinCont != myNbBounds)
1500 {// Standard_Failure::Raise("Curves are not joined");
9230430c 1501#ifdef PLATE_DEB
7fd59977 1502 cout<<"WARNING : Courbes non jointives a "<<myTol3d<<" pres"<<endl;
1503#endif
1504 nopt = 1;
1505 }
1506
1507 Standard_Real LenT=0;
1508 Standard_Integer Npt=0;
1509 Standard_Integer NTPoint=20*NTLinCont;
1510 Standard_Integer i ;
1511 for ( i=1;i<=NTLinCont;i++)
1512 LenT+=myLinCont->Value(i)->Length();
1513 for (i=1;i<=NTLinCont;i++)
1514 { Standard_Integer NbPoint= (Standard_Integer )( NTPoint*(myLinCont->Value(i)->Length())/LenT);
1515 if (NbPoint<10)
1516 NbPoint=10;
1517 Npt+=NbPoint;
1518 }
1519 // Tableau contenant le nuage de point pour le calcul du plan
1520 Handle(TColgp_HArray1OfPnt) Pts = new TColgp_HArray1OfPnt(1,20*NTLinCont+NTPntCont);
1521 Standard_Integer NbPoint=20;
1522 Standard_Real Uinit,Ufinal, Uif;
1523 for ( i=1;i<=NTLinCont;i++)
1524 { Uinit=myLinCont->Value(i)->FirstParameter();
1525 Ufinal=myLinCont->Value(i)->LastParameter();
1526 Uif=Ufinal-Uinit;
1527 if (mySense->Value(i) == 1)
1528 { Uinit = Ufinal;
1529 Uif = -Uif;
1530 }
1531 for (Standard_Integer j=0; j<NbPoint; j++)
1532 { // Nombre de point par courbe = 20
1533 // Repartition lineaire
1534 Standard_Real Inter=j*Uif/(NbPoint);
1535 gp_Pnt P;
1536 myLinCont->Value(i)->D0(Uinit+Inter,P);
1537 Pts->SetValue(Np++,P);
1538 }
1539 }
1540 for (i=1;i<=NTPntCont;i++)
1541 { gp_Pnt P;
1542 myPntCont->Value(i)->D0(P);
1543 Pts->SetValue(Np++,P);
1544 }
1545 if (!CourbeJoint)
1546 myNbBounds = 0;
1547 GeomPlate_BuildAveragePlane BAP( Pts, NbPoint*myNbBounds, myTol3d/1000, popt, nopt );
1548 if (!BAP.IsPlane()) Standard_Failure::Raise("the initial surface is not a plane.");
1549 Standard_Real u1,u2,v1,v2;
1550 BAP.MinMaxBox(u1,u2,v1,v2);
1551 // On agrandit le bazar pour les projections
1552 Standard_Real du = u2 - u1;
1553 Standard_Real dv = v2 - v1;
1554 u1 -= du; u2 += du; v1 -= dv; v2 += dv;
1555 mySurfInit = new Geom_RectangularTrimmedSurface(BAP.Plane(),u1,u2,v1,v2);
1556 } //if (nopt != 3)
1557
1558 //Comparing metrics of curves and projected curves
1559 if (NTLinCont != 0 && myIsLinear)
1560 {
1561 Handle( Geom_Surface ) InitPlane =
1562 (Handle( Geom_RectangularTrimmedSurface )::DownCast(mySurfInit))->BasisSurface();
1563
1564 Standard_Real Ratio, R1 = 2., R2 = 0.6; //R1 = 3, R2 = 0.5;//R1 = 1.4, R2 = 0.8; //R1 = 5., R2 = 0.2;
1565 Handle( GeomAdaptor_HSurface ) hsur =
1566 new GeomAdaptor_HSurface( InitPlane );
1567 Standard_Integer NbPoint = 20;
1568// gp_Pnt P;
1569// gp_Vec DerC, DerCproj, DU, DV;
1570// gp_Pnt2d P2d;
1571// gp_Vec2d DProj;
1572
1573
1574 for (Standard_Integer i = 1; i <= NTLinCont && myIsLinear; i++)
1575 {
1576 Standard_Real FirstPar = myLinCont->Value(i)->FirstParameter();
1577 Standard_Real LastPar = myLinCont->Value(i)->LastParameter();
1578 Standard_Real Uif = (LastPar - FirstPar)/(NbPoint);
1579
1580 Handle( Adaptor3d_HCurve ) Curve = myLinCont->Value(i)->Curve3d();
1581 ProjLib_CompProjectedCurve Projector( hsur, Curve, myTol3d, myTol3d );
1582 Handle( ProjLib_HCompProjectedCurve ) ProjCurve =
1583 new ProjLib_HCompProjectedCurve();
1584 ProjCurve->Set( Projector );
1585 Adaptor3d_CurveOnSurface AProj(ProjCurve, hsur);
1586
1587 gp_Pnt P;
1588 gp_Vec DerC, DerCproj;
1589 for (Standard_Integer j = 1 ; j < NbPoint && myIsLinear; j++)
1590 {
1591 Standard_Real Inter = FirstPar + j*Uif;
1592 Curve->D1(Inter, P, DerC);
1593 AProj.D1(Inter, P, DerCproj);
1594
1595 Standard_Real A1 = DerC.Magnitude();
1596 Standard_Real A2 = DerCproj.Magnitude();
1597 if (A2 <= 1.e-20) {
1598 Ratio = 1.e20;
1599 }
1600 else {
1601 Ratio = A1 / A2;
1602 }
1603 if (Ratio > R1 || Ratio < R2)
1604 {
1605 myIsLinear = Standard_False;
1606 break;
1607 }
1608 }
1609 }
9230430c 1610#if PLATE_DEB
7fd59977 1611 if (! myIsLinear)
1612 cout <<"Metrics are too different :"<< Ratio<<endl;
1613#endif
1614// myIsLinear = Standard_True; // !!
1615 } //comparing metrics of curves and projected curves
1616
1617
1618 if (! myIsLinear)
1619 {
1620 myPlanarSurfInit = mySurfInit;
1621#if DRAW
1622 if (Affich) {
1623 char* name = new char[100];
1624 sprintf(name,"planinit_%d",NbPlan+1);
1625 DrawTrSurf::Set(name, mySurfInit);
1626 }
1627#endif
1628 Standard_Real u1,v1,u2,v2;
1629 mySurfInit->Bounds(u1,v1,u2,v2);
1630 GeomAdaptor_Surface Surf(mySurfInit);
1631 myTolU = Surf.UResolution(myTol3d);
1632 myTolV = Surf.VResolution(myTol3d);
1633 myProj.Initialize(Surf,u1,v1,u2,v2,
1634 myTolU,myTolV);
1635
1636 //======================================================================
1637 // Projection des courbes
1638 //======================================================================
1639 Standard_Integer i;
1640 for (i = 1; i <= NTLinCont; i++)
1641 if (myLinCont->Value(i)->Curve2dOnSurf().IsNull())
1642 myLinCont->ChangeValue(i)->SetCurve2dOnSurf(ProjectCurve(myLinCont->Value(i)->Curve3d()));
1643
1644 //======================================================================
1645 // Projection des points
1646 //======================================================================
1647 for (i = 1; i<=NTPntCont; i++)
1648 { gp_Pnt P;
1649 myPntCont->Value(i)->D0(P);
1650 if (!myPntCont->Value(i)->HasPnt2dOnSurf())
1651 myPntCont->ChangeValue(i)->SetPnt2dOnSurf(ProjectPoint(P));
1652 }
1653
1654 //======================================================================
1655 // Nbre de point par courbe
1656 //======================================================================
1657 if ((NTLinCont !=0)&&(myNbPtsOnCur!=0))
1658 CalculNbPtsInit();
1659
1660 //======================================================================
1661 // Gestion des incompatibilites entre courbes
1662 //======================================================================
1663 Handle(GeomPlate_HArray1OfSequenceOfReal) PntInter;
1664 Handle(GeomPlate_HArray1OfSequenceOfReal) PntG1G1;
1665 if (NTLinCont != 0)
1666 {
1667 PntInter = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
1668 PntG1G1 = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
1669 Intersect(PntInter, PntG1G1);
1670 }
1671
1672 //====================================================================
1673 // Discretisation des courbes
1674 //====================================================================
1675 Discretise(PntInter, PntG1G1);
1676 //====================================================================
1677 //Preparation des points de contrainte pour plate
1678 //====================================================================
1679 LoadCurve( 0, 0 );
1680 if (myPntCont->Length() != 0)
1681 LoadPoint( 0, 0 );
1682 //====================================================================
1683 //Resolution de la surface
1684 //====================================================================
1685 myPlate.SolveTI(2, ComputeAnisotropie());
1686 if (!myPlate.IsDone())
1687 Standard_Failure::Raise("GeomPlate : abort calcul of Plate.");
1688
1689 myGeomPlateSurface = new GeomPlate_Surface( mySurfInit, myPlate );
1690
1691 GeomPlate_MakeApprox App(myGeomPlateSurface, myTol3d,
1692 1, 3,
1693 15*myTol3d,
1694 -1, GeomAbs_C0);
1695 mySurfInit = App.Surface();
1696
1697 mySurfInitIsGive = Standard_True;
1698 myPlate.Init(); // Reset
1699
1700 for (i=1;i<=NTLinCont;i++)
1701 {
1702 Handle( Geom2d_Curve ) NullCurve;
1703 NullCurve.Nullify();
1704 myLinCont->ChangeValue(i)->SetCurve2dOnSurf( NullCurve);
1705 }
1706 }
1707
1708#if DRAW
1709 if (Affich) {
1710 char* name = new char[100];
1711 sprintf(name,"surfinit_%d",++NbPlan);
1712 DrawTrSurf::Set(name, mySurfInit);
1713 }
1714#endif
1715}
1716
1717//---------------------------------------------------------
1718// fonction : Intersect
1719//---------------------------------------------------------
1720// Recherche les intersections entre les courbe 2d
1721// Si intersection et compatible( dans le cas G1-G1)
1722// enleve le point sur une des deux courbes
1723//---------------------------------------------------------
1724void GeomPlate_BuildPlateSurface::
1725Intersect(Handle(GeomPlate_HArray1OfSequenceOfReal)& PntInter,
1726 Handle(GeomPlate_HArray1OfSequenceOfReal)& PntG1G1)
1727{
1728 Standard_Integer NTLinCont = myLinCont->Length();
1729 Geom2dInt_GInter Intersection;
1730 Geom2dAdaptor_Curve Ci, Cj;
1731 IntRes2d_IntersectionPoint int2d;
1732 gp_Pnt P1,P2;
1733 gp_Pnt2d P2d;
1734 gp_Vec2d V2d;
1735// if (!mySurfInitIsGive)
1736 for (Standard_Integer i=1;i<=NTLinCont;i++)
1737 {
1738 //Standard_Real NbPnt_i=myLinCont->Value(i)->NbPoints();
1739 // Cherche l'intersection avec chaque courbe y compris elle meme.
1740 Ci.Load(myLinCont->Value(i)->Curve2dOnSurf());
1741 for(Standard_Integer j=i; j<=NTLinCont; j++)
1742 {
1743 Cj.Load(myLinCont->Value(j)->Curve2dOnSurf());
1744 if (i==j)
1745 Intersection.Perform(Ci, myTol2d*10, myTol2d*10);
1746 else
1747 Intersection.Perform(Ci, Cj, myTol2d*10, myTol2d*10);
1748
1749 if (!Intersection.IsEmpty())
1750 { // il existe une intersection
1751 Standard_Integer nbpt = Intersection.NbPoints();
1752 // nombre de point d'intersection
1753 for (Standard_Integer k = 1; k <= nbpt; k++)
1754 { int2d = Intersection.Point(k);
1755 myLinCont->Value(i)->D0(int2d.ParamOnFirst(),P1);
1756 myLinCont->Value(j)->D0(int2d.ParamOnSecond(),P2);
9230430c 1757#if PLATE_DEB
7fd59977 1758 if (Affich> 1)
1759 {
1760 cout << " Intersection "<< k << " entre " << i
1761 << " &" << j << endl;
1762 cout <<" Distance = "<< P1.Distance(P2) << endl;
1763 }
1764#endif
1765 if (P1.Distance( P2 ) < myTol3d)
1766 { // L'intersection 2d correspond a des points 3d proches
1767 // On note l'intervalle dans lequel il faudra enlever
1768 // un point pour eviter les doublons ce qui fait une
1769 // erreur dans plate. le point sur la courbe i est enleve
1770 // on conserve celui de la courbe j
1771 // la longueur de l'intervalle est une longueur 2d
1772 // correspondant en 3d a myTol3d
1773 Standard_Real tolint = Ci.Resolution(myTol3d);
1774 Ci.D1(int2d.ParamOnFirst(),P2d, V2d);
1775 Standard_Real aux = V2d.Magnitude();
1776 if (aux > 1.e-7)
1777 {
1778 aux = myTol3d/aux;
1779 if (aux > 100*tolint) tolint*=100;
1780 else tolint = aux;
1781 }
1782 else tolint*=100;
1783
1784 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() - tolint);
1785 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() + tolint);
1786 // Si G1-G1
1787 if ( (myLinCont->Value(i)->Order()==1)
1788 &&(myLinCont->Value(j)->Order()==1))
1789 { gp_Vec v11,v12,v13,v14,v15,v16,v21,v22,v23,v24,v25,v26;
1790 myLinCont->Value(i)->D2( int2d.ParamOnFirst(), P1, v11, v12, v13, v14, v15 );
1791 myLinCont->Value(j)->D2( int2d.ParamOnSecond(), P2, v21, v22, v23, v24, v25 );
1792 v16=v11^v12;v26=v21^v22;
1793 Standard_Real ant=v16.Angle(v26);
c6541a0c
D
1794 if (ant>(M_PI/2))
1795 ant= M_PI -ant;
7fd59977 1796 if ((Abs(v16*v15-v16*v25)>(myTol3d/1000))
1797 ||(Abs(ant)>myTol3d/1000))
1798 // Pas compatible ==> on enleve une zone en
1799 // contrainte G1 correspondant
1800 // a une tolerance 3D de 0.01
1801 { Standard_Real coin;
1802 Standard_Real Tol= 100 * myTol3d;
1803 Standard_Real A1;
1804 gp_Pnt2d P1,P2;
1805 gp_Vec2d V1,V2;
1806 myLinCont->Value(i)->Curve2dOnSurf()->D1( int2d.ParamOnFirst(), P1, V1);
1807 myLinCont->Value(j)->Curve2dOnSurf()->D1( int2d.ParamOnSecond(), P2, V2);
1808 A1 = V1.Angle(V2);
c6541a0c
D
1809 if (A1>(M_PI/2))
1810 A1= M_PI - A1;
1811 if (Abs(Abs(A1)-M_PI)<myTolAng) Tol = 100000 * myTol3d;
9230430c 1812#if PLATE_DEB
7fd59977 1813 if (Affich) cout <<"Angle entre Courbe "<<i<<","<<j
c6541a0c 1814 <<" "<<Abs(Abs(A1)-M_PI)<<endl;
7fd59977 1815#endif
1816
1817 coin = Ci.Resolution(Tol);
1818 Standard_Real Par1=int2d.ParamOnFirst()-coin,
1819 Par2=int2d.ParamOnFirst()+coin;
1820 // Stockage de l'intervalle pour la courbe i
1821 PntG1G1->ChangeValue(i).Append(Par1);
1822 PntG1G1->ChangeValue(i).Append(Par2);
1823 coin = Cj.Resolution(Tol);
1824 Par1=int2d.ParamOnSecond()-coin;
1825 Par2=int2d.ParamOnSecond()+coin;
1826 // Stockage de l'intervalle pour la courbe j
1827 PntG1G1->ChangeValue(j).Append(Par1);
1828 PntG1G1->ChangeValue(j).Append(Par2);
1829 }
1830 }
1831 // Si G0-G1
1832 if ((myLinCont->Value(i)->Order()==0 && myLinCont->Value(j)->Order()==1) ||
1833 (myLinCont->Value(i)->Order()==1 && myLinCont->Value(j)->Order()==0))
1834 {
1835 gp_Vec vec, vecU, vecV, N;
1836 if (myLinCont->Value(i)->Order() == 0)
1837 {
1838 Handle( Adaptor3d_HCurve ) theCurve = myLinCont->Value(i)->Curve3d();
1839 theCurve->D1( int2d.ParamOnFirst(), P1, vec );
1840 myLinCont->Value(j)->D1( int2d.ParamOnSecond(), P2, vecU, vecV );
1841 }
1842 else
1843 {
1844 Handle( Adaptor3d_HCurve ) theCurve = myLinCont->Value(j)->Curve3d();
1845 theCurve->D1( int2d.ParamOnSecond(), P2, vec );
1846 myLinCont->Value(i)->D1( int2d.ParamOnFirst(), P1, vecU, vecV );
1847 }
1848 N = vecU ^ vecV;
1849 Standard_Real Angle = vec.Angle( N );
c6541a0c 1850 Angle = Abs( M_PI/2-Angle );
7fd59977 1851 if (Angle > myTolAng/10.) //????????? //if (Abs( scal ) > myTol3d/100)
1852 { // Pas compatible ==> on enleve une zone en
1853 // contrainte G0 et G1 correspondant
1854 // a une tolerance 3D de 0.01
1855 Standard_Real coin;
1856 Standard_Real Tol= 100 * myTol3d;
1857 Standard_Real A1;
1858 gp_Pnt2d P1,P2;
1859 gp_Vec2d V1,V2;
1860 myLinCont->Value(i)->Curve2dOnSurf()->D1( int2d.ParamOnFirst(), P1, V1);
1861 myLinCont->Value(j)->Curve2dOnSurf()->D1( int2d.ParamOnSecond(), P2, V2);
1862 A1 = V1.Angle( V2 );
c6541a0c
D
1863 if (A1 > M_PI/2)
1864 A1= M_PI - A1;
1865 if (Abs(Abs(A1) - M_PI) < myTolAng) Tol = 100000 * myTol3d;
9230430c 1866#if PLATE_DEB
7fd59977 1867 if (Affich) cout <<"Angle entre Courbe "<<i<<","<<j
c6541a0c 1868 <<" "<<Abs(Abs(A1)-M_PI)<<endl;
7fd59977 1869#endif
1870 if (myLinCont->Value(i)->Order() == 1)
1871 {
1872 coin = Ci.Resolution(Tol);
1873 coin *= Angle / myTolAng * 10.;
9230430c 1874#if PLATE_DEB
7fd59977 1875 cout<<endl<<"coin = "<<coin<<endl;
9230430c 1876#endif
7fd59977 1877 Standard_Real Par1 = int2d.ParamOnFirst() - coin;
1878 Standard_Real Par2 = int2d.ParamOnFirst() + coin;
1879 // Stockage de l'intervalle pour la courbe i
1880 PntG1G1->ChangeValue(i).Append( Par1 );
1881 PntG1G1->ChangeValue(i).Append( Par2 );
1882 }
1883 else
1884 {
1885 coin = Cj.Resolution(Tol);
1886 coin *= Angle / myTolAng * 10.;
9230430c 1887#if PLATE_DEB
7fd59977 1888 cout<<endl<<"coin = "<<coin<<endl;
9230430c 1889#endif
7fd59977 1890 Standard_Real Par1 = int2d.ParamOnSecond() - coin;
1891 Standard_Real Par2 = int2d.ParamOnSecond() + coin;
1892 // Stockage de l'intervalle pour la courbe j
1893 PntG1G1->ChangeValue(j).Append( Par1 );
1894 PntG1G1->ChangeValue(j).Append( Par2 );
1895 }
1896 }
1897 }
1898 } //if (P1.Distance( P2 ) < myTol3d)
1899 else {
1900 //L'intersection 2d correspond a des points 3d eloigne
1901 // On note l'intervalle dans lequel il faudra enlever
1902 // des points pour eviter les doublons ce qui fait une
1903 // erreur dans plate. le point sur la courbe i est enleve
1904 // on conserve celui de la courbe j.
1905 // la longueur de l'intervalle est une longueur 2d
1906 // correspondant a la distance des points en 3d a myTol3d
1907 Standard_Real tolint, Dist;
1908 Dist = P1.Distance(P2);
1909 tolint = Ci.Resolution(Dist);
1910 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() - tolint);
1911 PntInter->ChangeValue(i).Append( int2d.ParamOnFirst() + tolint);
1912 if (j!=i)
1913 {
1914 tolint = Cj.Resolution(Dist);
1915 PntInter->ChangeValue(j).
1916 Append( int2d.ParamOnSecond() - tolint);
1917 PntInter->ChangeValue(j).
1918 Append( int2d.ParamOnSecond() + tolint);
1919 }
1920
9230430c 1921#ifdef PLATE_DEB
7fd59977 1922 cout<<"Attention: Deux points 3d ont la meme projection dist="
1923 <<Dist<<endl;
1924#endif
1925#ifdef DRAW
1926 if (Affich > 1)
1927 {
1928 Handle(Draw_Marker3D) mark = new (Draw_Marker3D)(P1, Draw_X, Draw_vert);
1929 char* name = new char[100];
1930 sprintf(name,"mark_%d",++NbMark);
1931 Draw::Set(name, mark);
1932 }
1933#endif
1934 }
1935 }
1936 }
1937 }
1938 }
1939}
1940
1941//---------------------------------------------------------
1942// fonction : Discretise
1943//---------------------------------------------------------
1944// Discretise les courbes selon le parametre de celle-ci
1945// le tableau des sequences Parcont contiendera tous les
1946// parametres des points sur les courbes
1947// Le champ myPlateCont contiendera les parametre des points envoye a plate
1948// il excluera les points en double et les zones imcompatibles.
1949// La premiere partie correspond a la verfication de la compatibilite
1950// et a la supprssion des points en double.
1951//---------------------------------------------------------
1952void GeomPlate_BuildPlateSurface::
1953Discretise(const Handle(GeomPlate_HArray1OfSequenceOfReal)& PntInter,
1954 const Handle(GeomPlate_HArray1OfSequenceOfReal)& PntG1G1)
1955{
1956 Standard_Integer NTLinCont = myLinCont->Length();
1957 Standard_Boolean ACR;
1958 Handle(Geom2d_Curve) C2d;
1959 Geom2dAdaptor_Curve AC2d;
1960// Handle(Adaptor_HCurve2d) HC2d;
1961 Handle(Law_Interpol) acrlaw = new (Law_Interpol) ();
1962 myPlateCont = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
1963 myParCont = new GeomPlate_HArray1OfSequenceOfReal(1,NTLinCont);
1964
1965
1966 //===========================================================================
1967 // Construction du tableau contenant les parametres des points de contrainte
1968 //===========================================================================
1969 Standard_Real Uinit, Ufinal, Length2d=0, Inter;
1970 Standard_Real CurLength;
1971 Standard_Integer NbPnt_i, NbPtInter, NbPtG1G1;
1972 Handle(GeomPlate_CurveConstraint) LinCont;
1973
1974 for (Standard_Integer i=1;i<=NTLinCont;i++) {
1975 LinCont = myLinCont->Value(i);
1976 Uinit=LinCont->FirstParameter();
1977 Ufinal=LinCont->LastParameter();
1978// HC2d=LinCont->ProjectedCurve();
1979// if(HC2d.IsNull())
1980// ACR = (!HC2d.IsNull() || !C2d.IsNull());
1981 C2d= LinCont->Curve2dOnSurf();
1982 ACR = (!C2d.IsNull());
1983 if (ACR) {
1984 // On Construit une loi proche de l'abscisse curviligne
1985 if(!C2d.IsNull()) AC2d.Load(C2d);
1986// AC2d.Load(LinCont->Curve2dOnSurf());
1987 Standard_Integer ii, Nbint = 20;
1988 Standard_Real U;
1989 TColgp_Array1OfPnt2d tabP2d(1, Nbint+1);
1990 tabP2d(1).SetY(Uinit);
1991 tabP2d(1).SetX(0.);
1992 tabP2d(Nbint+1).SetY(Ufinal);
1993/* if (!HC2d.IsNull())
1994
1995 Length2d = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, Ufinal);
1996 else*/
1997 Length2d = GCPnts_AbscissaPoint::Length(AC2d, Uinit, Ufinal);
1998
1999 tabP2d(Nbint+1).SetX(Length2d);
2000 for (ii = 2; ii<= Nbint; ii++) {
c6541a0c 2001 U = Uinit + (Ufinal-Uinit)*((1-Cos((ii-1)*M_PI/(Nbint)))/2);
7fd59977 2002 tabP2d(ii).SetY(U);
2003/* if (!HC2d.IsNull()) {
2004 Standard_Real L = GCPnts_AbscissaPoint::Length(HC2d->Curve2d(), Uinit, U);
2005 tabP2d(ii).SetX(L);
2006 }
2007 else*/
2008 tabP2d(ii).SetX(GCPnts_AbscissaPoint::Length(AC2d, Uinit, U));
2009 }
2010 acrlaw->Set(tabP2d);
2011 }
2012
2013 NbPnt_i = (Standard_Integer )( LinCont->NbPoints());
2014 NbPtInter= PntInter->Value(i).Length();
2015 NbPtG1G1= PntG1G1->Value(i).Length();
2016
9230430c 2017#if PLATE_DEB
7fd59977 2018 if (Affich > 1) {
2019 cout << "Courbe : " << i << endl;
2020 cout << " NbPnt, NbPtInter, NbPtG1G1 :" << NbPnt_i << ", "
2021 << NbPtInter << ", " << NbPtG1G1 << endl;
2022 }
2023#endif
2024 for (Standard_Integer j=1; j<=NbPnt_i; j++)
2025 { // repartition des points en cosinus selon l'ACR 2d
2026 // Afin d'eviter les points d'acumulation dans le 2d
c6541a0c 2027 //Inter=Uinit+(Uif)*((-cos(M_PI*((j-1)/(NbPnt_i-1)))+1)/2);
7fd59977 2028 if (j==NbPnt_i)
2029 Inter=Ufinal;//pour parer au bug sur sun
2030 else if (ACR) {
c6541a0c 2031 CurLength = Length2d*(1-Cos((j-1)*M_PI/(NbPnt_i-1)))/2;
7fd59977 2032 Inter = acrlaw->Value(CurLength);
2033 }
2034 else {
c6541a0c 2035 Inter=Uinit+(Ufinal-Uinit)*((1-Cos((j-1)*M_PI/(NbPnt_i-1)))/2);
7fd59977 2036 }
2037 myParCont->ChangeValue(i).Append(Inter);// on ajoute le point
2038 if (NbPtInter!=0)
2039 { for(Standard_Integer l=1;l<=NbPtInter;l+=2)
2040 //on cherche si le point Inter est dans l'intervalle
2041 //PntInter[i] PntInter[i+1]
2042 //auquelle cas il ne faudrait pas le stocker (pb de doublons)
2043 { if ((Inter>PntInter->Value(i).Value(l))
2044 &&(Inter<PntInter->Value(i).Value(l+1)))
2045 { l=NbPtInter+2;
2046 // pour sortir de la boucle sans stocker le point
2047 }
2048 else
2049 { if (l+1>=NbPtInter)
2050 // on a parcouru tout le tableau : Le point
2051 // n'appartient pas a un interval point commun
2052 if (NbPtG1G1!=0)
2053 // est qu'il existe un intervalle incompatible
2054 for(Standard_Integer k=1;k<=NbPtG1G1;k+=2)
2055 { if ((Inter>PntG1G1->Value(i).Value(k))
2056 &&(Inter<PntG1G1->Value(i).Value(k+1)))
2057 { k=NbPtG1G1+2; // pour sortir de la boucle
2058 // Ajouter les points de contrainte G0
2059 gp_Pnt P3d,PP,Pdif;
2060 gp_Pnt2d P2d;
2061
2062 AC2d.D0(Inter, P2d);
2063 LinCont->D0(Inter,P3d);
2064 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2065 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2066 -PP.Coord(2)+P3d.Coord(2),
2067 -PP.Coord(3)+P3d.Coord(3));
2068 Plate_PinpointConstraint PC(P2d.XY(),
2069 Pdif.XYZ(),0,0);
2070 myPlate.Load(PC);
2071
2072 }
2073 else // le point n'appartient pas a un interval G1
2074 { if (k+1>=NbPtG1G1)
2075 { myPlateCont->ChangeValue(i).Append(Inter);
2076 // on ajoute le point
2077 }
2078 }
2079 }
2080 else
2081 { myPlateCont->ChangeValue(i).Append(Inter);
2082 // on ajoute le point
2083 }
2084 }
2085 }
2086 }
2087 else
2088 { if (NbPtG1G1!=0) // est qu'il existe un intervalle incompatible
2089 for(Standard_Integer k=1;k<=NbPtG1G1;k+=2)
2090 { if ((Inter>PntG1G1->Value(i).Value(k))
2091 &&(Inter<PntG1G1->Value(i).Value(k+1)))
2092 { k=NbPtG1G1+2; // pour sortir de la boucle
2093 // Ajouter les points de contrainte G0
2094 gp_Pnt P3d,PP,Pdif;
2095 gp_Pnt2d P2d;
2096
2097 AC2d.D0(Inter, P2d);
2098 LinCont->D0(Inter,P3d);
2099 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2100 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2101 -PP.Coord(2)+P3d.Coord(2),
2102 -PP.Coord(3)+P3d.Coord(3));
2103 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2104 myPlate.Load(PC);
2105
2106 }
2107 else // le point n'appartient pas a un intervalle G1
2108 { if (k+1>=NbPtG1G1)
2109 { myPlateCont->ChangeValue(i).Append(Inter);
2110 // on ajoute le point
2111 }
2112 }
2113 }
2114 else
2115 { if ( ( (!mySurfInitIsGive)
2116 &&(Geom2dAdaptor_Curve(LinCont->Curve2dOnSurf()).GetType()!=GeomAbs_Circle))
2117 || ( (j>1) &&(j<NbPnt_i))) //on enleve les extremites
2118 myPlateCont->ChangeValue(i).Append(Inter);// on ajoute le point
2119 }
2120 }
2121 }
2122 }
2123}
2124//---------------------------------------------------------
2125// fonction : CalculNbPtsInit
2126//---------------------------------------------------------
2127// Calcul du nombre de points par courbe en fonction de la longueur
2128// pour la premiere iteration
2129//---------------------------------------------------------
2130void GeomPlate_BuildPlateSurface::CalculNbPtsInit ()
2131{
2132 Standard_Real LenT=0;
2133 Standard_Integer NTLinCont=myLinCont->Length();
2134 Standard_Integer NTPoint=(Standard_Integer )( myNbPtsOnCur*NTLinCont);
2135 Standard_Integer i ;
2136
2137 for ( i=1;i<=NTLinCont;i++)
2138 LenT+=myLinCont->Value(i)->Length();
2139 for ( i=1;i<=NTLinCont;i++)
2140 { Standard_Integer Cont=myLinCont->Value(i)->Order();
2141 switch(Cont)
2142 { case 0 : // Cas G0 *1.2
2143 myLinCont->ChangeValue(i)->SetNbPoints(
2144 Standard_Integer(1.2*NTPoint*(myLinCont->Value(i)->Length())/LenT));
2145 break;
2146 case 1 : // Cas G1 *1
2147 myLinCont->ChangeValue(i)->SetNbPoints(
2148 Standard_Integer(NTPoint*(myLinCont->Value(i)->Length())/LenT));
2149 break;
2150 case 2 : // Cas G2 *0.7
2151 myLinCont->ChangeValue(i)->SetNbPoints(
2152 Standard_Integer(0.7*NTPoint*(myLinCont->Value(i)->Length())/LenT));
2153 break;
2154 }
2155 if (myLinCont->Value(i)->NbPoints()<3)
2156 myLinCont->ChangeValue(i)->SetNbPoints(3);
2157 }
2158}
2159//---------------------------------------------------------
2160// fonction : LoadCurve
2161//---------------------------------------------------------
2162// A partir du tableau myParCont on charge tous les points notes dans plate
2163//---------------------------------------------------------
2164void GeomPlate_BuildPlateSurface::LoadCurve(const Standard_Integer NbBoucle,
2165 const Standard_Integer OrderMax )
2166{
2167 gp_Pnt P3d,Pdif,PP;
2168 gp_Pnt2d P2d;
2169 Standard_Integer NTLinCont=myLinCont->Length(), i, j;
2170 Standard_Integer Tang, Nt;
2171
2172
2173 for (i=1; i<=NTLinCont; i++){
2174 Handle(GeomPlate_CurveConstraint) CC = myLinCont->Value(i);
2175 if (CC ->Order()!=-1) {
2176 Tang = Min(CC->Order(), OrderMax);
2177 Nt = myPlateCont->Value(i).Length();
2178 if (Tang!=-1)
2179 for (j=1; j<=Nt; j++) {// Chargement des points G0 sur les frontieres
2180 CC ->D0(myPlateCont->Value(i).Value(j),P3d);
2181 if (!CC->ProjectedCurve().IsNull())
2182 P2d = CC->ProjectedCurve()->Value(myPlateCont->Value(i).Value(j));
2183
2184 else {
2185 if (!CC->Curve2dOnSurf().IsNull())
2186 P2d = CC->Curve2dOnSurf()->Value(myPlateCont->Value(i).Value(j));
2187 else
2188 P2d = ProjectPoint(P3d);
2189 }
2190 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2191 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2192 -PP.Coord(2)+P3d.Coord(2),
2193 -PP.Coord(3)+P3d.Coord(3));
2194 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2195 myPlate.Load(PC);
2196
2197 // Chargement des points G1
2198 if (Tang==1) { // ==1
2199 gp_Vec V1,V2,V3,V4;
2200 CC->D1( myPlateCont->Value(i).Value(j), PP, V1, V2 );
2201 mySurfInit->D1( P2d.Coord(1), P2d.Coord(2), PP, V3, V4 );
2202
2203 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2204 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2205 if (! myFree)
2206 {
2207 Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final );
2208 myPlate.Load(GCC);
2209 }
2210 else if (NbBoucle == 1)
2211 {
2212 Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final );
2213 myPlate.Load( FreeGCC );
2214 }
2215 else {
2216 gp_Vec DU, DV, Normal, DerPlateU, DerPlateV;
2217
2218 Normal = V1^V2;
2219 //Normal.Normalize();
2220 Standard_Real norm = Normal.Magnitude();
2221 if (norm > 1.e-12) Normal /= norm;
2222 DerPlateU = myPrevPlate.EvaluateDerivative( P2d.XY(), 1, 0 );
2223 DerPlateV = myPrevPlate.EvaluateDerivative( P2d.XY(), 0, 1 );
2224
2225 DU.SetLinearForm(-(V3 + DerPlateU).Dot(Normal), Normal, DerPlateU);
2226 DV.SetLinearForm(-(V4 + DerPlateV).Dot(Normal), Normal, DerPlateV);
2227 Plate_PinpointConstraint PinU( P2d.XY(), DU.XYZ(), 1, 0 );
2228 Plate_PinpointConstraint PinV( P2d.XY(), DV.XYZ(), 0, 1 );
2229 myPlate.Load( PinU );
2230 myPlate.Load( PinV );
2231 }
2232 }
2233 // Chargement des points G2
2234 if (Tang==2) // ==2
2235 { gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
2236 CC->D2(myPlateCont->Value(i).Value(j),
2237 PP,V1,V2,V5,V6,V7);
2238 mySurfInit->D2(P2d.Coord(1),P2d.Coord(2),PP,V3,V4,V8,V9,V10);
2239
2240 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2241 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2242 Plate_D2 D2final (V5.XYZ(),V6.XYZ(),V7.XYZ());
2243 Plate_D2 D2init (V8.XYZ(),V9.XYZ(),V10.XYZ());
2244// if (! myFree)
2245// {
2246 Plate_GtoCConstraint GCC(P2d.XY(),D1init,D1final,D2init,D2final);
2247 myPlate.Load(GCC);
2248// }
2249// else // Good but too expansive
2250// {
2251// Plate_FreeGtoCConstraint FreeGCC( P2d.XY(),
2252// D1init, D1final, D2init, D2final );
2253// myPlate.Load( FreeGCC );
2254// }
2255
2256 }
2257 }
2258 }
2259 }
2260}
2261
2262
2263//---------------------------------------------------------
2264//fonction : LoadPoint
2265//---------------------------------------------------------
2266//void GeomPlate_BuildPlateSurface::LoadPoint(const Standard_Integer NbBoucle,
2267void GeomPlate_BuildPlateSurface::LoadPoint(const Standard_Integer ,
2268 const Standard_Integer OrderMax)
2269{
2270 gp_Pnt P3d,Pdif,PP;
2271 gp_Pnt2d P2d;
2272 Standard_Integer NTPntCont=myPntCont->Length();
2273 Standard_Integer Tang, i;
2274// gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
2275 gp_Vec V1,V2,V3,V4;
2276
2277 // Chargement des points de contraintes ponctuel
2278 for (i=1;i<=NTPntCont;i++) {
2279 myPntCont->Value(i)->D0(P3d);
2280 P2d = myPntCont->Value(i)->Pnt2dOnSurf();
2281 mySurfInit->D0(P2d.Coord(1),P2d.Coord(2),PP);
2282 Pdif.SetCoord(-PP.Coord(1)+P3d.Coord(1),
2283 -PP.Coord(2)+P3d.Coord(2),
2284 -PP.Coord(3)+P3d.Coord(3));
2285 Plate_PinpointConstraint PC(P2d.XY(),Pdif.XYZ(),0,0);
2286 myPlate.Load(PC);
2287 Tang = Min(myPntCont->Value(i)->Order(), OrderMax);
2288 if (Tang==1) {// ==1
2289 myPntCont->Value(i)->D1(PP,V1,V2);
2290 mySurfInit->D1(P2d.Coord(1),P2d.Coord(2),PP,V3,V4);
2291 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2292 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2293 if (! myFree)
2294 {
2295 Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final );
2296 myPlate.Load( GCC );
2297 }
2298 else
2299 {
2300 Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final );
2301 myPlate.Load( FreeGCC );
2302 }
2303 }
2304 // Chargement des points G2 GeomPlate_PlateG0Criterion
2305 if (Tang==2) // ==2
2306 { gp_Vec V1,V2,V3,V4,V5,V6,V7,V8,V9,V10;
2307 myPntCont->Value(i)->D2(PP,V1,V2,V5,V6,V7);
2308// gp_Vec Tv2 = V1^V2;
2309 mySurfInit->D2(P2d.Coord(1),P2d.Coord(2),PP,V3,V4,V8,V9,V10);
2310 Plate_D1 D1final(V1.XYZ(),V2.XYZ());
2311 Plate_D1 D1init(V3.XYZ(),V4.XYZ());
2312 Plate_D2 D2final (V5.XYZ(),V6.XYZ(),V7.XYZ());
2313 Plate_D2 D2init (V8.XYZ(),V9.XYZ(),V10.XYZ());
2314// if (! myFree)
2315// {
2316 Plate_GtoCConstraint GCC( P2d.XY(), D1init, D1final, D2init, D2final );
2317 myPlate.Load( GCC );
2318// }
2319// else // Good but too expansive
2320// {
2321// Plate_FreeGtoCConstraint FreeGCC( P2d.XY(), D1init, D1final, D2init//, D2final );
2322// myPlate.Load( FreeGCC );
2323// }
2324 }
2325 }
2326
2327}
2328//---------------------------------------------------------
2329//fonction : VerifSurface
2330//---------------------------------------------------------
2331Standard_Boolean GeomPlate_BuildPlateSurface::
2332VerifSurface(const Standard_Integer NbBoucle)
2333{
2334 //======================================================================
2335 // Calcul des erreurs
2336 //======================================================================
2337 Standard_Integer NTLinCont=myLinCont->Length();
2338 Standard_Boolean Result=Standard_True;
2339
2340 // variable pour les calculs d erreur
2341 myG0Error=0,myG1Error =0, myG2Error=0;
2342
2343 for (Standard_Integer i=1;i<=NTLinCont;i++) {
2344 Handle(GeomPlate_CurveConstraint) LinCont;
2345 LinCont = myLinCont->Value(i);
2346 if (LinCont->Order()!=-1) {
2347 Standard_Integer NbPts_i = myParCont->Value(i).Length();
2348 if (NbPts_i<3)
2349 NbPts_i=4;
2350 Handle(TColStd_HArray1OfReal) tdist =
2351 new TColStd_HArray1OfReal(1,NbPts_i-1);
2352 Handle(TColStd_HArray1OfReal) tang =
2353 new TColStd_HArray1OfReal(1,NbPts_i-1);
2354 Handle(TColStd_HArray1OfReal) tcourb =
2355 new TColStd_HArray1OfReal(1,NbPts_i-1);
7fd59977 2356
2357 EcartContraintesMil (i,tdist,tang,tcourb);
2358
2359 Standard_Real diffDistMax=0,SdiffDist=0;
2360 Standard_Real diffAngMax=0,SdiffAng=0;
2361 Standard_Integer NdiffDist=0,NdiffAng=0;
2362
2363
2364 for (Standard_Integer j=1;j<NbPts_i;j++)
2365 { if (tdist->Value(j)>myG0Error)
2366 myG0Error=tdist->Value(j);
2367 if (tang->Value(j)>myG1Error)
2368 myG1Error=tang->Value(j);
2369 if (tcourb->Value(j)>myG2Error)
2370 myG2Error=tcourb->Value(j);
2371 Standard_Real U;
2372 if (myParCont->Value(i).Length()>3)
2373 U=(myParCont->Value(i).Value(j)+myParCont->Value(i).Value(j+1))/2;
2374 else
2375 U=LinCont->FirstParameter()+
2376 ( LinCont->LastParameter()-LinCont->FirstParameter())*(j-1)/(NbPts_i-2);
2377 Standard_Real diffDist=tdist->Value(j)-LinCont->G0Criterion(U),diffAng;
2378 if (LinCont->Order()>0)
2379 diffAng=tang->Value(j)-LinCont->G1Criterion(U);
2380 else diffAng=0;
2381 // recherche de la variation de l'erreur maxi et calcul de la moyenne
2382 if (diffDist>0) {
2383 diffDist = diffDist/LinCont->G0Criterion(U);
2384 if (diffDist>diffDistMax)
2385 diffDistMax = diffDist;
2386 SdiffDist+=diffDist;
2387 NdiffDist++;
2388#if DRAW
2389 if ((Affich) && (NbBoucle == myNbIter)) {
2390 gp_Pnt P;
2391 gp_Pnt2d P2d;
2392 LinCont->D0(U,P);
2393 Handle(Draw_Marker3D) mark =
2394 new (Draw_Marker3D)(P, Draw_X, Draw_orange);
2395 char* name = new char[100];
2396
2397 sprintf(name,"mark_%d",++NbMark);
2398 Draw::Set(name, mark);
2399 if (!LinCont->ProjectedCurve().IsNull())
2400 P2d = LinCont->ProjectedCurve()->Value(U);
2401 else
2402 { if (!LinCont->Curve2dOnSurf().IsNull())
2403 P2d = LinCont->Curve2dOnSurf()->Value(U);
2404 else
2405 P2d = ProjectPoint(P);
2406 }
2407 sprintf(name,"mark2d_%d",++NbMark);
2408 Handle(Draw_Marker2D) mark2d =
2409 new (Draw_Marker2D)(P2d, Draw_X, Draw_orange);
2410 Draw::Set(name, mark2d);
2411 }
2412#endif
2413 }
2414 else
2415 if ((diffAng>0)&&(LinCont->Order()==1)) {
2416 diffAng=diffAng/myLinCont->Value(i)->G1Criterion(U);
2417 if (diffAng>diffAngMax)
2418 diffAngMax = diffAng;
2419 SdiffAng+=diffAng;
2420 NdiffAng++;
2421#if DRAW
2422 if ((Affich) && (NbBoucle == myNbIter)) {
2423 gp_Pnt P;
2424 LinCont->D0(U,P);
2425 Handle(Draw_Marker3D) mark =
2426 new Draw_Marker3D(P, Draw_X, Draw_or);
2427 char* name = new char[100];
2428 sprintf(name,"mark_%d",++NbMark);
2429 Draw::Set(name, mark);
2430 }
2431#endif
2432 }
2433 }
2434
2435 if (NdiffDist>0) {// au moins un point n'est pas acceptable en G0
2436 Standard_Real Coef;
2437 if(LinCont->Order()== 0)
2438 Coef = 0.6*Log(diffDistMax+7.4);
2439 //7.4 correspond au calcul du coef mini = 1.2 donc e^1.2/0.6
2440 else
2441 Coef = Log(diffDistMax+3.3);
2442 //3.3 correspond au calcul du coef mini = 1.2 donc e^1.2
2443 if (Coef>3)
2444 Coef=3;
2445 //experimentalement apres le coef devient mauvais pour les cas L
2446 if ((NbBoucle>1)&&(diffDistMax>2))
2447 { Coef=1.6;
2448 }
2449
2450 if (LinCont->NbPoints()>=Floor(LinCont->NbPoints()*Coef))
2451 Coef=2;// pour assurer une augmentation du nombre de points
2452
2453 LinCont->SetNbPoints(Standard_Integer(LinCont->NbPoints() * Coef));
2454 Result=Standard_False;
2455 }
2456 else
2457 if (NdiffAng>0) // au moins 1 points ne sont pas accepable en G1
2458 { Standard_Real Coef=1.5;
2459 if ((LinCont->NbPoints()+1)>=Floor(LinCont->NbPoints()*Coef))
2460 Coef=2;
2461
2462 LinCont->SetNbPoints(Standard_Integer(LinCont->NbPoints()*Coef )) ;
2463 Result=Standard_False;
2464 }
2465 }
2466 }
2467 if (!Result)
2468 {
2469 if (myFree && NbBoucle == 1)
2470 myPrevPlate = myPlate;
2471 myPlate.Init();
2472 }
2473 return Result;
2474}
2475
2476
2477
2478//---------------------------------------------------------
2479//fonction : VerifPoint
2480//---------------------------------------------------------
2481void GeomPlate_BuildPlateSurface::
2482 VerifPoints (Standard_Real& Dist,
2483 Standard_Real& Ang,
2484 Standard_Real& Curv) const
2485{ Standard_Integer NTPntCont=myPntCont->Length();
2486 gp_Pnt Pi, Pf;
2487 gp_Pnt2d P2d;
2488 gp_Vec v1i, v1f, v2i, v2f, v3i, v3f;
2489 Ang=0;Dist=0,Curv=0;
2490 Handle(GeomPlate_PointConstraint) PntCont;
2491 for (Standard_Integer i=1;i<=NTPntCont;i++) {
2492 PntCont = myPntCont->Value(i);
2493 switch (PntCont->Order())
2494 { case 0 :
2495 P2d = PntCont->Pnt2dOnSurf();
2496 PntCont->D0(Pi);
2497 myGeomPlateSurface->D0( P2d.Coord(1), P2d.Coord(2), Pf);
2498 Dist = Pf.Distance(Pi);
2499 break;
2500 case 1 :
2501 PntCont->D1(Pi, v1i, v2i);
2502 P2d = PntCont->Pnt2dOnSurf();
2503 myGeomPlateSurface->D1( P2d.Coord(1), P2d.Coord(2), Pf, v1f, v2f);
2504 Dist = Pf.Distance(Pi);
2505 v3i = v1i^v2i; v3f=v1f^v2f;
2506 Ang=v3f.Angle(v3i);
c6541a0c
D
2507 if (Ang>(M_PI/2))
2508 Ang = M_PI -Ang;
7fd59977 2509 break;
2510 case 2 :
2511 Handle(Geom_Surface) Splate;
2512 Splate = Handle(Geom_Surface)::DownCast(myGeomPlateSurface);
2513 LocalAnalysis_SurfaceContinuity CG2;
2514 P2d = PntCont->Pnt2dOnSurf();
2515 GeomLProp_SLProps Prop (Splate, P2d.Coord(1), P2d.Coord(2),
2516 2, 0.001);
2517 CG2.ComputeAnalysis(Prop, PntCont->LPropSurf(),GeomAbs_G2);
2518 Dist=CG2.C0Value();
2519 Ang=CG2.G1Angle();
2520 Curv=CG2.G2CurvatureGap();
2521 break;
2522 }
2523 }
2524}
2525
2526Standard_Real GeomPlate_BuildPlateSurface::ComputeAnisotropie() const
2527{
2528 if (myAnisotropie)
2529 {
2530 //Temporary
2531 return 1.0;
2532 }
2533 else return 1.0;
2534}
2535
2536Standard_Boolean GeomPlate_BuildPlateSurface::IsOrderG1() const
2537{
2538 Standard_Boolean result = Standard_True;
2539 for (Standard_Integer i = 1; i <= myLinCont->Length(); i++)
2540 if (myLinCont->Value(i)->Order() < 1)
2541 {
2542 result = Standard_False;
2543 break;
2544 }
2545 return result;
2546}