0024510: Remove unused local variables
[occt.git] / src / GeomFill / GeomFill_ConstrainedFilling.cxx
CommitLineData
b311480e 1// Created on: 1995-10-26
2// Created by: Laurent BOURESCHE
3// Copyright (c) 1995-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
973c2be1 8// This library is free software; you can redistribute it and / or modify it
9// under the terms of the GNU Lesser General Public 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 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>
45static Standard_Boolean dodraw = 0;
46static Standard_Real drawfac = 0.1;
47#endif
48#ifdef DEB
49Standard_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);
55Standard_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>
62static OSD_Chronometer totclock, parclock, appclock, cstclock;
63#endif
64
65static 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
96static 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
108static 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();
c6541a0c 169 Standard_Real ang = M_PI - tgi.Angle(tgn);
7fd59977 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}
180static 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;
c6541a0c 197 if(an >= 0.5*M_PI) { twist = Standard_True; an = M_PI-an; }
7fd59977 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 DEB
221 cout<<"pb coons cnd coin : "<<i<<" fact = "<<killfactor<<endl;
222#endif
223 }
224 }
225 }
226 }
227}
228static 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 DEB
252 parclock.Start();
253#endif
254 bound[i]->Reparametrize(0.,1.,fnul,lnul,fscal,lscal,rev[i]);
255#ifdef DEB
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
274class 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
291void 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
306GeomFill_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
320void 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 DEB
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 DEB
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 DEB
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
414void 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 DEB
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 DEB
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 DEB
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
495void 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
510void 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
526Handle(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
538Handle(Geom_BSplineSurface) GeomFill_ConstrainedFilling::Surface() const
539{
540 return surf;
541}
542
543
544//=======================================================================
545//function : Build
546//purpose :
547//=======================================================================
548
549void 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 DEB
565 appclock.Start();
566#endif
567 if(nbd3) PerformApprox();
568#ifdef DEB
569 appclock.Stop();
570#endif
571 }
572 appdone = Standard_True;
573#ifdef DEB
574 cstclock.Start();
575#endif
576 MatchKnots();
577 PerformS0();
578 PerformS1();
579 PerformSurface();
580#ifdef DEB
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
605void 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
667void 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
863void 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
917void 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
1112void 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
1134Standard_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
1163void 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
1182Standard_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
1230void 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
1289void 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){
c6541a0c 1326 Standard_Real alpha = Abs(M_PI/2.-Abs(vnor.Angle(vtg)));
7fd59977 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
1340void 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){
c6541a0c 1360 Standard_Real alpha = Abs(M_PI/2.-Abs(vbound.Angle(vapp)));
7fd59977 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) {
c6541a0c 1379 maxang = maxang*180./M_PI;
7fd59977 1380 cout<<"Angle max : "<<maxang<<" deg"<<endl;
1381 }
1382}
1383
1384
1385//=======================================================================
1386//function : CheckResult
1387//purpose :
1388//=======================================================================
1389
1390void 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];
96a95605 1424#ifdef DRAW
7fd59977 1425 Standard_Real ang[31];
1426 Standard_Boolean hasang[31];
96a95605 1427#endif
7fd59977 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]));
c6541a0c 1440 alpha = Min(alpha,Abs(M_PI-alpha));
7fd59977 1441 if(alpha > maxang) maxang = alpha;
96a95605 1442#ifdef DRAW
7fd59977 1443 ang[k] = alpha;
1444 hasang[k] = 1;
96a95605 1445#endif
7fd59977 1446 }
96a95605 1447#ifdef DRAW
7fd59977 1448 else hasang[k] = 0;
96a95605 1449#endif
7fd59977 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) {
c6541a0c 1459 Standard_Real angdeg = maxang*180./M_PI;
7fd59977 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