src/GeomPlate/GeomPlate_BuildPlateSurface.cxx

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