src/GeomPlate/GeomPlate_BuildPlateSurface.cxx

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