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