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