Warning, /include/opencascade/AppParCurves_ResolConstraint.gxx is written in an unsupported language. File is not indexed.
0001 // Copyright (c) 1995-1999 Matra Datavision
0002 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0003 //
0004 // This file is part of Open CASCADE Technology software library.
0005 //
0006 // This library is free software; you can redistribute it and/or modify it under
0007 // the terms of the GNU Lesser General Public License version 2.1 as published
0008 // by the Free Software Foundation, with special exception defined in the file
0009 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0010 // distribution for complete text of the license and disclaimer of any warranty.
0011 //
0012 // Alternatively, this file may be used under the terms of Open CASCADE
0013 // commercial license or contractual agreement.
0014
0015 // lpa, le 11/09/91
0016
0017
0018 // Approximation d un ensemble de points contraints (MultiLine) avec une
0019 // solution approchee (MultiCurve). L algorithme utilise est l algorithme
0020 // d Uzawa du package mathematique.
0021
0022 #define No_Standard_RangeError
0023 #define No_Standard_OutOfRange
0024
0025 #include <math_Vector.hxx>
0026 #include <math_Matrix.hxx>
0027 #include <AppParCurves_Constraint.hxx>
0028 #include <AppParCurves_ConstraintCouple.hxx>
0029 #include <AppParCurves_MultiPoint.hxx>
0030 #include <AppParCurves.hxx>
0031 #include <Standard_DimensionError.hxx>
0032 #include <math_Uzawa.hxx>
0033 #include <TColStd_Array1OfInteger.hxx>
0034 #include <TColStd_Array2OfInteger.hxx>
0035 #include <gp_Pnt.hxx>
0036 #include <gp_Pnt2d.hxx>
0037 #include <gp_Vec.hxx>
0038 #include <gp_Vec2d.hxx>
0039 #include <TColgp_Array1OfPnt.hxx>
0040 #include <TColgp_Array1OfPnt2d.hxx>
0041 #include <TColgp_Array1OfVec.hxx>
0042 #include <TColgp_Array1OfVec2d.hxx>
0043
0044
0045 AppParCurves_ResolConstraint::
0046 AppParCurves_ResolConstraint(const MultiLine& SSP,
0047 AppParCurves_MultiCurve& SCurv,
0048 const Standard_Integer FirstPoint,
0049 const Standard_Integer LastPoint,
0050 const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
0051 const math_Matrix& Bern,
0052 const math_Matrix& DerivativeBern,
0053 const Standard_Real Tolerance):
0054 Cont(1, NbConstraints(SSP, FirstPoint,
0055 LastPoint, TheConstraints),
0056 1, NbColumns(SSP, SCurv.NbPoles()-1), 0.0),
0057 DeCont(1, NbConstraints(SSP, FirstPoint,
0058 LastPoint, TheConstraints),
0059 1, NbColumns(SSP, SCurv.NbPoles()-1), 0.0),
0060 Secont(1, NbConstraints(SSP, FirstPoint,
0061 LastPoint, TheConstraints), 0.0),
0062 CTCinv(1, NbConstraints(SSP, FirstPoint,
0063 LastPoint, TheConstraints),
0064 1, NbConstraints(SSP, FirstPoint,
0065 LastPoint, TheConstraints)),
0066 Vardua(1, NbConstraints(SSP, FirstPoint,
0067 LastPoint, TheConstraints)),
0068 IPas(1, LastPoint-FirstPoint+1),
0069 ITan(1, LastPoint-FirstPoint+1),
0070 ICurv(1, LastPoint-FirstPoint+1)
0071 {
0072 Standard_Integer i, j, k, NbCu= SCurv.NbCurves();
0073 Standard_Integer Npt;
0074 Standard_Integer Inc3, IncSec, IncCol, IP = 0, CCol;
0075 Standard_Integer myindex, Def = SCurv.NbPoles()-1;
0076 Standard_Integer nb3d, nb2d, Npol= Def+1, Npol2 = 2*Npol;
0077 Standard_Real T1 = 0., T2 = 0., T3, Tmax, Daij;
0078 gp_Vec V;
0079 gp_Vec2d V2d;
0080 gp_Pnt Poi;
0081 gp_Pnt2d Poi2d;
0082 AppParCurves_ConstraintCouple mycouple;
0083 AppParCurves_Constraint FC = AppParCurves_NoConstraint,
0084 LC = AppParCurves_NoConstraint, Cons;
0085
0086
0087
0088 // Boucle de calcul du nombre de points de passage afin de dimensionner
0089 // les matrices.
0090 IncPass = 0;
0091 IncTan = 0;
0092 IncCurv = 0;
0093 for (i = TheConstraints->Lower(); i <= TheConstraints->Upper(); i++) {
0094 mycouple = TheConstraints->Value(i);
0095 myindex = mycouple.Index();
0096 Cons = mycouple.Constraint();
0097 if (myindex == FirstPoint) FC = Cons;
0098 if (myindex == LastPoint) LC = Cons;
0099 if (Cons >= 1) {
0100 IncPass++; // IncPass = nbre de points de passage.
0101 IPas(IncPass) = myindex;
0102 }
0103 if (Cons >= 2) {
0104 IncTan++; // IncTan= nbre de points de tangence.
0105 ITan(IncTan) = myindex;
0106 }
0107 if (Cons == 3) {
0108 IncCurv++; // IncCurv = nbre de pts de courbure.
0109 ICurv(IncCurv) = myindex;
0110 }
0111 }
0112
0113
0114 if (IncPass == 0) {
0115 Done = Standard_True;
0116 return;
0117 }
0118
0119 nb3d = ToolLine::NbP3d(SSP);
0120 nb2d = ToolLine::NbP2d(SSP);
0121 Standard_Integer mynb3d=nb3d, mynb2d=nb2d;
0122 if (nb3d == 0) mynb3d = 1;
0123 if (nb2d == 0) mynb2d = 1;
0124 CCol = nb3d* 3 + nb2d* 2;
0125
0126
0127 // Declaration et initialisation des matrices et vecteurs de contraintes:
0128 math_Matrix ContInit(1, IncPass, 1, Npol);
0129 math_Vector Start(1, CCol*Npol);
0130 TColStd_Array2OfInteger Ibont(1, NbCu, 1, IncTan);
0131
0132
0133 // Remplissage de Cont pour les points de passage:
0134 // =================================================
0135 for (i =1; i <= IncPass; i++) { // Cette partie ne depend que de Bernstein
0136 Npt = IPas(i);
0137 for (j = 1; j <= Npol; j++) {
0138 ContInit(i, j) = Bern(Npt, j);
0139 }
0140 }
0141 for (i = 1; i <= CCol; i++) {
0142 Cont.Set(IncPass*(i-1)+1, IncPass*i, Npol*(i-1)+1, Npol*i, ContInit);
0143 }
0144
0145
0146 // recuperation des vecteurs de depart pour Uzawa. Ce vecteur represente les
0147 // poles de SCurv.
0148 // Remplissage de secont et resolution.
0149
0150 TColgp_Array1OfVec tabV(1, mynb3d);
0151 TColgp_Array1OfVec2d tabV2d(1, mynb2d);
0152 TColgp_Array1OfPnt tabP(1, mynb3d);
0153 TColgp_Array1OfPnt2d tabP2d(1, mynb2d);
0154
0155 Inc3 = CCol*IncPass +1;
0156 IncCol = 0;
0157 IncSec = 0;
0158 for (k = 1; k <= NbCu; k++) {
0159 if (k <= nb3d) {
0160 for (i = 1; i <= IncTan; i++) {
0161 Npt = ITan(i);
0162 // choix du maximum de tangence pour exprimer la colinearite:
0163 ToolLine::Tangency(SSP, Npt, tabV);
0164 V = tabV(k);
0165 V.Coord(T1, T2, T3);
0166 Tmax = Abs(T1);
0167 Ibont(k, i) = 1;
0168 if (Abs(T2) > Tmax) {
0169 Tmax = Abs(T2);
0170 Ibont(k, i) = 2;
0171 }
0172 if (Abs(T3) > Tmax) {
0173 Tmax = Abs(T3);
0174 Ibont(k, i) = 3;
0175 }
0176 if (Ibont(k, i) == 3) {
0177 for (j = 1; j <= Npol; j++) {
0178 Daij = DerivativeBern(Npt, j);
0179 Cont(Inc3, j+ Npol + IncCol) = Daij*T3/Tmax;
0180 Cont(Inc3, j + Npol2 + IncCol) = -Daij *T2/Tmax;
0181 Cont(Inc3+1, j+ IncCol) = Daij* T3/Tmax;
0182 Cont(Inc3+1, j+ Npol2 + IncCol) = -Daij*T1/Tmax;
0183 }
0184 }
0185 else if (Ibont(k, i) == 1) {
0186 for (j = 1; j <= Npol; j++) {
0187 Daij = DerivativeBern(Npt, j);
0188 Cont(Inc3, j + IncCol) = Daij*T3/Tmax;
0189 Cont(Inc3, j + Npol2 + IncCol) = -Daij *T1/Tmax;
0190 Cont(Inc3+1, j+ IncCol) = Daij* T2/Tmax;
0191 Cont(Inc3+1, j+ Npol +IncCol) = -Daij*T1/Tmax;
0192 }
0193 }
0194 else if (Ibont(k, i) == 2) {
0195 for (j = 1; j <= Def+1; j++) {
0196 Daij = DerivativeBern(Npt, j);
0197 Cont(Inc3, j+ Npol + IncCol) = Daij*T3/Tmax;
0198 Cont(Inc3, j + Npol2 + IncCol) = -Daij *T2/Tmax;
0199 Cont(Inc3+1, j+ IncCol) = Daij* T2/Tmax;
0200 Cont(Inc3+1, j+ Npol + IncCol) = -Daij*T1/Tmax;
0201 }
0202 }
0203 Inc3 = Inc3 + 2;
0204 }
0205
0206 // Remplissage du second membre:
0207 for (i = 1; i <= IncPass; i++) {
0208 ToolLine::Value(SSP, IPas(i), tabP);
0209 Poi = tabP(k);
0210 Secont(i + IncSec) = Poi.X();
0211 Secont(i + IncPass + IncSec) = Poi.Y();
0212 Secont(i + 2*IncPass + IncSec) = Poi.Z();
0213 }
0214 IncSec = IncSec + 3*IncPass;
0215
0216 // Vecteur de depart:
0217 for (j =1; j <= Npol; j++) {
0218 Poi = SCurv.Value(j).Point(k);
0219 Start(j + IncCol) = Poi.X();
0220 Start(j+ Npol + IncCol) = Poi.Y();
0221 Start(j + Npol2 + IncCol) = Poi.Z();
0222 }
0223 IncCol = IncCol + 3*Npol;
0224 }
0225
0226
0227 else {
0228 for (i = 1; i <= IncTan; i++) {
0229 Npt = ITan(i);
0230 ToolLine::Tangency(SSP, Npt, tabV2d);
0231 V2d = tabV2d(k-nb3d);
0232 T1 = V2d.X();
0233 T2 = V2d.Y();
0234 Tmax = Abs(T1);
0235 Ibont(k, i) = 1;
0236 if (Abs(T2) > Tmax) {
0237 Tmax = Abs(T2);
0238 Ibont(k, i) = 2;
0239 }
0240 for (j = 1; j <= Npol; j++) {
0241 Daij = DerivativeBern(Npt, j);
0242 Cont(Inc3, j + IncCol) = Daij*T2;
0243 Cont(Inc3, j + Npol + IncCol) = -Daij*T1;
0244 }
0245 Inc3 = Inc3 +1;
0246 }
0247
0248 // Remplissage du second membre:
0249 for (i = 1; i <= IncPass; i++) {
0250 ToolLine::Value(SSP, IPas(i), tabP2d);
0251 Poi2d = tabP2d(i-nb3d);
0252 Secont(i + IncSec) = Poi2d.X();
0253 Secont(i + IncPass + IncSec) = Poi2d.Y();
0254 }
0255 IncSec = IncSec + 2*IncPass;
0256
0257 // Remplissage du vecteur de depart:
0258 for (j = 1; j <= Npol; j++) {
0259 Poi2d = SCurv.Value(j).Point2d(k);
0260 Start(j + IncCol) = Poi2d.X();
0261 Start(j + Npol + IncCol) = Poi2d.Y();
0262 }
0263 IncCol = IncCol + Npol2;
0264 }
0265 }
0266
0267 // Equations exprimant le meme rapport de tangence sur chaque courbe:
0268 // On prend les coordonnees les plus significatives.
0269
0270 --Inc3;
0271 for (i = 1; i <= IncTan; ++i) {
0272 IncCol = 0;
0273 Npt = ITan(i);
0274 for (k = 1; k <= NbCu-1; ++k) {
0275 ++Inc3;
0276 // Initialize first relation variable (T1)
0277 Standard_Integer addIndex_1 = 0, aVal = Ibont(k, i);
0278 switch (aVal)
0279 {
0280 case 1: //T1 ~ T1x
0281 {
0282 if (k <= nb3d) {
0283 ToolLine::Tangency(SSP, Npt, tabV);
0284 V = tabV(k);
0285 T1 = V.X();
0286 IP = 3 * Npol;
0287 }
0288 else {
0289 ToolLine::Tangency(SSP, Npt, tabV2d);
0290 V2d = tabV2d(k-nb3d);
0291 T1 = V2d.X();
0292 IP = 2 * Npol;
0293 }
0294 addIndex_1 = 0;
0295 break;
0296 }
0297 case 2: //T1 ~ T1y
0298 {
0299 if (k <= nb3d) {
0300 ToolLine::Tangency(SSP, Npt, tabV);
0301 V = tabV(k);
0302 IP = 3 * Npol;
0303 }
0304 else {
0305 ToolLine::Tangency(SSP, Npt, tabV2d);
0306 V2d = tabV2d(k-nb3d);
0307 T1 = V2d.Y();
0308 IP = 2 * Npol;
0309 }
0310 addIndex_1 = Npol;
0311 break;
0312 }
0313 case 3: //T1 ~ T1z
0314 {
0315 ToolLine::Tangency(SSP, Npt, tabV);
0316 V = tabV(k);
0317 T1 = V.Z();
0318 IP = 3 * Npol;
0319 addIndex_1 = 2 * Npol;
0320 break;
0321 }
0322 }
0323 // Initialize second relation variable (T2)
0324 Standard_Integer addIndex_2 = 0, aNextVal = Ibont(k + 1, i);
0325 switch (aNextVal)
0326 {
0327 case 1: //T2 ~ T2x
0328 {
0329 if ((k+1) <= nb3d) {
0330 ToolLine::Tangency(SSP, Npt, tabV);
0331 V = tabV(k+1);
0332 T2 = V.X();
0333 }
0334 else {
0335 ToolLine::Tangency(SSP, Npt, tabV2d);
0336 V2d = tabV2d(k+1-nb3d);
0337 T2 = V2d.X();
0338 }
0339 addIndex_2 = 0;
0340 break;
0341 }
0342 case 2: //T2 ~ T2y
0343 {
0344 if ((k+1) <= nb3d) {
0345 ToolLine::Tangency(SSP, Npt, tabV);
0346 V = tabV(k+1);
0347 T2 = V.Y();
0348 }
0349 else {
0350 ToolLine::Tangency(SSP, Npt, tabV2d);
0351 V2d = tabV2d(k+1-nb3d);
0352 T2 = V2d.Y();
0353 }
0354 addIndex_2 = Npol;
0355 break;
0356 }
0357 case 3: //T2 ~ T2z
0358 {
0359 ToolLine::Tangency(SSP, Npt, tabV);
0360 V = tabV(k+1);
0361 T2 = V.Z();
0362 addIndex_2 = 2 * Npol;
0363 break;
0364 }
0365 }
0366
0367 // Relations between T1 and T2:
0368 for (j = 1; j <= Npol; j++) {
0369 Daij = DerivativeBern(Npt, j);
0370 Cont(Inc3, j + IncCol + addIndex_1) = Daij*T2;
0371 Cont(Inc3, j + IP + IncCol + addIndex_2) = -Daij*T1;
0372 }
0373 IncCol += IP;
0374 }
0375 }
0376
0377
0378
0379 // Equations concernant la courbure:
0380
0381
0382 /* Inc3 = Inc3 +1;
0383 IncCol = 0;
0384 for (k = 1; k <= NbCu; k++) {
0385 for (i = 1; i <= IncCurv; i++) {
0386 Npt = ICurv(i);
0387 AppDef_MultiPointConstraint MP = SSP.Value(Npt);
0388 DDA = SecondDerivativeBern(Parameters(Npt));
0389 if (SSP.Value(1).Dimension(k) == 3) {
0390 C1 = MP.Curv(k).X();
0391 C2 = MP.Curv(k).Y();
0392 C3 = MP.Curv(k).Z();
0393 for (j = 1; j <= Npol; j++) {
0394 Daij = DerivativeBern(Npt, j);
0395 D2Aij = DDA(j);
0396 Cont(Inc3, j + IncCol) = D2Aij;
0397 Cont(Inc3, j + Npol2 + IncCol) = -C2*Daij;
0398 Cont(Inc3, j + Npol + IncCol) = C3*Daij;
0399
0400 Cont(Inc3+1, j + Npol + IncCol) = D2Aij;
0401 Cont(Inc3+1, j +IncCol) = -C3*Daij;
0402 Cont(Inc3+1, j + Npol2 + IncCol) = C1*Daij;
0403
0404 Cont(Inc3+2, j + Npol2+IncCol) = D2Aij;
0405 Cont(Inc3+2, j + Npol+IncCol) = -C1*Daij;
0406 Cont(Inc3+2, j + IncCol) = C2*Daij;
0407 }
0408 Inc3 = Inc3 + 3;
0409 }
0410 else { // Dimension 2:
0411 C1 = MP.Curv2d(k).X();
0412 C2 = MP.Curv2d(k).Y();
0413 for (j = 1; j <= Npol; j++) {
0414 Daij = DerivativeBern(Npt, j);
0415 D2Aij = DDA(j);
0416 Cont(Inc3, j + IncCol) = D2Aij*C1;
0417 Cont(Inc3+1, j + Npol + IncCol) = D2Aij*C2;
0418 }
0419 Inc3 = Inc3 + 2;
0420 }
0421 }
0422 }
0423
0424 */
0425 // Resolution par Uzawa:
0426
0427 math_Uzawa UzaResol(Cont, Secont, Start, Tolerance);
0428 if (!(UzaResol.IsDone())) {
0429 Done = Standard_False;
0430 return;
0431 }
0432 CTCinv = UzaResol.InverseCont();
0433 UzaResol.Duale(Vardua);
0434 for (i = 1; i <= CTCinv.RowNumber(); i++) {
0435 for (j = i; j <= CTCinv.RowNumber(); j++) {
0436 CTCinv(i, j) = CTCinv(j, i);
0437 }
0438 }
0439 Done = Standard_True;
0440 math_Vector VecPoles (1, CCol*Npol);
0441 VecPoles = UzaResol.Value();
0442
0443 Standard_Integer polinit = 1, polfin=Npol;
0444 if (FC >= 1) polinit = 2;
0445 if (LC >= 1) polfin = Npol-1;
0446
0447 for (i = polinit; i <= polfin; i++) {
0448 IncCol = 0;
0449 AppParCurves_MultiPoint MPol(nb3d, nb2d);
0450 for (k = 1; k <= NbCu; k++) {
0451 if (k <= nb3d) {
0452 gp_Pnt Pol(VecPoles(IncCol + i),
0453 VecPoles(IncCol + Npol +i),
0454 VecPoles(IncCol + 2*Npol + i));
0455 MPol.SetPoint(k, Pol);
0456 IncCol = IncCol + 3*Npol;
0457 }
0458 else {
0459 gp_Pnt2d Pol2d(VecPoles(IncCol + i),
0460 VecPoles(IncCol + Npol + i));
0461 MPol.SetPoint2d(k, Pol2d);
0462 IncCol = IncCol + 2*Npol;
0463 }
0464 }
0465 SCurv.SetValue(i, MPol);
0466 }
0467
0468 }
0469
0470
0471
0472 Standard_Boolean AppParCurves_ResolConstraint::IsDone () const {
0473 return Done;
0474 }
0475
0476
0477 Standard_Integer AppParCurves_ResolConstraint::
0478 NbConstraints(const MultiLine& SSP,
0479 const Standard_Integer,
0480 const Standard_Integer,
0481 const Handle(AppParCurves_HArray1OfConstraintCouple)&
0482 TheConstraints)
0483 const
0484 {
0485 // Boucle de calcul du nombre de points de passage afin de dimensionner
0486 // les matrices.
0487 Standard_Integer aIncPass, aIncTan, aIncCurv, aCCol;
0488 Standard_Integer i;
0489 AppParCurves_Constraint Cons;
0490
0491 aIncPass = 0;
0492 aIncTan = 0;
0493 aIncCurv = 0;
0494
0495 for (i = TheConstraints->Lower(); i <= TheConstraints->Upper(); i++) {
0496 Cons = (TheConstraints->Value(i)).Constraint();
0497 if (Cons >= 1) {
0498 aIncPass++; // IncPass = nbre de points de passage.
0499 }
0500 if (Cons >= 2) {
0501 aIncTan++; // IncTan= nbre de points de tangence.
0502 }
0503 if (Cons == 3) {
0504 aIncCurv++; // IncCurv = nbre de pts de courbure.
0505 }
0506 }
0507 Standard_Integer nb3d = ToolLine::NbP3d(SSP);
0508 Standard_Integer nb2d = ToolLine::NbP2d(SSP);
0509 aCCol = nb3d* 3 + nb2d* 2;
0510
0511 return aCCol*aIncPass + aIncTan*(aCCol-1) + 3*aIncCurv;
0512
0513 }
0514
0515
0516 Standard_Integer AppParCurves_ResolConstraint::NbColumns(const MultiLine& SSP,
0517 const Standard_Integer Deg)
0518 const
0519 {
0520 Standard_Integer nb3d = ToolLine::NbP3d(SSP);
0521 Standard_Integer nb2d = ToolLine::NbP2d(SSP);
0522 Standard_Integer CCol = nb3d* 3 + nb2d* 2;
0523
0524 return CCol*(Deg+1);
0525 }
0526
0527 const math_Matrix& AppParCurves_ResolConstraint::ConstraintMatrix() const
0528 {
0529 return Cont;
0530 }
0531
0532 const math_Matrix& AppParCurves_ResolConstraint::InverseMatrix() const
0533 {
0534 return CTCinv;
0535 }
0536
0537
0538 const math_Vector& AppParCurves_ResolConstraint::Duale() const
0539 {
0540 return Vardua;
0541 }
0542
0543
0544
0545 const math_Matrix& AppParCurves_ResolConstraint::
0546 ConstraintDerivative(const MultiLine& SSP,
0547 const math_Vector& Parameters,
0548 const Standard_Integer Deg,
0549 const math_Matrix& DA)
0550 {
0551 Standard_Integer i, j, k, nb2d, nb3d, CCol, Inc3, IncCol, Npt;
0552 Standard_Integer Npol, Npol2, NbCu =ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP);
0553 Standard_Integer IP;
0554 Standard_Real Daij;
0555 Npol = Deg+1; Npol2 = 2*Npol;
0556 TColStd_Array2OfInteger Ibont(1, NbCu, 1, IncTan);
0557 math_Matrix DeCInit(1, IncPass, 1, Npol);
0558 math_Vector DDA(1, Npol);
0559 gp_Vec V;
0560 gp_Vec2d V2d;
0561 Standard_Real T1, T2, T3, Tmax, DDaij;
0562 // Standard_Integer FirstP = IPas(1);
0563 nb3d = ToolLine::NbP3d(SSP);
0564 nb2d = ToolLine::NbP2d(SSP);
0565 Standard_Integer mynb3d = nb3d, mynb2d = nb2d;
0566 if (nb3d == 0) mynb3d = 1;
0567 if (nb2d == 0) mynb2d = 1;
0568 CCol = nb3d* 3 + nb2d* 2;
0569
0570 TColgp_Array1OfVec tabV(1, mynb3d);
0571 TColgp_Array1OfVec2d tabV2d(1, mynb2d);
0572 TColgp_Array1OfPnt tabP(1, mynb3d);
0573 TColgp_Array1OfPnt2d tabP2d(1, mynb2d);
0574
0575 for (i = 1; i <= DeCont.RowNumber(); i++)
0576 for (j = 1; j <= DeCont.ColNumber(); j++)
0577 DeCont(i, j) = 0.0;
0578
0579
0580 // Remplissage de DK pour les points de passages:
0581
0582 for (i =1; i <= IncPass; i++) {
0583 Npt = IPas(i);
0584 for (j = 1; j <= Npol; j++) DeCInit(i, j) = DA(Npt, j);
0585 }
0586 for (i = 1; i <= CCol; i++) {
0587 DeCont.Set(IncPass*(i-1)+1, IncPass*i, Npol*(i-1)+1, Npol*i, DeCInit);
0588 }
0589
0590
0591 // Pour les points de tangence:
0592
0593 Inc3 = CCol*IncPass +1;
0594 IncCol = 0;
0595
0596 for (k = 1; k <= NbCu; k++) {
0597 if (k <= nb3d) {
0598 for (i = 1; i <= IncTan; i++) {
0599 Npt = ITan(i);
0600 // MultiPoint MPoint = ToolLine::Value(SSP, Npt);
0601 // choix du maximum de tangence pour exprimer la colinearite:
0602 ToolLine::Tangency(SSP, Npt, tabV);
0603 V = tabV(k);
0604 V.Coord(T1, T2, T3);
0605 Tmax = Abs(T1);
0606 Ibont(k, i) = 1;
0607 if (Abs(T2) > Tmax) {
0608 Tmax = Abs(T2);
0609 Ibont(k, i) = 2;
0610 }
0611 if (Abs(T3) > Tmax) {
0612 Tmax = Abs(T3);
0613 Ibont(k, i) = 3;
0614 }
0615 AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
0616 if (Ibont(k, i) == 3) {
0617 for (j = 1; j <= Npol; j++) {
0618 DDaij = DDA(j);
0619 DeCont(Inc3, j+ Npol + IncCol) = DDaij*T3/Tmax;
0620 DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T2/Tmax;
0621 DeCont(Inc3+1, j+ IncCol) = DDaij* T3/Tmax;
0622 DeCont(Inc3+1, j+ Npol2 + IncCol) = -DDaij*T1/Tmax;
0623 }
0624 }
0625 else if (Ibont(k, i) == 1) {
0626 for (j = 1; j <= Npol; j++) {
0627 DDaij = DDA(j);
0628 DeCont(Inc3, j + IncCol) = DDaij*T3/Tmax;
0629 DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T1/Tmax;
0630 DeCont(Inc3+1, j+ IncCol) = DDaij* T2/Tmax;
0631 DeCont(Inc3+1, j+ Npol +IncCol) = -DDaij*T1/Tmax;
0632 }
0633 }
0634 else if (Ibont(k, i) == 2) {
0635 for (j = 1; j <= Npol; j++) {
0636 DDaij = DDA(j);
0637 DeCont(Inc3, j+ Npol + IncCol) = DDaij*T3/Tmax;
0638 DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T2/Tmax;
0639 DeCont(Inc3+1, j+ IncCol) = DDaij* T2/Tmax;
0640 DeCont(Inc3+1, j+ Npol + IncCol) = -DDaij*T1/Tmax;
0641 }
0642 }
0643 Inc3 = Inc3 + 2;
0644 }
0645 IncCol = IncCol + 3*Npol;
0646 }
0647 else {
0648 for (i = 1; i <= IncTan; i++) {
0649 Npt = ITan(i);
0650 AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
0651 ToolLine::Tangency(SSP, Npt, tabV2d);
0652 V2d = tabV2d(k);
0653 V2d.Coord(T1, T2);
0654 Tmax = Abs(T1);
0655 Ibont(k, i) = 1;
0656 if (Abs(T2) > Tmax) {
0657 Tmax = Abs(T2);
0658 Ibont(k, i) = 2;
0659 }
0660 for (j = 1; j <= Npol; j++) {
0661 DDaij = DDA(j);
0662 DeCont(Inc3, j + IncCol) = DDaij*T2;
0663 DeCont(Inc3, j + Npol + IncCol) = -DDaij*T1;
0664 }
0665 Inc3 = Inc3 +1;
0666 }
0667 }
0668 }
0669
0670 // Equations exprimant le meme rapport de tangence sur chaque courbe:
0671 // On prend les coordonnees les plus significatives.
0672
0673 Inc3 = Inc3 -1;
0674 for (i =1; i <= IncTan; i++) {
0675 IncCol = 0;
0676 Npt = ITan(i);
0677 AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
0678 // MultiPoint MP = ToolLine::Value(SSP, Npt);
0679 for (k = 1; k <= NbCu-1; k++) {
0680 Inc3 = Inc3 + 1;
0681 if (Ibont(k, i) == 1) {
0682 if (k <= nb3d) {
0683 ToolLine::Tangency(SSP, Npt, tabV);
0684 V = tabV(k);
0685 T1 = V.X();
0686 IP = 3*Npol;
0687 }
0688 else {
0689 ToolLine::Tangency(SSP, Npt, tabV2d);
0690 V2d = tabV2d(k);
0691 T1 = V2d.X();
0692 IP = 2*Npol;
0693 }
0694 if (Ibont(k+1, i) == 1) { // Relations entre T1x et T2x:
0695 if ((k+1) <= nb3d) {
0696 ToolLine::Tangency(SSP, Npt, tabV);
0697 V = tabV(k+1);
0698 T2 = V.X();
0699 }
0700 else {
0701 ToolLine::Tangency(SSP, Npt, tabV2d);
0702 V2d = tabV2d(k+1);
0703 T2 = V2d.X();
0704 }
0705 for (j = 1; j <= Npol; j++) {
0706 Daij = DDA(j);
0707 Cont(Inc3, j + IncCol) = Daij*T2;
0708 Cont(Inc3, j + IP + IncCol) = -Daij*T1;
0709 }
0710 IncCol = IncCol + IP;
0711 }
0712 else if (Ibont(k+1, i) == 2) { // Relations entre T1x et T2y:
0713 if ((k+1) <= nb3d) {
0714 ToolLine::Tangency(SSP, Npt, tabV);
0715 V = tabV(k+1);
0716 T2 = V.Y();
0717 }
0718 else {
0719 ToolLine::Tangency(SSP, Npt, tabV2d);
0720 V2d = tabV2d(k+1);
0721 T2 = V2d.Y();
0722 }
0723 for (j = 1; j <= Npol; j++) {
0724 Daij = DDA(j);
0725 Cont(Inc3, j + IncCol) = Daij*T2;
0726 Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
0727 }
0728 IncCol = IncCol + IP;
0729 }
0730 else if (Ibont(k+1,i) == 3) { // Relations entre T1x et T2z:
0731 ToolLine::Tangency(SSP, Npt, tabV);
0732 V = tabV(k+1);
0733 T2 = V.Z();
0734 for (j = 1; j <= Npol; j++) {
0735 Daij = DDA(j);
0736 Cont(Inc3, j + IncCol) = Daij*T2;
0737 Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
0738 }
0739 IncCol = IncCol + IP;
0740 }
0741 }
0742 else if (Ibont(k,i) == 2) {
0743 if (k <= nb3d) {
0744 ToolLine::Tangency(SSP, Npt, tabV);
0745 V = tabV(k);
0746 T1 = V.Y();
0747 IP = 3*Npol;
0748 }
0749 else {
0750 ToolLine::Tangency(SSP, Npt, tabV2d);
0751 V2d = tabV2d(k);
0752 T1 = V2d.Y();
0753 IP = 2*Npol;
0754 }
0755 if (Ibont(k+1, i) == 1) { // Relations entre T1y et T2x:
0756 if ((k+1) <= nb3d) {
0757 ToolLine::Tangency(SSP, Npt, tabV);
0758 V = tabV(k+1);
0759 T2 = V.X();
0760 }
0761 else {
0762 ToolLine::Tangency(SSP, Npt, tabV2d);
0763 V2d = tabV2d(k+1);
0764 T2 = V2d.X();
0765 }
0766 for (j = 1; j <= Npol; j++) {
0767 Daij = DDA( j);
0768 Cont(Inc3, j + Npol + IncCol) = Daij*T2;
0769 Cont(Inc3, j + IP + IncCol) = -Daij*T1;
0770 }
0771 IncCol = IncCol + IP;
0772
0773 }
0774 else if (Ibont(k+1, i) == 2) { // Relations entre T1y et T2y:
0775 if ((k+1) <= nb3d) {
0776 ToolLine::Tangency(SSP, Npt, tabV);
0777 V = tabV(k+1);
0778 T2 = V.Y();
0779 }
0780 else {
0781 ToolLine::Tangency(SSP, Npt, tabV2d);
0782 V2d = tabV2d(k+1);
0783 T2 = V2d.Y();
0784 }
0785 for (j = 1; j <= Npol; j++) {
0786 Daij = DDA(j);
0787 Cont(Inc3, j + Npol + IncCol) = Daij*T2;
0788 Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
0789 }
0790 IncCol = IncCol + IP;
0791
0792 }
0793 else if (Ibont(k+1,i)== 3) { // Relations entre T1y et T2z:
0794 ToolLine::Tangency(SSP, Npt, tabV);
0795 V = tabV(k+1);
0796 T2 = V.Z();
0797 for (j = 1; j <= Npol; j++) {
0798 Daij = DDA(j);
0799 Cont(Inc3, j + Npol +IncCol) = Daij*T2;
0800 Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
0801 }
0802 IncCol = IncCol + IP;
0803 }
0804 }
0805
0806 else {
0807 ToolLine::Tangency(SSP, Npt, tabV);
0808 V = tabV(k);
0809 T1 = V.Z();
0810 IP = 3*Npol;
0811 if (Ibont(k+1, i) == 1) { // Relations entre T1z et T2x:
0812 if ((k+1) <= nb3d) {
0813 ToolLine::Tangency(SSP, Npt, tabV);
0814 V = tabV(k+1);
0815 T2 = V.X();
0816 }
0817 else {
0818 ToolLine::Tangency(SSP, Npt, tabV2d);
0819 V2d = tabV2d(k+1);
0820 T2 = V2d.X();
0821 }
0822 for (j = 1; j <= Npol; j++) {
0823 Daij = DDA(j);
0824 Cont(Inc3, j + 2*Npol +IncCol) = Daij*T2;
0825 Cont(Inc3, j + IP + IncCol) = -Daij*T1;
0826 }
0827 IncCol = IncCol + IP;
0828 }
0829
0830 else if (Ibont(k+1, i) == 2) { // Relations entre T1z et T2y:
0831 if ((k+1) <= nb3d) {
0832 ToolLine::Tangency(SSP, Npt, tabV);
0833 V = tabV(k+1);
0834 T2 = V.Y();
0835 }
0836 else {
0837 ToolLine::Tangency(SSP, Npt, tabV2d);
0838 V2d = tabV2d(k+1);
0839 T2 = V2d.Y();
0840 }
0841 for (j = 1; j <= Npol; j++) {
0842 Daij = DDA(j);
0843 Cont(Inc3, j + 2*Npol +IncCol) = Daij*T2;
0844 Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
0845 }
0846 IncCol = IncCol + IP;
0847 }
0848
0849 else if (Ibont(k+1,i)==3) { // Relations entre T1z et T2z:
0850 ToolLine::Tangency(SSP, Npt, tabV);
0851 V = tabV(k+1);
0852 T2 = V.Z();
0853 for (j = 1; j <= Npol; j++) {
0854 Daij = DDA(j);
0855 Cont(Inc3, j + 2*Npol + IncCol) = Daij*T2;
0856 Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
0857 }
0858 IncCol = IncCol + IP;
0859 }
0860 }
0861 }
0862 }
0863
0864 return DeCont;
0865 }
