OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
feitspline.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 "feitspline.h"
36 #include "floatarray.h"
37 #include "floatmatrix.h"
38 #include "mathfem.h"
39 #include "iga.h"
40 
41 namespace oofem {
42 // optimized version of A4.4 for d=1
43 #define OPTIMIZED_VERSION_A4dot4
44 
45 
47 {
48  for ( int i = 0; i <= numberOfControlPoints [ 0 ]; i++ ) {
49  for ( int j = 0; j < nsd; j++ ) {
50  delete [] localIndexKnotVector [ i ] [ j ];
51  }
52 
53  delete [] localIndexKnotVector [ i ];
54  }
55 
56  delete [] localIndexKnotVector;
57 
58  delete [] openLocalKnotVector;
59 }
60 
61 
63 {
64  IRResultType result; // Required by IR_GIVE_FIELD macro
65 
67 
68  IntArray localIndexKnotVector_tmp;
69  int *indexKnotVec, indexKnotVal;
70  int pos, p;
71 
72  InputFieldType IFT_localIndexKnotVector [ 3 ] = {
76  };
77  int max_deg = 0;
78  for ( int i = 0; i < nsd; i++ ) {
79  if ( degree [ i ] > max_deg ) {
80  max_deg = degree [ i ];
81  }
82  }
83 
84  openLocalKnotVector = new double [ 3 * max_deg + 2 ];
85 
87  for ( int i = 0; i < totalNumberOfControlPoints; i++ ) {
88  localIndexKnotVector [ i ] = new int * [ nsd ];
89  }
90 
91  for ( int n = 0; n < nsd; n++ ) {
92  localIndexKnotVector_tmp.clear();
93  IR_GIVE_FIELD(ir, localIndexKnotVector_tmp, IFT_localIndexKnotVector [ n ]);
94  if ( localIndexKnotVector_tmp.giveSize() != totalNumberOfControlPoints * ( degree [ n ] + 2 ) ) {
95  OOFEM_WARNING("invalid size of knot vector %s", IFT_localIndexKnotVector [ n ]);
96  return IRRT_BAD_FORMAT;
97  }
98 
99  pos = 0;
100  for ( int i = 0; i < totalNumberOfControlPoints; i++ ) {
101  indexKnotVec = localIndexKnotVector [ i ] [ n ] = new int [ degree [ n ] + 2 ];
102 
103  p = 0;
104  for ( int j = 0; j < degree [ n ] + 2; j++ ) {
105  indexKnotVec [ p++ ] = localIndexKnotVector_tmp(pos++);
106  }
107 
108  // check for monotonicity of local index knot vector with multiplicity
109  indexKnotVal = indexKnotVec [ 0 ];
110  for ( int j = 1; j < degree [ n ] + 2; j++ ) {
111  if ( indexKnotVal > indexKnotVec [ j ] ) {
112  OOFEM_WARNING("local index knot vector %s of control point %d is not monotonic",
113  IFT_localIndexKnotVector [ n ], i + 1);
114  return IRRT_BAD_FORMAT;
115  }
116 
117  /* this is only for the case when TSpline = NURBS
118  * if(indexKnotVal+1 < indexKnotVec[j])
119  * OOFEM_ERROR("local index knot vector %s of control point %d is not continuous",
120  * IFT_localIndexKnotVectorString[n], i+1);
121  */
122  indexKnotVal = indexKnotVec [ j ];
123  }
124 
125  // check for nondegeneracy of local index knot vector
126  if ( indexKnotVal == indexKnotVec [ 0 ] ) {
127  OOFEM_WARNING("local index knot vector %s of control point %d is degenerated",
128  IFT_localIndexKnotVector [ n ], i + 1);
129  return IRRT_BAD_FORMAT;
130  }
131 
132  // check for range of local index knot vector
133  if ( indexKnotVec [ 0 ] <= 0 || indexKnotVal > knotValues [ n ].giveSize() ) {
134  OOFEM_WARNING("local index knot vector %s of control point %d out of range",
135  IFT_localIndexKnotVector [ n ], i + 1);
136  return IRRT_BAD_FORMAT;
137  }
138  }
139  }
140 
141  return IRRT_OK;
142 }
143 
144 
145 void TSplineInterpolation :: evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
146 {
148  FloatArray N(nsd);
149  IntArray span(nsd);
150  IntArray mask;
151  double sum = 0.0, val;
152  int count;
153 
154  if ( nsd != 2 ) {
155  OOFEM_ERROR("implemented for nsd = %d", nsd);
156  }
157 
158  if ( gw->knotSpan ) {
159  span = * gw->knotSpan;
160  } else {
161  for ( int i = 0; i < nsd; i++ ) {
162  span(i) = this->findSpan(numberOfControlPoints [ i ], degree [ i ], lcoords(i), knotVector [ i ]);
163  }
164  }
165 
166  // identify which basis functions are nonzero
167  giveKnotSpanBasisFuncMask(span, mask);
168  count = mask.giveSize();
169  answer.resize(count);
170 
171  if ( nsd == 2 ) {
172  for ( int k = 0; k < count; k++ ) {
173  for ( int i = 0; i < nsd; i++ ) {
174  N(i) = this->basisFunction(lcoords(i), degree [ i ], * giveKnotValues(i + 1), localIndexKnotVector [ mask(k) - 1 ] [ i ]);
175  }
176 
177  answer(k) = val = N(0) * N(1) * cellgeo.giveVertexCoordinates( mask(k) )->at(3); // Nu*Nv*w
178  sum += val;
179  }
180  }
181 
182  while ( count ) {
183  answer.at(count--) /= sum;
184  }
185 }
186 
187 
188 double TSplineInterpolation :: evaldNdx(FloatMatrix &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
189 {
191  const FloatArray *vertexCoordsPtr;
192  FloatMatrix jacobian(nsd, nsd);
193  FloatArray temp(nsd);
194  IntArray span(nsd);
195  IntArray mask;
196  double Jacob = 0., product, w, xw, yw, weight;
197  int count;
198  std :: vector< FloatArray > tmp_ders(nsd);
199  std :: vector< FloatMatrix > ders(nsd);
200 
201  /*
202  * IntArray Bin(2,2); // binomial coefficients from 0 to d=1
203  * // Bin(n,k)=(n above k)=n!/k!(n-k)! for n>=k>=0
204  * // lower triangle corresponds to Pascal triangle
205  * // according to A4.4 it seems that only coefficients in lower triangle except the first column are used
206  */
207  if ( nsd != 2 ) {
208  OOFEM_ERROR("not implemented for nsd = %d", nsd);
209  }
210 
211  if ( gw->knotSpan ) {
212  span = * gw->knotSpan;
213  } else {
214  for ( int i = 0; i < nsd; i++ ) {
215  span(i) = this->findSpan(numberOfControlPoints [ i ], degree [ i ], lcoords(i), knotVector [ i ]);
216  }
217  }
218 
219  // identify which basis functions are nonzero
220  giveKnotSpanBasisFuncMask(span, mask);
221  count = mask.giveSize();
222  answer.resize(count, nsd);
223 
224  for ( int i = 0; i < nsd; i++ ) {
225  ders [ i ].resize(2, count);
226  }
227 
228  if ( nsd == 2 ) {
229  FloatMatrix Aders [ 2 ]; // derivatives in each coordinate direction on BSpline
230  //FloatMatrix Sders[nsd]; // derivatives in each coordinate direction on TSpline
231  FloatMatrix wders; // derivatives in w direction on BSpline
232 
233  // resizing to (2,2) has nothing common with nsd
234  // it is related to the fact that 0th and 1st derivatives are computed
235  for ( int i = 0; i < nsd; i++ ) {
236  Aders [ i ].resize(2, 2);
237  Aders [ i ].zero();
238  //Sders[i].resize(2,2);
239  }
240 
241  wders.resize(2, 2);
242  wders.zero();
243 
244  for ( int k = 0; k < count; k++ ) {
245  for ( int i = 0; i < nsd; i++ ) {
246  // it would be simpler if I could pass k-th column of ders[i] directly to dersBasisFunction HUHU array
247  this->dersBasisFunction(1, lcoords(i), degree [ i ], * giveKnotValues(i + 1), localIndexKnotVector [ mask(k) - 1 ] [ i ], tmp_ders [ i ]);
248  ders [ i ](0, k) = tmp_ders [ i ](0);
249  ders [ i ](1, k) = tmp_ders [ i ](1);
250  }
251 
252  // calculation of jacobian matrix in similar fashion as A4.4
253  // calculate values and derivatives of nonrational Bspline surface with weights at first (Aders, wders)
254  vertexCoordsPtr = cellgeo.giveVertexCoordinates( mask(k) );
255  w = vertexCoordsPtr->at(3);
256  xw = vertexCoordsPtr->at(1) * w;
257  yw = vertexCoordsPtr->at(2) * w;
258 
259  product = tmp_ders [ 0 ](0) * tmp_ders [ 1 ](0); // Nu*Nv
260  Aders [ 0 ](0, 0) += product * xw; // x=sum Nu*Nv*x*w
261  Aders [ 1 ](0, 0) += product * yw; // x=sum Nu*Nv*y*w
262  wders(0, 0) += product * w; // w=sum Nu*Nv*w
263 
264  product = tmp_ders [ 0 ](1) * tmp_ders [ 1 ](0); // dNu/du*Nv
265  Aders [ 0 ](1, 0) += product * xw; // dx/du=sum dNu/du*Nv*x*w
266  Aders [ 1 ](1, 0) += product * yw; // dy/du=sum dNu/du*Nv*y*w
267  wders(1, 0) += product * w; // dw/du=sum dNu/du*Nv*w
268 
269  product = tmp_ders [ 0 ](0) * tmp_ders [ 1 ](1); // Nu*dNv/dv
270  Aders [ 0 ](0, 1) += product * xw; // dx/dv=sum Nu*dNv/dv*x*w
271  Aders [ 1 ](0, 1) += product * yw; // dy/dv=sum Nu*dNv/dv*y*w
272  wders(0, 1) += product * w; // dw/dv=sum Nu*dNv/dv*w
273  }
274 
275  weight = wders(0, 0);
276 
277  // optimized version of A4.4 for d=1, binomial coefficients ignored
278  /*
279  * Sders[0](0,0) = Aders[0](0,0) / weight;
280  * Sders[1](0,0) = Aders[1](0,0) / weight;
281  * Sders[0](0,1) = (Aders[0](0,1)-wders(0,1)*Sders[0](0,0)) / weight;
282  * Sders[1](0,1) = (Aders[1](0,1)-wders(0,1)*Sders[1](0,0)) / weight;
283  * Sders[0](1,0) = (Aders[0](1,0)-wders(1,0)*Sders[0](0,0)) / weight;
284  * Sders[1](1,0) = (Aders[1](1,0)-wders(1,0)*Sders[1](0,0)) / weight;
285  *
286  * jacobian(0,0) = Sders[0](1,0); // dx/du
287  * jacobian(0,1) = Sders[1](1,0); // dy/du
288  * jacobian(1,0) = Sders[0](0,1); // dx/dv
289  * jacobian(1,1) = Sders[1](0,1); // dy/dv
290  */
291 
292  temp(0) = Aders [ 0 ](0, 0) / weight;
293  temp(1) = Aders [ 1 ](0, 0) / weight;
294  jacobian(1, 0) = ( Aders [ 0 ](0, 1) - wders(0, 1) * temp(0) ) / weight; // dx/dv
295  jacobian(1, 1) = ( Aders [ 1 ](0, 1) - wders(0, 1) * temp(1) ) / weight; // dy/dv
296  jacobian(0, 0) = ( Aders [ 0 ](1, 0) - wders(1, 0) * temp(0) ) / weight; // dx/du
297  jacobian(0, 1) = ( Aders [ 1 ](1, 0) - wders(1, 0) * temp(1) ) / weight; // dy/du
298 
299  Jacob = jacobian.giveDeterminant();
300 
301  //calculation of derivatives of TSpline basis functions with respect to local parameters
302  product = Jacob * weight * weight;
303 
304  for ( int k = 0; k < count; k++ ) {
305  w = cellgeo.giveVertexCoordinates( mask(k) )->at(3);
306  // dNu/du*Nv*w*sum(Nv*Nu*w) - Nu*Nv*w*sum(dNu/du*Nv*w)
307  temp(0) = ders [ 0 ](1, k) * ders [ 1 ](0, k) * w * weight - ders [ 0 ](0, k) * ders [ 1 ](0, k) * w * wders(1, 0);
308  // Nu*dNv/dv*w*sum(Nv*Nu*w) - Nu*Nv*w*sum(Nu*dNv/dv*w)
309  temp(1) = ders [ 0 ](0, k) * ders [ 1 ](1, k) * w * weight - ders [ 0 ](0, k) * ders [ 1 ](0, k) * w * wders(0, 1);
310 
311  answer(k, 0) = ( jacobian(1, 1) * temp(0) - jacobian(0, 1) * temp(1) ) / product;
312  answer(k, 1) = ( -jacobian(1, 0) * temp(0) + jacobian(0, 0) * temp(1) ) / product;
313  }
314  }
315 
316  return Jacob;
317 }
318 
319 
320 void TSplineInterpolation :: local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
321 {
322  /* Based on SurfacePoint A4.3 implementation*/
324  const FloatArray *vertexCoordsPtr;
325  FloatArray N(nsd);
326  IntArray span(nsd);
327  IntArray mask;
328  double w, xw, yw, product, weight = 0.0;
329  int count;
330 
331  if ( nsd != 2 ) {
332  OOFEM_ERROR("not implemented for nsd = %d", nsd);
333  }
334 
335  if ( gw->knotSpan ) {
336  span = * gw->knotSpan;
337  } else {
338  for ( int i = 0; i < nsd; i++ ) {
339  span(i) = this->findSpan(numberOfControlPoints [ i ], degree [ i ], lcoords(i), knotVector [ i ]);
340  }
341  }
342 
343  // identify which basis functions are nonzero
344  giveKnotSpanBasisFuncMask(span, mask);
345  count = mask.giveSize();
346 
347  answer.resize(nsd);
348  answer.zero();
349 
350  if ( nsd == 2 ) {
351  for ( int k = 0; k < count; k++ ) {
352  for ( int i = 0; i < nsd; i++ ) {
353  N(i) = this->basisFunction(lcoords(i), degree [ i ], * giveKnotValues(i + 1), localIndexKnotVector [ mask(k) - 1 ] [ i ]);
354  }
355 
356  vertexCoordsPtr = cellgeo.giveVertexCoordinates( mask(k) );
357  w = vertexCoordsPtr->at(3);
358  xw = vertexCoordsPtr->at(1) * w;
359  yw = vertexCoordsPtr->at(2) * w;
360 
361  product = N(0) * N(1); // Nu*Nv
362  answer(0) += product * xw; // x=sum Nu*Nv*x*w
363  answer(1) += product * yw; // y=sum Nu*Nv*y*w
364  weight += product * w; // w=sum Nu*Nv*w
365  }
366  }
367 
368  answer.times(1.0 / weight);
369 }
370 
371 
373 {
374  //
375  // Based on Algorithm A4.4 (p. 137) for d=1
376  //
378  const FloatArray *vertexCoordsPtr;
379  FloatArray temp(nsd);
380  IntArray span(nsd);
381  IntArray mask;
382  double w, xw, yw, product, weight;
383  int count;
384  std :: vector< FloatArray > ders(nsd);
385  jacobian.resize(nsd, nsd);
386 
387  /*
388  * IntArray Bin(2,2); // binomial coefficients from 0 to d=1
389  * // Bin(n,k)=(n above k)=n!/k!(n-k)! for n>=k>=0
390  * // lower triangle corresponds to Pascal triangle
391  * // according to A4.4 it seems that only coefficients in lower triangle except the first column are used
392  */
393  if ( nsd != 2 ) {
394  OOFEM_ERROR("not implemented for nsd = %d", nsd);
395  }
396 
397  if ( gw->knotSpan ) {
398  span = * gw->knotSpan;
399  } else {
400  for ( int i = 0; i < nsd; i++ ) {
401  span(i) = this->findSpan(numberOfControlPoints [ i ], degree [ i ], lcoords(i), knotVector [ i ]);
402  }
403  }
404 
405  // identify which basis functions are nonzero
406  giveKnotSpanBasisFuncMask(span, mask);
407  count = mask.giveSize();
408 
409  if ( nsd == 2 ) {
410  FloatMatrix Aders [ 2 ]; // derivatives in each coordinate direction on BSpline
411  //FloatMatrix Sders[nsd]; // derivatives in each coordinate direction on TSpline
412  FloatMatrix wders; // derivatives in w direction on BSpline
413 
414  // resizing to (2,2) has nothing common with nsd
415  // it is related to the fact that 0th and 1st derivatives are computed
416  for ( int i = 0; i < nsd; i++ ) {
417  Aders [ i ].resize(2, 2);
418  Aders [ i ].zero();
419  //Sders[i].resize(2,2);
420  }
421 
422  wders.resize(2, 2);
423  wders.zero();
424 
425  for ( int k = 0; k < count; k++ ) {
426  for ( int i = 0; i < nsd; i++ ) {
427  this->dersBasisFunction(1, lcoords(i), degree [ i ], * giveKnotValues(i + 1), localIndexKnotVector [ mask(k) - 1 ] [ i ], ders [ i ]);
428  }
429 
430  // calculation of jacobian matrix in similar fashion as A4.4
431  // calculate values and derivatives of nonrational Bspline surface with weights at first (Aders, wders)
432  vertexCoordsPtr = cellgeo.giveVertexCoordinates( mask(k) );
433  w = vertexCoordsPtr->at(3);
434  xw = vertexCoordsPtr->at(1) * w;
435  yw = vertexCoordsPtr->at(2) * w;
436 
437  product = ders [ 0 ](0) * ders [ 1 ](0); // Nu*Nv
438  Aders [ 0 ](0, 0) += product * xw; // x=sum Nu*Nv*x*w
439  Aders [ 1 ](0, 0) += product * yw; // x=sum Nu*Nv*y*w
440  wders(0, 0) += product * w; // w=sum Nu*Nv*w
441 
442  product = ders [ 0 ](1) * ders [ 1 ](0); // dNu/du*Nv
443  Aders [ 0 ](1, 0) += product * xw; // dx/du=sum dNu/du*Nv*x*w
444  Aders [ 1 ](1, 0) += product * yw; // dy/du=sum dNu/du*Nv*y*w
445  wders(1, 0) += product * w; // dw/du=sum dNu/du*Nv*w
446 
447  product = ders [ 0 ](0) * ders [ 1 ](1); // Nu*dNv/dv
448  Aders [ 0 ](0, 1) += product * xw; // dx/dv=sum Nu*dNv/dv*x*w
449  Aders [ 1 ](0, 1) += product * yw; // dy/dv=sum Nu*dNv/dv*y*w
450  wders(0, 1) += product * w; // dw/dv=sum Nu*dNv/dv*w
451  }
452 
453  weight = wders(0, 0);
454 
455  // optimized version of A4.4 for d=1, binomial coefficients ignored
456 #if 0
457  Sders[0](0,0) = Aders[0](0,0) / weight;
458  Sders[1](0,0) = Aders[1](0,0) / weight;
459  Sders[0](0,1) = (Aders[0](0,1)-wders(0,1)*Sders[0](0,0)) / weight;
460  Sders[1](0,1) = (Aders[1](0,1)-wders(0,1)*Sders[1](0,0)) / weight;
461  Sders[0](1,0) = (Aders[0](1,0)-wders(1,0)*Sders[0](0,0)) / weight;
462  Sders[1](1,0) = (Aders[1](1,0)-wders(1,0)*Sders[1](0,0)) / weight;
463 
464  jacobian(0,0) = Sders[0](1,0); // dx/du
465  jacobian(0,1) = Sders[1](1,0); // dy/du
466  jacobian(1,0) = Sders[0](0,1); // dx/dv
467  jacobian(1,1) = Sders[1](0,1); // dy/dv
468 #endif
469 
470  temp(0) = Aders [ 0 ](0, 0) / weight;
471  temp(1) = Aders [ 1 ](0, 0) / weight;
472  jacobian(1, 0) = ( Aders [ 0 ](0, 1) - wders(0, 1) * temp(0) ) / weight; // dx/dv
473  jacobian(1, 1) = ( Aders [ 1 ](0, 1) - wders(0, 1) * temp(1) ) / weight; // dy/dv
474  jacobian(0, 0) = ( Aders [ 0 ](1, 0) - wders(1, 0) * temp(0) ) / weight; // dx/du
475  jacobian(0, 1) = ( Aders [ 1 ](1, 0) - wders(1, 0) * temp(1) ) / weight; // dy/du
476  }
477 }
478 
479 
480 // knotSpan corresponds to knot span in terms of BSpline;
481 // it should not matter which part of IGAIntegrationElement if covering more spans is addressed;
482 // this implementation relies on the fact that IGAIntegrationElements are those subsets
483 // of T-mesh cells on which there are fully (this means not only partially) nonzero relevant basis functions
484 
486 {
487  FloatArray knotStart(nsd), knotEnd(nsd);
488 
489  // resize the mask initially to the size corresponding to BSpline case
490  // but there may be more nonzero basis functions
491  if ( nsd == 2 ) {
492  mask.preallocate( ( degree [ 0 ] + 1 ) * ( degree [ 1 ] + 1 ) );
493  } else {
494  OOFEM_ERROR("not implemented for nsd = %d", nsd);
495  }
496 
497  // get starting and ending knots
498  for ( int j = 0; j < nsd; j++ ) {
499  knotStart(j) = knotVector [ j ] [ knotSpan(j) ];
500  knotEnd(j) = knotVector [ j ] [ knotSpan(j) + 1 ];
501  }
502 
503  // for each control point check
504  for ( int i = 0; i < totalNumberOfControlPoints; i++ ) {
505  // whether local knot vector overlaps the given knot span
506  int nonzero = 1;
507  for ( int j = 0; j < nsd; j++ ) {
508  if ( ( knotEnd(j) <= knotValues [ j ].at(localIndexKnotVector [ i ] [ j ] [ 0 ]) ) ||
509  ( knotStart(j) >= knotValues [ j ].at(localIndexKnotVector [ i ] [ j ] [ degree [ j ] + 1 ]) ) ) {
510  nonzero = 0;
511  break;
512  }
513  }
514 
515  if ( nonzero ) {
516  mask.followedBy(i + 1, 4);
517  }
518  }
519 
520  return 1;
521 }
522 
523 
524 
525 // knotSpan corresponds to knot span in terms of BSpline;
526 // it should not matter which part of IGAIntegrationElement if covering more spans is addressed;
527 // this implementation relies on the fact that IGAIntegrationElements are those subsets
528 // of T-mesh cells on which there are fully (this means not only partially) nonzero relevant basis functions
529 
531 {
532  int answer = 0;
533  FloatArray knotStart(nsd), knotEnd(nsd);
534 
535  // get starting and ending knots
536  for ( int j = 0; j < nsd; j++ ) {
537  knotStart(j) = knotVector [ j ] [ knotSpan(j) ];
538  knotEnd(j) = knotVector [ j ] [ knotSpan(j) + 1 ];
539  }
540 
541  // for each control point check
542  for ( int i = 0; i < totalNumberOfControlPoints; i++ ) {
543  answer++;
544  // whether local knot vector overlaps the given knot span
545  for ( int j = 0; j < nsd; j++ ) {
546  if ( ( knotEnd(j) <= knotValues [ j ].at(localIndexKnotVector [ i ] [ j ] [ 0 ]) ) ||
547  ( knotStart(j) >= knotValues [ j ].at(localIndexKnotVector [ i ] [ j ] [ degree [ j ] + 1 ]) ) ) {
548  answer--;
549  break;
550  }
551  }
552  }
553 
554  return answer;
555 }
556 
557 
558 
559 // starKnotSpan and endKnotSpan correspond to knot span in terms of BSpline;
560 // should the number of non-zero basis function be calculated for single knot span
561 // starKnotSpan and endKnotSpan are equal;
562 // some of the basis function may not cover the whole knot span interval !!!
563 
564 int TSplineInterpolation :: giveKnotSpanBasisFuncMask(const IntArray &startKnotSpan, const IntArray &endKnotSpan, IntArray &mask)
565 {
566  int nonzero;
567  FloatArray knotStart(nsd), knotEnd(nsd);
568 
569  // resize the mask initially to the size corresponding to BSpline case
570  // but there may be more nonzero basis functions
571  if ( nsd == 2 ) {
572  mask.preallocate( ( degree [ 0 ] + 1 ) * ( degree [ 1 ] + 1 ) );
573  } else {
574  OOFEM_ERROR("not implemented for nsd = %d", nsd);
575  }
576 
577  // get starting and ending knots
578  for ( int j = 0; j < nsd; j++ ) {
579  knotStart(j) = knotVector [ j ] [ startKnotSpan(j) ];
580  knotEnd(j) = knotVector [ j ] [ endKnotSpan(j) + 1 ];
581  }
582 
583  // for each control point check
584  for ( int i = 0; i < totalNumberOfControlPoints; i++ ) {
585  // whether local knot vector overlaps at least partially the knot span interval
586  nonzero = 1;
587  for ( int j = 0; j < nsd; j++ ) {
588  if ( ( knotEnd(j) <= knotValues [ j ].at(localIndexKnotVector [ i ] [ j ] [ 0 ]) ) ||
589  ( knotStart(j) >= knotValues [ j ].at(localIndexKnotVector [ i ] [ j ] [ degree [ j ] + 1 ]) ) ) {
590  nonzero = 0;
591  break;
592  }
593  }
594 
595  if ( nonzero ) {
596  mask.followedBy(i + 1, 4);
597  }
598  }
599 
600  return 1;
601 }
602 
603 
604 // starKnotSpan and endKnotSpan correspond to knot apan in terms of BSpline;
605 // should the number of non-zero basis function be calculated for single knot span
606 // starKnotSpan and endKnotSpan are equal
607 // some of the basis function may not cover the whole knot span interval !!!
608 
610 {
611  int answer = 0;
612  FloatArray knotStart(nsd), knotEnd(nsd);
613 
614  // get starting and ending knots
615  for ( int j = 0; j < nsd; j++ ) {
616  knotStart(j) = knotVector [ j ] [ startKnotSpan(j) ];
617  knotEnd(j) = knotVector [ j ] [ endKnotSpan(j) + 1 ];
618  }
619 
620  // for each control point check
621  for ( int i = 0; i < totalNumberOfControlPoints; i++ ) {
622  answer++;
623  // whether local knot vector overlaps at least partially the knot span interval
624  for ( int j = 0; j < nsd; j++ ) {
625  if ( ( knotEnd(j) <= knotValues [ j ].at(localIndexKnotVector [ i ] [ j ] [ 0 ]) ) ||
626  ( knotStart(j) >= knotValues [ j ].at(localIndexKnotVector [ i ] [ j ] [ degree [ j ] + 1 ]) ) ) {
627  answer--;
628  break;
629  }
630  }
631  }
632 
633  return answer;
634 }
635 
636 
637 
638 // call corresponding BSpline methods for open local knot vector
639 
640 double TSplineInterpolation :: basisFunction(double u, int p, const FloatArray &U, const int *I)
641 {
642  int span, prepend, append;
643  FloatArray N;
644 
645  createLocalKnotVector(p, U, I, & prepend, & append);
646  span = BSplineInterpolation :: findSpan(prepend + append, p, u, openLocalKnotVector);
648 
649  // extract the middle basis function
650  // this corresponds to index p-span, however prepended knotspans must be considered
651  return N(p - span + prepend);
652 }
653 
654 
655 // call corresponding BSpline methods for open local knot vector
656 
657 void TSplineInterpolation :: dersBasisFunction(int n, double u, int p, const FloatArray &U, const int *I, FloatArray &ders)
658 {
659  int span, prepend, append;
660  FloatMatrix Ders;
661 
662  createLocalKnotVector(p, U, I, & prepend, & append);
663  span = BSplineInterpolation :: findSpan(prepend + append, p, u, openLocalKnotVector);
665 
666  // extract the middle basis function and its derivatives
667  // this corresponds to index p-span, however prepended knotspans must be considered
668  ders.resize(n + 1);
669  for ( int i = 0; i <= n; i++ ) {
670  ders(i) = Ders(i, p - span + prepend);
671  }
672 }
673 
674 
675 void TSplineInterpolation :: createLocalKnotVector(int p, const FloatArray &U, const int *I, int *prepend, int *append)
676 {
677  int j = 0, index_first = I [ 0 ], index_last = I [ p + 1 ], mult_first = 1, mult_last = 1;
678  double first = U.at(index_first), last = U.at(index_last);
679 
680  for ( int i = 1; i < p + 1; i++ ) {
681  if ( I [ i ] != index_first ) {
682  break;
683  }
684 
685  mult_first++;
686  }
687 
688  for ( int i = p; i > 0; i-- ) {
689  if ( I [ i ] != index_last ) {
690  break;
691  }
692 
693  mult_last++;
694  }
695 
696  * prepend = p + 1 - mult_first;
697  * append = p + 1 - mult_last;
698 
699  // prepend first knot (once more)
700  for ( int i = 0; i <= * prepend; i++ ) {
701  openLocalKnotVector [ j++ ] = first;
702  }
703 
704  // copy middle of knot vector (without first and last)
705  for ( int i = 1; i <= p; i++ ) {
706  openLocalKnotVector [ j++ ] = U.at(I [ i ]);
707  }
708 
709  // append last knot (once more)
710  for ( int i = 0; i <= * append; i++ ) {
711  openLocalKnotVector [ j++ ] = last;
712  }
713 }
714 } // end namespace oofem
double giveDeterminant() const
Returns the trace (sum of diagonal components) of the receiver.
Definition: floatmatrix.C:1408
const IntArray * knotSpan
Definition: iga.h:60
#define _IFT_TSplineInterpolation_localIndexKnotVectorV
Definition: feitspline.h:53
void createLocalKnotVector(int p, const FloatArray &U, const int *I, int *prepend, int *append)
Creates local open knot vector.
Definition: feitspline.C:675
virtual const FloatArray * giveKnotValues(int dim)
Returns the knot values of the receiver.
Definition: feibspline.h:187
void dersBasisFunction(int n, double u, int p, const FloatArray &U, const int *I, FloatArray &ders)
Computes the middle basis function and it derivatives on local knot vector at u.
Definition: feitspline.C:657
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
virtual const FloatArray * giveVertexCoordinates(int i) const =0
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
int * degree
Degree in each direction.
Definition: feibspline.h:66
Class representing a general abstraction for cell geometry.
Definition: feinterpol.h:62
virtual int giveNumberOfKnotSpanBasisFunctions(const IntArray &knotSpan)
Returns the number of nonzero basis functions at individual knot span,.
Definition: feitspline.C:530
virtual void local2global(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates global coordinates from given local ones.
Definition: feitspline.C:320
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: feitspline.C:188
const char * InputFieldType
Identifier of fields in input records.
Definition: inputrecord.h:52
double * openLocalKnotVector
Temporary open local knot vector to enable use of BSpline algorithms (common for all directions) [3*m...
Definition: feitspline.h:70
Class implementing an array of integers.
Definition: intarray.h:61
#define _IFT_TSplineInterpolation_localIndexKnotVectorW
Definition: feitspline.h:54
virtual void evalN(FloatArray &answer, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Evaluates the array of interpolation functions (shape functions) at given point.
Definition: feitspline.C:145
#define OOFEM_ERROR(...)
Definition: error.h:61
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: feibspline.C:59
void clear()
Clears the array (zero size).
Definition: intarray.h:177
double ** knotVector
Knot vectors [nsd][knot_vector_size].
Definition: feibspline.h:78
FloatArray * knotValues
Knot values [nsd].
Definition: feibspline.h:68
#define N(p, q)
Definition: mdm.C:367
virtual void giveJacobianMatrixAt(FloatMatrix &jacobianMatrix, const FloatArray &lcoords, const FEICellGeometry &cellgeo)
Gives the jacobian matrix at the local coordinates.
Definition: feitspline.C:372
Class representing vector of real numbers.
Definition: floatarray.h:82
IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: feitspline.C:62
double basisFunction(double u, int p, const FloatArray &U, const int *I)
Evaluates the middle basis function on local knot vector at u.
Definition: feitspline.C:640
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
int *** localIndexKnotVector
Local index knot vector of the dimensions [totalNumberOfControlPoints][nsd][degree+2].
Definition: feitspline.h:65
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
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
Class representing the general Input Record.
Definition: inputrecord.h:101
void followedBy(const IntArray &b, int allocChunk=0)
Appends array b at the end of receiver.
Definition: intarray.C:145
void zero()
Zeroes all coefficients of receiver.
Definition: floatarray.C:658
void times(double s)
Multiplies receiver with scalar.
Definition: floatarray.C:818
void zero()
Zeroes all coefficient of receiver.
Definition: floatmatrix.C:1326
Geometry wrapper for IGA elements.
Definition: iga.h:57
int * numberOfControlPoints
numberOfControlPoints[nsd]
Definition: feibspline.h:76
int giveSize() const
Definition: intarray.h:203
void preallocate(int futureSize)
Preallocates receiver to given futureSize if larger then allocatedSize.
Definition: intarray.C:130
the oofem namespace is to define a context or scope in which all oofem names are defined.
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
#define _IFT_TSplineInterpolation_localIndexKnotVectorU
Definition: feitspline.h:52
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
int nsd
Number of spatial directions.
Definition: feibspline.h:64
virtual int giveKnotSpanBasisFuncMask(const IntArray &knotSpan, IntArray &mask)
Returns indices (zero based) of nonzero basis functions for given knot span.
Definition: feitspline.C:485
#define OOFEM_WARNING(...)
Definition: error.h:62
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