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