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