src/GeomFill/GeomFill_ConstrainedFilling.cxx

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