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