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