OOFEM 3.0
Loading...
Searching...
No Matches
feinurbs.C
Go to the documentation of this file.
1/*
2 *
3 * ##### ##### ###### ###### ### ###
4 * ## ## ## ## ## ## ## ### ##
5 * ## ## ## ## #### #### ## # ##
6 * ## ## ## ## ## ## ## ##
7 * ## ## ## ## ## ## ## ##
8 * ##### ##### ## ###### ## ##
9 *
10 *
11 * OOFEM : Object Oriented Finite Element Code
12 *
13 * Copyright (C) 1993 - 2025 Borek Patzak
14 *
15 *
16 *
17 * Czech Technical University, Faculty of Civil Engineering,
18 * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19 *
20 * This library is free software; you can redistribute it and/or
21 * modify it under the terms of the GNU Lesser General Public
22 * License as published by the Free Software Foundation; either
23 * version 2.1 of the License, or (at your option) any later version.
24 *
25 * This program is distributed in the hope that it will be useful,
26 * but WITHOUT ANY WARRANTY; without even the implied warranty of
27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28 * Lesser General Public License for more details.
29 *
30 * You should have received a copy of the GNU Lesser General Public
31 * License along with this library; if not, write to the Free Software
32 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33 */
34
35#include "feinurbs.h"
36#include "floatarray.h"
37#include "floatmatrix.h"
38#include "iga.h"
39#include "parametermanager.h"
40#include "paramkey.h"
41
42namespace oofem {
43// optimized version of A4.4 for d=1
44#define OPTIMIZED_VERSION_A4dot4
45
46ParamKey NURBSInterpolation::IPK_NURBSInterpolation_weights("weights");
47
48void
49NURBSInterpolation :: initializeFrom(InputRecord &ir, ParameterManager&pm, int elnum, int priority)
50{
51 BSplineInterpolation :: initializeFrom(ir, pm, elnum, priority);
52 PM_UPDATE_PARAMETER(weights, pm, ir, elnum, IPK_NURBSInterpolation_weights, priority);
53}
54
55void
56NURBSInterpolation :: postInitialize(ParameterManager &pm, int elnum)
57{
58 BSplineInterpolation :: postInitialize(pm, elnum);
59 PM_ELEMENT_ERROR_IFNOTSET(pm, elnum, IPK_NURBSInterpolation_weights);
60}
61
62void NURBSInterpolation :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
63{
64 const FEIIGAElementGeometryWrapper &gw = static_cast< const FEIIGAElementGeometryWrapper& >(cellgeo);
65 IntArray span(nsd);
66 double sum = 0.0, val;
67 int count, c = 1;
68 std :: vector< FloatArray >N;
69
70 if ( gw.knotSpan ) {
71 span = * gw.knotSpan;
72 } else {
73 for ( int i = 0; i < nsd; i++ ) {
74 span[i] = this->findSpan(numberOfControlPoints [ i ], degree [ i ], lcoords[i], knotVector [ i ]);
75 }
76 }
77
78 for ( int i = 0; i < nsd; i++ ) {
79 this->basisFuns(N [ i ], span[i], lcoords[i], degree [ i ], knotVector [ i ]);
80 }
81
82 count = giveNumberOfKnotSpanBasisFunctions(span);
83 answer.resize(count);
84
85 if ( nsd == 1 ) {
86 int uind = span[0] - degree [ 0 ];
87 int ind = uind + 1;
88 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
89 answer.at(c++) = val = N [ 0 ][k] * weights.at(ind + k); // Nu*w
90 sum += val;
91 }
92 } else if ( nsd == 2 ) {
93 int uind = span[0] - degree [ 0 ];
94 int vind = span[1] - degree [ 1 ];
95 int ind = vind * numberOfControlPoints [ 0 ] + uind + 1;
96 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
97 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
98 answer.at(c++) = val = N [ 0 ][k] * N [ 1 ][l] * weights.at(ind + k); // Nu*Nv*w
99 sum += val;
100 }
101
102 ind += numberOfControlPoints [ 0 ];
103 }
104 } else if ( nsd == 3 ) {
105 int uind = span[0] - degree [ 0 ];
106 int vind = span[1] - degree [ 1 ];
107 int tind = span[2] - degree [ 2 ];
108 int ind = tind * numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ] + vind * numberOfControlPoints [ 0 ] + uind + 1;
109 for ( int m = 0; m <= degree [ 2 ]; m++ ) {
110 int indx = ind;
111 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
112 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
113 answer.at(c++) = val = N [ 0 ][k] * N [ 1 ][l] * N [ 2 ][m] * weights.at(ind + k); // Nu*Nv*Nt*w
114 sum += val;
115 }
116
117 ind += numberOfControlPoints [ 0 ];
118 }
119
120 ind = indx + numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ];
121 }
122 } else {
123 OOFEM_ERROR("not implemented for nsd = %d", nsd);
124 }
125
126 answer.times(1./sum);
127}
128
129
130double NURBSInterpolation :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
131{
132 const FEIIGAElementGeometryWrapper &gw = static_cast< const FEIIGAElementGeometryWrapper& >(cellgeo);
133 FloatMatrix jacobian(nsd, nsd);
134 IntArray span(nsd);
135 double Jacob = 0.;
136 int count;
137 std :: vector< FloatArray > N(nsd);
138 std :: vector< FloatMatrix > ders(nsd);
139
140 if ( gw.knotSpan ) {
141 span = * gw.knotSpan;
142 } else {
143 for ( int i = 0; i < nsd; i++ ) {
144 span[i] = this->findSpan(numberOfControlPoints [ i ], degree [ i ], lcoords[i], knotVector [ i ]);
145 }
146 }
147
148 for ( int i = 0; i < nsd; i++ ) {
149 this->dersBasisFuns(1, lcoords[i], span[i], degree [ i ], knotVector [ i ], ders [ i ]);
150 }
151
152 count = giveNumberOfKnotSpanBasisFunctions(span);
153 answer.resize(count, nsd);
154
155#if 0 // code according NURBS book (too general allowing higher derivatives)
156 if ( nsd == 2 ) {
157 FloatArray tmp1(nsd + 1), tmp2(nsd + 1); // allow for weight
158
159 FloatMatrix Aders [ 2 ]; // derivatives in each coordinate direction on BSpline
160 #ifndef OPTIMIZED_VERSION_A4dot4
161 FloatMatrix Sders [ 2 ]; // derivatives in each coordinate direction on NURBS
162 #endif
163 FloatMatrix wders; // derivatives in w direction on BSpline
164 /*
165 * IntArray Bin(2,2); // binomial coefficients from 0 to d=1
166 * // Bin(n,k)=(n above k)=n!/k!(n-k)! for n>=k>=0
167 * // lower triangle corresponds to Pascal triangle
168 * // according to A4.4 it seems that only coefficients in lower triangle except the first column are used
169 */
170 // resizing to (2,2) has nothing common with nsd
171 // it is related to the fact that 0th and 1st derivatives are computed in each direction
172 for ( int i = 0; i < nsd; i++ ) {
173 Aders [ i ].resize(2, 2);
174 Aders [ i ].zero();
175 #ifndef OPTIMIZED_VERSION_A4dot4
176 Sders [ i ].resize(2, 2);
177 #endif
178 }
179
180 wders.resize(2, 2);
181 wders.zero();
182
183 // calculation of jacobian matrix according to A4.4
184 // calculate values and derivatives of nonrational Bspline surface with weights at first (Aders, wders)
185 int uind = span[0] - degree [ 0 ];
186 int vind = span[1] - degree [ 1 ];
187 int ind = vind * numberOfControlPoints [ 0 ] + uind + 1;
188 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
189 tmp1.zero();
190 tmp2.zero();
191 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
192 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
193 double w = this->weights.at(ind + k);
194
195 tmp1[0] += ders [ 0 ](0, k) * vertexCoords[0] * w; // sum(Nu*x*w)
196 tmp1[1] += ders [ 0 ](0, k) * vertexCoords[1] * w; // sum(Nu*y*w)
197 tmp1[2] += ders [ 0 ](0, k) * w; // sum(Nu*w)
198
199 tmp2[0] += ders [ 0 ](1, k) * vertexCoords[0] * w; // sum(dNu/du*x*w)
200 tmp2[1] += ders [ 0 ](1, k) * vertexCoords[1] * w; // sum(dNu/du*y*w)
201 tmp2[2] += ders [ 0 ](1, k) * w; // sum(dNu/du*w)
202 }
203
204 ind += numberOfControlPoints [ 0 ];
205
206 Aders [ 0 ](0, 0) += ders [ 1 ](0, l) * tmp1[0]; // xw=sum(Nv*sum(Nu*x*w))
207 Aders [ 1 ](0, 0) += ders [ 1 ](0, l) * tmp1[1]; // yw=sum(Nv*sum(Nu*y*w))
208 wders(0, 0) += ders [ 1 ](0, l) * tmp1[2]; // w=sum(Nv*sum(Nu*w))
209
210 Aders [ 0 ](0, 1) += ders [ 1 ](1, l) * tmp1[0]; // dxw/dv=sum(dNv/dv*sum(Nu*x*w))
211 Aders [ 1 ](0, 1) += ders [ 1 ](1, l) * tmp1[1]; // dyw/dv=sum(dNv/dv*sum(Nu*y*w))
212 wders(0, 1) += ders [ 1 ](1, l) * tmp1[2]; // dw/dv=sum(dNv/dv*sum(Nu*w))
213
214 Aders [ 0 ](1, 0) += ders [ 1 ](0, l) * tmp2[0]; // dxw/du=sum(Nv*sum(dNu/du*x*w))
215 Aders [ 1 ](1, 0) += ders [ 1 ](0, l) * tmp2[1]; // dyw/du=sum(Nv*sum(dNu/du*y*w))
216 wders(1, 0) += ders [ 1 ](0, l) * tmp2[2]; // dw/du=sum(Nv*sum(dNu/du*w))
217 }
218
219 double weight = wders(0, 0);
220
221 #ifndef OPTIMIZED_VERSION_A4dot4
222 const int d = 1;
223 // calculate values and derivatives of NURBS surface (A4.4)
224 // since all entries in Pascal triangle up to d=1 are 1, binomial coefficients are ignored
225 for ( int k = 0; k <= d; k++ ) {
226 for ( int l = 0; l <= d - k; l++ ) {
227 tmp1[0] = Aders [ 0 ](k, l);
228 tmp1[1] = Aders [ 1 ](k, l);
229 for ( j = 1; j <= l; j++ ) {
230 tmp1[0] -= wders(0, j) * Sders [ 0 ](k, l - j); // *Bin(l,j)
231 tmp1[1] -= wders(0, j) * Sders [ 1 ](k, l - j); // *Bin(l,j)
232 }
233
234 for ( int i = 1; i <= k; i++ ) {
235 tmp1[0] -= wders(i, 0) * Sders [ 0 ](k - i, l); // *Bin(k,i)
236 tmp1[1] -= wders(i, 0) * Sders [ 1 ](k - i, l); // *Bin(k,i)
237 tmp2.zero();
238 for ( int j = 1; j <= l; j++ ) {
239 tmp2[0] += wders(i, j) * Sders [ 0 ](k - i, l - j); // *Bin(l,j)
240 tmp2[1] += wders(i, j) * Sders [ 1 ](k - i, l - j); // *Bin(l,j)
241 }
242
243 tmp1[0] -= tmp2[0]; // *Bin(k,i)
244 tmp1[1] -= tmp2[1]; // *Bin(k,i)
245 }
246
247 Sders [ 0 ](k, l) = tmp1[0] / weight;
248 Sders [ 1 ](k, l) = tmp1[1] / weight;
249 }
250 }
251
252 jacobian(0, 0) = Sders [ 0 ](1, 0); // dx/du
253 jacobian(0, 1) = Sders [ 1 ](1, 0); // dy/du
254 jacobian(1, 0) = Sders [ 0 ](0, 1); // dx/dv
255 jacobian(1, 1) = Sders [ 1 ](0, 1); // dy/dv
256 #else
257 // optimized version of A4.4 for d=1, binomial coefficients ignored
258
259 /*
260 * // k=0 l=0 loop
261 * Sders[0](0,0) = Aders[0](0,0) / weight;
262 * Sders[1](0,0) = Aders[1](0,0) / weight;
263 * // k=1 l=0 loop
264 * Sders[0](1,0) = (Aders[0](1,0)-wders(1,0)*Sders[0](0,0)) / weight;
265 * Sders[1](1,0) = (Aders[1](1,0)-wders(1,0)*Sders[1](0,0)) / weight;
266 * // k=0 l=1 loop
267 * Sders[0](0,1) = (Aders[0](0,1)-wders(0,1)*Sders[0](0,0)) / weight;
268 * Sders[1](0,1) = (Aders[1](0,1)-wders(0,1)*Sders[1](0,0)) / weight;
269 *
270 * jacobian(0,0) = Sders[0](1,0); // dx/du
271 * jacobian(0,1) = Sders[1](1,0); // dy/du
272 * jacobian(1,0) = Sders[0](0,1); // dx/dv
273 * jacobian(1,1) = Sders[1](0,1); // dy/dv
274 */
275
276 // k=0 l=0 loop
277 tmp1[0] = Aders [ 0 ](0, 0) / weight;
278 tmp1[1] = Aders [ 1 ](0, 0) / weight;
279 // k=1 l=0 loop
280 jacobian(0, 0) = ( Aders [ 0 ](1, 0) - wders(1, 0) * tmp1[0] ) / weight; // dx/du
281 jacobian(0, 1) = ( Aders [ 1 ](1, 0) - wders(1, 0) * tmp1[1] ) / weight; // dy/du
282 // k=0 l=1 loop
283 jacobian(1, 0) = ( Aders [ 0 ](0, 1) - wders(0, 1) * tmp1[0] ) / weight; // dx/dv
284 jacobian(1, 1) = ( Aders [ 1 ](0, 1) - wders(0, 1) * tmp1[1] ) / weight; // dy/dv
285 #endif
286
287 Jacob = jacobian.giveDeterminant();
288
289 //calculation of derivatives of NURBS basis functions with respect to local parameters is not covered by NURBS book
290 double product = Jacob * weight * weight;
291 int cnt = 0;
292 ind = vind * numberOfControlPoints [ 0 ] + uind + 1;
293 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
294 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
295 double w = weights.at(ind + k);
296 // dNu/du*Nv*w*sum(Nv*Nu*w) - Nu*Nv*w*sum(dNu/du*Nv*w)
297 tmp1[0] = ders [ 0 ](1, k) * ders [ 1 ](0, l) * w * weight - ders [ 0 ](0, k) * ders [ 1 ](0, l) * w * wders(1, 0);
298 // Nu*dNv/dv*w*sum(Nv*Nu*w) - Nu*Nv*w*sum(Nu*dNv/dv*w)
299 tmp1[1] = ders [ 0 ](0, k) * ders [ 1 ](1, l) * w * weight - ders [ 0 ](0, k) * ders [ 1 ](0, l) * w * wders(0, 1);
300
301 answer(cnt, 0) = ( +jacobian(1, 1) * tmp1[0] - jacobian(0, 1) * tmp1[1] ) / product;
302 answer(cnt, 1) = ( -jacobian(1, 0) * tmp1[0] + jacobian(0, 0) * tmp1[1] ) / product;
303 cnt++;
304 }
305
306 ind += numberOfControlPoints [ 0 ];
307 }
308 } else {
309 OOFEM_ERROR("not implemented for nsd = %d", nsd);
310 }
311
312#else
313 std :: vector< FloatArray > Aders(nsd);
314 FloatArray wders; // 0th and 1st derivatives in w direction on BSpline
315
316 for ( int i = 0; i < nsd; i++ ) {
317 Aders [ i ].resize(nsd + 1);
318 Aders [ i ].zero();
319 }
320
321 wders.resize(nsd + 1);
322 wders.zero();
323
324 if ( nsd == 1 ) {
325 // calculate values and derivatives of nonrational Bspline curve with weights at first (Aders, wders)
326 int uind = span[0] - degree [ 0 ];
327 int ind = uind + 1;
328 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
329 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
330 double w = this->weights.at(ind + k);
331
332 Aders [ 0 ][0] += ders [ 0 ](0, k) * vertexCoords[0] * w; // xw=sum(Nu*x*w)
333 wders[0] += ders [ 0 ](0, k) * w; // w=sum(Nu*w)
334
335 Aders [ 0 ][1] += ders [ 0 ](1, k) * vertexCoords[0] * w; // dxw/du=sum(dNu/du*x*w)
336 wders[1] += ders [ 0 ](1, k) * w; // dw/du=sum(dNu/du*w)
337 }
338
339 double weight = wders[0];
340
341 // calculation of jacobian matrix according to Eq 4.7
342 jacobian(0, 0) = ( Aders [ 0 ][1] - wders[1] * Aders [ 0 ][0] / weight ) / weight; // dx/du
343
344 Jacob = jacobian.giveDeterminant();
345
346 //calculation of derivatives of NURBS basis functions with respect to local parameters is not covered by NURBS book
347 double product = Jacob * weight * weight;
348 int cnt = 0;
349 ind = uind + 1;
350 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
351 double w = this->weights.at(ind + k);
352 // [dNu/du*w*sum(Nu*w) - Nu*w*sum(dNu/du*w)] / [J*sum(Nu*w)^2]
353 answer(cnt, 0) = ders [ 0 ](1, k) * w * weight - ders [ 0 ](0, k) * w * wders[1] / product;
354 cnt++;
355 }
356 } else if ( nsd == 2 ) {
357 FloatArray tmp1(nsd + 1), tmp2(nsd + 1); // allow for weight
358
359 // calculate values and derivatives of nonrational Bspline surface with weights at first (Aders, wders)
360 int uind = span[0] - degree [ 0 ];
361 int vind = span[1] - degree [ 1 ];
362 int ind = vind * numberOfControlPoints [ 0 ] + uind + 1;
363 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
364 tmp1.zero();
365 tmp2.zero();
366 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
367 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
368 double w = this->weights.at(ind + k);
369
370 tmp1[0] += ders [ 0 ](0, k) * vertexCoords[0] * w; // sum(Nu*x*w)
371 tmp1[1] += ders [ 0 ](0, k) * vertexCoords[1] * w; // sum(Nu*y*w)
372 tmp1[2] += ders [ 0 ](0, k) * w; // sum(Nu*w)
373
374 tmp2[0] += ders [ 0 ](1, k) * vertexCoords[0] * w; // sum(dNu/du*x*w)
375 tmp2[1] += ders [ 0 ](1, k) * vertexCoords[1] * w; // sum(dNu/du*y*w)
376 tmp2[2] += ders [ 0 ](1, k) * w; // sum(dNu/du*w)
377 }
378
379 ind += numberOfControlPoints [ 0 ];
380
381 Aders [ 0 ][0] += ders [ 1 ](0, l) * tmp1[0]; // xw=sum(Nv*sum(Nu*x*w)
382 Aders [ 1 ][0] += ders [ 1 ](0, l) * tmp1[1]; // yw=sum(Nv*sum(Nu*y*w)
383 wders[0] += ders [ 1 ](0, l) * tmp1[2]; // w=sum(Nv*sum(Nu*w)
384
385 Aders [ 0 ][1] += ders [ 1 ](0, l) * tmp2[0]; // dxw/du=sum(Nv*sum(dNu/du*x*w)
386 Aders [ 1 ][1] += ders [ 1 ](0, l) * tmp2[1]; // dyw/du=sum(Nv*sum(dNu/du*y*w)
387 wders[1] += ders [ 1 ](0, l) * tmp2[2]; // dw/du=sum(Nv*sum(dNu/du*w)
388
389 Aders [ 0 ][2] += ders [ 1 ](1, l) * tmp1[0]; // dxw/dv=sum(dNv/dv*sum(Nu*x*w)
390 Aders [ 1 ][2] += ders [ 1 ](1, l) * tmp1[1]; // dyw/dv=sum(dNv/dv*sum(Nu*y*w)
391 wders[2] += ders [ 1 ](1, l) * tmp1[2]; // dw/dv=sum(dNv/dv*sum(Nu*w)
392 }
393
394 double weight = wders[0];
395
396 // calculation of jacobian matrix according to Eq 4.19
397 tmp1[0] = Aders [ 0 ][0] / weight;
398 tmp1[1] = Aders [ 1 ][0] / weight;
399 jacobian(0, 0) = ( Aders [ 0 ][1] - wders[1] * tmp1[0] ) / weight; // dx/du
400 jacobian(0, 1) = ( Aders [ 1 ][1] - wders[1] * tmp1[1] ) / weight; // dy/du
401 jacobian(1, 0) = ( Aders [ 0 ][2] - wders[2] * tmp1[0] ) / weight; // dx/dv
402 jacobian(1, 1) = ( Aders [ 1 ][2] - wders[2] * tmp1[1] ) / weight; // dy/dv
403
404 Jacob = jacobian.giveDeterminant();
405
406 //calculation of derivatives of NURBS basis functions with respect to local parameters is not covered by NURBS book
407 double product = Jacob * weight * weight;
408 int cnt = 0;
409 ind = vind * numberOfControlPoints [ 0 ] + uind + 1;
410 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
411 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
412 double w = this->weights.at(ind + k);
413 // dNu/du*Nv*w*sum(Nu*Nv*w) - Nu*Nv*w*sum(dNu/du*Nv*w)
414 tmp1[0] = ders [ 0 ](1, k) * ders [ 1 ](0, l) * w * weight - ders [ 0 ](0, k) * ders [ 1 ](0, l) * w * wders[1];
415 // Nu*dNv/dv*w*sum(Nu*Nv*w) - Nu*Nv*w*sum(Nu*dNv/dv*w)
416 tmp1[1] = ders [ 0 ](0, k) * ders [ 1 ](1, l) * w * weight - ders [ 0 ](0, k) * ders [ 1 ](0, l) * w * wders[2];
417
418 answer(cnt, 0) = ( +jacobian(1, 1) * tmp1[0] - jacobian(0, 1) * tmp1[1] ) / product;
419 answer(cnt, 1) = ( -jacobian(1, 0) * tmp1[0] + jacobian(0, 0) * tmp1[1] ) / product;
420 cnt++;
421 }
422
423 ind += numberOfControlPoints [ 0 ];
424 }
425 } else if ( nsd == 3 ) {
426 FloatArray tmp1(nsd + 1), tmp2(nsd + 1); // allow for weight
427 FloatArray temp1(nsd + 1), temp2(nsd + 1), temp3(nsd + 1); // allow for weight
428
429 // calculate values and derivatives of nonrational Bspline solid with weights at first (Aders, wders)
430 int uind = span[0] - degree [ 0 ];
431 int vind = span[1] - degree [ 1 ];
432 int tind = span[2] - degree [ 2 ];
433 int ind = tind * numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ] + vind * numberOfControlPoints [ 0 ] + uind + 1;
434 for ( int m = 0; m <= degree [ 2 ]; m++ ) {
435 temp1.zero();
436 temp2.zero();
437 temp3.zero();
438 int indx = ind;
439 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
440 tmp1.zero();
441 tmp2.zero();
442 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
443 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
444 double w = this->weights.at(ind + k);
445
446 tmp1[0] += ders [ 0 ](0, k) * vertexCoords[0] * w; // sum(Nu*x*w)
447 tmp1[1] += ders [ 0 ](0, k) * vertexCoords[1] * w; // sum(Nu*y*w)
448 tmp1[2] += ders [ 0 ](0, k) * vertexCoords[2] * w; // sum(Nu*z*w)
449 tmp1[3] += ders [ 0 ](0, k) * w; // sum(Nu*w)
450
451 tmp2[0] += ders [ 0 ](1, k) * vertexCoords[0] * w; // sum(dNu/du*x*w)
452 tmp2[1] += ders [ 0 ](1, k) * vertexCoords[1] * w; // sum(dNu/du*y*w)
453 tmp2[2] += ders [ 0 ](1, k) * vertexCoords[2] * w; // sum(dNu/du*z*w)
454 tmp2[3] += ders [ 0 ](1, k) * w; // sum(dNu/du*w)
455 }
456
457 ind += numberOfControlPoints [ 0 ];
458
459 temp1[0] += ders [ 1 ](0, l) * tmp1[0]; // sum(Nv*sum(Nu*x*w))
460 temp1[1] += ders [ 1 ](0, l) * tmp1[1]; // sum(Nv*sum(Nu*y*w))
461 temp1[2] += ders [ 1 ](0, l) * tmp1[2]; // sum(Nv*sum(Nu*z*w))
462 temp1[3] += ders [ 1 ](0, l) * tmp1[3]; // sum(Nv*sum(Nu*w))
463
464 temp2[0] += ders [ 1 ](0, l) * tmp2[0]; // sum(Nv*sum(dNu/du*x*w))
465 temp2[1] += ders [ 1 ](0, l) * tmp2[1]; // sum(Nv*sum(dNu/du*y*w))
466 temp2[2] += ders [ 1 ](0, l) * tmp2[2]; // sum(Nv*sum(dNu/du*z*w))
467 temp2[3] += ders [ 1 ](0, l) * tmp2[3]; // sum(Nv*sum(dNu/du*w))
468
469 temp3[0] += ders [ 1 ](1, l) * tmp1[0]; // sum(dNv/dv*sum(Nu*x*w))
470 temp3[1] += ders [ 1 ](1, l) * tmp1[1]; // sum(dNv/dv*sum(Nu*y*w))
471 temp3[2] += ders [ 1 ](1, l) * tmp1[2]; // sum(dNv/dv*sum(Nu*z*w))
472 temp3[3] += ders [ 1 ](1, l) * tmp1[3]; // sum(dNv/dv*sum(Nu*w))
473 }
474
475 ind = indx + numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ];
476
477 Aders [ 0 ][0] += ders [ 2 ](0, m) * temp1[0]; // x=sum(Nt*sum(Nv*sum(Nu*x*w)))
478 Aders [ 1 ][0] += ders [ 2 ](0, m) * temp1[1]; // y=sum(Nt*sum(Nv*sum(Nu*y*w)))
479 Aders [ 2 ][0] += ders [ 2 ](0, m) * temp1[2]; // y=sum(Nt*sum(Nv*sum(Nu*y*w)))
480 wders[0] += ders [ 2 ](0, m) * temp1[3]; // w=sum(Nt*sum(Nv*sum(Nu*w)))
481
482 Aders [ 0 ][1] += ders [ 2 ](0, m) * temp2[0]; // dx/du=sum(Nt*sum(Nv*sum(dNu/du*x*w)))
483 Aders [ 1 ][1] += ders [ 2 ](0, m) * temp2[1]; // dy/du=sum(Nt*sum(Nv*sum(dNu/du*y*w)))
484 Aders [ 2 ][1] += ders [ 2 ](0, m) * temp2[2]; // dy/du=sum(Nt*sum(Nv*sum(dNu/du*y*w)))
485 wders[1] += ders [ 2 ](0, m) * temp2[3]; // dw/du=sum(Nt*sum(Nv*sum(dNu/du*w)))
486
487 Aders [ 0 ][2] += ders [ 2 ](0, m) * temp3[0]; // dx/dv=sum(Nt*sum(dNv/dv*sum(Nu*x*w)))
488 Aders [ 1 ][2] += ders [ 2 ](0, m) * temp3[1]; // dy/dv=sum(Nt*sum(dNv/dv*sum(Nu*y*w)))
489 Aders [ 2 ][2] += ders [ 2 ](0, m) * temp3[2]; // dy/dv=sum(Nt*sum(dNv/dv*sum(Nu*y*w)))
490 wders[2] += ders [ 2 ](0, m) * temp3[3]; // dw/dv=sum(Nt*sum(dNv/dv*sum(Nu*w)))
491
492 Aders [ 0 ][3] += ders [ 2 ](1, m) * temp1[0]; // dx/dt=sum(dNt/dt*sum(Nv*sum(Nu*x*w)))
493 Aders [ 1 ][3] += ders [ 2 ](1, m) * temp1[1]; // dy/dt=sum(dNt/dt*sum(Nv*sum(Nu*y*w)))
494 Aders [ 2 ][3] += ders [ 2 ](1, m) * temp1[2]; // dy/dt=sum(dNt/dt*sum(Nv*sum(Nu*y*w)))
495 wders[3] += ders [ 2 ](1, m) * temp1[3]; // dw/dt=sum(dNt/dt*sum(Nv*sum(Nu*w)))
496 }
497
498 double weight = wders[0];
499
500 // calculation of jacobian matrix
501 tmp1[0] = Aders [ 0 ][0] / weight;
502 tmp1[1] = Aders [ 1 ][0] / weight;
503 tmp1[2] = Aders [ 2 ][0] / weight;
504 jacobian(0, 0) = ( Aders [ 0 ][1] - wders[1] * tmp1[0] ) / weight; // dx/du
505 jacobian(0, 1) = ( Aders [ 1 ][1] - wders[1] * tmp1[1] ) / weight; // dy/du
506 jacobian(0, 2) = ( Aders [ 2 ][1] - wders[1] * tmp1[2] ) / weight; // dz/du
507 jacobian(1, 0) = ( Aders [ 0 ][2] - wders[2] * tmp1[0] ) / weight; // dx/dv
508 jacobian(1, 1) = ( Aders [ 1 ][2] - wders[2] * tmp1[1] ) / weight; // dy/dv
509 jacobian(1, 2) = ( Aders [ 2 ][2] - wders[2] * tmp1[2] ) / weight; // dz/dv
510 jacobian(2, 0) = ( Aders [ 0 ][3] - wders[3] * tmp1[0] ) / weight; // dx/dt
511 jacobian(2, 1) = ( Aders [ 1 ][3] - wders[3] * tmp1[1] ) / weight; // dy/dt
512 jacobian(2, 2) = ( Aders [ 2 ][3] - wders[3] * tmp1[2] ) / weight; // dz/dt
513
514 Jacob = jacobian.giveDeterminant();
515
516 //calculation of derivatives of NURBS basis functions with respect to local parameters is not covered by NURBS book
517 double product = Jacob * weight * weight;
518 int cnt = 0;
519 ind = tind * numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ] + vind * numberOfControlPoints [ 0 ] + uind + 1;
520 for ( int m = 0; m <= degree [ 2 ]; m++ ) {
521 int indx = ind;
522 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
523 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
524 double w = this->weights.at(ind + k);
525 // dNu/du*Nv*Nt*w*sum(Nu*Nv*Nt*w) - Nu*Nv*Nt*w*sum(dNu/du*Nv*Nt*w)
526 tmp1[0] = ders [ 0 ](1, k) * ders [ 1 ](0, l) * ders [ 2 ](0, m) * w * weight - ders [ 0 ](0, k) * ders [ 1 ](0, l) * ders [ 2 ](0, m) * w * wders[1];
527 // Nu*dNv/dv*Nt*w*sum(Nu*Nv*Nt*w) - Nu*Nv*Nt*w*sum(Nu*dNv/dv*Nt*w)
528 tmp1[1] = ders [ 0 ](0, k) * ders [ 1 ](1, l) * ders [ 2 ](0, m) * w * weight - ders [ 0 ](0, k) * ders [ 1 ](0, l) * ders [ 2 ](0, m) * w * wders[2];
529 // Nu*Nv*dNt/dt*w*sum(Nu*Nv*Nt*w) - Nu*Nv*Nt*w*sum(Nu*Nv*dNt/dt*w)
530 tmp1[2] = ders [ 0 ](0, k) * ders [ 1 ](0, l) * ders [ 2 ](1, m) * w * weight - ders [ 0 ](0, k) * ders [ 1 ](0, l) * ders [ 2 ](0, m) * w * wders[3];
531
532 answer(cnt, 0) = ( ( jacobian(1, 1) * jacobian(2, 2) - jacobian(1, 2) * jacobian(2, 1) ) * tmp1[0] +
533 ( jacobian(0, 2) * jacobian(2, 1) - jacobian(0, 1) * jacobian(2, 2) ) * tmp1[1] +
534 ( jacobian(0, 1) * jacobian(1, 2) - jacobian(0, 2) * jacobian(1, 1) ) * tmp1[2] ) / product; // dN/dx
535 answer(cnt, 1) = ( ( jacobian(1, 2) * jacobian(2, 0) - jacobian(1, 0) * jacobian(2, 2) ) * tmp1[0] +
536 ( jacobian(0, 0) * jacobian(2, 2) - jacobian(0, 2) * jacobian(2, 0) ) * tmp1[1] +
537 ( jacobian(0, 2) * jacobian(1, 0) - jacobian(0, 0) * jacobian(1, 2) ) * tmp1[2] ) / product; // dN/dy
538 answer(cnt, 2) = ( ( jacobian(1, 0) * jacobian(2, 1) - jacobian(1, 1) * jacobian(2, 0) ) * tmp1[0] +
539 ( jacobian(0, 1) * jacobian(2, 0) - jacobian(0, 0) * jacobian(2, 1) ) * tmp1[1] +
540 ( jacobian(0, 0) * jacobian(1, 1) - jacobian(0, 1) * jacobian(1, 0) ) * tmp1[2] ) / product; // dN/dz
541 cnt++;
542 }
543
544 ind += numberOfControlPoints [ 0 ];
545 }
546
547 ind = indx + numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ];
548 }
549 } else {
550 OOFEM_ERROR("not implemented for nsd = %d", nsd);
551 }
552
553#endif
554
555 return Jacob;
556}
557
558
559void NURBSInterpolation :: local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
560{
561 /* Based on SurfacePoint A4.3 implementation*/
562 const FEIIGAElementGeometryWrapper &gw = static_cast< const FEIIGAElementGeometryWrapper& >(cellgeo);
563 IntArray span(nsd);
564 double weight = 0.0;
565 std :: vector< FloatArray > N(nsd);
566
567 if ( gw.knotSpan ) {
568 span = * gw.knotSpan;
569 } else {
570 for ( int i = 0; i < nsd; i++ ) {
571 span[i] = this->findSpan(numberOfControlPoints [ i ], degree [ i ], lcoords[i], knotVector [ i ]);
572 }
573 }
574
575 for ( int i = 0; i < nsd; i++ ) {
576 this->basisFuns(N [ i ], span[i], lcoords[i], degree [ i ], knotVector [ i ]);
577 }
578
579 answer.resize(nsd);
580 answer.zero();
581
582 if ( nsd == 1 ) {
583 int uind = span[0] - degree [ 0 ];
584 int ind = uind + 1;
585 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
586 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
587 double w = this->weights.at(ind + k);
588 answer[0] += N [ 0 ][k] * vertexCoords[0] * w; // xw=sum(Nu*x*w)
589 weight += N [ 0 ][k] * w; // w=sum(Nu*w)
590 }
591 } else if ( nsd == 2 ) {
592 FloatArray tmp(nsd + 1); // allow for weight
593
594 int uind = span[0] - degree [ 0 ];
595 int vind = span[1] - degree [ 1 ];
596 int ind = vind * numberOfControlPoints [ 0 ] + uind + 1;
597 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
598 tmp.zero();
599 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
600 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
601 double w = this->weights.at(ind + k);
602 tmp[0] += N [ 0 ][k] * vertexCoords[0] * w; // sum(Nu*x*w)
603 tmp[1] += N [ 0 ][k] * vertexCoords[1] * w; // sum(Nu*y*w)
604 tmp[2] += N [ 0 ][k] * w; // sum(Nu*w)
605 }
606
607 ind += numberOfControlPoints [ 0 ];
608 answer[0] += N [ 1 ][l] * tmp[0]; // xw=sum(Nv*Nu*x*w)
609 answer[1] += N [ 1 ][l] * tmp[1]; // yw=sum(Nv*Nu*y*w)
610 weight += N [ 1 ][l] * tmp[2]; // w=sum(Nv*Nu*w)
611 }
612 } else if ( nsd == 3 ) {
613 FloatArray tmp(nsd + 1), temp(nsd + 1); // allow for weight
614
615 int uind = span[0] - degree [ 0 ];
616 int vind = span[1] - degree [ 1 ];
617 int tind = span[2] - degree [ 2 ];
618 int ind = tind * numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ] + vind * numberOfControlPoints [ 0 ] + uind + 1;
619 for ( int m = 0; m <= degree [ 2 ]; m++ ) {
620 temp.zero();
621 int indx = ind;
622 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
623 tmp.zero();
624 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
625 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
626 double w = this->weights.at(ind + k);
627 tmp[0] += N [ 0 ][k] * vertexCoords[0] * w; // sum(Nu*x*w)
628 tmp[1] += N [ 0 ][k] * vertexCoords[1] * w; // sum(Nu*y*w)
629 tmp[2] += N [ 0 ][k] * vertexCoords[2] * w; // sum(Nu*z*w)
630 tmp[3] += N [ 0 ][k] * w; // sum(Nu*w)
631 }
632
633 ind += numberOfControlPoints [ 0 ];
634
635 temp[0] += N [ 1 ][l] * tmp[0]; // sum(Nv*Nu*x*w)
636 temp[1] += N [ 1 ][l] * tmp[1]; // sum(Nv*Nu*y*w)
637 temp[2] += N [ 1 ][l] * tmp[2]; // sum(Nv*Nu*z*w)
638 temp[3] += N [ 1 ][l] * tmp[3]; // sum(Nv*Nu*w)
639 }
640
641 ind = indx + numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ];
642
643 answer[0] += N [ 2 ][m] * temp[0]; // xw=sum(Nv*Nu*Nt*x*w)
644 answer[1] += N [ 2 ][m] * temp[1]; // yw=sum(Nv*Nu*Nt*y*w)
645 answer[2] += N [ 2 ][m] * temp[2]; // zw=sum(Nv*Nu*Nt*z*w)
646 weight += N [ 2 ][m] * temp[3]; // w=sum(Nv*Nu*Nt*w)
647 }
648 } else {
649 OOFEM_ERROR("lnot implemented for nsd = %d", nsd);
650 }
651
652 answer.times(1.0 / weight);
653}
654
655
656void NURBSInterpolation :: giveJacobianMatrixAt(FloatMatrix &jacobian, const FloatArray &lcoords, const FEICellGeometry &cellgeo) const
657{
658 //
659 // Based on Algorithm A4.4 (p. 137) for d=1
660 //
661 const FEIIGAElementGeometryWrapper &gw = static_cast< const FEIIGAElementGeometryWrapper& >(cellgeo);
662 IntArray span(nsd);
663 std :: vector< FloatMatrix > ders(nsd);
664 jacobian.resize(nsd, nsd);
665
666 if ( gw.knotSpan ) {
667 span = * gw.knotSpan;
668 } else {
669 for ( int i = 0; i < nsd; i++ ) {
670 span[i] = this->findSpan(numberOfControlPoints [ i ], degree [ i ], lcoords[i], knotVector [ i ]);
671 }
672 }
673
674 for ( int i = 0; i < nsd; i++ ) {
675 this->dersBasisFuns(1, lcoords[i], span[i], degree [ i ], knotVector [ i ], ders [ i ]);
676 }
677
678#if 0 // code according NURBS book (too general allowing higher derivatives)
679 if ( nsd == 2 ) {
680 FloatArray tmp1(nsd + 1), tmp2(nsd + 1); // allow for weight
681
682 FloatMatrix Aders [ 2 ]; // derivatives in each coordinate direction on BSpline
683 #ifndef OPTIMIZED_VERSION_A4dot4
684 FloatMatrix Sders [ 2 ]; // derivatives in each coordinate direction on NURBS
685 #endif
686 FloatMatrix wders; // derivatives in w direction on BSpline
687 /*
688 * IntArray Bin(2,2); // binomial coefficients from 0 to d=1
689 * // Bin(n,k)=(n above k)=n!/k!(n-k)! for n>=k>=0
690 * // lower triangle corresponds to Pascal triangle
691 * // according to A4.4 it seems that only coefficients in lower triangle except the first column are used
692 */
693 // resizing to (2,2) has nothing common with nsd
694 // it is related to the fact that 0th and 1st derivatives are computed
695 for ( int i = 0; i < nsd; i++ ) {
696 Aders [ i ].resize(2, 2);
697 Aders [ i ].zero();
698 #ifndef OPTIMIZED_VERSION_A4dot4
699 Sders [ i ].resize(2, 2);
700 #endif
701 }
702
703 wders.resize(2, 2);
704 wders.zero();
705
706 // calculation of jacobian matrix according to A4.4
707 // calculate values and derivatives of nonrational Bspline surface with weights at first (Aders, wders)
708 int uind = span[0] - degree [ 0 ];
709 int vind = span[1] - degree [ 1 ];
710 int ind = vind * numberOfControlPoints [ 0 ] + uind + 1;
711 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
712 tmp1.zero();
713 tmp2.zero();
714 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
715 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
716 double w = this->weights.at(ind + k);
717
718 tmp1[0] += ders [ 0 ](0, k) * vertexCoords[0] * w; // sum(Nu*x*w)
719 tmp1[1] += ders [ 0 ](0, k) * vertexCoords[1] * w; // sum(Nu*y*w)
720 tmp1[2] += ders [ 0 ](0, k) * w; // sum(Nu*w)
721
722 tmp2[0] += ders [ 0 ](1, k) * vertexCoords[0] * w; // sum(dNu/du*x*w)
723 tmp2[1] += ders [ 0 ](1, k) * vertexCoords[1] * w; // sum(dNu/du*y*w)
724 tmp2[2] += ders [ 0 ](1, k) * w; // sum(dNu/du*w)
725 }
726
727 ind += numberOfControlPoints [ 0 ];
728
729 Aders [ 0 ](0, 0) += ders [ 1 ](0, l) * tmp1[0]; // xw=sum(Nv*sum(Nu*x*w))
730 Aders [ 1 ](0, 0) += ders [ 1 ](0, l) * tmp1[1]; // yw=sum(Nv*sum(Nu*y*w))
731 wders(0, 0) += ders [ 1 ](0, l) * tmp1[2]; // w=sum(Nv*sum(Nu*w))
732
733 Aders [ 0 ](0, 1) += ders [ 1 ](1, l) * tmp1[0]; // dxw/dv=sum(dNv/dv*sum(Nu*x*w))
734 Aders [ 1 ](0, 1) += ders [ 1 ](1, l) * tmp1[1]; // dyw/dv=sum(dNv/dv*sum(Nu*y*w))
735 wders(0, 1) += ders [ 1 ](1, l) * tmp1[2]; // dw/dv=sum(dNv/dv*sum(Nu*w))
736
737 Aders [ 0 ](1, 0) += ders [ 1 ](0, l) * tmp2[0]; // dxw/du=sum(Nv*sum(dNu/du*x*w))
738 Aders [ 1 ](1, 0) += ders [ 1 ](0, l) * tmp2[1]; // dyw/du=sum(Nv*sum(dNu/du*y*w))
739 wders(1, 0) += ders [ 1 ](0, l) * tmp2[2]; // dw/du=sum(Nv*sum(dNu/du*w))
740 }
741
742 double weight = wders(0, 0);
743
744 #ifndef OPTIMIZED_VERSION_A4dot4
745 const int d = 1;
746 // calculate values and derivatives of NURBS surface (A4.4)
747 // since all entries in Pascal triangle up to d=1 are 1, binomial coefficients are ignored
748 for ( int k = 0; k <= d; k++ ) {
749 for ( int l = 0; l <= d - k; l++ ) {
750 tmp1[0] = Aders [ 0 ](k, l);
751 tmp1[1] = Aders [ 1 ](k, l);
752 for ( j = 1; j <= l; j++ ) {
753 tmp1[0] -= wders(0, j) * Sders [ 0 ](k, l - j); // *Bin(l,j)
754 tmp1[1] -= wders(0, j) * Sders [ 1 ](k, l - j); // *Bin(l,j)
755 }
756
757 for ( int i = 1; i <= k; i++ ) {
758 tmp1[0] -= wders(i, 0) * Sders [ 0 ](k - i, l); // *Bin(k,i)
759 tmp1[1] -= wders(i, 0) * Sders [ 1 ](k - i, l); // *Bin(k,i)
760 tmp2.zero();
761 for ( int j = 1; j <= l; j++ ) {
762 tmp2[0] += wders(i, j) * Sders [ 0 ](k - i, l - j); // *Bin(l,j)
763 tmp2[1] += wders(i, j) * Sders [ 1 ](k - i, l - j); // *Bin(l,j)
764 }
765
766 tmp1[0] -= tmp2[0]; // *Bin(k,i)
767 tmp1[1] -= tmp2[1]; // *Bin(k,i)
768 }
769
770 Sders [ 0 ](k, l) = tmp1[0] / weight;
771 Sders [ 1 ](k, l) = tmp1[1] / weight;
772 }
773 }
774
775 jacobian(0, 0) = Sders [ 0 ](1, 0); // dx/du
776 jacobian(0, 1) = Sders [ 1 ](1, 0); // dy/du
777 jacobian(1, 0) = Sders [ 0 ](0, 1); // dx/dv
778 jacobian(1, 1) = Sders [ 1 ](0, 1); // dy/dv
779 #else
780 // optimized version of A4.4 for d=1, binomial coefficients ignored
781
782 /*
783 * // k=0 l=0 loop
784 * Sders[0](0,0) = Aders[0](0,0) / weight;
785 * Sders[1](0,0) = Aders[1](0,0) / weight;
786 * // k=1 l=0 loop
787 * Sders[0](1,0) = (Aders[0](1,0)-wders(1,0)*Sders[0](0,0)) / weight;
788 * Sders[1](1,0) = (Aders[1](1,0)-wders(1,0)*Sders[1](0,0)) / weight;
789 * // k=0 l=1 loop
790 * Sders[0](0,1) = (Aders[0](0,1)-wders(0,1)*Sders[0](0,0)) / weight;
791 * Sders[1](0,1) = (Aders[1](0,1)-wders(0,1)*Sders[1](0,0)) / weight;
792 *
793 * jacobian(0,0) = Sders[0](1,0); // dx/du
794 * jacobian(0,1) = Sders[1](1,0); // dy/du
795 * jacobian(1,0) = Sders[0](0,1); // dx/dv
796 * jacobian(1,1) = Sders[1](0,1); // dy/dv
797 */
798
799 // k=0 l=0 loop
800 tmp1[0] = Aders [ 0 ](0, 0) / weight;
801 tmp1[1] = Aders [ 1 ](0, 0) / weight;
802 // k=1 l=0 loop
803 jacobian(0, 0) = ( Aders [ 0 ](1, 0) - wders(1, 0) * tmp1[0] ) / weight; // dx/du
804 jacobian(0, 1) = ( Aders [ 1 ](1, 0) - wders(1, 0) * tmp1[1] ) / weight; // dy/du
805 // k=0 l=1 loop
806 jacobian(1, 0) = ( Aders [ 0 ](0, 1) - wders(0, 1) * tmp1[0] ) / weight; // dx/dv
807 jacobian(1, 1) = ( Aders [ 1 ](0, 1) - wders(0, 1) * tmp1[1] ) / weight; // dy/dv
808 #endif
809 } else {
810 OOFEM_ERROR("not implemented for nsd = %d", nsd);
811 }
812
813#else
814 std :: vector< FloatArray > Aders(nsd);
815 FloatArray wders; // 0th and 1st derivatives in w direction on BSpline
816
817 for ( int i = 0; i < nsd; i++ ) {
818 Aders [ i ].resize(nsd + 1);
819 Aders [ i ].zero();
820 }
821
822 wders.resize(nsd + 1);
823 wders.zero();
824
825 if ( nsd == 1 ) {
826 // calculate values and derivatives of nonrational Bspline curve with weights at first (Aders, wders)
827 int uind = span[0] - degree [ 0 ];
828 int ind = uind + 1;
829 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
830 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
831 double w = this->weights.at(ind + k);
832
833 Aders [ 0 ][0] += ders [ 0 ](0, k) * vertexCoords[0] * w; // xw=sum(Nu*x*w)
834 wders[0] += ders [ 0 ](0, k) * w; // w=sum(Nu*w)
835
836 Aders [ 0 ][1] += ders [ 0 ](1, k) * vertexCoords[0] * w; // dxw/du=sum(dNu/du*x*w)
837 wders[1] += ders [ 0 ](1, k) * w; // dw/du=sum(dNu/du*w)
838 }
839
840 double weight = wders[0];
841
842 // calculation of jacobian matrix according to Eq 4.7
843 jacobian(0, 0) = ( Aders [ 0 ][1] - wders[1] * Aders [ 0 ][0] / weight ) / weight; // dx/du
844 } else if ( nsd == 2 ) {
845 FloatArray tmp1(nsd + 1), tmp2(nsd + 1); // allow for weight
846
847 // calculate values and derivatives of nonrational Bspline surface with weights at first (Aders, wders)
848 int uind = span[0] - degree [ 0 ];
849 int vind = span[1] - degree [ 1 ];
850 int ind = vind * numberOfControlPoints [ 0 ] + uind + 1;
851 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
852 tmp1.zero();
853 tmp2.zero();
854 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
855 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
856 double w = this->weights.at(ind + k);
857
858 tmp1[0] += ders [ 0 ](0, k) * vertexCoords[0] * w; // sum(Nu*x*w)
859 tmp1[1] += ders [ 0 ](0, k) * vertexCoords[1] * w; // sum(Nu*y*w)
860 tmp1[2] += ders [ 0 ](0, k) * w; // sum(Nu*w)
861
862 tmp2[0] += ders [ 0 ](1, k) * vertexCoords[0] * w; // sum(dNu/du*x*w)
863 tmp2[1] += ders [ 0 ](1, k) * vertexCoords[1] * w; // sum(dNu/du*y*w)
864 tmp2[2] += ders [ 0 ](1, k) * w; // sum(dNu/du*w)
865 }
866
867 ind += numberOfControlPoints [ 0 ];
868
869 Aders [ 0 ][0] += ders [ 1 ](0, l) * tmp1[0]; // xw=sum(Nv*sum(Nu*x*w)
870 Aders [ 1 ][0] += ders [ 1 ](0, l) * tmp1[1]; // yw=sum(Nv*sum(Nu*y*w)
871 wders[0] += ders [ 1 ](0, l) * tmp1[2]; // w=sum(Nv*sum(Nu*w)
872
873 Aders [ 0 ][1] += ders [ 1 ](0, l) * tmp2[0]; // dxw/du=sum(Nv*sum(dNu/du*x*w)
874 Aders [ 1 ][1] += ders [ 1 ](0, l) * tmp2[1]; // dyw/du=sum(Nv*sum(dNu/du*y*w)
875 wders[1] += ders [ 1 ](0, l) * tmp2[2]; // dw/du=sum(Nv*sum(dNu/du*w)
876
877 Aders [ 0 ][2] += ders [ 1 ](1, l) * tmp1[0]; // dxw/dv=sum(dNv/dv*sum(Nu*x*w)
878 Aders [ 1 ][2] += ders [ 1 ](1, l) * tmp1[1]; // dyw/dv=sum(dNv/dv*sum(Nu*y*w)
879 wders[2] += ders [ 1 ](1, l) * tmp1[2]; // dw/dv=sum(dNv/dv*sum(Nu*w)
880 }
881
882 double weight = wders[0];
883
884 // calculation of jacobian matrix according to Eq 4.19
885 tmp1[0] = Aders [ 0 ][0] / weight;
886 tmp1[1] = Aders [ 1 ][0] / weight;
887 jacobian(0, 0) = ( Aders [ 0 ][1] - wders[1] * tmp1[0] ) / weight; // dx/du
888 jacobian(0, 1) = ( Aders [ 1 ][1] - wders[1] * tmp1[1] ) / weight; // dy/du
889 jacobian(1, 0) = ( Aders [ 0 ][2] - wders[2] * tmp1[0] ) / weight; // dx/dv
890 jacobian(1, 1) = ( Aders [ 1 ][2] - wders[2] * tmp1[1] ) / weight; // dy/dv
891 } else if ( nsd == 3 ) {
892 FloatArray tmp1(nsd + 1), tmp2(nsd + 1); // allow for weight
893 FloatArray temp1(nsd + 1), temp2(nsd + 1), temp3(nsd + 1); // allow for weight
894
895 // calculate values and derivatives of nonrational Bspline solid with weights at first (Aders, wders)
896 int uind = span[0] - degree [ 0 ];
897 int vind = span[1] - degree [ 1 ];
898 int tind = span[2] - degree [ 2 ];
899 int ind = tind * numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ] + vind * numberOfControlPoints [ 0 ] + uind + 1;
900 for ( int m = 0; m <= degree [ 2 ]; m++ ) {
901 temp1.zero();
902 temp2.zero();
903 temp3.zero();
904 int indx = ind;
905 for ( int l = 0; l <= degree [ 1 ]; l++ ) {
906 tmp1.zero();
907 tmp2.zero();
908 for ( int k = 0; k <= degree [ 0 ]; k++ ) {
909 const auto &vertexCoords = cellgeo.giveVertexCoordinates(ind + k);
910 double w = this->weights.at(ind + k);
911
912 tmp1[0] += ders [ 0 ](0, k) * vertexCoords[0] * w; // sum(Nu*x*w)
913 tmp1[1] += ders [ 0 ](0, k) * vertexCoords[1] * w; // sum(Nu*y*w)
914 tmp1[2] += ders [ 0 ](0, k) * vertexCoords[2] * w; // sum(Nu*z*w)
915 tmp1[3] += ders [ 0 ](0, k) * w; // sum(Nu*w)
916
917 tmp2[0] += ders [ 0 ](1, k) * vertexCoords[0] * w; // sum(dNu/du*x*w)
918 tmp2[1] += ders [ 0 ](1, k) * vertexCoords[1] * w; // sum(dNu/du*y*w)
919 tmp2[2] += ders [ 0 ](1, k) * vertexCoords[2] * w; // sum(dNu/du*z*w)
920 tmp2[3] += ders [ 0 ](1, k) * w; // sum(dNu/du*w)
921 }
922
923 ind += numberOfControlPoints [ 0 ];
924
925 temp1[0] += ders [ 1 ](0, l) * tmp1[0]; // sum(Nv*sum(Nu*x*w))
926 temp1[1] += ders [ 1 ](0, l) * tmp1[1]; // sum(Nv*sum(Nu*y*w))
927 temp1[2] += ders [ 1 ](0, l) * tmp1[2]; // sum(Nv*sum(Nu*z*w))
928 temp1[3] += ders [ 1 ](0, l) * tmp1[3]; // sum(Nv*sum(Nu*w))
929
930 temp2[0] += ders [ 1 ](0, l) * tmp2[0]; // sum(Nv*sum(dNu/du*x*w))
931 temp2[1] += ders [ 1 ](0, l) * tmp2[1]; // sum(Nv*sum(dNu/du*y*w))
932 temp2[2] += ders [ 1 ](0, l) * tmp2[2]; // sum(Nv*sum(dNu/du*z*w))
933 temp2[3] += ders [ 1 ](0, l) * tmp2[3]; // sum(Nv*sum(dNu/du*w))
934
935 temp3[0] += ders [ 1 ](1, l) * tmp1[0]; // sum(dNv/dv*sum(Nu*x*w))
936 temp3[1] += ders [ 1 ](1, l) * tmp1[1]; // sum(dNv/dv*sum(Nu*y*w))
937 temp3[2] += ders [ 1 ](1, l) * tmp1[2]; // sum(dNv/dv*sum(Nu*z*w))
938 temp3[3] += ders [ 1 ](1, l) * tmp1[3]; // sum(dNv/dv*sum(Nu*w))
939 }
940
941 ind = indx + numberOfControlPoints [ 0 ] * numberOfControlPoints [ 1 ];
942
943 Aders [ 0 ][0] += ders [ 2 ](0, m) * temp1[0]; // x=sum(Nt*sum(Nv*sum(Nu*x*w)))
944 Aders [ 1 ][0] += ders [ 2 ](0, m) * temp1[1]; // y=sum(Nt*sum(Nv*sum(Nu*y*w)))
945 Aders [ 2 ][0] += ders [ 2 ](0, m) * temp1[2]; // y=sum(Nt*sum(Nv*sum(Nu*y*w)))
946 wders[0] += ders [ 2 ](0, m) * temp1[3]; // w=sum(Nt*sum(Nv*sum(Nu*w)))
947
948 Aders [ 0 ][1] += ders [ 2 ](0, m) * temp2[0]; // dx/du=sum(Nt*sum(Nv*sum(dNu/du*x*w)))
949 Aders [ 1 ][1] += ders [ 2 ](0, m) * temp2[1]; // dy/du=sum(Nt*sum(Nv*sum(dNu/du*y*w)))
950 Aders [ 2 ][1] += ders [ 2 ](0, m) * temp2[2]; // dy/du=sum(Nt*sum(Nv*sum(dNu/du*y*w)))
951 wders[1] += ders [ 2 ](0, m) * temp2[3]; // dw/du=sum(Nt*sum(Nv*sum(dNu/du*w)))
952
953 Aders [ 0 ][2] += ders [ 2 ](0, m) * temp3[0]; // dx/dv=sum(Nt*sum(dNv/dv*sum(Nu*x*w)))
954 Aders [ 1 ][2] += ders [ 2 ](0, m) * temp3[1]; // dy/dv=sum(Nt*sum(dNv/dv*sum(Nu*y*w)))
955 Aders [ 2 ][2] += ders [ 2 ](0, m) * temp3[2]; // dy/dv=sum(Nt*sum(dNv/dv*sum(Nu*y*w)))
956 wders[2] += ders [ 2 ](0, m) * temp3[3]; // dw/dv=sum(Nt*sum(dNv/dv*sum(Nu*w)))
957
958 Aders [ 0 ][3] += ders [ 2 ](1, m) * temp1[0]; // dx/dt=sum(dNt/dt*sum(Nv*sum(Nu*x*w)))
959 Aders [ 1 ][3] += ders [ 2 ](1, m) * temp1[1]; // dy/dt=sum(dNt/dt*sum(Nv*sum(Nu*y*w)))
960 Aders [ 2 ][3] += ders [ 2 ](1, m) * temp1[2]; // dy/dt=sum(dNt/dt*sum(Nv*sum(Nu*y*w)))
961 wders[3] += ders [ 2 ](1, m) * temp1[3]; // dw/dt=sum(dNt/dt*sum(Nv*sum(Nu*w)))
962 }
963
964 double weight = wders[0];
965
966 // calculation of jacobian matrix
967 tmp1[0] = Aders [ 0 ][0] / weight;
968 tmp1[1] = Aders [ 1 ][0] / weight;
969 tmp1[2] = Aders [ 2 ][0] / weight;
970 jacobian(0, 0) = ( Aders [ 0 ][1] - wders[1] * tmp1[0] ) / weight; // dx/du
971 jacobian(0, 1) = ( Aders [ 1 ][1] - wders[1] * tmp1[1] ) / weight; // dy/du
972 jacobian(0, 2) = ( Aders [ 2 ][1] - wders[1] * tmp1[2] ) / weight; // dz/du
973 jacobian(1, 0) = ( Aders [ 0 ][2] - wders[2] * tmp1[0] ) / weight; // dx/dv
974 jacobian(1, 1) = ( Aders [ 1 ][2] - wders[2] * tmp1[1] ) / weight; // dy/dv
975 jacobian(1, 2) = ( Aders [ 2 ][2] - wders[2] * tmp1[2] ) / weight; // dz/dv
976 jacobian(2, 0) = ( Aders [ 0 ][3] - wders[3] * tmp1[0] ) / weight; // dx/dt
977 jacobian(2, 1) = ( Aders [ 1 ][3] - wders[3] * tmp1[1] ) / weight; // dy/dt
978 jacobian(2, 2) = ( Aders [ 2 ][3] - wders[3] * tmp1[2] ) / weight; // dz/dt
979 } else {
980 OOFEM_ERROR("not implemented for nsd = %d", nsd);
981 }
982
983#endif
984}
985} // end namespace oofem
N
#define N(a, b)
oofem::BSplineInterpolation::degree
std::array< int, 3 > degree
Degree in each direction.
Definition feibspline.h:68
oofem::BSplineInterpolation::knotVector
std::array< FloatArray, 3 > knotVector
Knot vectors [nsd][knot_vector_size].
Definition feibspline.h:80
oofem::BSplineInterpolation::numberOfControlPoints
std::array< int, 3 > numberOfControlPoints
numberOfControlPoints[nsd]
Definition feibspline.h:78
oofem::BSplineInterpolation::giveNumberOfKnotSpanBasisFunctions
int giveNumberOfKnotSpanBasisFunctions(const IntArray &knotSpan) const override
Definition feibspline.C:625
oofem::BSplineInterpolation::basisFuns
void basisFuns(FloatArray &N, int span, double u, int p, const FloatArray &U) const
Definition feibspline.C:644
oofem::BSplineInterpolation::findSpan
int findSpan(int n, int p, double u, const FloatArray &U) const
Definition feibspline.C:774
oofem::BSplineInterpolation::dersBasisFuns
void dersBasisFuns(int n, double u, int span, int p, const FloatArray &U, FloatMatrix &ders) const
Definition feibspline.C:675
oofem::BSplineInterpolation::nsd
int nsd
Number of spatial directions.
Definition feibspline.h:66
oofem::FEICellGeometry::giveVertexCoordinates
virtual const FloatArray giveVertexCoordinates(int i) const =0
oofem::FEIIGAElementGeometryWrapper::knotSpan
const IntArray * knotSpan
Definition iga.h:60
oofem::FloatArray::resize
void resize(Index s)
Definition floatarray.C:94
oofem::FloatArray::at
double & at(Index i)
Definition floatarray.h:202
oofem::FloatArray::zero
void zero()
Zeroes all coefficients of receiver.
Definition floatarray.C:683
oofem::FloatArray::times
void times(double s)
Definition floatarray.C:834
oofem::FloatMatrix::resize
void resize(Index rows, Index cols)
Definition floatmatrix.C:79
oofem::FloatMatrix::zero
void zero()
Zeroes all coefficient of receiver.
Definition floatmatrix.C:1064
oofem::FloatMatrix::giveDeterminant
double giveDeterminant() const
Definition floatmatrix.C:1085
oofem::NURBSInterpolation::IPK_NURBSInterpolation_weights
static ParamKey IPK_NURBSInterpolation_weights
Definition feinurbs.h:57
OOFEM_ERROR
#define OOFEM_ERROR(...)
Definition error.h:79
oofem::sum
double sum(const FloatArray &x)
Definition floatarray.C:1046
oofem::product
double product(const FloatArray &x)
Definition floatarray.C:1051
PM_ELEMENT_ERROR_IFNOTSET
#define PM_ELEMENT_ERROR_IFNOTSET(_pm, _componentnum, _paramkey)
Definition parametermanager.h:103
PM_UPDATE_PARAMETER
#define PM_UPDATE_PARAMETER(_val, _pm, _ir, _componentnum, _paramkey, _prio)
Definition parametermanager.h:56
This page is part of the OOFEM-3.0 documentation. Copyright Copyright (C) 1994-2025 Borek Patzak Bořek Patzák
Project e-mail: oofem@fsv.cvut.cz
Generated at for OOFEM by doxygen 1.15.0 written by Dimitri van Heesch, © 1997-2011