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