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