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