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