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