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