src/AppParCurves/AppParCurves_ResolConstraint.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   3 //
   4 // The content of this file is subject to the Open CASCADE Technology Public
   5 // License Version 6.5 (the "License"). You may not use the content of this file
   6 // except in compliance with the License. Please obtain a copy of the License
   7 // at http://www.opencascade.org and read it completely before using this file.
   8 //
   9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  11 //
  12 // The Original Code and all software distributed under the License is
  13 // distributed on an "AS IS" basis, without warranty of any kind, and the
  14 // Initial Developer hereby disclaims all such warranties, including without
  15 // limitation, any warranties of merchantability, fitness for a particular
  16 // purpose or non-infringement. Please see the License for the specific terms
  17 // and conditions governing the rights and limitations under the License.
  18
  19 // lpa, le 11/09/91
  20
  21
  22 // Approximation d un ensemble de points contraints (MultiLine) avec une
  23 // solution approchee (MultiCurve). L algorithme utilise est l algorithme
  24 // d Uzawa du package mathematique.
  25
  26 #define No_Standard_RangeError
  27 #define No_Standard_OutOfRange
  28
  29 #include <math_Vector.hxx>
  30 #include <math_Matrix.hxx>
  31 #include <AppParCurves_Constraint.hxx>
  32 #include <AppParCurves_ConstraintCouple.hxx>
  33 #include <AppParCurves_MultiPoint.hxx>
  34 #include <AppParCurves.hxx>
  35 #include <Standard_DimensionError.hxx>
  36 #include <math_Uzawa.hxx>
  37 #include <TColStd_Array1OfInteger.hxx>
  38 #include <TColStd_Array2OfInteger.hxx>
  39 #include <gp_Pnt.hxx>
  40 #include <gp_Pnt2d.hxx>
  41 #include <gp_Vec.hxx>
  42 #include <gp_Vec2d.hxx>
  43 #include <TColgp_Array1OfPnt.hxx>
  44 #include <TColgp_Array1OfPnt2d.hxx>
  45 #include <TColgp_Array1OfVec.hxx>
  46 #include <TColgp_Array1OfVec2d.hxx>
  47
  48
  49 AppParCurves_ResolConstraint::
  50   AppParCurves_ResolConstraint(const MultiLine& SSP,
  51           AppParCurves_MultiCurve& SCurv,
  52           const Standard_Integer FirstPoint,
  53           const Standard_Integer LastPoint,
  54           const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
  55           const math_Matrix& Bern,
  56           const math_Matrix& DerivativeBern,
  57           const Standard_Real Tolerance):
  58           Cont(1, NbConstraints(SSP, FirstPoint,
  59                       LastPoint, TheConstraints),
  60                         1, NbColumns(SSP, SCurv.NbPoles()-1), 0.0),
  61           DeCont(1, NbConstraints(SSP, FirstPoint,
  62                       LastPoint, TheConstraints),
  63                         1, NbColumns(SSP, SCurv.NbPoles()-1), 0.0),
  64           Secont(1, NbConstraints(SSP, FirstPoint,
  65                            LastPoint, TheConstraints), 0.0),
  66           CTCinv(1, NbConstraints(SSP, FirstPoint,
  67                            LastPoint, TheConstraints),
  68                           1, NbConstraints(SSP, FirstPoint,
  69                                    LastPoint, TheConstraints)),
  70           Vardua(1, NbConstraints(SSP, FirstPoint,
  71                            LastPoint, TheConstraints)),
  72           IPas(1, LastPoint-FirstPoint+1),
  73           ITan(1, LastPoint-FirstPoint+1),
  74           ICurv(1, LastPoint-FirstPoint+1)
  75  {
  76   Standard_Integer i, j, k, NbCu= SCurv.NbCurves();
  77   Standard_Integer Npt;
  78   Standard_Integer Inc3, IncSec, IncCol, IP = 0, CCol;
  79   Standard_Integer myindex, Def = SCurv.NbPoles()-1;
  80   Standard_Integer nb3d, nb2d, Npol= Def+1, Npol2 = 2*Npol;
  81   Standard_Real T1 = 0., T2 = 0., T3, Tmax, Daij;
  82   Standard_Boolean Ok;
  83   gp_Vec V;
  84   gp_Vec2d V2d;
  85   gp_Pnt Poi;
  86   gp_Pnt2d Poi2d;
  87   AppParCurves_ConstraintCouple mycouple;
  88   AppParCurves_Constraint FC = AppParCurves_NoConstraint,
  89                           LC = AppParCurves_NoConstraint, Cons;
  90
  91
  92
  93   // Boucle de calcul du nombre de points de passage afin de dimensionner
  94   // les matrices.
  95   IncPass = 0;
  96   IncTan = 0;
  97   IncCurv = 0;
  98   for (i = TheConstraints->Lower(); i <= TheConstraints->Upper(); i++) {
  99     mycouple = TheConstraints->Value(i);
 100     myindex = mycouple.Index();
 101     Cons = mycouple.Constraint();
 102     if (myindex == FirstPoint) FC = Cons;
 103     if (myindex == LastPoint)  LC = Cons;
 104     if (Cons >= 1) {
 105       IncPass++;                      // IncPass = nbre de points de passage.
 106       IPas(IncPass) = myindex;
 107     }
 108     if (Cons >= 2) {
 109       IncTan++;                       // IncTan= nbre de points de tangence.
 110       ITan(IncTan) = myindex;
 111     }
 112     if (Cons == 3) {
 113       IncCurv++;                      // IncCurv = nbre de pts de courbure.
 114       ICurv(IncCurv) = myindex;
 115     }
 116   }
 117
 118
 119   if (IncPass == 0) {
 120     Done = Standard_True;
 121     return;
 122   }
 123
 124   nb3d = ToolLine::NbP3d(SSP);
 125   nb2d = ToolLine::NbP2d(SSP);
 126   Standard_Integer mynb3d=nb3d, mynb2d=nb2d;
 127   if (nb3d == 0) mynb3d = 1;
 128   if (nb2d == 0) mynb2d = 1;
 129   CCol = nb3d* 3 + nb2d* 2;
 130
 131
 132   // Declaration et initialisation des matrices et vecteurs de contraintes:
 133   math_Matrix ContInit(1, IncPass, 1, Npol);
 134   math_Vector Start(1, CCol*Npol);
 135   TColStd_Array2OfInteger Ibont(1, NbCu, 1, IncTan);
 136
 137
 138   // Remplissage de Cont pour les points de passage:
 139   // =================================================
 140   for (i =1; i <= IncPass; i++) {   // Cette partie ne depend que de Bernstein
 141     Npt = IPas(i);
 142     for (j = 1; j <= Npol; j++) {
 143       ContInit(i, j) = Bern(Npt, j);
 144     }
 145   }
 146   for (i = 1; i <= CCol; i++) {
 147     Cont.Set(IncPass*(i-1)+1, IncPass*i, Npol*(i-1)+1, Npol*i, ContInit);
 148   }
 149
 150
 151 // recuperation des vecteurs de depart pour Uzawa. Ce vecteur represente les
 152 // poles de SCurv.
 153 // Remplissage de secont et resolution.
 154
 155   TColgp_Array1OfVec tabV(1, mynb3d);
 156   TColgp_Array1OfVec2d tabV2d(1, mynb2d);
 157   TColgp_Array1OfPnt tabP(1, mynb3d);
 158   TColgp_Array1OfPnt2d tabP2d(1, mynb2d);
 159
 160   Inc3 = CCol*IncPass +1;
 161   IncCol = 0;
 162   IncSec = 0;
 163   for (k = 1; k <= NbCu; k++) {
 164     if (k <= nb3d) {
 165       for (i = 1; i <= IncTan; i++) {
 166         Npt = ITan(i);
 167         // choix du maximum de tangence pour exprimer la colinearite:
 168         Ok = ToolLine::Tangency(SSP, Npt, tabV);
 169         V = tabV(k);
 170         V.Coord(T1, T2, T3);
 171         Tmax = Abs(T1);
 172         Ibont(k, i) = 1;
 173         if (Abs(T2) > Tmax) {
 174           Tmax = Abs(T2);
 175           Ibont(k, i) = 2;
 176         }
 177         if (Abs(T3) > Tmax) {
 178           Tmax = Abs(T3);
 179           Ibont(k, i) = 3;
 180         }
 181         if (Ibont(k, i) == 3) {
 182           for (j = 1; j <= Npol; j++) {
 183             Daij = DerivativeBern(Npt, j);
 184             Cont(Inc3, j+ Npol + IncCol) = Daij*T3/Tmax;
 185             Cont(Inc3, j + Npol2 + IncCol) = -Daij *T2/Tmax;
 186             Cont(Inc3+1, j+ IncCol) = Daij* T3/Tmax;
 187             Cont(Inc3+1, j+ Npol2 + IncCol) = -Daij*T1/Tmax;
 188           }
 189         }
 190         else if (Ibont(k, i) == 1) {
 191           for (j = 1; j <= Npol; j++) {
 192             Daij = DerivativeBern(Npt, j);
 193             Cont(Inc3, j + IncCol) = Daij*T3/Tmax;
 194             Cont(Inc3, j + Npol2 + IncCol) = -Daij *T1/Tmax;
 195             Cont(Inc3+1, j+ IncCol) = Daij* T2/Tmax;
 196             Cont(Inc3+1, j+ Npol +IncCol) = -Daij*T1/Tmax;
 197           }
 198         }
 199         else if (Ibont(k, i) == 2) {
 200           for (j = 1; j <= Def+1; j++) {
 201             Daij = DerivativeBern(Npt, j);
 202             Cont(Inc3, j+ Npol + IncCol) = Daij*T3/Tmax;
 203             Cont(Inc3, j + Npol2 + IncCol) = -Daij *T2/Tmax;
 204             Cont(Inc3+1, j+ IncCol) = Daij* T2/Tmax;
 205             Cont(Inc3+1, j+ Npol + IncCol) = -Daij*T1/Tmax;
 206           }
 207         }
 208         Inc3 = Inc3 + 2;
 209       }
 210
 211       // Remplissage du second membre:
 212       for (i = 1; i <= IncPass; i++) {
 213         ToolLine::Value(SSP, IPas(i), tabP);
 214         Poi = tabP(k);
 215         Secont(i + IncSec) = Poi.X();
 216         Secont(i + IncPass + IncSec) = Poi.Y();
 217         Secont(i + 2*IncPass + IncSec) = Poi.Z();
 218       }
 219       IncSec = IncSec + 3*IncPass;
 220
 221       // Vecteur de depart:
 222       for (j =1; j <= Npol; j++) {
 223         Poi = SCurv.Value(j).Point(k);
 224         Start(j + IncCol) = Poi.X();
 225         Start(j+ Npol + IncCol) = Poi.Y();
 226         Start(j + Npol2 + IncCol) = Poi.Z();
 227       }
 228       IncCol = IncCol + 3*Npol;
 229     }
 230
 231
 232     else {
 233       for (i = 1; i <= IncTan; i++) {
 234         Npt = ITan(i);
 235         Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 236         V2d = tabV2d(k-nb3d);
 237         T1 = V2d.X();
 238         T2 = V2d.Y();
 239         Tmax = Abs(T1);
 240         Ibont(k, i) = 1;
 241         if (Abs(T2) > Tmax) {
 242           Tmax = Abs(T2);
 243           Ibont(k, i) = 2;
 244         }
 245         for (j = 1; j <= Npol; j++) {
 246           Daij = DerivativeBern(Npt, j);
 247           Cont(Inc3, j + IncCol) = Daij*T2;
 248           Cont(Inc3, j + Npol + IncCol) = -Daij*T1;
 249         }
 250         Inc3 = Inc3 +1;
 251       }
 252
 253       // Remplissage du second membre:
 254       for (i = 1; i <= IncPass; i++) {
 255         ToolLine::Value(SSP, IPas(i), tabP2d);
 256         Poi2d = tabP2d(i-nb3d);
 257         Secont(i + IncSec) = Poi2d.X();
 258         Secont(i + IncPass + IncSec) = Poi2d.Y();
 259       }
 260       IncSec = IncSec + 2*IncPass;
 261
 262       // Remplissage du vecteur de depart:
 263       for (j = 1; j <= Npol; j++) {
 264         Poi2d = SCurv.Value(j).Point2d(k);
 265         Start(j + IncCol) = Poi2d.X();
 266         Start(j + Npol + IncCol) = Poi2d.Y();
 267       }
 268       IncCol = IncCol + Npol2;
 269     }
 270   }
 271
 272   // Equations exprimant le meme rapport de tangence sur chaque courbe:
 273   // On prend les coordonnees les plus significatives.
 274
 275   --Inc3;
 276   for (i = 1; i <= IncTan; ++i) {
 277     IncCol = 0;
 278     Npt = ITan(i);
 279     for (k = 1; k <= NbCu-1; ++k) {
 280       ++Inc3;
 281 // Initialize first relation variable (T1)
 282       Standard_Integer addIndex_1 = 0, aVal = Ibont(k, i);
 283       switch (aVal)
 284       {
 285         case 1: //T1 ~ T1x
 286         {
 287           if (k <= nb3d) {
 288             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 289             V = tabV(k);
 290             T1 = V.X();
 291             IP = 3 * Npol;
 292           }
 293           else {
 294             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 295             V2d = tabV2d(k-nb3d);
 296             T1 = V2d.X();
 297             IP = 2 * Npol;
 298           }
 299           addIndex_1 = 0;
 300           break;
 301         }
 302         case 2: //T1 ~ T1y
 303         {
 304           if (k <= nb3d) {
 305             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 306             V = tabV(k);
 307             IP = 3 * Npol;
 308           }
 309           else {
 310             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 311             V2d = tabV2d(k-nb3d);
 312             T1 = V2d.Y();
 313             IP = 2 * Npol;
 314           }
 315           addIndex_1 = Npol;
 316           break;
 317         }
 318         case 3: //T1 ~ T1z
 319         {
 320           Ok = ToolLine::Tangency(SSP, Npt, tabV);
 321           V = tabV(k);
 322           T1 = V.Z();
 323           IP = 3 * Npol;
 324           addIndex_1 = 2 * Npol;
 325           break;
 326         }
 327       }
 328       // Initialize second relation variable (T2)
 329       Standard_Integer addIndex_2 = 0, aNextVal = Ibont(k + 1, i);
 330       switch (aNextVal)
 331       {
 332         case 1: //T2 ~ T2x
 333         {
 334           if ((k+1) <= nb3d) {
 335             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 336             V = tabV(k+1);
 337             T2 = V.X();
 338           }
 339           else {
 340             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 341             V2d = tabV2d(k+1-nb3d);
 342             T2 = V2d.X();
 343           }
 344           addIndex_2 = 0;
 345           break;
 346         }
 347         case 2: //T2 ~ T2y
 348         {
 349           if ((k+1) <= nb3d) {
 350             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 351             V = tabV(k+1);
 352             T2 = V.Y();
 353           }
 354           else {
 355             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 356             V2d = tabV2d(k+1-nb3d);
 357             T2 = V2d.Y();
 358           }
 359           addIndex_2 = Npol;
 360           break;
 361         }
 362         case 3: //T2 ~ T2z
 363         {
 364           Ok = ToolLine::Tangency(SSP, Npt, tabV);
 365           V = tabV(k+1);
 366           T2 = V.Z();
 367           addIndex_2 = 2 * Npol;
 368           break;
 369         }
 370       }
 371
 372       // Relations between T1 and T2:
 373       for (j = 1; j <=  Npol; j++) {
 374         Daij = DerivativeBern(Npt, j);
 375         Cont(Inc3, j + IncCol + addIndex_1) = Daij*T2;
 376         Cont(Inc3, j + IP + IncCol + addIndex_2) = -Daij*T1;
 377       }
 378       IncCol += IP;
 379     }
 380   }
 381
 382
 383
 384   // Equations concernant la courbure:
 385
 386
 387 /*  Inc3 = Inc3 +1;
 388   IncCol = 0;
 389   for (k = 1; k <= NbCu; k++) {
 390     for (i = 1; i <= IncCurv; i++) {
 391       Npt = ICurv(i);
 392       AppDef_MultiPointConstraint MP = SSP.Value(Npt);
 393       DDA = SecondDerivativeBern(Parameters(Npt));
 394       if (SSP.Value(1).Dimension(k) == 3) {
 395         C1 = MP.Curv(k).X();
 396         C2 = MP.Curv(k).Y();
 397         C3 = MP.Curv(k).Z();
 398         for (j = 1; j <= Npol; j++) {
 399           Daij = DerivativeBern(Npt, j);
 400           D2Aij = DDA(j);
 401           Cont(Inc3, j + IncCol) = D2Aij;
 402           Cont(Inc3, j + Npol2 + IncCol) = -C2*Daij;
 403           Cont(Inc3, j + Npol + IncCol) = C3*Daij;
 404
 405           Cont(Inc3+1, j + Npol + IncCol) = D2Aij;
 406           Cont(Inc3+1, j +IncCol) = -C3*Daij;
 407           Cont(Inc3+1, j + Npol2 + IncCol) = C1*Daij;
 408
 409           Cont(Inc3+2, j + Npol2+IncCol) = D2Aij;
 410           Cont(Inc3+2, j + Npol+IncCol) = -C1*Daij;
 411           Cont(Inc3+2, j + IncCol) = C2*Daij;
 412         }
 413         Inc3 = Inc3 + 3;
 414       }
 415       else {        // Dimension 2:
 416         C1 = MP.Curv2d(k).X();
 417         C2 = MP.Curv2d(k).Y();
 418         for (j = 1; j <= Npol; j++) {
 419           Daij = DerivativeBern(Npt, j);
 420           D2Aij = DDA(j);
 421           Cont(Inc3, j + IncCol) = D2Aij*C1;
 422           Cont(Inc3+1, j + Npol + IncCol) = D2Aij*C2;
 423         }
 424         Inc3 = Inc3 + 2;
 425       }
 426     }
 427   }
 428
 429 */
 430   // Resolution par Uzawa:
 431
 432   math_Uzawa UzaResol(Cont, Secont, Start, Tolerance);
 433   if (!(UzaResol.IsDone())) {
 434     Done = Standard_False;
 435     return;
 436   }
 437   CTCinv = UzaResol.InverseCont();
 438   UzaResol.Duale(Vardua);
 439   for (i = 1; i <= CTCinv.RowNumber(); i++) {
 440     for (j = i; j <= CTCinv.RowNumber(); j++) {
 441       CTCinv(i, j) = CTCinv(j, i);
 442     }
 443   }
 444   Done = Standard_True;
 445   math_Vector VecPoles (1, CCol*Npol);
 446   VecPoles = UzaResol.Value();
 447
 448   Standard_Integer polinit = 1, polfin=Npol;
 449   if (FC >= 1) polinit = 2;
 450   if (LC >= 1) polfin = Npol-1;
 451
 452   for (i = polinit; i <= polfin; i++) {
 453     IncCol = 0;
 454     AppParCurves_MultiPoint MPol(nb3d, nb2d);
 455     for (k = 1; k <= NbCu; k++) {
 456       if (k <= nb3d) {
 457         gp_Pnt Pol(VecPoles(IncCol + i),
 458                    VecPoles(IncCol + Npol +i),
 459                    VecPoles(IncCol + 2*Npol + i));
 460         MPol.SetPoint(k, Pol);
 461         IncCol = IncCol + 3*Npol;
 462       }
 463       else {
 464         gp_Pnt2d Pol2d(VecPoles(IncCol + i),
 465                        VecPoles(IncCol + Npol + i));
 466         MPol.SetPoint2d(k, Pol2d);
 467         IncCol = IncCol + 2*Npol;
 468       }
 469     }
 470     SCurv.SetValue(i, MPol);
 471   }
 472
 473 }
 474
 475
 476
 477 Standard_Boolean AppParCurves_ResolConstraint::IsDone () const {
 478   return Done;
 479 }
 480
 481
 482 Standard_Integer AppParCurves_ResolConstraint::
 483   NbConstraints(const MultiLine& SSP,
 484                 const Standard_Integer,
 485                 const Standard_Integer,
 486                 const Handle(AppParCurves_HArray1OfConstraintCouple)&
 487                 TheConstraints)
 488 const
 489 {
 490   // Boucle de calcul du nombre de points de passage afin de dimensionner
 491   // les matrices.
 492   Standard_Integer aIncPass, aIncTan, aIncCurv, aCCol;
 493   Standard_Integer i;
 494   AppParCurves_Constraint Cons;
 495
 496   aIncPass = 0;
 497   aIncTan = 0;
 498   aIncCurv = 0;
 499
 500   for (i = TheConstraints->Lower(); i <= TheConstraints->Upper(); i++) {
 501     Cons = (TheConstraints->Value(i)).Constraint();
 502     if (Cons >= 1) {
 503       aIncPass++;                       // IncPass = nbre de points de passage.
 504     }
 505     if (Cons >= 2) {
 506       aIncTan++;                          // IncTan= nbre de points de tangence.
 507     }
 508     if (Cons == 3) {
 509       aIncCurv++;                      // IncCurv = nbre de pts de courbure.
 510     }
 511   }
 512   Standard_Integer nb3d = ToolLine::NbP3d(SSP);
 513   Standard_Integer nb2d = ToolLine::NbP2d(SSP);
 514   aCCol = nb3d* 3 + nb2d* 2;
 515
 516   return aCCol*aIncPass + aIncTan*(aCCol-1) + 3*aIncCurv;
 517
 518 }
 519
 520
 521 Standard_Integer AppParCurves_ResolConstraint::NbColumns(const MultiLine& SSP,
 522                                                          const Standard_Integer Deg)
 523 const
 524 {
 525   Standard_Integer nb3d = ToolLine::NbP3d(SSP);
 526   Standard_Integer nb2d = ToolLine::NbP2d(SSP);
 527   Standard_Integer CCol = nb3d* 3 + nb2d* 2;
 528
 529   return CCol*(Deg+1);
 530 }
 531
 532 const math_Matrix& AppParCurves_ResolConstraint::ConstraintMatrix() const
 533 {
 534   return Cont;
 535 }
 536
 537 const math_Matrix& AppParCurves_ResolConstraint::InverseMatrix() const
 538 {
 539   return CTCinv;
 540 }
 541
 542
 543 const math_Vector& AppParCurves_ResolConstraint::Duale() const
 544 {
 545   return Vardua;
 546 }
 547
 548
 549
 550 const math_Matrix& AppParCurves_ResolConstraint::
 551                    ConstraintDerivative(const MultiLine& SSP,
 552                                         const math_Vector& Parameters,
 553                                         const Standard_Integer Deg,
 554                                         const math_Matrix& DA)
 555 {
 556   Standard_Integer i, j, k, nb2d, nb3d, CCol, Inc3, IncCol, Npt;
 557   Standard_Integer Npol, Npol2, NbCu =ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP);
 558   Standard_Integer IP;
 559   Standard_Real Daij;
 560   Npol = Deg+1; Npol2 = 2*Npol;
 561   Standard_Boolean Ok;
 562   TColStd_Array2OfInteger Ibont(1, NbCu, 1, IncTan);
 563   math_Matrix DeCInit(1, IncPass, 1, Npol);
 564   math_Vector DDA(1, Npol);
 565   gp_Vec V;
 566   gp_Vec2d V2d;
 567   Standard_Real T1, T2, T3, Tmax, DDaij;
 568 //  Standard_Integer FirstP = IPas(1);
 569   nb3d = ToolLine::NbP3d(SSP);
 570   nb2d = ToolLine::NbP2d(SSP);
 571   Standard_Integer mynb3d = nb3d, mynb2d = nb2d;
 572   if (nb3d == 0) mynb3d = 1;
 573   if (nb2d == 0) mynb2d = 1;
 574   CCol = nb3d* 3 + nb2d* 2;
 575
 576   TColgp_Array1OfVec tabV(1, mynb3d);
 577   TColgp_Array1OfVec2d tabV2d(1, mynb2d);
 578   TColgp_Array1OfPnt tabP(1, mynb3d);
 579   TColgp_Array1OfPnt2d tabP2d(1, mynb2d);
 580
 581   for (i = 1; i <= DeCont.RowNumber(); i++)
 582     for (j = 1; j <= DeCont.ColNumber(); j++)
 583       DeCont(i, j) = 0.0;
 584
 585
 586   //  Remplissage de DK pour les points de passages:
 587
 588   for (i =1; i <= IncPass; i++) {
 589     Npt = IPas(i);
 590     for (j = 1; j <= Npol; j++) DeCInit(i, j) = DA(Npt, j);
 591   }
 592   for (i = 1; i <= CCol; i++) {
 593     DeCont.Set(IncPass*(i-1)+1, IncPass*i, Npol*(i-1)+1, Npol*i, DeCInit);
 594   }
 595
 596
 597   // Pour les points de tangence:
 598
 599   Inc3 = CCol*IncPass +1;
 600   IncCol = 0;
 601
 602   for (k = 1; k <= NbCu; k++) {
 603     if (k <= nb3d) {
 604       for (i = 1; i <= IncTan; i++) {
 605         Npt = ITan(i);
 606 //      MultiPoint MPoint = ToolLine::Value(SSP, Npt);
 607         // choix du maximum de tangence pour exprimer la colinearite:
 608         Ok = ToolLine::Tangency(SSP, Npt, tabV);
 609         V = tabV(k);
 610         V.Coord(T1, T2, T3);
 611         Tmax = Abs(T1);
 612         Ibont(k, i) = 1;
 613         if (Abs(T2) > Tmax) {
 614           Tmax = Abs(T2);
 615           Ibont(k, i) = 2;
 616         }
 617         if (Abs(T3) > Tmax) {
 618           Tmax = Abs(T3);
 619           Ibont(k, i) = 3;
 620         }
 621         AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
 622         if (Ibont(k, i) == 3) {
 623           for (j = 1; j <= Npol; j++) {
 624             DDaij = DDA(j);
 625             DeCont(Inc3, j+ Npol + IncCol) = DDaij*T3/Tmax;
 626             DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T2/Tmax;
 627             DeCont(Inc3+1, j+ IncCol) = DDaij* T3/Tmax;
 628             DeCont(Inc3+1, j+ Npol2 + IncCol) = -DDaij*T1/Tmax;
 629           }
 630         }
 631         else if (Ibont(k, i) == 1) {
 632           for (j = 1; j <= Npol; j++) {
 633             DDaij = DDA(j);
 634             DeCont(Inc3, j + IncCol) = DDaij*T3/Tmax;
 635             DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T1/Tmax;
 636             DeCont(Inc3+1, j+ IncCol) = DDaij* T2/Tmax;
 637             DeCont(Inc3+1, j+ Npol +IncCol) = -DDaij*T1/Tmax;
 638           }
 639         }
 640         else if (Ibont(k, i) == 2) {
 641           for (j = 1; j <= Npol; j++) {
 642             DDaij = DDA(j);
 643             DeCont(Inc3, j+ Npol + IncCol) = DDaij*T3/Tmax;
 644             DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T2/Tmax;
 645             DeCont(Inc3+1, j+ IncCol) = DDaij* T2/Tmax;
 646             DeCont(Inc3+1, j+ Npol + IncCol) = -DDaij*T1/Tmax;
 647           }
 648         }
 649         Inc3 = Inc3 + 2;
 650       }
 651       IncCol = IncCol + 3*Npol;
 652     }
 653     else {
 654       for (i = 1; i <= IncTan; i++) {
 655         Npt = ITan(i);
 656         AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
 657         Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 658         V2d = tabV2d(k);
 659         V2d.Coord(T1, T2);
 660         Tmax = Abs(T1);
 661         Ibont(k, i) = 1;
 662         if (Abs(T2) > Tmax) {
 663           Tmax = Abs(T2);
 664           Ibont(k, i) = 2;
 665         }
 666         for (j = 1; j <= Npol; j++) {
 667           DDaij = DDA(j);
 668           DeCont(Inc3, j + IncCol) = DDaij*T2;
 669           DeCont(Inc3, j + Npol + IncCol) = -DDaij*T1;
 670         }
 671         Inc3 = Inc3 +1;
 672       }
 673     }
 674   }
 675
 676   // Equations exprimant le meme rapport de tangence sur chaque courbe:
 677   // On prend les coordonnees les plus significatives.
 678
 679   Inc3 = Inc3 -1;
 680   for (i =1; i <= IncTan; i++) {
 681     IncCol = 0;
 682     Npt = ITan(i);
 683     AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
 684 //    MultiPoint MP = ToolLine::Value(SSP, Npt);
 685     for (k = 1; k <= NbCu-1; k++) {
 686       Inc3 = Inc3 + 1;
 687       if (Ibont(k, i) == 1) {
 688         if (k <= nb3d) {
 689           Ok = ToolLine::Tangency(SSP, Npt, tabV);
 690           V = tabV(k);
 691           T1 = V.X();
 692           IP = 3*Npol;
 693         }
 694         else {
 695           Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 696           V2d = tabV2d(k);
 697           T1 = V2d.X();
 698           IP = 2*Npol;
 699         }
 700         if (Ibont(k+1, i) == 1) {  // Relations entre T1x et T2x:
 701           if ((k+1) <= nb3d) {
 702             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 703             V = tabV(k+1);
 704             T2 = V.X();
 705           }
 706           else {
 707             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 708             V2d = tabV2d(k+1);
 709             T2 = V2d.X();
 710           }
 711           for (j = 1; j <=  Npol; j++) {
 712             Daij = DDA(j);
 713             Cont(Inc3, j + IncCol) = Daij*T2;
 714             Cont(Inc3, j + IP + IncCol) = -Daij*T1;
 715           }
 716           IncCol = IncCol + IP;
 717         }
 718         else if (Ibont(k+1, i) == 2) {  // Relations entre T1x et T2y:
 719           if ((k+1) <= nb3d) {
 720             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 721             V = tabV(k+1);
 722             T2 = V.Y();
 723           }
 724           else {
 725             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 726             V2d = tabV2d(k+1);
 727             T2 = V2d.Y();
 728           }
 729           for (j = 1; j <=  Npol; j++) {
 730             Daij = DDA(j);
 731             Cont(Inc3, j + IncCol) = Daij*T2;
 732             Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
 733           }
 734           IncCol = IncCol + IP;
 735         }
 736         else if (Ibont(k+1,i) == 3) {    // Relations entre T1x et T2z:
 737           Ok = ToolLine::Tangency(SSP, Npt, tabV);
 738           V = tabV(k+1);
 739           T2 = V.Z();
 740           for (j = 1; j <=  Npol; j++) {
 741             Daij = DDA(j);
 742             Cont(Inc3, j + IncCol) = Daij*T2;
 743             Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
 744           }
 745           IncCol = IncCol + IP;
 746         }
 747       }
 748       else if (Ibont(k,i) == 2) {
 749         if (k <= nb3d) {
 750           Ok = ToolLine::Tangency(SSP, Npt, tabV);
 751           V = tabV(k);
 752           T1 = V.Y();
 753           IP = 3*Npol;
 754         }
 755         else {
 756           Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 757           V2d = tabV2d(k);
 758           T1 = V2d.Y();
 759           IP = 2*Npol;
 760         }
 761         if (Ibont(k+1, i) == 1) {  // Relations entre T1y et T2x:
 762           if ((k+1) <= nb3d) {
 763             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 764             V = tabV(k+1);
 765             T2 = V.X();
 766           }
 767           else {
 768             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 769             V2d = tabV2d(k+1);
 770             T2 = V2d.X();
 771           }
 772           for (j = 1; j <=  Npol; j++) {
 773             Daij = DDA( j);
 774             Cont(Inc3, j + Npol + IncCol) = Daij*T2;
 775             Cont(Inc3, j + IP + IncCol) = -Daij*T1;
 776           }
 777           IncCol = IncCol + IP;
 778
 779         }
 780         else if (Ibont(k+1, i) == 2) {  // Relations entre T1y et T2y:
 781           if ((k+1) <= nb3d) {
 782             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 783             V = tabV(k+1);
 784             T2 = V.Y();
 785           }
 786           else {
 787             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 788             V2d = tabV2d(k+1);
 789             T2 = V2d.Y();
 790           }
 791           for (j = 1; j <=  Npol; j++) {
 792             Daij = DDA(j);
 793             Cont(Inc3, j + Npol + IncCol) = Daij*T2;
 794             Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
 795           }
 796           IncCol = IncCol + IP;
 797
 798         }
 799         else if (Ibont(k+1,i)== 3) {    // Relations entre T1y et T2z:
 800           Ok = ToolLine::Tangency(SSP, Npt, tabV);
 801           V = tabV(k+1);
 802           T2 = V.Z();
 803           for (j = 1; j <=  Npol; j++) {
 804             Daij = DDA(j);
 805             Cont(Inc3, j + Npol +IncCol) = Daij*T2;
 806             Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
 807           }
 808           IncCol = IncCol + IP;
 809         }
 810       }
 811
 812       else {
 813         Ok = ToolLine::Tangency(SSP, Npt, tabV);
 814         V = tabV(k);
 815         T1 = V.Z();
 816         IP = 3*Npol;
 817         if (Ibont(k+1, i) == 1) {  // Relations entre T1z et T2x:
 818           if ((k+1) <= nb3d) {
 819             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 820             V = tabV(k+1);
 821             T2 = V.X();
 822           }
 823           else {
 824             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 825             V2d = tabV2d(k+1);
 826             T2 = V2d.X();
 827           }
 828           for (j = 1; j <=  Npol; j++) {
 829             Daij = DDA(j);
 830             Cont(Inc3, j + 2*Npol +IncCol) = Daij*T2;
 831             Cont(Inc3, j + IP + IncCol) = -Daij*T1;
 832           }
 833           IncCol = IncCol + IP;
 834         }
 835
 836         else if (Ibont(k+1, i) == 2) {  // Relations entre T1z et T2y:
 837           if ((k+1) <= nb3d) {
 838             Ok = ToolLine::Tangency(SSP, Npt, tabV);
 839             V = tabV(k+1);
 840             T2 = V.Y();
 841           }
 842           else {
 843             Ok = ToolLine::Tangency(SSP, Npt, tabV2d);
 844             V2d = tabV2d(k+1);
 845             T2 = V2d.Y();
 846           }
 847           for (j = 1; j <=  Npol; j++) {
 848             Daij = DDA(j);
 849             Cont(Inc3, j + 2*Npol +IncCol) = Daij*T2;
 850             Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
 851           }
 852           IncCol = IncCol + IP;
 853         }
 854
 855         else if (Ibont(k+1,i)==3) {      // Relations entre T1z et T2z:
 856           Ok = ToolLine::Tangency(SSP, Npt, tabV);
 857           V = tabV(k+1);
 858           T2 = V.Z();
 859           for (j = 1; j <=  Npol; j++) {
 860             Daij = DDA(j);
 861             Cont(Inc3, j + 2*Npol + IncCol) = Daij*T2;
 862             Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
 863           }
 864           IncCol = IncCol + IP;
 865         }
 866       }
 867     }
 868   }
 869
 870   return DeCont;
 871 }