src/Plate/Plate_Plate.cxx

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