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