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