0024428: Implementation of LGPL license
[occt.git] / src / AppParCurves / AppParCurves_Variational_8.gxx
CommitLineData
b311480e 1// Created on: 1999-02-15
2// Created by: Igor FEOKTISTOV
3// Copyright (c) 1999-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
973c2be1 8// This library is free software; you can redistribute it and / or modify it
9// under the terms of the GNU Lesser General Public 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 16
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>
22
23void AppParCurves_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve,
24 const TColStd_Array1OfReal& Parameters,
25 const Standard_Real CBLONG,
26 FEmTool_Assembly& A ) const
27{
28
29 Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
30 NbElm = Curve->NbElements(),
31 NbDim = Curve->Dimension();
32
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);
37
38 Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1,
39 NTang3d, NTang2d,
40 Point, TypOfConstr,
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;
45
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;
53 NPass = 0;
54 NTang3d = 3 * NgPC1;
55 NTang2d = 2 * NgPC1;
56
57 TColStd_Array1OfReal& Intervals = Curve->Knots();
58
59 Standard_Real t, R1, R2;
60
61 Handle(PLib_Base) myBase = Curve->Base();
62 Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));
63 Standard_Integer Order = myHermitJacobi->NivConstr() + 1;
64
65 Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1;
66
67 A.NullifyConstraint();
68
69 ipnt = -1;
70 ityp = 0;
71 for(i = 1; i <= NbConstr; i++) {
72
73 ipnt += 2; ityp += 2;
74
75 Point = myTypConstraints->Value(ipnt);
76 TypOfConstr = myTypConstraints->Value(ityp);
77
78 t = Parameters(p0 + Point);
79
80 for(el = curel; el <= NbElm; ) {
81 if( t <= Intervals(++el) ) {
82 curel = el - 1;
83 break;
84 }
85 }
86
87
88 UFirst = Intervals(curel); ULast = Intervals(curel + 1);
89 coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.;
90
91 t = (t - c0) / coeff;
92
93 if(TypOfConstr == 0) {
94 myBase->D0(t, G0);
95 for(k = 1; k < Order; k++) {
96 mfact = Pow(coeff, k);
97 G0(k) *= mfact;
98 G0(k + Order) *= mfact;
99 }
100 }
101 else if(TypOfConstr == 1) {
102 myBase->D1(t, G0, G1);
103 for(k = 1; k < Order; k++) {
104 mfact = Pow(coeff, k);
105 G0(k) *= mfact;
106 G0(k + Order) *= mfact;
107 G1(k) *= mfact;
108 G1(k + Order) *= mfact;
109 }
110 mfact = 1./coeff;
111 for(k = 0; k <= MxDeg; k++) {
112 G1(k) *= mfact;
113 }
114 }
115 else {
116 myBase->D2(t, G0, G1, G2);
117 for(k = 1; k < Order; k++) {
118 mfact = Pow(coeff, k);
119 G0(k) *= mfact;
120 G0(k + Order) *= mfact;
121 G1(k) *= mfact;
122 G1(k + Order) *= mfact;
123 G2(k) *= mfact;
124 G2(k + Order) *= mfact;
125 }
126 mfact = 1. / coeff;
127 mfact1 = mfact / coeff;
128 for(k = 0; k <= MxDeg; k++) {
129 G1(k) *= mfact;
130 G2(k) *= mfact1;
131 }
132 }
133
134 NPass++;
135
136 j = NbDim * (Point - myFirstPoint);
137 Standard_Integer n0 = NPass;
138 curdim = 0;
139 for(pnt = 1; pnt <= myNbP3d; pnt++) {
140 IndexOfConstraint = n0;
141 for(k = 1; k <= 3; k++) {
142 curdim++;
143 A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
144 IndexOfConstraint += NgPC1;
145 }
146 j +=3;
147 n0 += Ng3d;
148 }
149
150 n0 = NPass + NBeg2d;
151 for(pnt = 1; pnt <= myNbP2d; pnt++) {
152 IndexOfConstraint = n0;
153 for(k = 1; k <= 2; k++) {
154 curdim++;
155 A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
156 IndexOfConstraint += NgPC1;
157 }
158 j +=2;
159 n0 += Ng2d;
160 }
161
162/* if(TypOfConstr == 1) {
163
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.);
170 }
171 IndexOfConstraint += Ng3d;
172 jt += 6;
173 }
174
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.);
179 }
180 IndexOfConstraint += Ng2d;
181 jt += 2;
182 }
183
184 NTang3d += 2;
185 NTang2d += 1;
186 } */
187 if(TypOfConstr == 1) {
188
189 NPass++;
190 n0 = NPass;
191 j = 2 * NbDim * (i - 1);
192 curdim = 0;
193 for(pnt = 1; pnt <= myNbP3d; pnt++) {
194 IndexOfConstraint = n0;
195 for(k = 1; k <= 3; k++) {
196 curdim++;
197 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
198 IndexOfConstraint += NgPC1;
199 }
200 n0 += Ng3d;
201 j += 6;
202 }
203
204 n0 = NPass + NBeg2d;
205 for(pnt = 1; pnt <= myNbP2d; pnt++) {
206 IndexOfConstraint = n0;
207 for(k = 1; k <= 2; k++) {
208 curdim++;
209 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
210 IndexOfConstraint += NgPC1;
211 }
212 n0 += Ng2d;
213 j += 4;
214 }
215 }
216 if(TypOfConstr == 2) {
217
218 NPass++;
219 n0 = NPass;
220 j = 2 * NbDim * (i - 1);
221 curdim = 0;
222 for(pnt = 1; pnt <= myNbP3d; pnt++) {
223 IndexOfConstraint = n0;
224 for(k = 1; k <= 3; k++) {
225 curdim++;
226 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
227 IndexOfConstraint += NgPC1;
228 }
229 n0 += Ng3d;
230 j += 6;
231 }
232
233 n0 = NPass + NBeg2d;
234 for(pnt = 1; pnt <= myNbP2d; pnt++) {
235 IndexOfConstraint = n0;
236 for(k = 1; k <= 2; k++) {
237 curdim++;
238 A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
239 IndexOfConstraint += NgPC1;
240 }
241 n0 += Ng2d;
242 j += 4;
243 }
244
245 j = 2 * NbDim * (i - 1) + 3;
246 jt = Ntheta * (i - 1);
247 IndexOfConstraint = NTang3d + 1;
248 curdim = 0;
249 for(pnt = 1; pnt <= myNbP3d; pnt++) {
250 R1 = 0.; R2 = 0.;
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);
254 }
255 R1 *= CBLONG * CBLONG;
256 R2 *= CBLONG * CBLONG;
257 for(k = 1; k <= 3; k++) {
258 curdim++;
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);
262 }
263 IndexOfConstraint += Ng3d;
264 j += 6;
265 jt += 6;
266 }
267
268 j--;
269 IndexOfConstraint = NBeg2d + NTang2d + 1;
270 for(pnt = 1; pnt <= myNbP2d; pnt++) {
271 R1 = 0.;
272 for(k = 1; k <= 2; k++) {
273 R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
274 }
275 R1 *= CBLONG * CBLONG;
276 for(k = 1; k <= 2; k++) {
277 curdim++;
278 if(k > 1) R1 = 0.;
279 A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
280 }
281 IndexOfConstraint += Ng2d;
282 j += 4;
283 jt += 2;
284 }
285
286 NTang3d += 2;
287 NTang2d += 1;
288 }
289
290 }
291
292}