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