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