1 // Created on: 1995-10-19
2 // Created by: Andre LIEUTIER
3 // Copyright (c) 1995-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
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
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.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <gp_XY.hxx>
18 #include <gp_XYZ.hxx>
19 #include <math_Gauss.hxx>
20 #include <math_Matrix.hxx>
21 #include <math_Vector.hxx>
22 #include <Plate_FreeGtoCConstraint.hxx>
23 #include <Plate_GlobalTranslationConstraint.hxx>
24 #include <Plate_GtoCConstraint.hxx>
25 #include <Plate_LinearScalarConstraint.hxx>
26 #include <Plate_LinearXYZConstraint.hxx>
27 #include <Plate_LineConstraint.hxx>
28 #include <Plate_PinpointConstraint.hxx>
29 #include <Plate_PlaneConstraint.hxx>
30 #include <Plate_Plate.hxx>
31 #include <Plate_SampledCurveConstraint.hxx>
32 #include <Standard_ErrorHandler.hxx>
33
34 //=======================================================================
35 //function : Plate_Plate
36 //purpose :
37 //=======================================================================
Plate_Plate()38 Plate_Plate::Plate_Plate()
39 : order(0), n_el(0), n_dim(0),
40 solution(0),points(0),deru(0),derv(0),
41 OK(Standard_False),maxConstraintOrder(0),
42 Uold (1.e20),
43 Vold (1.e20),
44 U2 (0.0),
45 R (0.0),
46 L (0.0)
47 {
48 PolynomialPartOnly = Standard_False;
49 memset (ddu, 0, sizeof (ddu));
50 memset (ddv, 0, sizeof (ddv));
51 }
52
53 //=======================================================================
54 //function : Plate_Plate
55 //purpose :
56 //=======================================================================
57
Plate_Plate(const Plate_Plate & Ref)58 Plate_Plate::Plate_Plate(const Plate_Plate& Ref)
59 : order(Ref.order),n_el(Ref.n_el),n_dim(Ref.n_dim),
60 solution(0),points(0),deru(0),derv(0),
61 OK (Ref.OK),
62 Uold (1.e20),
63 Vold (1.e20),
64 U2 (0.0),
65 R (0.0),
66 L (0.0)
67 {
68 Standard_Integer i;
69 if (Ref.OK) {
70 if (n_dim >0 && Ref.solution != 0) {
71 solution = new gp_XYZ[n_dim];
72 for(i=0; i<n_dim ;i++) {
73 Solution(i) = Ref.Solution(i);
74 }
75 }
76
77 if (n_el >0) {
78 if (Ref.points != 0) {
79 points = new gp_XY[n_el];
80 for(i=0; i<n_el;i++) {
81 Points(i) = Ref.Points(i);
82 }
83 }
84
85 if (Ref.deru != 0) {
86 deru = new Standard_Integer[n_el] ;
87 for (i = 0 ; i < n_el ; i++) {
88 Deru(i) = Ref.Deru(i);
89 }
90 }
91
92 if (Ref.derv != 0) {
93 derv = new Standard_Integer[n_el] ;
94 for (i = 0 ; i < n_el ; i++) {
95 Derv(i) = Ref.Derv(i);
96 }
97 }
98 }
99 }
100
101 myConstraints = Ref.myConstraints;
102 myLXYZConstraints = Ref.myLXYZConstraints;
103 myLScalarConstraints = Ref.myLScalarConstraints;
104 maxConstraintOrder = Ref.maxConstraintOrder;
105 PolynomialPartOnly = Ref.PolynomialPartOnly;
106 for (i=0; i<10;i++) {
107 ddu[i]=Ref.ddu[i];
108 ddv[i]=Ref.ddv[i];
109 }
110 }
111 //=======================================================================
112 //function : Copy
113 //purpose :
114 //=======================================================================
Copy(const Plate_Plate & Ref)115 Plate_Plate& Plate_Plate::Copy(const Plate_Plate& Ref)
116 {
117 Init();
118 order = Ref.order;
119 n_el = Ref.n_el;
120 n_dim = Ref.n_dim;
121 OK = Ref.OK;
122 Standard_Integer i;
123 if (Ref.OK) {
124 if (n_dim >0 && Ref.solution != 0) {
125 solution = new gp_XYZ[n_dim];
126 for(i=0; i<n_dim ;i++) {
127 Solution(i) = Ref.Solution(i);
128 }
129 }
130
131 if (n_el >0) {
132 if (Ref.points != 0) {
133 points = new gp_XY[n_el];
134 for(i=0; i<n_el;i++) {
135 Points(i) = Ref.Points(i);
136 }
137 }
138
139 if (Ref.deru != 0) {
140 deru = new Standard_Integer[n_el] ;
141 for (i = 0 ; i < n_el ; i++) {
142 Deru(i) = Ref.Deru(i);
143 }
144 }
145
146 if (Ref.derv != 0) {
147 derv = new Standard_Integer[n_el] ;
148 for (i = 0 ; i < n_el ; i++) {
149 Derv(i) = Ref.Derv(i);
150 }
151 }
152 }
153 }
154
155 myConstraints = Ref.myConstraints;
156 myLXYZConstraints = Ref.myLXYZConstraints;
157 myLScalarConstraints = Ref.myLScalarConstraints;
158 maxConstraintOrder = Ref.maxConstraintOrder;
159 PolynomialPartOnly = Ref.PolynomialPartOnly;
160
161 for (i=0; i<10;i++) {
162 ddu[i]=Ref.ddu[i];
163 ddv[i]=Ref.ddv[i];
164 }
165 return *this;
166 }
167 //=======================================================================
168 //function : Load
169 //purpose :
170 //=======================================================================
171
Load(const Plate_PinpointConstraint & PConst)172 void Plate_Plate::Load(const Plate_PinpointConstraint& PConst)
173 {
174 OK = Standard_False;
175 n_el++;
176 myConstraints.Append(PConst);
177 Standard_Integer OrdreConst = PConst.Idu() + PConst.Idv();
178 if(maxConstraintOrder<OrdreConst) maxConstraintOrder = OrdreConst;
179 }
180
Load(const Plate_LinearXYZConstraint & LXYZConst)181 void Plate_Plate::Load(const Plate_LinearXYZConstraint& LXYZConst)
182 {
183 OK = Standard_False;
184 n_el += LXYZConst.Coeff().RowLength();
185
186 myLXYZConstraints.Append(LXYZConst);
187 for(Standard_Integer j=1;j <= LXYZConst.GetPPC().Length() ; j++)
188 {
189 Standard_Integer OrdreConst = LXYZConst.GetPPC()(j).Idu() + LXYZConst.GetPPC()(j).Idv();
190 if(maxConstraintOrder<OrdreConst) maxConstraintOrder = OrdreConst;
191 }
192 }
193
Load(const Plate_LinearScalarConstraint & LScalarConst)194 void Plate_Plate::Load(const Plate_LinearScalarConstraint& LScalarConst)
195 {
196 OK = Standard_False;
197 n_el += LScalarConst.Coeff().RowLength();
198 myLScalarConstraints.Append(LScalarConst);
199 for(Standard_Integer j=1;j <= LScalarConst.GetPPC().Length() ; j++)
200 {
201 Standard_Integer OrdreConst = LScalarConst.GetPPC()(j).Idu() + LScalarConst.GetPPC()(j).Idv();
202 if(maxConstraintOrder<OrdreConst) maxConstraintOrder = OrdreConst;
203 }
204 }
205
Load(const Plate_LineConstraint & LConst)206 void Plate_Plate::Load(const Plate_LineConstraint& LConst)
207 {
208 Load(LConst.LSC());
209 }
210
Load(const Plate_PlaneConstraint & PConst)211 void Plate_Plate::Load(const Plate_PlaneConstraint& PConst)
212 {
213 Load(PConst.LSC());
214 }
215
Load(const Plate_SampledCurveConstraint & SCConst)216 void Plate_Plate::Load(const Plate_SampledCurveConstraint& SCConst)
217 {
218 Load(SCConst.LXYZC());
219 }
220
Load(const Plate_GtoCConstraint & GtoCConst)221 void Plate_Plate::Load(const Plate_GtoCConstraint& GtoCConst)
222 {
223 for(Standard_Integer i=0;i< GtoCConst.nb_PPC();i++)
224 Load(GtoCConst.GetPPC(i));
225 }
226
Load(const Plate_FreeGtoCConstraint & FGtoCConst)227 void Plate_Plate::Load(const Plate_FreeGtoCConstraint& FGtoCConst)
228 {
229 Standard_Integer i ;
230 for( i=0;i< FGtoCConst.nb_PPC();i++)
231 Load(FGtoCConst.GetPPC(i));
232 for(i=0;i< FGtoCConst.nb_LSC();i++)
233 Load(FGtoCConst.LSC(i));
234 }
235
Load(const Plate_GlobalTranslationConstraint & GTConst)236 void Plate_Plate::Load(const Plate_GlobalTranslationConstraint& GTConst)
237 {
238 Load(GTConst.LXYZC());
239 }
240
241 //=======================================================================
242 //function : SolveTI
243 //purpose : to solve the set of constraints
244 //=======================================================================
245
SolveTI(const Standard_Integer ord,const Standard_Real anisotropie,const Message_ProgressRange & theProgress)246 void Plate_Plate::SolveTI(const Standard_Integer ord,
247 const Standard_Real anisotropie,
248 const Message_ProgressRange& theProgress)
249 {
250 Standard_Integer IterationNumber=0;
251 OK = Standard_False;
252 order = ord;
253 if(ord <=1) return;
254 if(ord > 9) return;
255 if(n_el<1) return;
256 if(anisotropie < 1.e-6) return;
257 if(anisotropie > 1.e+6) return;
258
259 // computation of the bounding box of the 2d PPconstraints
260 Standard_Real xmin,xmax,ymin,ymax;
261 UVBox(xmin,xmax,ymin,ymax);
262
263 Standard_Real du = 0.5*(xmax - xmin);
264 if(anisotropie >1.) du *= anisotropie;
265 if(du < 1.e-10) return;
266 ddu[0] = 1;
267 Standard_Integer i ;
268 for( i=1;i<=9;i++) ddu[i] = ddu[i-1] / du;
269
270 Standard_Real dv = 0.5*(ymax - ymin);
271 if(anisotropie <1.) dv /= anisotropie;
272 if(dv < 1.e-10) return;
273 ddv[0] = 1;
274 for(i=1;i<=9;i++) ddv[i] = ddv[i-1] / dv;
275
276 if(myLScalarConstraints.IsEmpty())
277 {
278 if(myLXYZConstraints.IsEmpty())
279 SolveTI1(IterationNumber, theProgress);
280 else
281 SolveTI2(IterationNumber, theProgress);
282 }
283 else
284 SolveTI3(IterationNumber, theProgress);
285
286 }
287
288 //=======================================================================
289 //function : SolveTI1
290 //purpose : to solve the set of constraints in the easiest case,
291 // only PinPointConstraints are loaded
292 //=======================================================================
293
SolveTI1(const Standard_Integer IterationNumber,const Message_ProgressRange & theProgress)294 void Plate_Plate::SolveTI1(const Standard_Integer IterationNumber,
295 const Message_ProgressRange& theProgress)
296 {
297 // computation of square matrix members
298
299
300 n_dim = n_el + order*(order+1)/2;
301 math_Matrix mat(0, n_dim-1, 0, n_dim-1, 0.);
302
303 delete [] (gp_XY*)points;
304 points = new gp_XY[n_el];
305 Standard_Integer i ;
306 for( i=0; i<n_el;i++) Points(i) = myConstraints(i+1).Pnt2d();
307
308 delete [] (Standard_Integer*)deru;
309 deru = new Standard_Integer[n_el];
310 for(i=0; i<n_el;i++) Deru(i) = myConstraints(i+1).Idu();
311
312 delete [] (Standard_Integer*)derv;
313 derv = new Standard_Integer[n_el];
314 for(i=0; i<n_el;i++) Derv(i) = myConstraints(i+1).Idv();
315
316 for(i=0; i<n_el;i++) {
317 for(Standard_Integer j=0;j<i;j++) {
318 Standard_Real signe = 1;
319 if ( ((Deru(j)+Derv(j))%2) == 1) signe = -1;
320 Standard_Integer iu = Deru(i) + Deru(j);
321 Standard_Integer iv = Derv(i) + Derv(j);
322 mat(i,j) = signe * SolEm(Points(i) - Points(j),iu,iv);
323 }
324 }
325
326 i = n_el;
327 for(Standard_Integer iu = 0; iu< order; iu++) {
328 for(Standard_Integer iv =0; iu+iv < order; iv++) {
329 for(Standard_Integer j=0;j<n_el;j++) {
330 Standard_Integer idu = Deru(j);
331 Standard_Integer idv = Derv(j);
332 mat(i,j) = Polm (Points(j), iu, iv, idu, idv);
333 }
334 i++;
335 }
336 }
337
338 for(i=0;i<n_dim;i++) {
339 for(Standard_Integer j = i+1; j<n_dim;j++) {
340 mat(i,j) = mat(j,i);
341 }
342 }
343
344 // initialisation of the Gauss algorithm
345 Standard_Real pivot_max = 1.e-12;
346 OK = Standard_True;
347
348 Message_ProgressScope aScope (theProgress, "Plate_Plate::SolveTI1()", 10);
349 math_Gauss algo_gauss(mat,pivot_max, aScope.Next (7));
350
351 if (aScope.UserBreak())
352 {
353 OK = Standard_False;
354 return;
355 }
356
357 if(!algo_gauss.IsDone()) {
358 Standard_Integer nbm = order*(order+1)/2;
359 for(i=n_el;i<n_el+nbm;i++) {
360 mat(i,i) = 1.e-8;
361 }
362 pivot_max = 1.e-18;
363
364 math_Gauss thealgo(mat,pivot_max, aScope.Next (3));
365
366 if (aScope.UserBreak())
367 {
368 OK = Standard_False;
369 return;
370 }
371 algo_gauss = thealgo;
372 OK = algo_gauss.IsDone();
373 }
374
375 if (OK) {
376 // computation of the linear system solution for the X, Y and Z coordinates
377 math_Vector sec_member( 0, n_dim-1, 0.);
378 math_Vector sol(0,n_dim-1);
379
380 delete [] (gp_XYZ*) solution;
381 solution = new gp_XYZ[n_dim];
382
383 for(Standard_Integer icoor=1; icoor<=3;icoor++) {
384 for(i=0;i<n_el;i++) {
385 sec_member(i) = myConstraints(i+1).Value().Coord(icoor);
386 }
387 algo_gauss.Solve(sec_member, sol);
388 //alr iteration pour affiner la solution
389 {
390 math_Vector sol1(0,n_dim-1);
391 math_Vector sec_member1(0,n_dim-1);
392 for(i=1;i<=IterationNumber;i++)
393 {
394 sec_member1 = sec_member - mat*sol;
395 algo_gauss.Solve(sec_member1, sol1);
396 sol += sol1;
397 }
398 }
399 //finalr
400
401 for(i=0;i<n_dim;i++) {
402 Solution(i).SetCoord (icoor, sol(i));
403 }
404 }
405 }
406 }
407
408 //=======================================================================
409 //function : SolveTI2
410 //purpose : to solve the set of constraints in the medium case,
411 // LinearXYZ constraints are provided but no LinearScalar one
412 //=======================================================================
413
SolveTI2(const Standard_Integer IterationNumber,const Message_ProgressRange & theProgress)414 void Plate_Plate::SolveTI2(const Standard_Integer IterationNumber,
415 const Message_ProgressRange& theProgress)
416 {
417 // computation of square matrix members
418
419 Standard_Integer nCC1 = myConstraints.Length();
420 Standard_Integer nCC2 = 0;
421 Standard_Integer i ;
422 for( i = 1; i<= myLXYZConstraints.Length(); i++)
423 nCC2 += myLXYZConstraints(i).Coeff().ColLength();
424
425 Standard_Integer n_dimat = nCC1 + nCC2 + order*(order+1)/2;
426
427
428 delete [] (gp_XY*)points;
429 points = new gp_XY[n_el];
430 delete [] (Standard_Integer*)deru;
431 deru = new Standard_Integer[n_el];
432 delete [] (Standard_Integer*)derv;
433 derv = new Standard_Integer[n_el];
434
435
436 for(i=0; i< nCC1;i++)
437 {
438 Points(i) = myConstraints(i+1).Pnt2d();
439 Deru(i) = myConstraints(i+1).Idu();
440 Derv(i) = myConstraints(i+1).Idv();
441 }
442
443 Standard_Integer k = nCC1;
444 for( i = 1; i<= myLXYZConstraints.Length(); i++)
445 for(Standard_Integer j=1;j <= myLXYZConstraints(i).GetPPC().Length() ; j++)
446 {
447 Points(k) = myLXYZConstraints(i).GetPPC()(j).Pnt2d();
448 Deru(k) = myLXYZConstraints(i).GetPPC()(j).Idu();
449 Derv(k) = myLXYZConstraints(i).GetPPC()(j).Idv();
450 k++;
451 }
452
453 math_Matrix mat(0, n_dimat-1, 0, n_dimat-1, 0.);
454
455 fillXYZmatrix(mat,0,0,nCC1,nCC2);
456
457
458 // initialisation of the Gauss algorithm
459 Standard_Real pivot_max = 1.e-12;
460 OK = Standard_True; // ************ JHH
461
462 Message_ProgressScope aScope (theProgress, "Plate_Plate::SolveTI2()", 10);
463 math_Gauss algo_gauss(mat,pivot_max, aScope.Next (7));
464
465 if (aScope.UserBreak())
466 {
467 OK = Standard_False;
468 return;
469 }
470
471 if(!algo_gauss.IsDone()) {
472 for(i=nCC1+nCC2;i<n_dimat;i++) {
473 mat(i,i) = 1.e-8;
474 }
475 pivot_max = 1.e-18;
476
477 math_Gauss thealgo1(mat,pivot_max, aScope.Next (3));
478
479 if (aScope.UserBreak())
480 {
481 OK = Standard_False;
482 return;
483 }
484 algo_gauss = thealgo1;
485 OK = algo_gauss.IsDone();
486 }
487
488 if (OK) {
489 // computation of the linear system solution for the X, Y and Z coordinates
490 math_Vector sec_member( 0, n_dimat-1, 0.);
491 math_Vector sol(0,n_dimat-1);
492
493 delete [] (gp_XYZ*) solution;
494 n_dim = n_el+order*(order+1)/2;
495 solution = new gp_XYZ[n_dim];
496
497 for(Standard_Integer icoor=1; icoor<=3;icoor++) {
498 for(i=0;i<nCC1;i++) {
499 sec_member(i) = myConstraints(i+1).Value().Coord(icoor);
500 }
501
502 k = nCC1;
503 for(i = 1; i<= myLXYZConstraints.Length(); i++) {
504 for(Standard_Integer irow =1; irow <= myLXYZConstraints(i).Coeff().ColLength(); irow++) {
505 for(Standard_Integer icol=1; icol<=myLXYZConstraints(i).Coeff().RowLength();icol++)
506 sec_member(k) += myLXYZConstraints(i).Coeff()(irow,icol)
507 * myLXYZConstraints(i).GetPPC()(icol).Value().Coord(icoor);
508 k++;
509 }
510 }
511
512 algo_gauss.Solve(sec_member, sol);
513 //alr iteration pour affiner la solution
514 {
515 math_Vector sol1(0,n_dimat-1);
516 math_Vector sec_member1(0,n_dimat-1);
517 for(i=1;i<=IterationNumber;i++)
518 {
519 sec_member1 = sec_member - mat*sol;
520 algo_gauss.Solve(sec_member1, sol1);
521 sol += sol1;
522 }
523 }
524 //finalr
525
526 for(i=0;i<nCC1;i++) Solution(i).SetCoord (icoor, sol(i));
527
528 Standard_Integer kSolution = nCC1;
529 Standard_Integer ksol = nCC1;
530
531 for(i = 1; i<= myLXYZConstraints.Length(); i++) {
532 for(Standard_Integer icol=1; icol<=myLXYZConstraints(i).Coeff().RowLength();icol++){
533 Standard_Real vsol = 0;
534 for(Standard_Integer irow =1; irow <= myLXYZConstraints(i).Coeff().ColLength(); irow++)
535 vsol += myLXYZConstraints(i).Coeff()(irow,icol)*sol(ksol+irow-1);
536 Solution(kSolution).SetCoord (icoor, vsol);
537 kSolution++;
538 }
539 ksol += myLXYZConstraints(i).Coeff().ColLength();
540 }
541
542 for(i=0;i<order*(order+1)/2; i++) {
543 Solution(n_el+i).SetCoord (icoor, sol(ksol+i));
544 }
545 }
546 }
547 }
548
549 //=======================================================================
550 //function : SolveTI3
551 //purpose : to solve the set of constraints in the most general situation
552 //=======================================================================
553
SolveTI3(const Standard_Integer IterationNumber,const Message_ProgressRange & theProgress)554 void Plate_Plate::SolveTI3(const Standard_Integer IterationNumber,
555 const Message_ProgressRange& theProgress)
556 {
557 // computation of square matrix members
558
559 Standard_Integer nCC1 = myConstraints.Length();
560
561 Standard_Integer nCC2 = 0;
562 Standard_Integer i ;
563 for( i = 1; i<= myLXYZConstraints.Length(); i++)
564 nCC2 += myLXYZConstraints(i).Coeff().ColLength();
565
566 Standard_Integer nCC3 = 0;
567 for(i = 1; i<= myLScalarConstraints.Length(); i++)
568 nCC3 += myLScalarConstraints(i).Coeff().ColLength();
569
570 Standard_Integer nbm = order*(order+1)/2;
571 Standard_Integer n_dimsousmat = nCC1 + nCC2 + nbm ;
572 Standard_Integer n_dimat =3*n_dimsousmat + nCC3;
573
574
575 delete [] (gp_XY*)points;
576 points = new gp_XY[n_el];
577 delete [] (Standard_Integer*)deru;
578 deru = new Standard_Integer[n_el];
579 delete [] (Standard_Integer*)derv;
580 derv = new Standard_Integer[n_el];
581
582
583 for(i=0; i< nCC1;i++)
584 {
585 Points(i) = myConstraints(i+1).Pnt2d();
586 Deru(i) = myConstraints(i+1).Idu();
587 Derv(i) = myConstraints(i+1).Idv();
588 }
589
590 Standard_Integer k = nCC1;
591 for(i = 1; i<= myLXYZConstraints.Length(); i++)
592 for(Standard_Integer j=1;j <= myLXYZConstraints(i).GetPPC().Length() ; j++)
593 {
594 Points(k) = myLXYZConstraints(i).GetPPC()(j).Pnt2d();
595 Deru(k) = myLXYZConstraints(i).GetPPC()(j).Idu();
596 Derv(k) = myLXYZConstraints(i).GetPPC()(j).Idv();
597 k++;
598 }
599 Standard_Integer nPPC2 = k;
600 for(i = 1; i<= myLScalarConstraints.Length(); i++)
601 for(Standard_Integer j=1;j <= myLScalarConstraints(i).GetPPC().Length() ; j++)
602 {
603 Points(k) = myLScalarConstraints(i).GetPPC()(j).Pnt2d();
604 Deru(k) = myLScalarConstraints(i).GetPPC()(j).Idu();
605 Derv(k) = myLScalarConstraints(i).GetPPC()(j).Idv();
606 k++;
607 }
608
609 math_Matrix mat(0, n_dimat-1, 0, n_dimat-1, 0.);
610
611 fillXYZmatrix(mat,0,0,nCC1,nCC2);
612 fillXYZmatrix(mat,n_dimsousmat,n_dimsousmat,nCC1,nCC2);
613 fillXYZmatrix(mat,2*n_dimsousmat,2*n_dimsousmat,nCC1,nCC2);
614
615 k = 3*n_dimsousmat;
616 Standard_Integer kppc = nPPC2;
617 Standard_Integer j ;
618 for(i = 1; i<= myLScalarConstraints.Length(); i++) {
619 for( j=0;j<nCC1;j++){
620
621 math_Vector vmat(1,myLScalarConstraints(i).GetPPC().Length());
622
623 for(Standard_Integer ippc=1;ippc <= myLScalarConstraints(i).GetPPC().Length() ; ippc++) {
624 Standard_Real signe = 1;
625 if ( ((Deru(j)+Derv(j))%2) == 1) signe = -1;
626 Standard_Integer iu = Deru(kppc+ippc-1) + Deru(j);
627 Standard_Integer iv = Derv(kppc+ippc-1) + Derv(j);
628 vmat(ippc) = signe * SolEm(Points(kppc+ippc-1) - Points(j),iu,iv);
629 }
630
631 for(Standard_Integer irow=1;irow <= myLScalarConstraints(i).Coeff().ColLength() ; irow++)
632 for(Standard_Integer icol=1;icol <= myLScalarConstraints(i).Coeff().RowLength() ; icol++){
633 mat(k+irow-1,j) += myLScalarConstraints(i).Coeff()(irow,icol).X()*vmat(icol);
634 mat(k+irow-1,n_dimsousmat+j) += myLScalarConstraints(i).Coeff()(irow,icol).Y()*vmat(icol);
635 mat(k+irow-1,2*n_dimsousmat+j) += myLScalarConstraints(i).Coeff()(irow,icol).Z()*vmat(icol);
636 }
637 }
638
639 Standard_Integer k2 = nCC1;
640 Standard_Integer kppc2 = nCC1;
641 Standard_Integer i2 ;
642 for( i2 = 1; i2<=myLXYZConstraints.Length() ; i2++){
643
644 math_Matrix tmpmat(1,myLScalarConstraints(i).GetPPC().Length(),1,myLXYZConstraints(i2).GetPPC().Length() );
645
646 for(Standard_Integer ippc=1;ippc <= myLScalarConstraints(i).GetPPC().Length() ; ippc++)
647 for(Standard_Integer ippc2=1;ippc2 <= myLXYZConstraints(i2).GetPPC().Length() ; ippc2++){
648 Standard_Real signe = 1;
649 if ( ((Deru(kppc2+ippc2-1)+Derv(kppc2+ippc2-1))%2) == 1) signe = -1;
650 Standard_Integer iu = Deru(kppc+ippc-1) + Deru(kppc2+ippc2-1);
651 Standard_Integer iv = Derv(kppc+ippc-1) + Derv(kppc2+ippc2-1);
652 tmpmat(ippc,ippc2) = signe * SolEm(Points(kppc+ippc-1) - Points(kppc2+ippc2-1),iu,iv);
653 }
654
655 for(Standard_Integer irow=1;irow <= myLScalarConstraints(i).Coeff().ColLength() ; irow++)
656 for(Standard_Integer irow2=1;irow2 <= myLXYZConstraints(i2).Coeff().ColLength() ; irow2++)
657 for(Standard_Integer icol=1;icol <= myLScalarConstraints(i).Coeff().RowLength() ; icol++)
658 for(Standard_Integer icol2=1;icol2 <= myLXYZConstraints(i2).Coeff().RowLength() ; icol2++){
659 mat(k+irow-1,k2+irow2-1) +=
660 myLScalarConstraints(i).Coeff()(irow,icol).X()*myLXYZConstraints(i2).Coeff()(irow2,icol2)*tmpmat(icol,icol2);
661 mat(k+irow-1,n_dimsousmat+k2+irow2-1) +=
662 myLScalarConstraints(i).Coeff()(irow,icol).Y()*myLXYZConstraints(i2).Coeff()(irow2,icol2)*tmpmat(icol,icol2);
663 mat(k+irow-1,2*n_dimsousmat+k2+irow2-1) +=
664 myLScalarConstraints(i).Coeff()(irow,icol).Z()*myLXYZConstraints(i2).Coeff()(irow2,icol2)*tmpmat(icol,icol2);
665 }
666
667 k2 += myLXYZConstraints(i2).Coeff().ColLength();
668 kppc2 += myLXYZConstraints(i2).Coeff().RowLength();
669 }
670
671
672
673 j = nCC1+nCC2;
674 for(Standard_Integer iu = 0; iu< order; iu++)
675 for(Standard_Integer iv =0; iu+iv < order; iv++) {
676
677 math_Vector vmat(1,myLScalarConstraints(i).GetPPC().Length());
678 for(Standard_Integer ippc=1;ippc <= myLScalarConstraints(i).GetPPC().Length() ; ippc++){
679 Standard_Integer idu = Deru(kppc+ippc-1);
680 Standard_Integer idv = Derv(kppc+ippc-1);
681 vmat(ippc) = Polm (Points(kppc+ippc-1),iu,iv,idu,idv);
682 }
683
684 for(Standard_Integer irow=1;irow <= myLScalarConstraints(i).Coeff().ColLength() ; irow++)
685 for(Standard_Integer icol=1;icol <= myLScalarConstraints(i).Coeff().RowLength() ; icol++){
686 mat(k+irow-1,j) += myLScalarConstraints(i).Coeff()(irow,icol).X()*vmat(icol);
687 mat(k+irow-1,n_dimsousmat+j) += myLScalarConstraints(i).Coeff()(irow,icol).Y()*vmat(icol);
688 mat(k+irow-1,2*n_dimsousmat+j) += myLScalarConstraints(i).Coeff()(irow,icol).Z()*vmat(icol);
689 }
690
691 j++;
692 }
693
694
695 k2 = 3*n_dimsousmat;
696 kppc2 = nPPC2;
697 for(i2 = 1; i2<=i ; i2++){
698
699 math_Matrix tmpmat(1,myLScalarConstraints(i).GetPPC().Length(),1,myLScalarConstraints(i2).GetPPC().Length() );
700
701 for(Standard_Integer ippc=1;ippc <= myLScalarConstraints(i).GetPPC().Length() ; ippc++)
702 for(Standard_Integer ippc2=1;ippc2 <= myLScalarConstraints(i2).GetPPC().Length() ; ippc2++){
703 Standard_Real signe = 1;
704 if ( ((Deru(kppc2+ippc2-1)+Derv(kppc2+ippc2-1))%2) == 1) signe = -1;
705 Standard_Integer a_iu = Deru(kppc+ippc-1) + Deru(kppc2+ippc2-1);
706 Standard_Integer iv = Derv(kppc+ippc-1) + Derv(kppc2+ippc2-1);
707 tmpmat(ippc,ippc2) = signe * SolEm(Points(kppc+ippc-1) - Points(kppc2+ippc2-1),a_iu,iv);
708 }
709
710 for(Standard_Integer irow=1;irow <= myLScalarConstraints(i).Coeff().ColLength() ; irow++)
711 for(Standard_Integer irow2=1;irow2 <= myLScalarConstraints(i2).Coeff().ColLength() ; irow2++)
712 for(Standard_Integer icol=1;icol <= myLScalarConstraints(i).Coeff().RowLength() ; icol++)
713 for(Standard_Integer icol2=1;icol2 <= myLScalarConstraints(i2).Coeff().RowLength() ; icol2++){
714 mat(k+irow-1,k2+irow2-1) +=
715 myLScalarConstraints(i).Coeff()(irow,icol)*myLScalarConstraints(i2).Coeff()(irow2,icol2)*tmpmat(icol,icol2);
716 }
717
718 k2 += myLScalarConstraints(i2).Coeff().ColLength();
719 kppc2 += myLScalarConstraints(i2).Coeff().RowLength();
720 }
721
722 k += myLScalarConstraints(i).Coeff().ColLength();
723 kppc += myLScalarConstraints(i).Coeff().RowLength();
724 }
725
726 for( j=3*n_dimsousmat;j<n_dimat;j++)
727 for(i=0;i<j;i++)
728 mat(i,j)= mat(j,i);
729
730
731
732 // initialisation of the Gauss algorithm
733 Standard_Real pivot_max = 1.e-12;
734 OK = Standard_True; // ************ JHH
735
736 Message_ProgressScope aScope (theProgress, "Plate_Plate::SolveTI3()", 10);
737 math_Gauss algo_gauss(mat,pivot_max, aScope.Next (7));
738
739 if (aScope.UserBreak())
740 {
741 OK = Standard_False;
742 return;
743 }
744
745 if(!algo_gauss.IsDone()) {
746 for(i=nCC1+nCC2;i<nCC1+nCC2+nbm;i++) {
747 mat(i,i) = 1.e-8;
748 mat(n_dimsousmat+i,n_dimsousmat+i) = 1.e-8;
749 mat(2*n_dimsousmat+i,2*n_dimsousmat+i) = 1.e-8;
750 }
751 pivot_max = 1.e-18;
752
753 math_Gauss thealgo2(mat,pivot_max, aScope.Next (3));
754
755 if (aScope.UserBreak())
756 {
757 OK = Standard_False;
758 return;
759 }
760 algo_gauss = thealgo2;
761 OK = algo_gauss.IsDone();
762 }
763
764 if (OK) {
765 // computation of the linear system solution for the X, Y and Z coordinates
766 math_Vector sec_member( 0, n_dimat-1, 0.);
767 math_Vector sol(0,n_dimat-1);
768
769 delete [] (gp_XYZ*) solution;
770 n_dim = n_el+order*(order+1)/2;
771 solution = new gp_XYZ[n_dim];
772
773 Standard_Integer icoor ;
774 for( icoor=1; icoor<=3;icoor++){
775 for(i=0;i<nCC1;i++)
776 sec_member((icoor-1)*n_dimsousmat+i) = myConstraints(i+1).Value().Coord(icoor);
777
778
779 k = nCC1;
780 for(i = 1; i<= myLXYZConstraints.Length(); i++)
781 for(Standard_Integer irow =1; irow <= myLXYZConstraints(i).Coeff().ColLength(); irow++) {
782 for(Standard_Integer icol=1; icol<=myLXYZConstraints(i).Coeff().RowLength();icol++)
783 sec_member((icoor-1)*n_dimsousmat+k) += myLXYZConstraints(i).Coeff()(irow,icol)
784 * myLXYZConstraints(i).GetPPC()(icol).Value().Coord(icoor);
785 k++;
786 }
787 }
788 k = 3*n_dimsousmat;
789 for(i = 1; i<= myLScalarConstraints.Length(); i++)
790 for(Standard_Integer irow =1; irow <= myLScalarConstraints(i).Coeff().ColLength(); irow++) {
791 for(Standard_Integer icol=1; icol<=myLScalarConstraints(i).Coeff().RowLength();icol++)
792 sec_member(k) += myLScalarConstraints(i).Coeff()(irow,icol)
793 * myLScalarConstraints(i).GetPPC()(icol).Value();
794 k++;
795 }
796
797 algo_gauss.Solve(sec_member, sol);
798 // iteration to refine the solution
799 {
800 math_Vector sol1(0,n_dimat-1);
801 math_Vector sec_member1(0,n_dimat-1);
802 for(i=1;i<=IterationNumber;i++)
803 {
804 sec_member1 = sec_member - mat*sol;
805 algo_gauss.Solve(sec_member1, sol1);
806 sol += sol1;
807 }
808 }
809
810 for(icoor=1; icoor<=3;icoor++){
811 for(i=0;i<nCC1;i++) Solution(i).SetCoord (icoor, sol((icoor-1)*n_dimsousmat+i));
812
813 Standard_Integer kSolution = nCC1;
814 Standard_Integer ksol = nCC1;
815
816 for(i = 1; i<= myLXYZConstraints.Length(); i++) {
817 for(Standard_Integer icol=1; icol<=myLXYZConstraints(i).Coeff().RowLength();icol++){
818 Standard_Real vsol = 0;
819 for(Standard_Integer irow =1; irow <= myLXYZConstraints(i).Coeff().ColLength(); irow++)
820 vsol += myLXYZConstraints(i).Coeff()(irow,icol)*sol((icoor-1)*n_dimsousmat+ksol+irow-1);
821 Solution(kSolution).SetCoord (icoor, vsol);
822 kSolution++;
823 }
824 ksol += myLXYZConstraints(i).Coeff().ColLength();
825 }
826
827 ksol = nCC1+nCC2;
828 for(i=0;i<order*(order+1)/2; i++) {
829 Solution(n_el+i).SetCoord (icoor, sol((icoor-1)*n_dimsousmat+ksol+i));
830 }
831 }
832
833 Standard_Integer ksol = 3*n_dimsousmat;
834 Standard_Integer kSolution = nPPC2;
835 for(i = 1; i<= myLScalarConstraints.Length(); i++) {
836 for(Standard_Integer icol=1; icol<=myLScalarConstraints(i).Coeff().RowLength();icol++){
837 gp_XYZ Vsol(0.,0.,0.);
838 for(Standard_Integer irow =1; irow <= myLScalarConstraints(i).Coeff().ColLength(); irow++)
839 Vsol += myLScalarConstraints(i).Coeff()(irow,icol)*sol(ksol+irow-1);
840 Solution(kSolution) = Vsol;
841 kSolution++;
842 }
843 ksol += myLScalarConstraints(i).Coeff().ColLength();
844 }
845 }
846 }
847
848 //=======================================================================
849 //function : fillXYZmatrix
850 //purpose :
851 //=======================================================================
fillXYZmatrix(math_Matrix & mat,const Standard_Integer i0,const Standard_Integer j0,const Standard_Integer ncc1,const Standard_Integer ncc2) const852 void Plate_Plate::fillXYZmatrix(math_Matrix &mat,
853 const Standard_Integer i0,
854 const Standard_Integer j0,
855 const Standard_Integer ncc1,
856 const Standard_Integer ncc2) const
857 {
858 Standard_Integer i,j ;
859 for( i=0; i<ncc1;i++) {
860 for( j=0;j<i;j++) {
861 Standard_Real signe = 1;
862 if ( ((Deru(j)+Derv(j))%2) == 1) signe = -1;
863 Standard_Integer iu = Deru(i) + Deru(j);
864 Standard_Integer iv = Derv(i) + Derv(j);
865 mat(i0+i,j0+j) = signe * SolEm(Points(i) - Points(j),iu,iv);
866 }
867 }
868
869 Standard_Integer k = ncc1;
870 Standard_Integer kppc = ncc1;
871 for( i = 1; i<= myLXYZConstraints.Length(); i++){
872
873 for(Standard_Integer a_j=0; a_j < ncc1; a_j++){
874
875 math_Vector vmat(1,myLXYZConstraints(i).GetPPC().Length());
876
877 for(Standard_Integer ippc=1;ippc <= myLXYZConstraints(i).GetPPC().Length() ; ippc++) {
878 Standard_Real signe = 1;
879 if ( ((Deru(a_j)+Derv(a_j))%2) == 1) signe = -1;
880 Standard_Integer iu = Deru(kppc+ippc-1) + Deru(a_j);
881 Standard_Integer iv = Derv(kppc+ippc-1) + Derv(a_j);
882 vmat(ippc) = signe * SolEm(Points(kppc+ippc-1) - Points(a_j),iu,iv);
883 }
884
885 for(Standard_Integer irow=1;irow <= myLXYZConstraints(i).Coeff().ColLength() ; irow++)
886 for(Standard_Integer icol=1;icol <= myLXYZConstraints(i).Coeff().RowLength() ; icol++)
887 mat(i0+k+irow-1,j0+a_j) += myLXYZConstraints(i).Coeff()(irow,icol)*vmat(icol);
888 }
889
890 Standard_Integer k2 = ncc1;
891 Standard_Integer kppc2 = ncc1;
892 for(Standard_Integer i2 = 1; i2<= i; i2++){
893
894 math_Matrix tmpmat(1,myLXYZConstraints(i).GetPPC().Length(),1,myLXYZConstraints(i2).GetPPC().Length() );
895
896 for(Standard_Integer ippc=1;ippc <= myLXYZConstraints(i).GetPPC().Length() ; ippc++)
897 for(Standard_Integer ippc2=1;ippc2 <= myLXYZConstraints(i2).GetPPC().Length() ; ippc2++){
898 Standard_Real signe = 1;
899 if ( ((Deru(kppc2+ippc2-1)+Derv(kppc2+ippc2-1))%2) == 1) signe = -1;
900 Standard_Integer iu = Deru(kppc+ippc-1) + Deru(kppc2+ippc2-1);
901 Standard_Integer iv = Derv(kppc+ippc-1) + Derv(kppc2+ippc2-1);
902 tmpmat(ippc,ippc2) = signe * SolEm(Points(kppc+ippc-1) - Points(kppc2+ippc2-1),iu,iv);
903 }
904
905 for(Standard_Integer irow=1;irow <= myLXYZConstraints(i).Coeff().ColLength() ; irow++)
906 for(Standard_Integer irow2=1;irow2 <= myLXYZConstraints(i2).Coeff().ColLength() ; irow2++)
907 for(Standard_Integer icol=1;icol <= myLXYZConstraints(i).Coeff().RowLength() ; icol++)
908 for(Standard_Integer icol2=1;icol2 <= myLXYZConstraints(i2).Coeff().RowLength() ; icol2++)
909 mat(i0+k+irow-1,j0+k2+irow2-1) +=
910 myLXYZConstraints(i).Coeff()(irow,icol)*myLXYZConstraints(i2).Coeff()(irow2,icol2)*tmpmat(icol,icol2);
911
912
913 k2 += myLXYZConstraints(i2).Coeff().ColLength();
914 kppc2 += myLXYZConstraints(i2).Coeff().RowLength();
915 }
916
917 k += myLXYZConstraints(i).Coeff().ColLength();
918 kppc += myLXYZConstraints(i).Coeff().RowLength();
919 }
920
921
922
923
924 i = ncc1+ncc2;
925 for(Standard_Integer iu = 0; iu< order; iu++)
926 for(Standard_Integer iv =0; iu+iv < order; iv++) {
927 for(Standard_Integer a_j=0; a_j < ncc1; a_j++) {
928 Standard_Integer idu = Deru(a_j);
929 Standard_Integer idv = Derv(a_j);
930 mat(i0+i,j0+a_j) = Polm (Points(a_j), iu, iv, idu, idv);
931 }
932
933 Standard_Integer k2 = ncc1;
934 Standard_Integer kppc2 = ncc1;
935 for(Standard_Integer i2 = 1; i2<= myLXYZConstraints.Length(); i2++){
936 math_Vector vmat(1,myLXYZConstraints(i2).GetPPC().Length());
937 for(Standard_Integer ippc2=1;ippc2 <= myLXYZConstraints(i2).GetPPC().Length() ; ippc2++){
938 Standard_Integer idu = Deru(kppc2+ippc2-1);
939 Standard_Integer idv = Derv(kppc2+ippc2-1);
940 vmat(ippc2) = Polm (Points(kppc2+ippc2-1),iu,iv,idu,idv);
941 }
942
943 for(Standard_Integer irow2=1;irow2 <= myLXYZConstraints(i2).Coeff().ColLength() ; irow2++)
944 for(Standard_Integer icol2=1;icol2 <= myLXYZConstraints(i2).Coeff().RowLength() ; icol2++)
945 mat(i0+i,j0+k2+irow2-1) += myLXYZConstraints(i2).Coeff()(irow2,icol2)*vmat(icol2);
946
947 k2 += myLXYZConstraints(i2).Coeff().ColLength();
948 kppc2 += myLXYZConstraints(i2).Coeff().RowLength();
949 }
950
951 i++;
952 }
953
954 Standard_Integer n_dimat = ncc1 + ncc2 + order*(order+1)/2;
955
956 for(i=0;i<n_dimat;i++) {
957 for(Standard_Integer a_j = i+1; a_j < n_dimat; a_j++) {
958 mat(i0+i,j0+a_j) = mat(i0+a_j,j0+i);
959 }
960 }
961
962 }
963
964 //=======================================================================
965 //function : IsDone
966 //purpose :
967 //=======================================================================
968
IsDone() const969 Standard_Boolean Plate_Plate::IsDone() const
970 {
971 return OK;
972 }
973
974
975 //=======================================================================
976 //function : destroy
977 //purpose :
978 //=======================================================================
979
destroy()980 void Plate_Plate::destroy()
981 {
982 Init();
983 }
984
985 //=======================================================================
986 //function : Init
987 //purpose :
988 //=======================================================================
989
Init()990 void Plate_Plate::Init()
991 {
992 myConstraints.Clear();
993 myLXYZConstraints.Clear();
994 myLScalarConstraints.Clear();
995
996
997 delete [] (gp_XYZ*)solution;
998 solution = 0;
999
1000 delete [] (gp_XY*)points;
1001 points = 0;
1002
1003 delete [] (Standard_Integer*)deru;
1004 deru = 0;
1005
1006 delete [] (Standard_Integer*)derv;
1007 derv = 0;
1008
1009 order = 0;
1010 n_el = 0;
1011 n_dim = 0;
1012 OK = Standard_True;
1013 maxConstraintOrder=0;
1014 }
1015
1016 //=======================================================================
1017 //function : Evaluate
1018 //purpose :
1019 //=======================================================================
1020
Evaluate(const gp_XY & point2d) const1021 gp_XYZ Plate_Plate::Evaluate(const gp_XY& point2d) const
1022 {
1023 if(solution == 0) return gp_XYZ(0,0,0);
1024 if(!OK) return gp_XYZ(0,0,0);
1025
1026 gp_XYZ valeur(0,0,0);
1027
1028 if(!PolynomialPartOnly)
1029 {
1030 for(Standard_Integer i=0; i<n_el;i++)
1031 {
1032 Standard_Real signe = 1;
1033 if ( ((Deru(i)+Derv(i))%2) == 1) signe = -1;
1034 valeur += Solution(i) * (signe*SolEm(point2d - Points(i), Deru(i), Derv(i))) ;
1035 }
1036 }
1037 Standard_Integer i = n_el;
1038 for(Standard_Integer idu = 0; idu< order; idu++)
1039 for(Standard_Integer idv =0; idu+idv < order; idv++)
1040 {
1041 valeur += Solution(i) * Polm( point2d, idu,idv,0,0);
1042 i++;
1043 }
1044 return valeur;
1045 }
1046
1047 //=======================================================================
1048 //function : EvaluateDerivative
1049 //purpose :
1050 //=======================================================================
1051
EvaluateDerivative(const gp_XY & point2d,const Standard_Integer iu,const Standard_Integer iv) const1052 gp_XYZ Plate_Plate::EvaluateDerivative(const gp_XY& point2d, const Standard_Integer iu, const Standard_Integer iv) const
1053 {
1054 if(solution == 0) return gp_XYZ(0,0,0);
1055 if(!OK) return gp_XYZ(0,0,0);
1056
1057 gp_XYZ valeur(0,0,0);
1058 if(!PolynomialPartOnly)
1059 {
1060 for(Standard_Integer i=0; i<n_el;i++)
1061 {
1062 Standard_Real signe = 1;
1063 if ( ((Deru(i)+Derv(i))%2) == 1) signe = -1;
1064 valeur += Solution(i) * (signe*SolEm(point2d - Points(i), Deru(i)+iu, Derv(i)+iv)) ;
1065 }
1066 }
1067 Standard_Integer i = n_el;
1068 for(Standard_Integer idu = 0; idu< order; idu++)
1069 for(Standard_Integer idv =0; idu+idv < order; idv++)
1070 {
1071 valeur += Solution(i) * Polm( point2d, idu,idv, iu,iv);
1072 i++;
1073 }
1074 return valeur;
1075 }
1076 //=======================================================================
1077 //function : Plate_Plate::CoefPol
1078 //purpose : give back the array of power basis coefficient of
1079 // the polynomial part of the Plate function
1080 //=======================================================================
1081
CoefPol(Handle (TColgp_HArray2OfXYZ)& Coefs) const1082 void Plate_Plate::CoefPol(Handle(TColgp_HArray2OfXYZ)& Coefs) const
1083 {
1084 Coefs = new TColgp_HArray2OfXYZ(0,order-1,0,order-1,gp_XYZ(0.,0.,0.));
1085 Standard_Integer i = n_el;
1086 for(Standard_Integer iu = 0; iu< order; iu++)
1087 for(Standard_Integer iv =0; iu+iv < order; iv++)
1088 {
1089 Coefs->ChangeValue(iu,iv) = Solution(i)*ddu[iu]*ddv[iv];
1090 //Coefs->ChangeValue(idu,idv) = Solution(i);
1091 // it is necessary to reset this line if one remove factors in method Polm.
1092 i++;
1093 }
1094
1095 }
1096 //=======================================================================
1097 //function : Plate_Plate::Continuity
1098 //purpose : give back the continuity order of the Plate function
1099 //=======================================================================
1100
Continuity() const1101 Standard_Integer Plate_Plate::Continuity() const
1102 {
1103 return 2*order - 3 - maxConstraintOrder;
1104 }
1105
1106 //=======================================================================
1107 //function : Plate_Plate::SolEm
1108 //purpose : compute the (iu,iv)th derivative of the fundamental solution
1109 // of Laplcian at the power order
1110 //=======================================================================
1111
1112
SolEm(const gp_XY & point2d,const Standard_Integer iu,const Standard_Integer iv) const1113 Standard_Real Plate_Plate::SolEm(const gp_XY& point2d, const Standard_Integer iu, const Standard_Integer iv) const
1114 {
1115 Plate_Plate* aThis = const_cast<Plate_Plate*>(this);
1116 Standard_Real U,V;
1117 Standard_Integer IU,IV;
1118
1119 if(iv>iu)
1120 {
1121 // SolEm is symmetric in (u<->v) : we swap u and v if iv>iu
1122 // to avoid some code
1123 IU = iv;
1124 IV = iu;
1125 U = point2d.Y() *ddv[1];
1126 V = point2d.X() *ddu[1];
1127 }
1128 else
1129 {
1130 IU = iu;
1131 IV = iv;
1132 U = point2d.X() *ddu[1];
1133 V = point2d.Y() *ddv[1];
1134 }
1135
1136 if((U==Uold)&&(V==Vold) )
1137 {
1138 if (R<1.e-20) return 0;
1139 }
1140 else
1141 {
1142 aThis->Uold = U;
1143 aThis->Vold = V;
1144 aThis->U2 = U*U;
1145 aThis->R = U2+V*V;
1146 if (R<1.e-20) return 0;
1147 aThis->L = log(R);
1148 }
1149 Standard_Real DUV = 0;
1150
1151 Standard_Integer m = order;
1152 Standard_Integer mm1 = m-1;
1153 Standard_Real &r = aThis->R;
1154
1155
1156 //Standard_Real pr = pow(R, mm1 - IU - IV);
1157 // this expression takes a lot of time
1158 //(does not take into account a small integer value of the exponent)
1159 //
1160
1161 Standard_Integer expo = mm1 - IU - IV;
1162 Standard_Real pr;
1163 if(expo<0)
1164 {
1165 pr = R;
1166 for(Standard_Integer i=1;i<-expo;i++) pr *= R;
1167 pr = 1./pr;
1168 }
1169 else if(expo>0)
1170 {
1171 pr = R;
1172 for(Standard_Integer i=1;i<expo;i++) pr *= R;
1173 }
1174 else pr = 1.;
1175
1176
1177 switch (IU)
1178 {
1179 case 0:
1180 switch (IV)
1181 {
1182 case 0:
1183 {
1184 DUV = pr*L;
1185 }
1186 break;
1187
1188 default:
1189 break;
1190 }
1191 break;
1192
1193 case 1:
1194 switch (IV)
1195 {
1196 case 0:
1197 {
1198 DUV = 2*pr*U*(1+L*mm1);
1199 }
1200 break;
1201
1202 case 1:
1203 {
1204 Standard_Real m2 = m*m;
1205 //DUV = 4*pr*U*V*(-3+2*L+2*m-3*L*m+L*m2);
1206 DUV = 4*pr*U*V*((2*m-3)+(m2-3*m+2)*L);
1207 }
1208 break;
1209
1210 default:
1211 break;
1212 }
1213 break;
1214
1215 case 2:
1216 switch (IV)
1217 {
1218 case 0:
1219 {
1220 Standard_Real m2 = m*m;
1221 DUV = 2*pr*(R-L*R+L*m*R-6*U2+4*L*U2+4*m*U2-6*L*m*U2+2*L*m2*U2);
1222 }
1223 break;
1224
1225 case 1:
1226 {
1227 Standard_Real m2 = m*m;
1228 Standard_Real m3 = m2*m;
1229 DUV = -3*R+2*L*R+2*m*R-3*L*m*R+L*m2*R+22*U2-12*L*U2-24*m*U2+22*L*m*U2+6*m2*U2-12*L*m2*U2+2*L*m3*U2;
1230 DUV = DUV * 4* pr*V;
1231 }
1232 break;
1233
1234 case 2:
1235 {
1236 Standard_Real m2 = m*m;
1237 Standard_Real m3 = m2*m;
1238 Standard_Real m4 = m2*m2;
1239 Standard_Real V2 = V*V;
1240 Standard_Real R2 = R*R;
1241 DUV = -3*R2+2*L*R2+2*m*R2-3*L*m*R2+L*m2*R2+22*R*U2-12*L*R*U2-24*m*R*U2+22*L*m*R*U2+6*m2*R*U2-12*L*m2*R*U2;
1242 DUV += 2*L*m3*R*U2+22*R*V2-12*L*R*V2-24*m*R*V2+22*L*m*R*V2+6*m2*R*V2-12*L*m2*R*V2+2*L*m3*R*V2-200*U2*V2+96*L*U2*V2;
1243 DUV += 280*m*U2*V2-200*L*m*U2*V2-120*m2*U2*V2+140*L*m2*U2*V2+16*m3*U2*V2-40*L*m3*U2*V2+4*L*m4*U2*V2;
1244 DUV = 4*pr*DUV;
1245 }
1246 break;
1247
1248 default:
1249 break;
1250 }
1251 break;
1252
1253 case 3:
1254 switch (IV)
1255 {
1256 case 0:
1257 {
1258 Standard_Real m2 = m*m;
1259 Standard_Real m3 = m2*m;
1260 DUV = -9*R+6*L*R+6*m*R-9*L*m*R+3*L*m2*R+22*U2-12*L*U2-24*m*U2+22*L*m*U2+6*m2*U2-12*L*m2*U2+2*L*m3*U2;
1261 DUV = DUV * 4* pr*U;
1262 }
1263 break;
1264
1265 case 1:
1266 {
1267 Standard_Real m2 = m*m;
1268 Standard_Real m3 = m2*m;
1269 Standard_Real m4 = m2*m2;
1270 DUV = 33*R-18*L*R-36*m*R+33*L*m*R+9*m2*R-18*L*m2*R+3*L*m3*R-100*U2+48*L*U2+140*m*U2-100*L*m*U2-60*m2*U2+70*L*m2*U2;
1271 DUV += 8*m3*U2-20*L*m3*U2+2*L*m4*U2;
1272 DUV = 8*pr*U*V*DUV;
1273 }
1274 break;
1275
1276 case 2:
1277 {
1278 Standard_Real m2 = m*m;
1279 Standard_Real m3 = m2*m;
1280 Standard_Real m4 = m2*m2;
1281 Standard_Real m5 = m4*m;
1282 Standard_Real ru2 = R*U2;
1283 Standard_Real v2 = V*V;
1284 Standard_Real rv2 = R*v2;
1285 Standard_Real u2v2 = v2*U2;
1286 Standard_Real r2 = r*r;
1287
1288 // copy-paste the mathematics
1289 DUV =
1290 -100*ru2 + 48*L*ru2 + 140*m*ru2 - 100*L*m*ru2 - 60*m2*ru2 + 70*L*m2*ru2 + 8*m3*ru2 -
1291 20*L*m3*ru2 + 2*L*m4*ru2 - 300*rv2 + 144*L*rv2 + 420*m*rv2 - 300*L*m*rv2 - 180*m2*rv2 + 210*L*m2*rv2 +
1292 24*m3*rv2 - 60*L*m3*rv2 + 6*L*m4*rv2 + 33*r2 - 18*L*r2 - 36*m*r2 + 33*L*m*r2 + 9*m2*r2 - 18*L*m2*r2 +
1293 3*L*m3*r2 + 1096*u2v2 - 480*L*u2v2 - 1800*m*u2v2 + 1096*L*m*u2v2 + 1020*m2*u2v2 - 900*L*m2*u2v2 -
1294 240*m3*u2v2 + 340*L*m3*u2v2 + 20*m4*u2v2 - 60*L*m4*u2v2 + 4*L*m5*u2v2;
1295
1296 DUV = 8*pr*U*DUV;
1297 }
1298 break;
1299
1300 case 3:
1301 {
1302 Standard_Real m2 = m*m;
1303 Standard_Real m3 = m2*m;
1304 Standard_Real m4 = m2*m2;
1305 Standard_Real m5 = m3*m2;
1306 Standard_Real m6 = m3*m3;
1307 Standard_Real ru2 = r*U2;
1308 Standard_Real v2 = V*V;
1309 Standard_Real rv2 = R*v2;
1310 Standard_Real u2v2 = v2*U2;
1311 Standard_Real r2 = r*r;
1312
1313 // copy-paste the mathematics
1314 DUV =
1315 1644*ru2 - 720*L*ru2 - 2700*m*ru2 + 1644*L*m*ru2 + 1530*m2*ru2 - 1350*L*m2*ru2 -
1316 360*m3*ru2 + 510*L*m3*ru2 + 30*m4*ru2 - 90*L*m4*ru2 + 6*L*m5*ru2 + 1644*rv2 - 720*L*rv2 - 2700*m*rv2 +
1317 1644*L*m*rv2 + 1530*m2*rv2 - 1350*L*m2*rv2 - 360*m3*rv2 + 510*L*m3*rv2 + 30*m4*rv2 - 90*L*m4*rv2 +
1318 6*L*m5*rv2 - 450*r2 + 216*L*r2 + 630*m*r2 - 450*L*m*r2 - 270*m2*r2 + 315*L*m2*r2 + 36*m3*r2 - 90*L*m3*r2 +
1319 9*L*m4*r2 - 7056*u2v2 + 2880*L*u2v2 + 12992*m*u2v2 - 7056*L*m*u2v2 - 8820*m2*u2v2 + 6496*L*m2*u2v2 +
1320 2800*m3*u2v2 - 2940*L*m3*u2v2 - 420*m4*u2v2 + 700*L*m4*u2v2 + 24*m5*u2v2 - 84*L*m5*u2v2 + 4*L*m6*u2v2;
1321
1322 DUV = 16*pr*U*V*DUV;
1323 }
1324 break;
1325
1326 default:
1327 break;
1328 }
1329 break;
1330
1331 case 4:
1332 switch (IV)
1333 {
1334 case 0:
1335 {
1336 Standard_Real m2 = m*m;
1337 Standard_Real m3 = m2*m;
1338 Standard_Real m4 = m2*m2;
1339 Standard_Real U4 = U2*U2;
1340 Standard_Real R2 = R*R;
1341 DUV = -9*R2+6*L*R2+6*m*R2-9*L*m*R2+3*L*m2*R2+132*R*U2-72*L*R*U2-144*m*R*U2+132*L*m*R*U2+36*m2*R*U2-72*L*m2*R*U2;
1342 DUV += 12*L*m3*R*U2-200*U4+96*L*U4+280*m*U4-200*L*m*U4-120*m2*U4+140*L*m2*U4+16*m3*U4-40*L*m3*U4+4*L*m4*U4;
1343 DUV = 4*pr*DUV;
1344 }
1345 break;
1346
1347 case 1:
1348 {
1349 Standard_Real m2 = m*m;
1350 Standard_Real m3 = m2*m;
1351 Standard_Real m4 = m2*m2;
1352 Standard_Real m5 = m2*m3;
1353 Standard_Real u4 = U2*U2;
1354 Standard_Real ru2 = R*U2;
1355 Standard_Real r2 = R*R;
1356
1357 // copy-paste the mathematics
1358 DUV =
1359 -600*ru2 + 288*L*ru2 + 840*m*ru2 - 600*L*m*ru2 - 360*m2*ru2 + 420*L*m2*ru2 + 48*m3*ru2 -
1360 120*L*m3*ru2 + 12*L*m4*ru2 + 33*r2 - 18*L*r2 - 36*m*r2 + 33*L*m*r2 + 9*m2*r2 - 18*L*m2*r2 + 3*L*m3*r2 +
1361 1096*u4 - 480*L*u4 - 1800*m*u4 + 1096*L*m*u4 + 1020*m2*u4 - 900*L*m2*u4 - 240*m3*u4 + 340*L*m3*u4 + 20*m4*u4 -
1362 60*L*m4*u4 + 4*L*m5*u4;
1363
1364 DUV = 8*pr*V*DUV;
1365 }
1366 break;
1367
1368 case 2:
1369 {
1370 Standard_Real m2 = m*m;
1371 Standard_Real m3 = m2*m;
1372 Standard_Real m4 = m2*m2;
1373 Standard_Real m5 = m2*m3;
1374 Standard_Real m6 = m3*m3;
1375 Standard_Real u4 = U2*U2;
1376 Standard_Real r2 = r*r;
1377 Standard_Real r3 = r2*r;
1378 Standard_Real v2 = V*V;
1379 Standard_Real u2v2 = v2*U2;
1380 Standard_Real ru2v2 = R*u2v2;
1381 Standard_Real u4v2 = u4*v2;
1382 Standard_Real r2u2 = r2*U2;
1383 Standard_Real ru4 = r*u4;
1384 Standard_Real r2v2 = r2*v2;
1385
1386 // copy-paste the mathematics
1387 DUV =
1388 6576*ru2v2 - 2880*L*ru2v2 - 10800*m*ru2v2 + 6576*L*m*ru2v2 + 6120*m2*ru2v2 - 5400*L*m2*ru2v2 -
1389 1440*m3*ru2v2 + 2040*L*m3*ru2v2 + 120*m4*ru2v2 - 360*L*m4*ru2v2 + 24*L*m5*ru2v2 + 1096*ru4 - 480*L*ru4 -
1390 1800*m*ru4 + 1096*L*m*ru4 + 1020*m2*ru4 - 900*L*m2*ru4 - 240*m3*ru4 + 340*L*m3*ru4 + 20*m4*ru4 - 60*L*m4*ru4 +
1391 4*L*m5*ru4 - 600*r2u2 + 288*L*r2u2 + 840*m*r2u2 - 600*L*m*r2u2 - 360*m2*r2u2 + 420*L*m2*r2u2 + 48*m3*r2u2 -
1392 120*L*m3*r2u2 + 12*L*m4*r2u2 - 300*r2v2 + 144*L*r2v2 + 420*m*r2v2 - 300*L*m*r2v2 - 180*m2*r2v2 + 210*L*m2*r2v2 +
1393 24*m3*r2v2 - 60*L*m3*r2v2 + 6*L*m4*r2v2 + 33*r3 - 18*L*r3 - 36*m*r3 + 33*L*m*r3 + 9*m2*r3 - 18*L*m2*r3 +
1394 3*L*m3*r3 - 14112*u4v2 + 5760*L*u4v2 + 25984*m*u4v2 - 14112*L*m*u4v2 - 17640*m2*u4v2 + 12992*L*m2*u4v2 +
1395 5600*m3*u4v2 - 5880*L*m3*u4v2 - 840*m4*u4v2 + 1400*L*m4*u4v2 + 48*m5*u4v2 - 168*L*m5*u4v2 + 8*L*m6*u4v2;
1396
1397 DUV = 8*pr*DUV;
1398 }
1399 break;
1400
1401 case 3:
1402 {
1403 Standard_Real m2 = m*m;
1404 Standard_Real m3 = m2*m;
1405 Standard_Real m4 = m2*m2;
1406 Standard_Real m5 = m2*m3;
1407 Standard_Real m6 = m3*m3;
1408 Standard_Real m7 = m3*m4;
1409 Standard_Real u4 = U2*U2;
1410 Standard_Real r2 = r*r;
1411 Standard_Real r3 = r2*r;
1412 Standard_Real v2 = V*V;
1413 Standard_Real u2v2 = v2*U2;
1414 Standard_Real ru2v2 = R*u2v2;
1415 Standard_Real u4v2 = u4*v2;
1416 Standard_Real r2u2 = r2*U2;
1417 Standard_Real r2v2 = r2*v2;
1418 Standard_Real ru4 = r*u4;
1419
1420 // copy-paste the mathematics
1421 DUV =
1422 -42336*ru2v2 + 17280*L*ru2v2 + 77952*m*ru2v2 - 42336*L*m*ru2v2 - 52920*m2*ru2v2 +
1423 38976*L*m2*ru2v2 + 16800*m3*ru2v2 - 17640*L*m3*ru2v2 - 2520*m4*ru2v2 + 4200*L*m4*ru2v2 + 144*m5*ru2v2 -
1424 504*L*m5*ru2v2 + 24*L*m6*ru2v2 - 21168*ru4 + 8640*L*ru4 + 38976*m*ru4 - 21168*L*m*ru4 - 26460*m2*ru4 +
1425 19488*L*m2*ru4 + 8400*m3*ru4 - 8820*L*m3*ru4 - 1260*m4*ru4 + 2100*L*m4*ru4 + 72*m5*ru4 - 252*L*m5*ru4 +
1426 12*L*m6*ru4 + 9864*r2u2 - 4320*L*r2u2 - 16200*m*r2u2 + 9864*L*m*r2u2 + 9180*m2*r2u2 - 8100*L*m2*r2u2 -
1427 2160*m3*r2u2 + 3060*L*m3*r2u2 + 180*m4*r2u2 - 540*L*m4*r2u2 + 36*L*m5*r2u2 + 1644*r2v2 - 720*L*r2v2 -
1428 2700*m*r2v2 + 1644*L*m*r2v2 + 1530*m2*r2v2 - 1350*L*m2*r2v2 - 360*m3*r2v2 + 510*L*m3*r2v2 + 30*m4*r2v2 -
1429 90*L*m4*r2v2 + 6*L*m5*r2v2 - 450*r3 + 216*L*r3 + 630*m*r3 - 450*L*m*r3 - 270*m2*r3 + 315*L*m2*r3 + 36*m3*r3 -
1430 90*L*m3*r3 + 9*L*m4*r3 + 104544*u4v2 - 40320*L*u4v2 - 210112*m*u4v2 + 104544*L*m*u4v2 + 162456*m2*u4v2 -
1431 105056*L*m2*u4v2 - 62720*m3*u4v2 + 54152*L*m3*u4v2 + 12880*m4*u4v2 - 15680*L*m4*u4v2 - 1344*m5*u4v2 +
1432 2576*L*m5*u4v2 + 56*m6*u4v2 - 224*L*m6*u4v2 + 8*L*m7*u4v2;
1433
1434 DUV = 16*pr*V*DUV;
1435 }
1436 break;
1437
1438 default:
1439 break;
1440 }
1441 break;
1442
1443 case 5:
1444 switch (IV)
1445 {
1446 case 0:
1447 {
1448 Standard_Real m2 = m*m;
1449 Standard_Real m3 = m2*m;
1450 Standard_Real m4 = m2*m2;
1451 Standard_Real m5 = m2*m3;
1452 Standard_Real u4 = U2*U2;
1453 Standard_Real r2 = R*R;
1454 Standard_Real ru2 = R*U2;
1455
1456 // copy-paste the mathematics
1457 DUV =
1458 -1000*ru2 + 480*L*ru2 + 1400*m*ru2 - 1000*L*m*ru2 - 600*m2*ru2 + 700*L*m2*ru2 + 80*m3*ru2 -
1459 200*L*m3*ru2 + 20*L*m4*ru2 + 165*r2 - 90*L*r2 - 180*m*r2 + 165*L*m*r2 + 45*m2*r2 - 90*L*m2*r2 + 15*L*m3*r2 +
1460 1096*u4 - 480*L*u4 - 1800*m*u4 + 1096*L*m*u4 + 1020*m2*u4 - 900*L*m2*u4 - 240*m3*u4 + 340*L*m3*u4 + 20*m4*u4 -
1461 60*L*m4*u4 + 4*L*m5*u4;
1462
1463 DUV = 8*pr*U*DUV;
1464 }
1465 break;
1466
1467 case 1:
1468 {
1469 Standard_Real m2 = m*m;
1470 Standard_Real m3 = m2*m;
1471 Standard_Real m4 = m2*m2;
1472 Standard_Real m5 = m2*m3;
1473 Standard_Real m6 = m3*m3;
1474 Standard_Real u4 = U2*U2;
1475 Standard_Real ru2 = r*U2;
1476 Standard_Real r2 = r*r;
1477
1478
1479 // copy-paste the mathematics
1480 DUV =
1481 5480*ru2 - 2400*L*ru2 - 9000*m*ru2 + 5480*L*m*ru2 + 5100*m2*ru2 - 4500*L*m2*ru2 - 1200*m3*ru2 +
1482 1700*L*m3*ru2 + 100*m4*ru2 - 300*L*m4*ru2 + 20*L*m5*ru2 - 750*r2 + 360*L*r2 + 1050*m*r2 - 750*L*m*r2 -
1483 450*m2*r2 + 525*L*m2*r2 + 60*m3*r2 - 150*L*m3*r2 + 15*L*m4*r2 - 7056*u4 + 2880*L*u4 + 12992*m*u4 - 7056*L*m*u4 -
1484 8820*m2*u4 + 6496*L*m2*u4 + 2800*m3*u4 - 2940*L*m3*u4 - 420*m4*u4 + 700*L*m4*u4 + 24*m5*u4 - 84*L*m5*u4 +
1485 4*L*m6*u4;
1486
1487 DUV = 16*pr*U*V*DUV;
1488 }
1489 break;
1490
1491 case 2:
1492 {
1493 Standard_Real m2 = m*m;
1494 Standard_Real m3 = m2*m;
1495 Standard_Real m4 = m2*m2;
1496 Standard_Real m5 = m2*m3;
1497 Standard_Real m6 = m3*m3;
1498 Standard_Real m7 = m3*m4;
1499 Standard_Real u4 = U2*U2;
1500 Standard_Real r2 = r*r;
1501 Standard_Real r3 = r2*r;
1502 Standard_Real v2 = V*V;
1503 Standard_Real u2v2 = v2*U2;
1504 Standard_Real ru2v2 = R*u2v2;
1505 Standard_Real u4v2 = u4*v2;
1506 Standard_Real r2u2 = r2*U2;
1507 Standard_Real r2v2 = r2*v2;
1508 Standard_Real ru4 = r*u4;
1509
1510 // copy-paste the mathematics
1511 DUV =
1512
1513 -70560*ru2v2 + 28800*L*ru2v2 + 129920*m*ru2v2 - 70560*L*m*ru2v2 - 88200*m2*ru2v2 +
1514 64960*L*m2*ru2v2 + 28000*m3*ru2v2 - 29400*L*m3*ru2v2 - 4200*m4*ru2v2 + 7000*L*m4*ru2v2 + 240*m5*ru2v2 -
1515 840*L*m5*ru2v2 + 40*L*m6*ru2v2 - 7056*ru4 + 2880*L*ru4 + 12992*m*ru4 - 7056*L*m*ru4 - 8820*m2*ru4 +
1516 6496*L*m2*ru4 + 2800*m3*ru4 - 2940*L*m3*ru4 - 420*m4*ru4 + 700*L*m4*ru4 + 24*m5*ru4 - 84*L*m5*ru4 + 4*L*m6*ru4 +
1517 5480*r2u2 - 2400*L*r2u2 - 9000*m*r2u2 + 5480*L*m*r2u2 + 5100*m2*r2u2 - 4500*L*m2*r2u2 - 1200*m3*r2u2 +
1518 1700*L*m3*r2u2 + 100*m4*r2u2 - 300*L*m4*r2u2 + 20*L*m5*r2u2 + 8220*r2v2 - 3600*L*r2v2 - 13500*m*r2v2 +
1519 8220*L*m*r2v2 + 7650*m2*r2v2 - 6750*L*m2*r2v2 - 1800*m3*r2v2 + 2550*L*m3*r2v2 + 150*m4*r2v2 - 450*L*m4*r2v2 +
1520 30*L*m5*r2v2 - 750*r3 + 360*L*r3 + 1050*m*r3 - 750*L*m*r3 - 450*m2*r3 + 525*L*m2*r3 + 60*m3*r3 - 150*L*m3*r3 +
1521 15*L*m4*r3 + 104544*u4v2 - 40320*L*u4v2 - 210112*m*u4v2 + 104544*L*m*u4v2 + 162456*m2*u4v2 - 105056*L*m2*u4v2 -
1522 62720*m3*u4v2 + 54152*L*m3*u4v2 + 12880*m4*u4v2 - 15680*L*m4*u4v2 - 1344*m5*u4v2 + 2576*L*m5*u4v2 + 56*m6*u4v2 -
1523 224*L*m6*u4v2 + 8*L*m7*u4v2;
1524
1525 DUV = 16*pr*U*DUV;
1526 }
1527 break;
1528
1529 default:
1530 break;
1531 }
1532 break;
1533
1534 case 6:
1535 switch (IV)
1536 {
1537 case 0:
1538 {
1539 Standard_Real m2 = m*m;
1540 Standard_Real m3 = m2*m;
1541 Standard_Real m4 = m2*m2;
1542 Standard_Real m5 = m2*m3;
1543 Standard_Real m6 = m3*m3;
1544 Standard_Real u4 = U2*U2;
1545 Standard_Real u6 = U2*u4;
1546 Standard_Real r2 = r*r;
1547 Standard_Real r3 = r2*r;
1548 Standard_Real r2u2 = r2*U2;
1549 Standard_Real ru4 = r*u4;
1550
1551 // copy-paste the mathematics
1552 DUV =
1553 16440*ru4 - 7200*L*ru4 - 27000*m*ru4 + 16440*L*m*ru4 + 15300*m2*ru4 - 13500*L*m2*ru4 -
1554 3600*m3*ru4 + 5100*L*m3*ru4 + 300*m4*ru4 - 900*L*m4*ru4 + 60*L*m5*ru4 - 4500*r2u2 + 2160*L*r2u2 + 6300*m*r2u2 -
1555 4500*L*m*r2u2 - 2700*m2*r2u2 + 3150*L*m2*r2u2 + 360*m3*r2u2 - 900*L*m3*r2u2 + 90*L*m4*r2u2 + 165*r3 - 90*L*r3 -
1556 180*m*r3 + 165*L*m*r3 + 45*m2*r3 - 90*L*m2*r3 + 15*L*m3*r3 - 14112*u6 + 5760*L*u6 + 25984*m*u6 - 14112*L*m*u6 -
1557 17640*m2*u6 + 12992*L*m2*u6 + 5600*m3*u6 - 5880*L*m3*u6 - 840*m4*u6 + 1400*L*m4*u6 + 48*m5*u6 - 168*L*m5*u6 +
1558 8*L*m6*u6;
1559
1560 DUV = 8*pr*DUV;
1561 }
1562 break;
1563
1564 default:
1565 break;
1566 }
1567 break;
1568
1569 default:
1570 break;
1571 }
1572
1573 return DUV * ddu[iu]*ddv[iv];
1574
1575 }
1576
1577
1578 //=======================================================================
1579 //function : UVBox
1580 //purpose :
1581 //=======================================================================
1582
UVBox(Standard_Real & UMin,Standard_Real & UMax,Standard_Real & VMin,Standard_Real & VMax) const1583 void Plate_Plate::UVBox(Standard_Real& UMin, Standard_Real& UMax,
1584 Standard_Real& VMin, Standard_Real& VMax) const
1585 {
1586 Standard_Integer i ;
1587 const Standard_Real Bmin = 1.e-3;
1588 UMin = myConstraints(1).Pnt2d().X();
1589 UMax = UMin;
1590 VMin = myConstraints(1).Pnt2d().Y();
1591 VMax = VMin;
1592
1593 for( i=2; i<=myConstraints.Length();i++)
1594 {
1595 Standard_Real x = myConstraints(i).Pnt2d().X();
1596 if(x<UMin) UMin = x;
1597 if(x>UMax) UMax = x;
1598 Standard_Real y = myConstraints(i).Pnt2d().Y();
1599 if(y<VMin) VMin = y;
1600 if(y>VMax) VMax = y;
1601 }
1602
1603 for(i=1; i<=myLXYZConstraints.Length();i++)
1604 for(Standard_Integer j=1;j<= myLXYZConstraints(i).GetPPC().Length(); j++)
1605 {
1606 Standard_Real x = myLXYZConstraints(i).GetPPC()(j).Pnt2d().X();
1607 if(x<UMin) UMin = x;
1608 if(x>UMax) UMax = x;
1609 Standard_Real y = myLXYZConstraints(i).GetPPC()(j).Pnt2d().Y();
1610 if(y<VMin) VMin = y;
1611 if(y>VMax) VMax = y;
1612 }
1613
1614 for(i=1; i<=myLScalarConstraints.Length();i++)
1615 for(Standard_Integer j=1;j<= myLScalarConstraints(i).GetPPC().Length(); j++)
1616 {
1617 Standard_Real x = myLScalarConstraints(i).GetPPC()(j).Pnt2d().X();
1618 if(x<UMin) UMin = x;
1619 if(x>UMax) UMax = x;
1620 Standard_Real y = myLScalarConstraints(i).GetPPC()(j).Pnt2d().Y();
1621 if(y<VMin) VMin = y;
1622 if(y>VMax) VMax = y;
1623 }
1624
1625
1626 if(UMax-UMin < Bmin)
1627 {
1628 Standard_Real UM = 0.5*(UMin+UMax);
1629 UMin = UM - 0.5*Bmin;
1630 UMax = UM + 0.5*Bmin;
1631 }
1632 if(VMax-VMin < Bmin)
1633 {
1634 Standard_Real VM = 0.5*(VMin+VMax);
1635 VMin = VM - 0.5*Bmin;
1636 VMax = VM + 0.5*Bmin;
1637 }
1638 }
1639
1640 //=======================================================================
1641 //function : UVConstraints
1642 //purpose :
1643 //=======================================================================
1644
UVConstraints(TColgp_SequenceOfXY & Seq) const1645 void Plate_Plate::UVConstraints(TColgp_SequenceOfXY& Seq) const
1646 {
1647 for (Standard_Integer i=1;i<=myConstraints.Length();i++) {
1648 if ((myConstraints.Value(i).Idu()==0) && (myConstraints.Value(i).Idv()==0))
1649 Seq.Append((myConstraints.Value(i)).Pnt2d());
1650 }
1651 }
1652 //=======================================================================
1653
SetPolynomialPartOnly(const Standard_Boolean PPOnly)1654 void Plate_Plate::SetPolynomialPartOnly(const Standard_Boolean PPOnly)
1655 {
1656 PolynomialPartOnly = PPOnly;
1657 }
1658