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