src/GeomFill/GeomFill_ConstrainedFilling.cxx

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