libMesh
fe_bernstein_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 <>
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 <>
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  libmesh_error_msg("We'll never get here!");
383  return 0.;
384 }
385 
386 
387 
388 template <>
390  const Order,
391  const unsigned int,
392  const unsigned int,
393  const Point &)
394 {
395  libmesh_error_msg("Bernstein polynomials require the element type \nbecause edge orientation is needed.");
396  return 0.;
397 }
398 
399 
400 
401 template <>
403  const Order order,
404  const unsigned int i,
405  const unsigned int j,
406  const Point & p)
407 {
408  libmesh_assert(elem);
409 
410  const ElemType type = elem->type();
411 
412  const Order totalorder = static_cast<Order>(order + elem->p_level());
413 
414  switch (type)
415  {
416  // Hierarchic shape functions on the quadrilateral.
417  case QUAD4:
418  case QUAD9:
419  {
420  // Compute quad shape functions as a tensor-product
421  const Real xi = p(0);
422  const Real eta = p(1);
423 
424  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
425 
426  unsigned int i0, i1;
427 
428  // Vertex DoFs
429  if (i == 0)
430  { i0 = 0; i1 = 0; }
431  else if (i == 1)
432  { i0 = 1; i1 = 0; }
433  else if (i == 2)
434  { i0 = 1; i1 = 1; }
435  else if (i == 3)
436  { i0 = 0; i1 = 1; }
437 
438 
439  // Edge DoFs
440  else if (i < totalorder + 3u)
441  { i0 = i - 2; i1 = 0; }
442  else if (i < 2u*totalorder + 2)
443  { i0 = 1; i1 = i - totalorder - 1; }
444  else if (i < 3u*totalorder + 1)
445  { i0 = i - 2u*totalorder; i1 = 1; }
446  else if (i < 4u*totalorder)
447  { i0 = 0; i1 = i - 3u*totalorder + 1; }
448  // Interior DoFs
449  else
450  {
451  unsigned int basisnum = i - 4*totalorder;
452  i0 = square_number_column[basisnum] + 2;
453  i1 = square_number_row[basisnum] + 2;
454  }
455 
456 
457  // Flip odd degree of freedom values if necessary
458  // to keep continuity on sides
459  if ((i>= 4 && i<= 4+ totalorder-2u) && elem->point(0) > elem->point(1)) i0=totalorder+2-i0;
460  else if ((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem->point(1) > elem->point(2)) i1=totalorder+2-i1;
461  else if ((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->point(3) > elem->point(2)) i0=totalorder+2-i0;
462  else if ((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->point(0) > elem->point(3)) i1=totalorder+2-i1;
463 
464  switch (j)
465  {
466  // d()/dxi
467  case 0:
468  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
469  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta));
470 
471  // d()/deta
472  case 1:
473  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)*
474  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
475 
476  default:
477  libmesh_error_msg("Invalid shape function derivative j = " << j);
478  }
479  }
480 
481  // Bernstein shape functions on the 8-noded quadrilateral
482  // is handled separately.
483  case QUAD8:
484  case QUADSHELL8:
485  {
486  libmesh_assert_less (totalorder, 3);
487 
488  const Real xi = p(0);
489  const Real eta = p(1);
490 
491  libmesh_assert_less (i, 8);
492 
493  // 0 1 2 3 4 5 6 7 8
494  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
495  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
496  static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
497  switch (j)
498  {
499  // d()/dxi
500  case 0:
501  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
502  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)
503  +scal[i]*
504  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)*
505  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[8], eta));
506 
507  // d()/deta
508  case 1:
509  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
510  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)
511  +scal[i]*
512  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[8], xi)*
513  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta));
514 
515  default:
516  libmesh_error_msg("Invalid shape function derivative j = " << j);
517  }
518  }
519 
520  case TRI3:
521  case TRISHELL3:
522  libmesh_assert_less (totalorder, 2);
523  libmesh_fallthrough();
524  case TRI6:
525  {
526  // I have been lazy here and am using finite differences
527  // to compute the derivatives!
528  const Real eps = 1.e-6;
529 
530  switch (j)
531  {
532  // d()/dxi
533  case 0:
534  {
535  const Point pp(p(0)+eps, p(1));
536  const Point pm(p(0)-eps, p(1));
537 
538  return (FE<2,BERNSTEIN>::shape(elem, totalorder, i, pp) -
539  FE<2,BERNSTEIN>::shape(elem, totalorder, i, pm))/2./eps;
540  }
541 
542  // d()/deta
543  case 1:
544  {
545  const Point pp(p(0), p(1)+eps);
546  const Point pm(p(0), p(1)-eps);
547 
548  return (FE<2,BERNSTEIN>::shape(elem, totalorder, i, pp) -
549  FE<2,BERNSTEIN>::shape(elem, totalorder, i, pm))/2./eps;
550  }
551 
552 
553  default:
554  libmesh_error_msg("Invalid shape function derivative j = " << j);
555  }
556  }
557 
558  default:
559  libmesh_error_msg("ERROR: Unsupported element type = " << type);
560  }
561 
562  libmesh_error_msg("We'll never get here!");
563  return 0.;
564 }
565 
566 
567 
568 template <>
570  const Order,
571  const unsigned int,
572  const unsigned int,
573  const Point &)
574 {
575  static bool warning_given = false;
576 
577  if (!warning_given)
578  libMesh::err << "Second derivatives for Bernstein elements "
579  << "are not yet implemented!"
580  << std::endl;
581 
582  warning_given = true;
583  return 0.;
584 }
585 
586 
587 
588 template <>
590  const Order,
591  const unsigned int,
592  const unsigned int,
593  const Point &)
594 {
595  static bool warning_given = false;
596 
597  if (!warning_given)
598  libMesh::err << "Second derivatives for Bernstein elements "
599  << "are not yet implemented!"
600  << std::endl;
601 
602  warning_given = true;
603  return 0.;
604 }
605 
606 } // namespace libMesh
607 
608 
609 #endif// LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
OStreamProxy err
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
const unsigned char square_number_column[]
A specific instantiation of the FEBase class.
Definition: fe.h:89
T pow(const T &x)
Definition: utility.h:195
PetscErrorCode Vec x
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
const unsigned char square_number_row[]
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)