fe_bernstein_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // Local includes
20 #include "libmesh/libmesh_config.h"
21 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
22 
23 #include "libmesh/fe.h"
24 #include "libmesh/elem.h"
25 #include "libmesh/number_lookups.h"
26 #include "libmesh/utility.h"
27 
28 
29 namespace libMesh
30 {
31 
32 
33 template <>
34 Real FE<2,BERNSTEIN>::shape(const ElemType,
35  const Order,
36  const unsigned int,
37  const Point &)
38 {
39  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
40  return 0.;
41 }
42 
43 
44 
45 template <>
46 Real FE<2,BERNSTEIN>::shape(const Elem * elem,
47  const Order order,
48  const unsigned int i,
49  const Point & p)
50 {
51  libmesh_assert(elem);
52 
53  const ElemType type = elem->type();
54 
55  const Order totalorder = static_cast<Order>(order + elem->p_level());
56 
57  // Declare that we are using our own special power function
58  // from the Utility namespace. This saves typing later.
59  using Utility::pow;
60 
61  switch (type)
62  {
63  // Hierarchic shape functions on the quadrilateral.
64  case QUAD4:
65  case QUADSHELL4:
66  case QUAD9:
67  {
68  // Compute quad shape functions as a tensor-product
69  const Real xi = p(0);
70  const Real eta = p(1);
71 
72  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
73 
74  // Example i, i0, i1 values for totalorder = 5:
75  // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
76  // static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 3, 2, 4, 4, 4, 3, 2, 5, 5, 5, 5, 4, 3, 2};
77  // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 3, 3, 2, 3, 4, 4, 4, 2, 3, 4, 5, 5, 5, 5};
78 
79  unsigned int i0, i1;
80 
81  // Vertex DoFs
82  if (i == 0)
83  { i0 = 0; i1 = 0; }
84  else if (i == 1)
85  { i0 = 1; i1 = 0; }
86  else if (i == 2)
87  { i0 = 1; i1 = 1; }
88  else if (i == 3)
89  { i0 = 0; i1 = 1; }
90 
91 
92  // Edge DoFs
93  else if (i < totalorder + 3u)
94  { i0 = i - 2; i1 = 0; }
95  else if (i < 2u*totalorder + 2)
96  { i0 = 1; i1 = i - totalorder - 1; }
97  else if (i < 3u*totalorder + 1)
98  { i0 = i - 2u*totalorder; i1 = 1; }
99  else if (i < 4u*totalorder)
100  { i0 = 0; i1 = i - 3u*totalorder + 1; }
101  // Interior DoFs. Use Roy's number look up
102  else
103  {
104  unsigned int basisnum = i - 4*totalorder;
105  i0 = square_number_column[basisnum] + 2;
106  i1 = square_number_row[basisnum] + 2;
107  }
108 
109 
110  // Flip odd degree of freedom values if necessary
111  // to keep continuity on sides.
112  if ((i>= 4 && i<= 4+ totalorder-2u) && elem->point(0) > elem->point(1)) i0=totalorder+2-i0;//
113  else if ((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem->point(1) > elem->point(2)) i1=totalorder+2-i1;
114  else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->point(3) > elem->point(2)) i0=totalorder+2-i0;
115  else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->point(0) > elem->point(3)) i1=totalorder+2-i1;
116 
117 
118  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)*
119  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta));
120  }
121  // handle serendipity QUAD8 element separately
122  case QUAD8:
123  case QUADSHELL8:
124  {
125  libmesh_assert_less (totalorder, 3);
126 
127  const Real xi = p(0);
128  const Real eta = p(1);
129 
130  libmesh_assert_less (i, 8);
131 
132  // 0 1 2 3 4 5 6 7 8
133  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
134  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
135  static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
136 
137  //B_t,i0(i)|xi * B_s,i1(i)|eta + scal(i) * B_t,2|xi * B_t,2|eta
138  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
139  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)
140  +scal[i]*
141  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[8], xi)*
142  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[8], eta));
143 
144  }
145 
146  case TRI3:
147  case TRISHELL3:
148  libmesh_assert_less (totalorder, 2);
149  libmesh_fallthrough();
150  case TRI6:
151  switch (totalorder)
152  {
153  case FIRST:
154  {
155  const Real x=p(0);
156  const Real y=p(1);
157  const Real r=1.-x-y;
158 
159  libmesh_assert_less (i, 3);
160 
161  switch(i)
162  {
163  case 0: return r; //f0,0,1
164  case 1: return x; //f0,1,1
165  case 2: return y; //f1,0,1
166 
167  default:
168  libmesh_error_msg("Invalid shape function index i = " << i);
169  }
170  }
171  case SECOND:
172  {
173  const Real x=p(0);
174  const Real y=p(1);
175  const Real r=1.-x-y;
176 
177  libmesh_assert_less (i, 6);
178 
179  switch(i)
180  {
181  case 0: return r*r;
182  case 1: return x*x;
183  case 2: return y*y;
184 
185  case 3: return 2.*x*r;
186  case 4: return 2.*x*y;
187  case 5: return 2.*r*y;
188 
189  default:
190  libmesh_error_msg("Invalid shape function index i = " << i);
191  }
192  }
193  case THIRD:
194  {
195  const Real x=p(0);
196  const Real y=p(1);
197  const Real r=1.-x-y;
198  libmesh_assert_less (i, 10);
199 
200  unsigned int shape=i;
201 
202 
203  if ((i==3||i==4) && elem->point(0) > elem->point(1)) shape=7-i;
204  if ((i==5||i==6) && elem->point(1) > elem->point(2)) shape=11-i;
205  if ((i==7||i==8) && elem->point(0) > elem->point(2)) shape=15-i;
206 
207  switch(shape)
208  {
209  case 0: return r*r*r;
210  case 1: return x*x*x;
211  case 2: return y*y*y;
212 
213  case 3: return 3.*x*r*r;
214  case 4: return 3.*x*x*r;
215 
216  case 5: return 3.*y*x*x;
217  case 6: return 3.*y*y*x;
218 
219  case 7: return 3.*y*r*r;
220  case 8: return 3.*y*y*r;
221 
222  case 9: return 6.*x*y*r;
223 
224  default:
225  libmesh_error_msg("Invalid shape function index shape = " << shape);
226  }
227  }
228  case FOURTH:
229  {
230  const Real x=p(0);
231  const Real y=p(1);
232  const Real r=1-x-y;
233  unsigned int shape=i;
234 
235  libmesh_assert_less (i, 15);
236 
237  if ((i==3||i== 5) && elem->point(0) > elem->point(1)) shape=8-i;
238  if ((i==6||i== 8) && elem->point(1) > elem->point(2)) shape=14-i;
239  if ((i==9||i==11) && elem->point(0) > elem->point(2)) shape=20-i;
240 
241 
242  switch(shape)
243  {
244  // point functions
245  case 0: return r*r*r*r;
246  case 1: return x*x*x*x;
247  case 2: return y*y*y*y;
248 
249  // edge functions
250  case 3: return 4.*x*r*r*r;
251  case 4: return 6.*x*x*r*r;
252  case 5: return 4.*x*x*x*r;
253 
254  case 6: return 4.*y*x*x*x;
255  case 7: return 6.*y*y*x*x;
256  case 8: return 4.*y*y*y*x;
257 
258  case 9: return 4.*y*r*r*r;
259  case 10: return 6.*y*y*r*r;
260  case 11: return 4.*y*y*y*r;
261 
262  // inner functions
263  case 12: return 12.*x*y*r*r;
264  case 13: return 12.*x*x*y*r;
265  case 14: return 12.*x*y*y*r;
266 
267  default:
268  libmesh_error_msg("Invalid shape function index shape = " << shape);
269  }
270  }
271  case FIFTH:
272  {
273  const Real x=p(0);
274  const Real y=p(1);
275  const Real r=1-x-y;
276  unsigned int shape=i;
277 
278  libmesh_assert_less (i, 21);
279 
280  if ((i>= 3&&i<= 6) && elem->point(0) > elem->point(1)) shape=9-i;
281  if ((i>= 7&&i<=10) && elem->point(1) > elem->point(2)) shape=17-i;
282  if ((i>=11&&i<=14) && elem->point(0) > elem->point(2)) shape=25-i;
283 
284  switch(shape)
285  {
286  //point functions
287  case 0: return pow<5>(r);
288  case 1: return pow<5>(x);
289  case 2: return pow<5>(y);
290 
291  //edge functions
292  case 3: return 5.*x *pow<4>(r);
293  case 4: return 10.*pow<2>(x)*pow<3>(r);
294  case 5: return 10.*pow<3>(x)*pow<2>(r);
295  case 6: return 5.*pow<4>(x)*r;
296 
297  case 7: return 5.*y *pow<4>(x);
298  case 8: return 10.*pow<2>(y)*pow<3>(x);
299  case 9: return 10.*pow<3>(y)*pow<2>(x);
300  case 10: return 5.*pow<4>(y)*x;
301 
302  case 11: return 5.*y *pow<4>(r);
303  case 12: return 10.*pow<2>(y)*pow<3>(r);
304  case 13: return 10.*pow<3>(y)*pow<2>(r);
305  case 14: return 5.*pow<4>(y)*r;
306 
307  //inner functions
308  case 15: return 20.*x*y*pow<3>(r);
309  case 16: return 30.*x*pow<2>(y)*pow<2>(r);
310  case 17: return 30.*pow<2>(x)*y*pow<2>(r);
311  case 18: return 20.*x*pow<3>(y)*r;
312  case 19: return 20.*pow<3>(x)*y*r;
313  case 20: return 30.*pow<2>(x)*pow<2>(y)*r;
314 
315  default:
316  libmesh_error_msg("Invalid shape function index shape = " << shape);
317  }
318  }
319  case SIXTH:
320  {
321  const Real x=p(0);
322  const Real y=p(1);
323  const Real r=1-x-y;
324  unsigned int shape=i;
325 
326  libmesh_assert_less (i, 28);
327 
328  if ((i>= 3&&i<= 7) && elem->point(0) > elem->point(1)) shape=10-i;
329  if ((i>= 8&&i<=12) && elem->point(1) > elem->point(2)) shape=20-i;
330  if ((i>=13&&i<=17) && elem->point(0) > elem->point(2)) shape=30-i;
331 
332  switch(shape)
333  {
334  //point functions
335  case 0: return pow<6>(r);
336  case 1: return pow<6>(x);
337  case 2: return pow<6>(y);
338 
339  //edge functions
340  case 3: return 6.*x *pow<5>(r);
341  case 4: return 15.*pow<2>(x)*pow<4>(r);
342  case 5: return 20.*pow<3>(x)*pow<3>(r);
343  case 6: return 15.*pow<4>(x)*pow<2>(r);
344  case 7: return 6.*pow<5>(x)*r;
345 
346  case 8: return 6.*y *pow<5>(x);
347  case 9: return 15.*pow<2>(y)*pow<4>(x);
348  case 10: return 20.*pow<3>(y)*pow<3>(x);
349  case 11: return 15.*pow<4>(y)*pow<2>(x);
350  case 12: return 6.*pow<5>(y)*x;
351 
352  case 13: return 6.*y *pow<5>(r);
353  case 14: return 15.*pow<2>(y)*pow<4>(r);
354  case 15: return 20.*pow<3>(y)*pow<3>(r);
355  case 16: return 15.*pow<4>(y)*pow<2>(r);
356  case 17: return 6.*pow<5>(y)*r;
357 
358  //inner functions
359  case 18: return 30.*x*y*pow<4>(r);
360  case 19: return 60.*x*pow<2>(y)*pow<3>(r);
361  case 20: return 60.* pow<2>(x)*y*pow<3>(r);
362  case 21: return 60.*x*pow<3>(y)*pow<2>(r);
363  case 22: return 60.*pow<3>(x)*y*pow<2>(r);
364  case 23: return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r);
365  case 24: return 30.*x*pow<4>(y)*r;
366  case 25: return 60.*pow<2>(x)*pow<3>(y)*r;
367  case 26: return 60.*pow<3>(x)*pow<2>(y)*r;
368  case 27: return 30.*pow<4>(x)*y*r;
369 
370  default:
371  libmesh_error_msg("Invalid shape function index shape = " << shape);
372  } // switch shape
373  } // case TRI6
374  default:
375  libmesh_error_msg("Invalid totalorder = " << totalorder);
376  } // switch order
377 
378  default:
379  libmesh_error_msg("ERROR: Unsupported element type = " << type);
380  } // switch type
381 }
382 
383 
384 
385 template <>
386 Real FE<2,BERNSTEIN>::shape_deriv(const ElemType,
387  const Order,
388  const unsigned int,
389  const unsigned int,
390  const Point &)
391 {
392  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
393  return 0.;
394 }
395 
396 
397 
398 template <>
399 Real FE<2,BERNSTEIN>::shape_deriv(const Elem * elem,
400  const Order order,
401  const unsigned int i,
402  const unsigned int j,
403  const Point & p)
404 {
405  libmesh_assert(elem);
406 
407  const ElemType type = elem->type();
408 
409  const Order totalorder = static_cast<Order>(order + elem->p_level());
410 
411  switch (type)
412  {
413  // Hierarchic shape functions on the quadrilateral.
414  case QUAD4:
415  case QUAD9:
416  {
417  // Compute quad shape functions as a tensor-product
418  const Real xi = p(0);
419  const Real eta = p(1);
420 
421  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
422 
423  unsigned int i0, i1;
424 
425  // Vertex DoFs
426  if (i == 0)
427  { i0 = 0; i1 = 0; }
428  else if (i == 1)
429  { i0 = 1; i1 = 0; }
430  else if (i == 2)
431  { i0 = 1; i1 = 1; }
432  else if (i == 3)
433  { i0 = 0; i1 = 1; }
434 
435 
436  // Edge DoFs
437  else if (i < totalorder + 3u)
438  { i0 = i - 2; i1 = 0; }
439  else if (i < 2u*totalorder + 2)
440  { i0 = 1; i1 = i - totalorder - 1; }
441  else if (i < 3u*totalorder + 1)
442  { i0 = i - 2u*totalorder; i1 = 1; }
443  else if (i < 4u*totalorder)
444  { i0 = 0; i1 = i - 3u*totalorder + 1; }
445  // Interior DoFs
446  else
447  {
448  unsigned int basisnum = i - 4*totalorder;
449  i0 = square_number_column[basisnum] + 2;
450  i1 = square_number_row[basisnum] + 2;
451  }
452 
453 
454  // Flip odd degree of freedom values if necessary
455  // to keep continuity on sides
456  if ((i>= 4 && i<= 4+ totalorder-2u) && elem->point(0) > elem->point(1)) i0=totalorder+2-i0;
457  else if ((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem->point(1) > elem->point(2)) i1=totalorder+2-i1;
458  else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->point(3) > elem->point(2)) i0=totalorder+2-i0;
459  else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->point(0) > elem->point(3)) i1=totalorder+2-i1;
460 
461  switch (j)
462  {
463  // d()/dxi
464  case 0:
465  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
466  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta));
467 
468  // d()/deta
469  case 1:
470  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)*
471  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
472 
473  default:
474  libmesh_error_msg("Invalid shape function derivative j = " << j);
475  }
476  }
477 
478  // Bernstein shape functions on the 8-noded quadrilateral
479  // is handled separately.
480  case QUAD8:
481  case QUADSHELL8:
482  {
483  libmesh_assert_less (totalorder, 3);
484 
485  const Real xi = p(0);
486  const Real eta = p(1);
487 
488  libmesh_assert_less (i, 8);
489 
490  // 0 1 2 3 4 5 6 7 8
491  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
492  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
493  static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
494  switch (j)
495  {
496  // d()/dxi
497  case 0:
498  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
499  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)
500  +scal[i]*
501  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)*
502  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[8], eta));
503 
504  // d()/deta
505  case 1:
506  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
507  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)
508  +scal[i]*
509  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[8], xi)*
510  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta));
511 
512  default:
513  libmesh_error_msg("Invalid shape function derivative j = " << j);
514  }
515  }
516 
517  case TRI3:
518  case TRISHELL3:
519  libmesh_assert_less (totalorder, 2);
520  libmesh_fallthrough();
521  case TRI6:
522  {
523  // I have been lazy here and am using finite differences
524  // to compute the derivatives!
525  const Real eps = 1.e-6;
526 
527  switch (j)
528  {
529  // d()/dxi
530  case 0:
531  {
532  const Point pp(p(0)+eps, p(1));
533  const Point pm(p(0)-eps, p(1));
534 
535  return (FE<2,BERNSTEIN>::shape(elem, totalorder, i, pp) -
536  FE<2,BERNSTEIN>::shape(elem, totalorder, i, pm))/2./eps;
537  }
538 
539  // d()/deta
540  case 1:
541  {
542  const Point pp(p(0), p(1)+eps);
543  const Point pm(p(0), p(1)-eps);
544 
545  return (FE<2,BERNSTEIN>::shape(elem, totalorder, i, pp) -
546  FE<2,BERNSTEIN>::shape(elem, totalorder, i, pm))/2./eps;
547  }
548 
549 
550  default:
551  libmesh_error_msg("Invalid shape function derivative j = " << j);
552  }
553  }
554 
555  default:
556  libmesh_error_msg("ERROR: Unsupported element type = " << type);
557  }
558 }
559 
560 
561 
562 template <>
563 Real FE<2,BERNSTEIN>::shape_second_deriv(const ElemType,
564  const Order,
565  const unsigned int,
566  const unsigned int,
567  const Point &)
568 {
569  static bool warning_given = false;
570 
571  if (!warning_given)
572  libMesh::err << "Second derivatives for Bernstein elements "
573  << "are not yet implemented!"
574  << std::endl;
575 
576  warning_given = true;
577  return 0.;
578 }
579 
580 
581 
582 template <>
583 Real FE<2,BERNSTEIN>::shape_second_deriv(const Elem *,
584  const Order,
585  const unsigned int,
586  const unsigned int,
587  const Point &)
588 {
589  static bool warning_given = false;
590 
591  if (!warning_given)
592  libMesh::err << "Second derivatives for Bernstein elements "
593  << "are not yet implemented!"
594  << std::endl;
595 
596  warning_given = true;
597  return 0.;
598 }
599 
600 } // namespace libMesh
601 
602 
603 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2555
libMesh::square_number_column
const unsigned char square_number_column[]
Definition: number_lookups.C:56
libMesh::FE
Template class which generates the different FE families and orders.
Definition: fe.h:89
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:198
libMesh::err
OStreamProxy err(std::cerr)
Definition: libmesh_common.h:229
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::square_number_row
const unsigned char square_number_row[]
Definition: number_lookups.C:69
libmesh_config.h
number_lookups.h
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1892
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
generated by doxygen