src/GeomFill/GeomFill_ConstrainedFilling.cxx

   1 // Created on: 1995-10-26
   2 // Created by: Laurent BOURESCHE
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 //  Modified by skv - Fri Jun 18 12:52:54 2004 OCC6129
  18
  19 #include <GeomFill_ConstrainedFilling.ixx>
  20
  21 #include <Standard_Failure.hxx>
  22 #include <Standard_NotImplemented.hxx>
  23 #include <TColStd_HArray1OfReal.hxx>
  24 #include <TColgp_Array1OfPnt.hxx>
  25 #include <gp_XYZ.hxx>
  26 #include <PLib.hxx>
  27 #include <BSplCLib.hxx>
  28 #include <AdvApprox_ApproxAFunction.hxx>
  29 #include <Law.hxx>
  30 #include <Law_Linear.hxx>
  31 #include <Law_BSpline.hxx>
  32 #include <Law_BSpFunc.hxx>
  33 #include <GeomFill_DegeneratedBound.hxx>
  34 #include <GeomFill_TgtOnCoons.hxx>
  35
  36 #ifdef DRAW
  37 // Pour le dessin.
  38 #include <Draw_Appli.hxx>
  39 #include <Draw_Display.hxx>
  40 #include <Draw.hxx>
  41 #include <Draw_Segment3D.hxx>
  42 #include <Draw_Segment2D.hxx>
  43 #include <Draw_Marker2D.hxx>
  44 #include <Draw_ColorKind.hxx>
  45 #include <Draw_MarkerShape.hxx>
  46 static Standard_Boolean dodraw = 0;
  47 static Standard_Real drawfac = 0.1;
  48 #endif
  49 #ifdef OCCT_DEBUG
  50 Standard_IMPORT void Law_draw1dcurve(const TColStd_Array1OfReal&    pol,
  51                             const TColStd_Array1OfReal&    knots,
  52                             const TColStd_Array1OfInteger& mults,
  53                             const Standard_Integer         deg,
  54                             const gp_Vec2d&                tra,
  55                             const Standard_Real            scal);
  56 Standard_IMPORT void Law_draw1dcurve(const Handle(Law_BSpline)&     bs,
  57                             const gp_Vec2d&                tra,
  58                             const Standard_Real            scal);
  59
  60
  61 // Pour les mesures.
  62 #include <OSD_Chronometer.hxx>
  63 static OSD_Chronometer totclock, parclock, appclock, cstclock;
  64 #endif
  65
  66 static Standard_Integer inqadd(const Standard_Real    d1,
  67                                const Standard_Real    d2,
  68                                Standard_Real*         k,
  69                                Standard_Integer*      m,
  70                                const Standard_Integer deg,
  71                                const Standard_Real    tolk)
  72 {
  73   Standard_Integer nbadd = 0;
  74   m[0] = m[1] = deg - 2;
  75   if (d1 != 1. && d2 != 1.){
  76     if(Abs(d1+d2-1.) < tolk) {
  77       k[0] = 0.5 * (d1 + 1. - d2);
  78       nbadd = 1;
  79     }
  80     else{
  81       nbadd = 2;
  82       k[0] = Min(d1,1. - d2);
  83       k[1] = Max(d1,1. - d2);
  84     }
  85   }
  86   else if (d1 != 1.) {
  87     k[0] = d1;
  88     nbadd = 1;
  89   }
  90   else if (d2 != 1.) {
  91     k[0] = d2;
  92     nbadd = 1;
  93   }
  94   return nbadd;
  95 }
  96
  97 static Handle(Law_Linear) mklin(const Handle(Law_Function)& func)
  98 {
  99   Handle(Law_Linear) fu = Handle(Law_Linear)::DownCast(func);
 100   if(fu.IsNull()) {
 101     fu = new Law_Linear();
 102     Standard_Real d,f;
 103     func->Bounds(d,f);
 104     fu->Set(d,func->Value(d),f,func->Value(f));
 105   }
 106   return fu;
 107 }
 108
 109 static void sortbounds(const Standard_Integer     nb,
 110                        Handle(GeomFill_Boundary)* bound,
 111                        Standard_Boolean*          rev,
 112                        GeomFill_CornerState*      stat)
 113 {
 114   // trier les bords (facon bourinos),
 115   // flaguer ceux a renverser,
 116   // flaguer les baillements au coins.
 117   Standard_Integer i,j;
 118   Handle(GeomFill_Boundary) temp;
 119   rev[0] = 0;
 120   gp_Pnt pf,pl;
 121   gp_Pnt qf,ql;
 122   for (i = 0; i < nb-1; i++){
 123     if(!rev[i]) bound[i]->Points(pf,pl);
 124     else bound[i]->Points(pl,pf);
 125     for (j = i+1; j <= nb-1; j++){
 126       bound[j]->Points(qf,ql);
 127       //  Modified by skv - Fri Jun 18 12:52:54 2004 OCC6129 Begin
 128       Standard_Real df = qf.Distance(pl);
 129       Standard_Real dl = ql.Distance(pl);
 130       if (df<dl) {
 131         if(df < stat[i+1].Gap()){
 132           temp = bound[i+1];
 133           bound[i+1] = bound[j];
 134           bound[j] = temp;
 135           stat[i+1].Gap(df);
 136           rev[i+1] = Standard_False;
 137         }
 138       } else {
 139         if(dl < stat[i+1].Gap()){
 140           temp = bound[i+1];
 141           bound[i+1] = bound[j];
 142           bound[j] = temp;
 143           stat[i+1].Gap(dl);
 144           rev[i+1] = Standard_True;
 145         }
 146       }
 147       //  Modified by skv - Fri Jun 18 12:52:54 2004 OCC6129 End
 148     }
 149   }
 150   if(!rev[nb-1]) bound[nb-1]->Points(pf,pl);
 151   else bound[nb-1]->Points(pl,pf);
 152   bound[0]->Points(qf,ql);
 153   stat[0].Gap(pl.Distance(qf));
 154
 155   // flaguer les angles entre tangentes au coins et entre les normales au
 156   // coins pour les bords contraints.
 157   gp_Pnt pbid;
 158   gp_Vec tgi, nori, tgn, norn;
 159   Standard_Real fi, fn, li, ln;
 160   for (i = 0; i < nb; i++){
 161     Standard_Integer next = (i+1)%nb;
 162     if(!rev[i]) bound[i]->Bounds(fi,li);
 163     else bound[i]->Bounds(li,fi);
 164     bound[i]->D1(li,pbid,tgi);
 165     if(rev[i]) tgi.Reverse();
 166     if(!rev[next]) bound[next]->Bounds(fn,ln);
 167     else bound[next]->Bounds(ln,fn);
 168     bound[next]->D1(fn,pbid,tgn);
 169     if(rev[next]) tgn.Reverse();
 170     Standard_Real ang = M_PI - tgi.Angle(tgn);
 171     stat[next].TgtAng(ang);
 172     if(bound[i]->HasNormals() && bound[next]->HasNormals()){
 173       stat[next].Constraint();
 174       nori = bound[i]->Norm(li);
 175       norn = bound[next]->Norm(fn);
 176       ang = nori.Angle(norn);
 177       stat[next].NorAng(ang);
 178     }
 179   }
 180 }
 181 static void coonscnd(const Standard_Integer     nb,
 182                      Handle(GeomFill_Boundary)* bound,
 183                      Standard_Boolean*          rev,
 184                      GeomFill_CornerState*      stat,
 185 //                   Handle(GeomFill_TgtField)* tga,
 186                      Handle(GeomFill_TgtField)* ,
 187                      Standard_Real*             mintg)
 188 {
 189   Standard_Real fact_normalization = 100.;
 190   Standard_Integer i;
 191   // Pour chaque coin contraint, on controle les bounds adjascents.
 192   for(i = 0; i < nb; i++){
 193     if(stat[i].HasConstraint()){
 194       Standard_Integer ip = (i-1+nb)%nb;
 195       Standard_Real tolang = Min(bound[ip]->Tolang(),bound[i]->Tolang());
 196       Standard_Real an = stat[i].NorAng();
 197       Standard_Boolean twist = Standard_False;
 198       if(an >= 0.5*M_PI) { twist = Standard_True; an = M_PI-an; }
 199       if(an > tolang) stat[i].DoKill(0.);
 200       else{
 201         Standard_Real fact = 0.5*27./4;
 202         tolang *= (Min(mintg[ip],mintg[i])*fact*fact_normalization);
 203         gp_Vec tgp, dnorp, tgi, dnori, vbid;
 204         gp_Pnt pbid;
 205         Standard_Real fp,lp,fi,li;
 206         if(!rev[ip]) bound[ip]->Bounds(fp,lp);
 207         else bound[ip]->Bounds(lp,fp);
 208         bound[ip]->D1(lp,pbid,tgp);
 209         bound[ip]->D1Norm(lp,vbid,dnorp);
 210         if(!rev[i]) bound[i]->Bounds(fi,li);
 211         else bound[i]->Bounds(li,fi);
 212         bound[i]->D1(fi,pbid,tgi);
 213         bound[i]->D1Norm(fi,vbid,dnori);
 214         Standard_Real scal1 = tgp.Dot(dnori);
 215         Standard_Real scal2 = tgi.Dot(dnorp);
 216         if(!twist) scal2 *= -1.;
 217         scal1 = Abs(scal1+scal2);
 218         if(scal1 > tolang) {
 219           Standard_Real killfactor = tolang/scal1;
 220           stat[i].DoKill(killfactor);
 221 #ifdef OCCT_DEBUG
 222           cout<<"pb coons cnd coin : "<<i<<" fact = "<<killfactor<<endl;
 223 #endif
 224         }
 225       }
 226     }
 227   }
 228 }
 229 static void killcorners(const Standard_Integer     nb,
 230                         Handle(GeomFill_Boundary)* bound,
 231                         Standard_Boolean*          rev,
 232                         Standard_Boolean*          nrev,
 233                         GeomFill_CornerState*      stat,
 234                         Handle(GeomFill_TgtField)* tga)
 235 {
 236   Standard_Integer i;
 237   // Pour chaque  bound, on  controle l etat  des extremites  et on flingue
 238   // eventuellement le champ tangent et les derivees du bound.
 239   for(i = 0; i < nb; i++){
 240     Standard_Integer inext = (i+1)%nb;
 241     Standard_Boolean fnul, lnul;
 242     Standard_Real fscal, lscal;
 243     if(!nrev[i]){
 244       fnul = stat[i].IsToKill(fscal);
 245       lnul = stat[inext].IsToKill(lscal);
 246     }
 247     else{
 248       lnul = stat[i].IsToKill(lscal);
 249       fnul = stat[inext].IsToKill(fscal);
 250     }
 251     if(fnul || lnul){
 252 #ifdef OCCT_DEBUG
 253       parclock.Start();
 254 #endif
 255       bound[i]->Reparametrize(0.,1.,fnul,lnul,fscal,lscal,rev[i]);
 256 #ifdef OCCT_DEBUG
 257       parclock.Stop();
 258 #endif
 259       if(bound[i]->HasNormals() && tga[i]->IsScalable()) {
 260         Handle(Law_BSpline) bs = Law::ScaleCub(0.,1.,fnul,lnul,fscal,lscal);
 261         tga[i]->Scale(bs);
 262 #ifdef DRAW
 263         if(dodraw) Law_draw1dcurve(bs,gp_Vec2d(1.,0.),1.);
 264 #endif
 265       }
 266     }
 267   }
 268 }
 269
 270 //=======================================================================
 271 //class : GeomFill_ConstrainedFilling_Eval
 272 //purpose: The evaluator for curve approximation
 273 //=======================================================================
 274
 275 class GeomFill_ConstrainedFilling_Eval : public AdvApprox_EvaluatorFunction
 276 {
 277  public:
 278   GeomFill_ConstrainedFilling_Eval (GeomFill_ConstrainedFilling& theTool)
 279     : curfil(theTool) {}
 280
 281   virtual void Evaluate (Standard_Integer *Dimension,
 282                          Standard_Real     StartEnd[2],
 283                          Standard_Real    *Parameter,
 284                          Standard_Integer *DerivativeRequest,
 285                          Standard_Real    *Result, // [Dimension]
 286                          Standard_Integer *ErrorCode);
 287
 288  private:
 289   GeomFill_ConstrainedFilling& curfil;
 290 };
 291
 292 void GeomFill_ConstrainedFilling_Eval::Evaluate (Standard_Integer *,/*Dimension*/
 293                                                  Standard_Real     /*StartEnd*/[2],
 294                                                  Standard_Real    *Parameter,
 295                                                  Standard_Integer *DerivativeRequest,
 296                                                  Standard_Real    *Result,// [Dimension]
 297                                                  Standard_Integer *ErrorCode)
 298 {
 299   *ErrorCode = curfil.Eval(*Parameter, *DerivativeRequest, Result[0]);
 300 }
 301
 302 //=======================================================================
 303 //function : GeomFill_ConstrainedFilling
 304 //purpose  :
 305 //=======================================================================
 306
 307 GeomFill_ConstrainedFilling::GeomFill_ConstrainedFilling
 308 (const Standard_Integer MaxDeg,
 309  const Standard_Integer MaxSeg) :
 310  degmax(MaxDeg),segmax(MaxSeg),appdone(Standard_False)
 311 {
 312   dom[0] = dom[1] = dom[2] = dom[3] = 1.;
 313 }
 314
 315
 316 //=======================================================================
 317 //function : Init
 318 //purpose  :
 319 //=======================================================================
 320
 321 void GeomFill_ConstrainedFilling::Init(const Handle(GeomFill_Boundary)& B1,
 322                                        const Handle(GeomFill_Boundary)& B2,
 323                                        const Handle(GeomFill_Boundary)& B3,
 324                                        const Standard_Boolean           NoCheck)
 325 {
 326 #ifdef OCCT_DEBUG
 327   totclock.Reset();
 328   parclock.Reset();
 329   appclock.Reset();
 330   cstclock.Reset();
 331   totclock.Start();
 332 #endif
 333   Standard_Boolean rev[3];
 334   rev[0] = rev[1] = rev[2] = Standard_False;
 335   Handle(GeomFill_Boundary) bound[3];
 336   bound[0] = B1; bound[1] = B2; bound[2] = B3;
 337   Standard_Integer i;
 338   sortbounds(3,bound,rev,stcor);
 339
 340   // on reoriente.
 341   rev[2] = !rev[2];
 342
 343   // on reparamettre tout le monde entre 0. et 1.
 344 #ifdef OCCT_DEBUG
 345   parclock.Start();
 346 #endif
 347   for (i = 0; i <= 2; i++){
 348     bound[i]->Reparametrize(0.,1.,0,0,1.,1.,rev[i]);
 349   }
 350 #ifdef OCCT_DEBUG
 351   parclock.Stop();
 352 #endif
 353
 354   // On cree le carreau algorithmique (u,(1-u)) et les champs tangents
 355   // 1er jus.
 356   // On cree donc le bord manquant.
 357   gp_Pnt p1 = bound[1]->Value(1.);
 358   gp_Pnt p2 = bound[2]->Value(1.);
 359   gp_Pnt ppp(0.5*(p1.XYZ()+p2.XYZ()));
 360   Standard_Real t3 = Max(bound[1]->Tol3d(),bound[2]->Tol3d());
 361   Handle(GeomFill_DegeneratedBound)
 362     DB = new GeomFill_DegeneratedBound(ppp,0.,1.,t3,10.);
 363
 364   ptch = new GeomFill_CoonsAlgPatch(bound[0],bound[1],DB,bound[2]);
 365
 366   Handle(GeomFill_TgtField) ttgalg[3];
 367   if(bound[0]->HasNormals())
 368     ttgalg[0] = tgalg[0] = new GeomFill_TgtOnCoons(ptch,0);
 369   if(bound[1]->HasNormals())
 370     ttgalg[1] = tgalg[1] = new GeomFill_TgtOnCoons(ptch,1);
 371   if(bound[2]->HasNormals())
 372     ttgalg[2] = tgalg[3] = new GeomFill_TgtOnCoons(ptch,3);
 373
 374   for (i = 0; i <= 3; i++){
 375     mig[i] = 1.;
 376     if(!tgalg[i].IsNull()) MinTgte(i);
 377   }
 378
 379   if(!NoCheck){
 380     // On  verifie enfin les conditions  de compatibilites sur les derivees
 381     // aux coins maintenant qu on a quelque chose a quoi les comparer.
 382     Standard_Boolean nrev[3];
 383     nrev[0] = nrev[1] = 0;
 384     nrev[2] = 1;
 385     mig[2] = mig[3];
 386     coonscnd(3,bound,nrev,stcor,ttgalg,mig);
 387     killcorners(3,bound,rev,nrev,stcor,ttgalg);
 388   }
 389   // on remet les coins en place (on duplique la pointe).
 390   stcor[3] = stcor[2];
 391
 392   for (i = 0; i <= 3; i++){
 393     mig[i] = 1.;
 394     if(!tgalg[i].IsNull()) {
 395       if(!CheckTgte(i)) {
 396         Handle(Law_Function) fu1,fu2;
 397         ptch->Func(fu1,fu2);
 398         fu1 = Law::MixBnd(*((Handle(Law_Linear)*) &fu1));
 399         fu2 = Law::MixBnd(*((Handle(Law_Linear)*) &fu2));
 400         ptch->Func(fu1,fu2);
 401         break;
 402       }
 403     }
 404   }
 405
 406   Build();
 407 }
 408
 409
 410 //=======================================================================
 411 //function : Init
 412 //purpose  :
 413 //=======================================================================
 414
 415 void GeomFill_ConstrainedFilling::Init(const Handle(GeomFill_Boundary)& B1,
 416                                        const Handle(GeomFill_Boundary)& B2,
 417                                        const Handle(GeomFill_Boundary)& B3,
 418                                        const Handle(GeomFill_Boundary)& B4,
 419                                        const Standard_Boolean           NoCheck)
 420 {
 421 #ifdef OCCT_DEBUG
 422   totclock.Reset();
 423   parclock.Reset();
 424   appclock.Reset();
 425   cstclock.Reset();
 426   totclock.Start();
 427 #endif
 428   Standard_Boolean rev[4];
 429   rev[0] = rev[1] = rev[2] = rev[3] = Standard_False;
 430   Handle(GeomFill_Boundary) bound[4];
 431   bound[0] = B1; bound[1] = B2; bound[2] = B3; bound[3] = B4;
 432   Standard_Integer i;
 433   sortbounds(4,bound,rev,stcor);
 434
 435   // on reoriente.
 436   rev[2] = !rev[2];
 437   rev[3] = !rev[3];
 438
 439   // on reparamettre tout le monde entre 0. et 1.
 440 #ifdef OCCT_DEBUG
 441   parclock.Start();
 442 #endif
 443   for (i = 0; i <= 3; i++){
 444     bound[i]->Reparametrize(0.,1.,0,0,1.,1.,rev[i]);
 445   }
 446 #ifdef OCCT_DEBUG
 447   parclock.Stop();
 448 #endif
 449
 450   // On cree le carreau algorithmique (u,(1-u)) et les champs tangents
 451   // 1er jus.
 452   ptch = new GeomFill_CoonsAlgPatch(bound[0],bound[1],bound[2],bound[3]);
 453   for (i = 0; i <= 3; i++){
 454     if(bound[i]->HasNormals()) tgalg[i] = new GeomFill_TgtOnCoons(ptch,i);
 455   }
 456   // on calcule le min de chacun des champs tangents pour l evaluation
 457   // des tolerances.
 458   for (i = 0; i <= 3; i++){
 459     mig[i] = 1.;
 460     if(!tgalg[i].IsNull()) MinTgte(i);
 461   }
 462
 463   if(!NoCheck){
 464     // On  verifie enfin les conditions  de compatibilites sur les derivees
 465     // aux coins maintenant qu on a quelque chose a quoi les comparer.
 466     Standard_Boolean nrev[4];
 467     nrev[0] = nrev[1] = 0;
 468     nrev[2] = nrev[3] = 1;
 469     coonscnd(4,bound,nrev,stcor,tgalg,mig);
 470     killcorners(4,bound,rev,nrev,stcor,tgalg);
 471   }
 472   // On verifie les champs tangents ne changent pas de direction.
 473   for (i = 0; i <= 3; i++){
 474     mig[i] = 1.;
 475     if(!tgalg[i].IsNull()) {
 476       if(!CheckTgte(i)) {
 477         Handle(Law_Function) fu1,fu2;
 478         ptch->Func(fu1,fu2);
 479         Handle(Law_Function) ffu1 = Law::MixBnd(*((Handle(Law_Linear)*) &fu1));
 480         Handle(Law_Function) ffu2 = Law::MixBnd(*((Handle(Law_Linear)*) &fu2));
 481         ptch->SetFunc(ffu1,ffu2);
 482         break;
 483       }
 484     }
 485   }
 486
 487   Build();
 488 }
 489
 490
 491 //=======================================================================
 492 //function : SetDomain
 493 //purpose  :
 494 //=======================================================================
 495
 496 void GeomFill_ConstrainedFilling::SetDomain
 497 (const Standard_Real l, const Handle(GeomFill_BoundWithSurf)& B)
 498 {
 499   if(B == ptch->Bound(0)) dom[0] = Min(1.,Abs(l));
 500   else if(B == ptch->Bound(1)) dom[1] = Min(1.,Abs(l));
 501   else if(B == ptch->Bound(2)) dom[2] = Min(1.,Abs(l));
 502   else if(B == ptch->Bound(3)) dom[3] = Min(1.,Abs(l));
 503 }
 504
 505
 506 //=======================================================================
 507 //function : ReBuild
 508 //purpose  :
 509 //=======================================================================
 510
 511 void GeomFill_ConstrainedFilling::ReBuild()
 512 {
 513   if(!appdone) Standard_Failure::Raise
 514     ("GeomFill_ConstrainedFilling::ReBuild Approx non faite");
 515   MatchKnots();
 516   PerformS0();
 517   PerformS1();
 518   PerformSurface();
 519 }
 520
 521
 522 //=======================================================================
 523 //function : Boundary
 524 //purpose  :
 525 //=======================================================================
 526
 527 Handle(GeomFill_Boundary) GeomFill_ConstrainedFilling::Boundary
 528 (const Standard_Integer I) const
 529 {
 530   return ptch->Bound(I);
 531 }
 532
 533
 534 //=======================================================================
 535 //function : Surface
 536 //purpose  :
 537 //=======================================================================
 538
 539 Handle(Geom_BSplineSurface) GeomFill_ConstrainedFilling::Surface() const
 540 {
 541   return surf;
 542 }
 543
 544
 545 //=======================================================================
 546 //function : Build
 547 //purpose  :
 548 //=======================================================================
 549
 550 void GeomFill_ConstrainedFilling::Build()
 551 {
 552   for (Standard_Integer count = 0; count < 2; count++){
 553     ibound[0] = count; ibound[1] = count+2;
 554     ctr[0] = ctr[1] = nbd3 = 0;
 555     Standard_Integer ii ;
 556     for ( ii = 0; ii < 2; ii++){
 557       if (ptch->Bound(ibound[ii])->HasNormals()) {
 558         ctr[ii] = 2;
 559       }
 560       else if (!ptch->Bound(ibound[ii])->IsDegenerated()){
 561         ctr[ii] = 1;
 562       }
 563       nbd3 += ctr[ii];
 564     }
 565 #ifdef OCCT_DEBUG
 566     appclock.Start();
 567 #endif
 568     if(nbd3) PerformApprox();
 569 #ifdef OCCT_DEBUG
 570     appclock.Stop();
 571 #endif
 572   }
 573   appdone = Standard_True;
 574 #ifdef OCCT_DEBUG
 575   cstclock.Start();
 576 #endif
 577   MatchKnots();
 578   PerformS0();
 579   PerformS1();
 580   PerformSurface();
 581 #ifdef OCCT_DEBUG
 582   cstclock.Stop();
 583   totclock.Stop();
 584   Standard_Real tottime, apptime, partime, csttime;
 585   totclock.Show(tottime);
 586   parclock.Show(partime);
 587   appclock.Show(apptime);
 588   cstclock.Show(csttime);
 589   cout<<"temp total : "<<tottime<<" secondes"<<endl;
 590   cout<<endl;
 591   cout<<"dont"<<endl;
 592   cout<<endl;
 593   cout<<"reparametrage         : "<<partime<<" secondes"<<endl;
 594   cout<<"approximation         : "<<apptime<<" secondes"<<endl;
 595   cout<<"construction formelle : "<<csttime<<" secondes"<<endl;
 596   cout<<endl;
 597 #endif
 598 }
 599
 600
 601 //=======================================================================
 602 //function : PerformApprox
 603 //purpose  :
 604 //=======================================================================
 605
 606 void GeomFill_ConstrainedFilling::PerformApprox()
 607 {
 608   Standard_Integer ii ;
 609   Handle(TColStd_HArray1OfReal) tol3d, tol2d, tol1d;
 610   if(nbd3) tol3d = new TColStd_HArray1OfReal(1,nbd3);
 611   Standard_Integer i3d = 0;
 612   for( ii = 0; ii <= 1; ii++){
 613     if (ctr[ii]) {tol3d->SetValue((++i3d),ptch->Bound(ibound[ii])->Tol3d());}
 614     if(ctr[ii] == 2){
 615       tol3d->SetValue(++i3d,0.5* mig[ibound[ii]] * ptch->Bound(ibound[ii])->Tolang());
 616     }
 617   }
 618   Standard_Real f,l;
 619   ptch->Bound(ibound[0])->Bounds(f,l);
 620
 621   GeomFill_ConstrainedFilling_Eval ev (*this);
 622   AdvApprox_ApproxAFunction  app(0,
 623                                  0,
 624                                  nbd3,
 625                                  tol1d,
 626                                  tol2d,
 627                                  tol3d,
 628                                  f,
 629                                  l,
 630                                  GeomAbs_C1,
 631                                  degmax,
 632                                  segmax,
 633                                  ev);
 634
 635   if (app.IsDone() || app.HasResult()){
 636     Standard_Integer imk = Min(ibound[0],ibound[1]);
 637     Standard_Integer nbpol = app.NbPoles();
 638     degree[imk] = app.Degree();
 639     mults[imk] = app.Multiplicities();
 640     knots[imk] = app.Knots();
 641     i3d = 0;
 642     for(ii = 0; ii <= 1; ii++){
 643       curvpol[ibound[ii]] = new TColgp_HArray1OfPnt(1,nbpol);
 644       TColgp_Array1OfPnt& cp = curvpol[ibound[ii]]->ChangeArray1();
 645       if (ctr[ii]){
 646         app.Poles(++i3d,cp);
 647       }
 648       else{
 649         gp_Pnt ppp =  ptch->Bound(ibound[ii])->Value(0.5*(f+l));
 650         for(Standard_Integer ij = 1; ij <= nbpol; ij++){
 651           cp(ij) = ppp;
 652         }
 653       }
 654       if(ctr[ii] == 2){
 655         tgtepol[ibound[ii]] = new TColgp_HArray1OfPnt(1,nbpol);
 656         app.Poles(++i3d,tgtepol[ibound[ii]]->ChangeArray1());
 657       }
 658     }
 659   }
 660 }
 661
 662
 663 //=======================================================================
 664 //function : MatchKnots
 665 //purpose  :
 666 //=======================================================================
 667
 668 void GeomFill_ConstrainedFilling::MatchKnots()
 669 {
 670   // on n insere rien si les domaines valent 1.
 671   Standard_Integer i, j, l;
 672   Standard_Integer ind[4];
 673   nm[0] = mults[0]; nm[1] = mults[1];
 674   nk[0] = knots[0]; nk[1] = knots[1];
 675   ind[0] = nk[1]->Length(); ind[2] = 1;
 676   ind[1] = 1; ind[3] = nk[0]->Length();
 677   ncpol[0] = curvpol[0]; ncpol[1] = curvpol[1];
 678   ncpol[2] = curvpol[2]; ncpol[3] = curvpol[3];
 679   ntpol[0] = tgtepol[0]; ntpol[1] = tgtepol[1];
 680   ntpol[2] = tgtepol[2]; ntpol[3] = tgtepol[3];
 681   Standard_Real kadd[2];
 682   Standard_Integer madd[2];
 683   Standard_Real tolk = 1./Max(10,2*knots[1]->Array1().Length());
 684   Standard_Integer nbadd = inqadd(dom[0],dom[2],kadd,madd,degree[1],tolk);
 685   if(nbadd){
 686     TColStd_Array1OfReal addk(kadd[0],1,nbadd);
 687     TColStd_Array1OfInteger addm(madd[0],1,nbadd);
 688     Standard_Integer nbnp, nbnk;
 689     if(BSplCLib::PrepareInsertKnots(degree[1],0,
 690                                     knots[1]->Array1(),
 691                                     mults[1]->Array1(),
 692                                     addk,addm,nbnp,nbnk,tolk,0)){
 693       nm[1] = new TColStd_HArray1OfInteger(1,nbnk);
 694       nk[1] = new TColStd_HArray1OfReal(1,nbnk);
 695       ncpol[1] = new TColgp_HArray1OfPnt(1,nbnp);
 696       ncpol[3] = new TColgp_HArray1OfPnt(1,nbnp);
 697       BSplCLib::InsertKnots(degree[1],0,
 698                             curvpol[1]->Array1(),PLib::NoWeights(),
 699                             knots[1]->Array1(),mults[1]->Array1(),
 700                             addk,addm,
 701                             ncpol[1]->ChangeArray1(),PLib::NoWeights(),
 702                             nk[1]->ChangeArray1(),nm[1]->ChangeArray1(),
 703                             tolk,0);
 704
 705       BSplCLib::InsertKnots(degree[1],0,
 706                             curvpol[3]->Array1(),PLib::NoWeights(),
 707                             knots[1]->Array1(),mults[1]->Array1(),
 708                             addk,addm,
 709                             ncpol[3]->ChangeArray1(),PLib::NoWeights(),
 710                             nk[1]->ChangeArray1(),nm[1]->ChangeArray1(),
 711                             tolk,0);
 712       if(!tgtepol[1].IsNull()){
 713         ntpol[1] = new TColgp_HArray1OfPnt(1,nbnp);
 714         BSplCLib::InsertKnots(degree[1],0,
 715                               tgtepol[1]->Array1(),PLib::NoWeights(),
 716                               knots[1]->Array1(),mults[1]->Array1(),
 717                               addk,addm,
 718                               ntpol[1]->ChangeArray1(),PLib::NoWeights(),
 719                               nk[1]->ChangeArray1(),nm[1]->ChangeArray1(),
 720                               tolk,0);
 721       }
 722       if(!tgtepol[3].IsNull()){
 723         ntpol[3] = new TColgp_HArray1OfPnt(1,nbnp);
 724         BSplCLib::InsertKnots(degree[1],0,
 725                               tgtepol[3]->Array1(),PLib::NoWeights(),
 726                               knots[1]->Array1(),mults[1]->Array1(),
 727                               addk,addm,
 728                               ntpol[3]->ChangeArray1(),PLib::NoWeights(),
 729                               nk[1]->ChangeArray1(),nm[1]->ChangeArray1(),
 730                               tolk,0);
 731       }
 732     }
 733     if(dom[0] != 1.) {
 734       for(i = 2; i <= nbnk; i++){
 735         if(Abs(dom[0]-nm[1]->Value(i)) < tolk){
 736           ind[0] = i;
 737           break;
 738         }
 739       }
 740     }
 741     if(dom[2] != 1.) {
 742       for(i = 1; i < nbnk; i++){
 743         if(Abs(1.-dom[2]-nm[1]->Value(i)) < tolk){
 744           ind[2] = i;
 745           break;
 746         }
 747       }
 748     }
 749   }
 750   tolk = 1./Max(10.,2.*knots[0]->Array1().Length());
 751   nbadd = inqadd(dom[1],dom[3],kadd,madd,degree[0],tolk);
 752   if(nbadd){
 753     TColStd_Array1OfReal addk(kadd[0],1,nbadd);
 754     TColStd_Array1OfInteger addm(madd[0],1,nbadd);
 755     Standard_Integer nbnp, nbnk;
 756     if(BSplCLib::PrepareInsertKnots(degree[0],0,
 757                                     knots[0]->Array1(),
 758                                     mults[0]->Array1(),
 759                                     addk,addm,nbnp,nbnk,tolk,0)){
 760       nm[0] = new TColStd_HArray1OfInteger(1,nbnk);
 761       nk[0] = new TColStd_HArray1OfReal(1,nbnk);
 762       ncpol[0] = new TColgp_HArray1OfPnt(1,nbnp);
 763       ncpol[2] = new TColgp_HArray1OfPnt(1,nbnp);
 764       BSplCLib::InsertKnots(degree[0],0,
 765                             curvpol[0]->Array1(),PLib::NoWeights(),
 766                             knots[0]->Array1(),mults[0]->Array1(),
 767                             addk,addm,
 768                             ncpol[0]->ChangeArray1(),PLib::NoWeights(),
 769                             nk[0]->ChangeArray1(),nm[0]->ChangeArray1(),
 770                             tolk,0);
 771
 772       BSplCLib::InsertKnots(degree[0],0,
 773                             curvpol[2]->Array1(),PLib::NoWeights(),
 774                             knots[0]->Array1(),mults[0]->Array1(),
 775                             addk,addm,
 776                             ncpol[2]->ChangeArray1(),PLib::NoWeights(),
 777                             nk[0]->ChangeArray1(),nm[0]->ChangeArray1(),
 778                             tolk,0);
 779       if(!tgtepol[0].IsNull()){
 780         ntpol[0] = new TColgp_HArray1OfPnt(1,nbnp);
 781         BSplCLib::InsertKnots(degree[0],0,
 782                               tgtepol[0]->Array1(),PLib::NoWeights(),
 783                               knots[0]->Array1(),mults[0]->Array1(),
 784                               addk,addm,
 785                               ntpol[0]->ChangeArray1(),PLib::NoWeights(),
 786                               nk[0]->ChangeArray1(),nm[0]->ChangeArray1(),
 787                               tolk,0);
 788       }
 789       if(!tgtepol[2].IsNull()){
 790         ntpol[2] = new TColgp_HArray1OfPnt(1,nbnp);
 791         BSplCLib::InsertKnots(degree[0],0,
 792                               tgtepol[2]->Array1(),PLib::NoWeights(),
 793                               knots[0]->Array1(),mults[0]->Array1(),
 794                               addk,addm,
 795                               ntpol[2]->ChangeArray1(),PLib::NoWeights(),
 796                               nk[0]->ChangeArray1(),nm[0]->ChangeArray1(),
 797                               tolk,0);
 798       }
 799     }
 800     if(dom[1] != 1.) {
 801       for(i = 2; i <= nbnk; i++){
 802         if(Abs(dom[1]-nm[0]->Value(i)) < tolk){
 803           ind[1] = i;
 804           break;
 805         }
 806       }
 807     }
 808     if(dom[3] != 1.) {
 809       for(i = 1; i < nbnk; i++){
 810         if(Abs(1.-dom[3]-nm[0]->Value(i)) < tolk){
 811           ind[3] = i;
 812           break;
 813         }
 814       }
 815     }
 816   }
 817   Handle(Law_Linear) fu = mklin(ptch->Func(0));
 818   ab[0] = Law::MixBnd(degree[1],nk[1]->Array1(),nm[1]->Array1(),fu);
 819   fu = mklin(ptch->Func(1));
 820   ab[1] = Law::MixBnd(degree[0],nk[0]->Array1(),nm[0]->Array1(),fu);
 821
 822   for(i = 0; i<2; i++){
 823     l = ab[i]->Length();
 824     ab[i+2] = new TColStd_HArray1OfReal(1,l);
 825     for(j = 1; j <= l; j++){
 826       ab[i+2]->SetValue(j,1.-ab[i]->Value(j));
 827     }
 828   }
 829   pq[0] = Law::MixTgt(degree[1],nk[1]->Array1(),nm[1]->Array1(),1,ind[0]);
 830   pq[2] = Law::MixTgt(degree[1],nk[1]->Array1(),nm[1]->Array1(),0,ind[2]);
 831
 832   pq[1] = Law::MixTgt(degree[0],nk[0]->Array1(),nm[0]->Array1(),0,ind[1]);
 833   pq[3] = Law::MixTgt(degree[0],nk[0]->Array1(),nm[0]->Array1(),1,ind[3]);
 834
 835 #ifdef DRAW
 836   if(dodraw){
 837     gp_Vec2d tra(0.,0.);
 838     Standard_Real scal = 1.;
 839     Law_draw1dcurve(ab[0]->Array1(),nk[1]->Array1(),nm[1]->Array1(),degree[1],tra,scal);
 840     tra.SetCoord(1.,0.);
 841     Law_draw1dcurve(ab[1]->Array1(),nk[0]->Array1(),nm[0]->Array1(),degree[0],tra,scal);
 842     tra.SetCoord(0.,1.);
 843     Law_draw1dcurve(ab[2]->Array1(),nk[1]->Array1(),nm[1]->Array1(),degree[1],tra,scal);
 844     tra.SetCoord(1.,1.);
 845     Law_draw1dcurve(ab[3]->Array1(),nk[0]->Array1(),nm[0]->Array1(),degree[0],tra,scal);
 846     tra.SetCoord(0.,0.);
 847     Law_draw1dcurve(pq[0]->Array1(),nk[1]->Array1(),nm[1]->Array1(),degree[1],tra,scal);
 848     tra.SetCoord(0.,1.);
 849     Law_draw1dcurve(pq[2]->Array1(),nk[1]->Array1(),nm[1]->Array1(),degree[1],tra,scal);
 850     tra.SetCoord(1.,0.);
 851     Law_draw1dcurve(pq[1]->Array1(),nk[0]->Array1(),nm[0]->Array1(),degree[0],tra,scal);
 852     tra.SetCoord(1.,1.);
 853     Law_draw1dcurve(pq[3]->Array1(),nk[0]->Array1(),nm[0]->Array1(),degree[0],tra,scal);
 854   }
 855 #endif
 856 }
 857
 858
 859 //=======================================================================
 860 //function : PerformS0
 861 //purpose  :
 862 //=======================================================================
 863
 864 void GeomFill_ConstrainedFilling::PerformS0()
 865 {
 866   // On construit les poles de S0 par combinaison des poles des bords,
 867   // des poles des fonctions ab, des points c selon la formule :
 868   // S0(i,j) = ab[0](j)*ncpol[0](i) + ab[1](i)*ncpol[1](j)
 869   //         + ab[2](j)*ncpol[2](i) + ab[3](i)*ncpol[3](j)
 870   //         - ab[3](i)*ab[0](j)*c[0] - ab[0](j)*ab[1](i)*c[1]
 871   //         - ab[1](i)*ab[2](j)*c[2] - ab[2](j)*ab[3](i)*c[3]
 872
 873   Standard_Integer i, j;
 874   Standard_Integer ni = ncpol[0]->Length();
 875   Standard_Integer nj = ncpol[1]->Length();
 876   S0 = new TColgp_HArray2OfPnt(1,ni,1,nj);
 877   TColgp_Array2OfPnt& ss0 = S0->ChangeArray2();
 878   const gp_XYZ& c0 = ptch->Corner(0).Coord();
 879   const gp_XYZ& c1 = ptch->Corner(1).Coord();
 880   const gp_XYZ& c2 = ptch->Corner(2).Coord();
 881   const gp_XYZ& c3 = ptch->Corner(3).Coord();
 882   for (i = 1; i <= ni; i++){
 883     Standard_Real ab1 = ab[1]->Value(i);
 884     Standard_Real ab3 = ab[3]->Value(i);
 885     const gp_XYZ& b0 = ncpol[0]->Value(i).Coord();
 886     const gp_XYZ& b2 = ncpol[2]->Value(i).Coord();
 887     for (j = 1; j <= nj; j++){
 888       Standard_Real ab0 = ab[0]->Value(j);
 889       Standard_Real ab2 = ab[2]->Value(j);
 890       const gp_XYZ& b1 = ncpol[1]->Value(j).Coord();
 891       const gp_XYZ& b3 = ncpol[3]->Value(j).Coord();
 892       gp_XYZ polij = b0.Multiplied(ab0);
 893       gp_XYZ temp = b1.Multiplied(ab1);
 894       polij.Add(temp);
 895       temp = b2.Multiplied(ab2);
 896       polij.Add(temp);
 897       temp = b3.Multiplied(ab3);
 898       polij.Add(temp);
 899       temp = c0.Multiplied(-ab3*ab0);
 900       polij.Add(temp);
 901       temp = c1.Multiplied(-ab0*ab1);
 902       polij.Add(temp);
 903       temp = c2.Multiplied(-ab1*ab2);
 904       polij.Add(temp);
 905       temp = c3.Multiplied(-ab2*ab3);
 906       polij.Add(temp);
 907       ss0(i,j).SetXYZ(polij);
 908     }
 909   }
 910 }
 911
 912
 913 //=======================================================================
 914 //function : PerformS1
 915 //purpose  :
 916 //=======================================================================
 917
 918 void GeomFill_ConstrainedFilling::PerformS1()
 919 {
 920   // on construit en temporaire les poles des champs tangents
 921   // definis par :
 922   // tgte[ibound](u) - d/dv (S0(u,vbound)) pour ibound = 0 ou 2
 923   // tgte[ibound](v) - d/du (S0(ubound,v)) pour ibound = 1 ou 3
 924   // sur les bords ou tgte est defini.
 925   gp_XYZ* nt[4];
 926   const TColgp_Array2OfPnt& ss0 = S0->Array2();
 927   Standard_Integer l, i, j, k;
 928   Standard_Integer ni = ss0.ColLength();
 929   Standard_Integer nj = ss0.RowLength();
 930   for(i = 0; i <= 3; i++){
 931     if(ntpol[i].IsNull()) nt[i] = 0;
 932     else {
 933       Standard_Real z=0;
 934       Standard_Integer nbp = ntpol[i]->Length();
 935       Standard_Integer i1=0,i2=0,j1=0,j2=0;
 936       Standard_Boolean inci=0;
 937       nt[i] = new gp_XYZ[nbp];
 938       switch(i){
 939       case 0 :
 940         z = - degree[1]/(nk[1]->Value(2) - nk[1]->Value(1));
 941         inci = Standard_True;
 942         i1 = 1; i2 = 1; j1 = 1; j2 = 2;
 943         break;
 944       case 1 :
 945         l = nk[0]->Length();
 946         z = - degree[0]/(nk[0]->Value(l) - nk[0]->Value(l-1));
 947         inci = Standard_False;
 948         i1 = ni-1; i2 = ni; j1 = 1; j2 = 1;
 949         break;
 950       case 2 :
 951         l = nk[1]->Length();
 952         z = - degree[1]/(nk[1]->Value(l) - nk[1]->Value(l-1));
 953         inci = Standard_True;
 954         i1 = 1; i2 = 1; j1 = nj-1; j2 = nj;
 955         break;
 956       case 3 :
 957         z = - degree[0]/(nk[0]->Value(2) - nk[0]->Value(1));
 958         inci = Standard_False;
 959         i1 = 1; i2 = 2; j1 = 1; j2 = 1;
 960         break;
 961       }
 962       for(k = 0; k < nbp; k++){
 963         nt[i][k] = S0->Value(i1,j1).XYZ();
 964         nt[i][k].Multiply(-1.);
 965         nt[i][k].Add(S0->Value(i2,j2).XYZ());
 966         nt[i][k].Multiply(z);
 967         nt[i][k].Add(ntpol[i]->Value(k+1).XYZ());
 968         if(inci) { i1++; i2++; }
 969         else { j1++; j2++; }
 970       }
 971     }
 972   }
 973   // on calcul les termes correctifs pour le melange.
 974   Standard_Real coef0 = degree[0]/(nk[0]->Value(2) - nk[0]->Value(1));
 975   Standard_Real coef1 = degree[1]/(nk[1]->Value(2) - nk[1]->Value(1));
 976   gp_XYZ vtemp, vtemp0, vtemp1;
 977   if(nt[0] && nt[3]){
 978     vtemp0 = nt[0][0].Multiplied(-1.);
 979     vtemp0.Add(nt[0][1]);
 980     vtemp0.Multiply(coef0);
 981     vtemp1 = nt[3][0].Multiplied(-1.);
 982     vtemp1.Add(nt[3][1]);
 983     vtemp1.Multiply(coef1);
 984     vtemp = vtemp0.Added(vtemp1);
 985     vtemp.Multiply(0.5);
 986     v[0].SetXYZ(vtemp);
 987   }
 988
 989   Standard_Integer ln0 = nk[0]->Length(), lp0 = ncpol[0]->Length();
 990   coef0 = degree[0]/(nk[0]->Value(ln0) - nk[0]->Value(ln0 - 1));
 991   coef1 = degree[1]/(nk[1]->Value(2) - nk[1]->Value(1));
 992   if(nt[0] && nt[1]){
 993     vtemp0 = nt[0][lp0 - 2].Multiplied(-1.);
 994     vtemp0.Add(nt[0][lp0 - 1]);
 995     vtemp0.Multiply(coef0);
 996     vtemp1 = nt[1][0].Multiplied(-1.);
 997     vtemp1.Add(nt[1][1]);
 998     vtemp1.Multiply(coef1);
 999     vtemp = vtemp0.Added(vtemp1);
1000     vtemp.Multiply(0.5);
1001     v[1].SetXYZ(vtemp);
1002   }
1003   ln0 = nk[0]->Length(); lp0 = ncpol[0]->Length();
1004   Standard_Integer ln1 = nk[1]->Length(), lp1 = ncpol[1]->Length();
1005   coef0 = degree[0]/(nk[0]->Value(ln0) - nk[0]->Value(ln0 - 1));
1006   coef1 = degree[1]/(nk[1]->Value(ln1) - nk[1]->Value(ln1 - 1));
1007   if(nt[1] && nt[2]){
1008     vtemp0 = nt[2][lp0 - 2].Multiplied(-1.);
1009     vtemp0.Add(nt[2][lp0 - 1]);
1010     vtemp0.Multiply(coef0);
1011     vtemp1 = nt[1][lp1 - 2].Multiplied(-1.);
1012     vtemp1.Add(nt[1][lp1 - 1]);
1013     vtemp1.Multiply(coef1);
1014     vtemp = vtemp0.Added(vtemp1);
1015     vtemp.Multiply(0.5);
1016     v[2].SetXYZ(vtemp);
1017   }
1018   ln1 = nk[1]->Length(); lp1 = ncpol[1]->Length();
1019   coef0 = degree[0]/(nk[0]->Value(2) - nk[0]->Value(1));
1020   coef1 = degree[1]/(nk[1]->Value(ln1) - nk[1]->Value(ln1 - 1));
1021   if(nt[2] && nt[3]){
1022     vtemp0 = nt[2][0].Multiplied(-1.);
1023     vtemp0.Add(nt[2][1]);
1024     vtemp0.Multiply(coef0);
1025     vtemp1 = nt[3][lp1 - 2].Multiplied(-1.);
1026     vtemp1.Add(nt[3][lp1 - 1]);
1027     vtemp1.Multiply(coef1);
1028     vtemp = vtemp0.Added(vtemp1);
1029     vtemp.Multiply(0.5);
1030     v[3].SetXYZ(vtemp);
1031   }
1032
1033   // On construit les poles de S1 par combinaison des poles des
1034   // champs tangents, des poles des fonctions pq, des duv au coins
1035   // selon la formule :
1036   // S1(i,j) = pq[0](j)*ntpol[0](i) + pq[1](i)*ntpol[1](j)
1037   //         + pq[2](j)*ntpol[2](i) + pq[3](i)*ntpol[3](j)
1038   //         - pq[3](i)*pq[0](j)*v[0] - pq[0](j)*pq[1](i)*v[1]
1039   //         - pq[1](i)*pq[2](j)*v[2] - pq[2](j)*pq[3](i)*v[3]
1040   S1 = new TColgp_HArray2OfPnt(1,ni,1,nj);
1041   TColgp_Array2OfPnt& ss1 = S1->ChangeArray2();
1042   const gp_XYZ& v0 = v[0].XYZ();
1043   const gp_XYZ& v1 = v[1].XYZ();
1044   const gp_XYZ& v2 = v[2].XYZ();
1045   const gp_XYZ& v3 = v[3].XYZ();
1046
1047   for (i = 1; i <= ni; i++){
1048     Standard_Real pq1=0, pq3=0;
1049     if(nt[1]) pq1 = -pq[1]->Value(i);
1050     if(nt[3]) pq3 = pq[3]->Value(i);
1051     gp_XYZ t0, t2;
1052     if(nt[0]) t0 = nt[0][i-1];
1053     if(nt[2]) t2 = nt[2][i-1];
1054     for (j = 1; j <= nj; j++){
1055       Standard_Real pq0=0, pq2=0;
1056       if(nt[0]) pq0 = pq[0]->Value(j);
1057       if(nt[2]) pq2 = -pq[2]->Value(j);
1058       gp_XYZ t1, t3;
1059       if(nt[1]) t1 = nt[1][j-1];
1060       if(nt[3]) t3 = nt[3][j-1];
1061
1062       gp_XYZ tpolij(0.,0.,0.), temp;
1063       if(nt[0]) {
1064         temp = t0.Multiplied(pq0);
1065         tpolij.Add(temp);
1066       }
1067       if(nt[1]) {
1068         temp = t1.Multiplied(pq1);
1069         tpolij.Add(temp);
1070       }
1071       if(nt[2]){
1072         temp = t2.Multiplied(pq2);
1073         tpolij.Add(temp);
1074       }
1075       if(nt[3]){
1076         temp = t3.Multiplied(pq3);
1077         tpolij.Add(temp);
1078       }
1079       if(nt[3] && nt[0]){
1080         temp = v0.Multiplied(-pq3*pq0);
1081         tpolij.Add(temp);
1082       }
1083       if(nt[0] && nt[1]){
1084         temp = v1.Multiplied(-pq0*pq1);
1085         tpolij.Add(temp);
1086       }
1087       if(nt[1] && nt[2]){
1088         temp = v2.Multiplied(-pq1*pq2);
1089         tpolij.Add(temp);
1090       }
1091       if(nt[2] && nt[3]){
1092         temp = v3.Multiplied(-pq2*pq3);
1093         tpolij.Add(temp);
1094       }
1095       ss1(i,j).SetXYZ(tpolij);
1096     }
1097   }
1098
1099   // Un petit menage
1100   for(i = 0; i <= 3; i++){
1101     if(nt[i]){
1102       delete[] nt[i];
1103     }
1104   }
1105 }
1106
1107
1108 //=======================================================================
1109 //function : PerformSurface
1110 //purpose  :
1111 //=======================================================================
1112
1113 void GeomFill_ConstrainedFilling::PerformSurface()
1114 {
1115   Standard_Integer ni = S0->ColLength(), nj = S0->RowLength(),i,j;
1116   TColgp_Array2OfPnt temp(1,ni,1,nj);
1117   const TColgp_Array2OfPnt& t0 = S0->Array2();
1118   const TColgp_Array2OfPnt& t1 = S1->Array2();
1119   for(i = 1; i <= ni; i++){
1120     for(j = 1; j <= nj; j++){
1121       temp(i,j).SetXYZ(t0(i,j).XYZ().Added(t1(i,j).XYZ()));
1122     }
1123   }
1124   surf = new Geom_BSplineSurface(temp,
1125                                  nk[0]->Array1(),nk[1]->Array1(),
1126                                  nm[0]->Array1(),nm[1]->Array1(),
1127                                  degree[0],degree[1]);
1128 }
1129
1130 //=======================================================================
1131 //function : CheckTgte
1132 //purpose  :
1133 //=======================================================================
1134
1135 Standard_Boolean GeomFill_ConstrainedFilling::CheckTgte(const Standard_Integer I)
1136 {
1137   Handle(GeomFill_Boundary) bou = ptch->Bound(I);
1138   if(!bou->HasNormals()) return Standard_True;
1139   // On prend 13 points le long du bord et on verifie que le triedre
1140   // forme par la tangente a la courbe la normale et la tangente du
1141   // peigne ne change pas d orientation.
1142   Standard_Real ll = 1./12., pmix=0;
1143   for (Standard_Integer iu = 0; iu < 13; iu++){
1144     Standard_Real uu = iu * ll;
1145     gp_Pnt pbid;
1146     gp_Vec tgte;
1147     bou->D1(uu,pbid,tgte);
1148     gp_Vec norm = bou->Norm(uu);
1149     gp_Vec vfield = tgalg[I]->Value(uu);
1150     if(iu == 0) pmix = vfield.Dot(tgte.Crossed(norm));
1151     else {
1152       Standard_Real pmixcur = vfield.Dot(tgte.Crossed(norm));
1153       if(pmix*pmixcur < 0.) return Standard_False;
1154     }
1155   }
1156   return Standard_True;
1157 }
1158
1159 //=======================================================================
1160 //function : MinTgte
1161 //purpose  :
1162 //=======================================================================
1163
1164 void GeomFill_ConstrainedFilling::MinTgte(const Standard_Integer I)
1165 {
1166   if(!ptch->Bound(I)->HasNormals()) return;
1167   Standard_Real minmag = RealLast();
1168   Standard_Real ll = 0.02;
1169   for (Standard_Integer iu = 0; iu <= 30; iu++){
1170     Standard_Real uu = 0.2 + iu * ll;
1171     gp_Vec vv = tgalg[I]->Value(uu);
1172     Standard_Real temp = vv.SquareMagnitude();
1173     if(temp < minmag) minmag = temp;
1174   }
1175   mig[I] = sqrt(minmag);
1176 }
1177
1178 //=======================================================================
1179 //function : Eval
1180 //purpose  :
1181 //=======================================================================
1182
1183 Standard_Integer GeomFill_ConstrainedFilling::Eval(const Standard_Real W,
1184                                                    const Standard_Integer Ord,
1185                                                    Standard_Real& Result)const
1186 {
1187   Standard_Real* res = &Result;
1188   Standard_Integer jmp = (3 * ctr[0]);
1189   switch(Ord){
1190   case 0 :
1191     if(ctr[0]){
1192       ptch->Bound(ibound[0])->Value(W).Coord(res[0],res[1],res[2]);
1193     }
1194     if(ctr[0] == 2){
1195       tgalg[ibound[0]]->Value(W).Coord(res[3],res[4],res[5]);
1196     }
1197     if(ctr[1]){
1198       ptch->Bound(ibound[1])->Value(W).Coord(res[jmp],res[jmp+1],res[jmp+2]);
1199     }
1200     if(ctr[1] == 2){
1201       tgalg[ibound[1]]->Value(W).Coord(res[jmp+3],res[jmp+4],res[jmp+5]);
1202     }
1203     break;
1204   case 1 :
1205     gp_Pnt pt;
1206     gp_Vec vt;
1207     if(ctr[0]){
1208       ptch->Bound(ibound[0])->D1(W,pt,vt);
1209       vt.Coord(res[0],res[1],res[2]);
1210     }
1211     if(ctr[0] == 2){
1212       tgalg[ibound[0]]->D1(W).Coord(res[3],res[4],res[5]);
1213     }
1214     if(ctr[1]){
1215       ptch->Bound(ibound[1])->D1(W,pt,vt);
1216       vt.Coord(res[jmp],res[jmp+1],res[jmp+2]);
1217     }
1218     if(ctr[1] == 2){
1219       tgalg[ibound[1]]->D1(W).Coord(res[jmp+3],res[jmp+4],res[jmp+5]);
1220     }
1221     break;
1222   }
1223   return 0;
1224 }
1225
1226 //=======================================================================
1227 //function : CheckCoonsAlgPatch
1228 //purpose  :
1229 //=======================================================================
1230
1231 void GeomFill_ConstrainedFilling::CheckCoonsAlgPatch(const Standard_Integer I)
1232 {
1233   Standard_Integer nbp = 30;
1234   Standard_Real uu=0,duu=0,vv=0,dvv=0,ww=0,dww=0,u1,u2,v1,v2;
1235   surf->Bounds(u1,u2,v1,v2);
1236   Standard_Boolean enu = Standard_False;
1237   switch(I){
1238   case 0:
1239     uu = ww = u1;
1240     vv = v1;
1241     duu = dww = (u2 - u1)/nbp;
1242     dvv = 0.;
1243     break;
1244   case 1:
1245     vv = ww = v1;
1246     uu = u2;
1247     dvv = dww = (v2 - v1)/nbp;
1248     duu = 0.;
1249     enu = Standard_True;
1250     break;
1251   case 2:
1252     uu = ww = u1;
1253     vv = v2;
1254     duu = dww = (u2 - u1)/nbp;
1255     dvv = 0.;
1256     break;
1257   case 3:
1258     vv = ww = v1;
1259     uu = u1;
1260     dvv = dww = (v2 - v1)/nbp;
1261     duu = 0.;
1262     enu = Standard_True;
1263     break;
1264   }
1265   gp_Pnt pbound;
1266   gp_Vec vptch;
1267   Handle(GeomFill_Boundary) bou = ptch->Bound(I);
1268   for (Standard_Integer k = 0; k <= nbp; k++){
1269     pbound = bou->Value(ww);
1270     if(enu) vptch = ptch->D1U(uu,vv);
1271     else vptch = ptch->D1V(uu,vv);
1272 #ifdef DRAW
1273     gp_Pnt pp;
1274     Handle(Draw_Segment3D) seg;
1275     pp = pbound.Translated(vptch);
1276     seg = new Draw_Segment3D(pbound,pp,Draw_jaune);
1277     dout << seg;
1278 #endif
1279     uu += duu;
1280     vv += dvv;
1281     ww += dww;
1282   }
1283 }
1284
1285 //=======================================================================
1286 //function : CheckTgteField
1287 //purpose  :
1288 //=======================================================================
1289
1290 void GeomFill_ConstrainedFilling::CheckTgteField(const Standard_Integer I)
1291 {
1292   if(tgalg[I].IsNull()) return;
1293 #ifdef DRAW
1294   gp_Pnt p1,p2;
1295 #else
1296   gp_Pnt p1;
1297 #endif
1298   gp_Vec d1;
1299   Standard_Boolean caplisse = 0;
1300   Standard_Real maxang = 0.,pmix=0,pmixcur;
1301   Handle(GeomFill_Boundary) bou =  ptch->Bound(I);
1302   for (Standard_Integer iu = 0; iu <= 30; iu++){
1303     Standard_Real uu = iu/30.;
1304     bou->D1(uu,p1,d1);
1305     gp_Vec vtg = tgalg[I]->Value(uu);
1306     gp_Vec vnor = bou->Norm(uu);
1307     gp_Vec vcros = d1.Crossed(vnor);
1308     vcros.Normalize();
1309     if(iu == 0) pmix = vtg.Dot(vcros);
1310     else {
1311       pmixcur = vtg.Dot(vcros);
1312       if(pmix*pmixcur < 0.) caplisse = 1;
1313     }
1314 #ifdef DRAW
1315     Handle(Draw_Segment3D) seg;
1316     p2 = p1.Translated(vtg);
1317     seg = new Draw_Segment3D(p1,p2,Draw_blanc);
1318     dout << seg;
1319     p2 = p1.Translated(vnor);
1320     seg = new Draw_Segment3D(p1,p2,Draw_rouge);
1321     dout << seg;
1322     p2 = p1.Translated(vcros);
1323     seg = new Draw_Segment3D(p1,p2,Draw_jaune);
1324     dout << seg;
1325 #endif
1326     if(vnor.Magnitude() > 1.e-15 && vtg.Magnitude() > 1.e-15){
1327       Standard_Real alpha = Abs(M_PI/2.-Abs(vnor.Angle(vtg)));
1328       if(Abs(alpha) > maxang) maxang = Abs(alpha);
1329     }
1330   }
1331   cout<<"KAlgo angle max sur bord "<<I<<" : "<<maxang<<endl;
1332   if(caplisse) cout<<"sur bord "<<I<<" le champ tangent change de cote!"<<endl;
1333 }
1334
1335
1336 //=======================================================================
1337 //function : CheckApprox
1338 //purpose  :
1339 //=======================================================================
1340
1341 void GeomFill_ConstrainedFilling::CheckApprox(const Standard_Integer I)
1342 {
1343   Standard_Boolean donor = !tgalg[I].IsNull();
1344   Standard_Integer nbp = 30;
1345   Standard_Real maxang = 0., maxdist = 0.;
1346   gp_Pnt pbound, papp, pbid;
1347   gp_Vec vbound, vapp;
1348   Handle(GeomFill_Boundary) bou = ptch->Bound(I);
1349   for (Standard_Integer iu = 0; iu <= nbp; iu++){
1350     Standard_Real uu = iu;
1351     uu /= nbp;
1352     pbound = bou->Value(uu);
1353     BSplCLib::D0(uu,0,degree[I%2],0,ncpol[I]->Array1(),BSplCLib::NoWeights(),
1354                  nk[I%2]->Array1(),nm[I%2]->Array1(),papp);
1355     if(donor) {
1356       BSplCLib::D0(uu,0,degree[I%2],0,ntpol[I]->Array1(),BSplCLib::NoWeights(),
1357                    nk[I%2]->Array1(),nm[I%2]->Array1(),pbid);
1358       vapp.SetXYZ(pbid.XYZ());
1359       vbound = bou->Norm(uu);
1360       if(vapp.Magnitude() > 1.e-15 && vbound.Magnitude() > 1.e-15){
1361         Standard_Real alpha = Abs(M_PI/2.-Abs(vbound.Angle(vapp)));
1362         if(Abs(alpha) > maxang) maxang = Abs(alpha);
1363       }
1364 #ifdef DRAW
1365       Handle(Draw_Segment3D) seg;
1366       gp_Pnt pp;
1367       pp = pbound.Translated(vbound);
1368       seg = new Draw_Segment3D(pbound,pp,Draw_blanc);
1369       dout << seg;
1370       pp = papp.Translated(vapp);
1371       seg = new Draw_Segment3D(papp,pp,Draw_rouge);
1372       dout << seg;
1373 #endif
1374     }
1375     if(papp.Distance(pbound) > maxdist) maxdist = papp.Distance(pbound);
1376   }
1377   cout<<"Controle approx/contrainte sur bord "<<I<<" : "<<endl;
1378   cout<<"Distance max : "<<maxdist<<endl;
1379   if (donor) {
1380     maxang = maxang*180./M_PI;
1381     cout<<"Angle max    : "<<maxang<<" deg"<<endl;
1382   }
1383 }
1384
1385
1386 //=======================================================================
1387 //function : CheckResult
1388 //purpose  :
1389 //=======================================================================
1390
1391 void GeomFill_ConstrainedFilling::CheckResult(const Standard_Integer I)
1392 {
1393   Standard_Boolean donor = !tgalg[I].IsNull();
1394   Standard_Real maxang = 0., maxdist = 0.;
1395   Standard_Real uu=0,duu=0,vv=0,dvv=0,ww=0,dww=0,u1,u2,v1,v2;
1396   surf->Bounds(u1,u2,v1,v2);
1397   switch(I){
1398   case 0:
1399     uu = ww = u1;
1400     vv = v1;
1401     duu = dww = (u2 - u1)/30;
1402     dvv = 0.;
1403     break;
1404   case 1:
1405     vv = ww = v1;
1406     uu = u2;
1407     dvv = dww = (v2 - v1)/30;
1408     duu = 0.;
1409     break;
1410   case 2:
1411     uu = ww = u1;
1412     vv = v2;
1413     duu = dww = (u2 - u1)/30;
1414     dvv = 0.;
1415     break;
1416   case 3:
1417     vv = ww = v1;
1418     uu = u1;
1419     dvv = dww = (v2 - v1)/30;
1420     duu = 0.;
1421     break;
1422   }
1423   gp_Pnt pbound[31],pres[31];
1424   gp_Vec vbound[31],vres[31];
1425 #ifdef DRAW
1426   Standard_Real ang[31];
1427   Standard_Boolean hasang[31];
1428 #endif
1429   Handle(GeomFill_Boundary) bou = ptch->Bound(I);
1430   Standard_Integer k ;
1431   for ( k = 0; k <= 30; k++){
1432     pbound[k] = bou->Value(ww);
1433     if(!donor) surf->D0(uu,vv,pres[k]);
1434     else{
1435       vbound[k] = bou->Norm(ww);
1436       gp_Vec V1,V2;
1437       surf->D1(uu,vv,pres[k],V1,V2);
1438       vres[k] = V1.Crossed(V2);
1439       if(vres[k].Magnitude() > 1.e-15 && vbound[k].Magnitude() > 1.e-15){
1440         Standard_Real alpha = Abs(vres[k].Angle(vbound[k]));
1441         alpha = Min(alpha,Abs(M_PI-alpha));
1442         if(alpha > maxang) maxang = alpha;
1443 #ifdef DRAW
1444         ang[k] = alpha;
1445         hasang[k] = 1;
1446 #endif
1447       }
1448 #ifdef DRAW
1449       else hasang[k] = 0;
1450 #endif
1451     }
1452     if(pres[k].Distance(pbound[k]) > maxdist) maxdist = pres[k].Distance(pbound[k]);
1453     uu += duu;
1454     vv += dvv;
1455     ww += dww;
1456   }
1457   cout<<"Controle resultat/contrainte sur bord "<<I<<" : "<<endl;
1458   cout<<"Distance max : "<<maxdist<<endl;
1459   if (donor) {
1460     Standard_Real angdeg = maxang*180./M_PI;
1461     cout<<"Angle max    : "<<angdeg<<" deg"<<endl;
1462   }
1463 #ifdef DRAW
1464   Standard_Boolean scale = maxang>1.e-10;
1465   for (k = 0; k <= 30; k++){
1466     if(hasang[k]){
1467       gp_Pnt pp;
1468       Handle(Draw_Segment3D) seg;
1469       vbound[k].Normalize();
1470       if(scale) vbound[k].Multiply(1.+3.*ang[k]/maxang);
1471       vbound[k].Multiply(drawfac);
1472       pp = pbound[k].Translated(vbound[k]);
1473       seg = new Draw_Segment3D(pbound[k],pp,Draw_blanc);
1474       dout << seg;
1475       vres[k].Normalize();
1476       if(scale) vres[k].Multiply(1.+3.*ang[k]/maxang);
1477       vres[k].Multiply(drawfac);
1478       pp = pres[k].Translated(vres[k]);
1479       seg = new Draw_Segment3D(pres[k],pp,Draw_rouge);
1480       dout << seg;
1481     }
1482   }
1483 #endif
1484 }
1485