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