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