src/Plate/Plate_Plate.cxx

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