1 // Created on: 1999-02-15
2 // Created by: Igor FEOKTISTOV
3 // Copyright (c) 1999-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 #include <PLib_Base.hxx>
18 #include <PLib_JacobiPolynomial.hxx>
19 #include <PLib_HermitJacobi.hxx>
20 #include <TColStd_Array1OfReal.hxx>
21 #include <math_Vector.hxx>
23 void AppParCurves_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve,
24 const TColStd_Array1OfReal& Parameters,
25 const Standard_Real CBLONG,
26 FEmTool_Assembly& A ) const
29 Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
30 NbElm = Curve->NbElements(),
31 NbDim = Curve->Dimension();
33 TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg);
34 math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg),
35 V1((Standard_Real*)&G1(0), 0, MxDeg),
36 V2((Standard_Real*)&G2(0), 0, MxDeg);
38 Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1,
41 p0 = Parameters.Lower() - myFirstPoint,
42 curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim,
43 jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d;
44 Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
46 // Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints;
47 // Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints;
48 Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints;
49 Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints;
50 NBeg2d = Ng3d * myNbP3d;
51 // NgPC1 = NbConstr + myNbCurvPoints;
52 NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints;
57 TColStd_Array1OfReal& Intervals = Curve->Knots();
59 Standard_Real t, R1, R2;
61 Handle(PLib_Base) myBase = Curve->Base();
62 Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));
63 Standard_Integer Order = myHermitJacobi->NivConstr() + 1;
65 Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1;
67 A.NullifyConstraint();
71 for(i = 1; i <= NbConstr; i++) {
75 Point = myTypConstraints->Value(ipnt);
76 TypOfConstr = myTypConstraints->Value(ityp);
78 t = Parameters(p0 + Point);
80 for(el = curel; el <= NbElm; ) {
81 if( t <= Intervals(++el) ) {
88 UFirst = Intervals(curel); ULast = Intervals(curel + 1);
89 coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.;
93 if(TypOfConstr == 0) {
95 for(k = 1; k < Order; k++) {
96 mfact = Pow(coeff, k);
98 G0(k + Order) *= mfact;
101 else if(TypOfConstr == 1) {
102 myBase->D1(t, G0, G1);
103 for(k = 1; k < Order; k++) {
104 mfact = Pow(coeff, k);
106 G0(k + Order) *= mfact;
108 G1(k + Order) *= mfact;
111 for(k = 0; k <= MxDeg; k++) {
116 myBase->D2(t, G0, G1, G2);
117 for(k = 1; k < Order; k++) {
118 mfact = Pow(coeff, k);
120 G0(k + Order) *= mfact;
122 G1(k + Order) *= mfact;
124 G2(k + Order) *= mfact;
127 mfact1 = mfact / coeff;
128 for(k = 0; k <= MxDeg; k++) {
136 j = NbDim * (Point - myFirstPoint);
137 Standard_Integer n0 = NPass;
139 for(pnt = 1; pnt <= myNbP3d; pnt++) {
140 IndexOfConstraint = n0;
141 for(k = 1; k <= 3; k++) {
143 A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
144 IndexOfConstraint += NgPC1;
151 for(pnt = 1; pnt <= myNbP2d; pnt++) {
152 IndexOfConstraint = n0;
153 for(k = 1; k <= 2; k++) {
155 A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
156 IndexOfConstraint += NgPC1;
162 /* if(TypOfConstr == 1) {
164 IndexOfConstraint = NTang3d + 1;
165 jt = Ntheta * (i - 1);
166 for(pnt = 1; pnt <= myNbP3d; pnt++) {
167 for(k = 1; k <= 3; k++) {
168 A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
169 A.AddConstraint(IndexOfConstraint + 1, curel, k, myTtheta->Value(jt + 3 + k) * V1, 0.);
171 IndexOfConstraint += Ng3d;
175 IndexOfConstraint = NBeg2d + NTang2d + 1;
176 for(pnt = 1; pnt <= myNbP2d; pnt++) {
177 for(k = 1; k <= 2; k++) {
178 A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
180 IndexOfConstraint += Ng2d;
187 if(TypOfConstr == 1) {
191 j = 2 * NbDim * (i - 1);
193 for(pnt = 1; pnt <= myNbP3d; pnt++) {
194 IndexOfConstraint = n0;
195 for(k = 1; k <= 3; k++) {
197 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
198 IndexOfConstraint += NgPC1;
205 for(pnt = 1; pnt <= myNbP2d; pnt++) {
206 IndexOfConstraint = n0;
207 for(k = 1; k <= 2; k++) {
209 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
210 IndexOfConstraint += NgPC1;
216 if(TypOfConstr == 2) {
220 j = 2 * NbDim * (i - 1);
222 for(pnt = 1; pnt <= myNbP3d; pnt++) {
223 IndexOfConstraint = n0;
224 for(k = 1; k <= 3; k++) {
226 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
227 IndexOfConstraint += NgPC1;
234 for(pnt = 1; pnt <= myNbP2d; pnt++) {
235 IndexOfConstraint = n0;
236 for(k = 1; k <= 2; k++) {
238 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
239 IndexOfConstraint += NgPC1;
245 j = 2 * NbDim * (i - 1) + 3;
246 jt = Ntheta * (i - 1);
247 IndexOfConstraint = NTang3d + 1;
249 for(pnt = 1; pnt <= myNbP3d; pnt++) {
251 for(k = 1; k <= 3; k++) {
252 R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
253 R2 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + 3 + k);
255 R1 *= CBLONG * CBLONG;
256 R2 *= CBLONG * CBLONG;
257 for(k = 1; k <= 3; k++) {
259 if(k > 1) R1 = R2 = 0.;
260 A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
261 A.AddConstraint(IndexOfConstraint + 1, curel, curdim, myTfthet->Value(jt + 3 + k) * V2, R2);
263 IndexOfConstraint += Ng3d;
269 IndexOfConstraint = NBeg2d + NTang2d + 1;
270 for(pnt = 1; pnt <= myNbP2d; pnt++) {
272 for(k = 1; k <= 2; k++) {
273 R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
275 R1 *= CBLONG * CBLONG;
276 for(k = 1; k <= 2; k++) {
279 A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
281 IndexOfConstraint += Ng2d;