fe_bernstein_shape_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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::err << "Bernstein polynomials require the element type\n"
40  << "because edge orientation is needed."
41  << std::endl;
42 
43  libmesh_error();
44  return 0.;
45 }
46 
47 
48 
49 template <>
51  const Order order,
52  const unsigned int i,
53  const Point& p)
54 {
55  libmesh_assert(elem);
56 
57  const ElemType type = elem->type();
58 
59  const Order totalorder = static_cast<Order>(order + elem->p_level());
60 
61  // Declare that we are using our own special power function
62  // from the Utility namespace. This saves typing later.
63  using Utility::pow;
64 
65  switch (type)
66  {
67  // Hierarchic shape functions on the quadrilateral.
68  case QUAD4:
69  case QUAD9:
70  {
71  // Compute quad shape functions as a tensor-product
72  const Real xi = p(0);
73  const Real eta = p(1);
74 
75  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
76 
77 // Example i, i0, i1 values for totalorder = 5:
78 // 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
79 // 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};
80 // 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};
81 
82  unsigned int i0, i1;
83 
84  // Vertex DoFs
85  if (i == 0)
86  { i0 = 0; i1 = 0; }
87  else if (i == 1)
88  { i0 = 1; i1 = 0; }
89  else if (i == 2)
90  { i0 = 1; i1 = 1; }
91  else if (i == 3)
92  { i0 = 0; i1 = 1; }
93 
94 
95  // Edge DoFs
96  else if (i < totalorder + 3u)
97  { i0 = i - 2; i1 = 0; }
98  else if (i < 2u*totalorder + 2)
99  { i0 = 1; i1 = i - totalorder - 1; }
100  else if (i < 3u*totalorder + 1)
101  { i0 = i - 2u*totalorder; i1 = 1; }
102  else if (i < 4u*totalorder)
103  { i0 = 0; i1 = i - 3u*totalorder + 1; }
104  // Interior DoFs. Use Roy's number look up
105  else
106  {
107  unsigned int basisnum = i - 4*totalorder;
108  i0 = square_number_column[basisnum] + 2;
109  i1 = square_number_row[basisnum] + 2;
110  }
111 
112 
113  // Flip odd degree of freedom values if necessary
114  // to keep continuity on sides.
115  if ((i>= 4 && i<= 4+ totalorder-2u) && elem->point(0) > elem->point(1)) i0=totalorder+2-i0; //
116  else if((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem->point(1) > elem->point(2)) i1=totalorder+2-i1;
117  else if((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->point(3) > elem->point(2)) i0=totalorder+2-i0;
118  else if((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->point(0) > elem->point(3)) i1=totalorder+2-i1;
119 
120 
121  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)*
122  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta));
123  }
124  // handle serendipity QUAD8 element separately
125  case QUAD8:
126  {
127  libmesh_assert_less (totalorder, 3);
128 
129  const Real xi = p(0);
130  const Real eta = p(1);
131 
132  libmesh_assert_less (i, 8);
133 
134  // 0 1 2 3 4 5 6 7 8
135  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
136  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
137  static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
138 
139  //B_t,i0(i)|xi * B_s,i1(i)|eta + scal(i) * B_t,2|xi * B_t,2|eta
140  return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
141  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)
142  +scal[i]*
143  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[8], xi)*
144  FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[8], eta));
145 
146  }
147 
148  case TRI3:
149  libmesh_assert_less (totalorder, 2);
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: libmesh_error(); return 0;
168  }
169  }
170  case SECOND:
171  {
172  const Real x=p(0);
173  const Real y=p(1);
174  const Real r=1.-x-y;
175 
176  libmesh_assert_less (i, 6);
177 
178  switch(i)
179  {
180  case 0: return r*r;
181  case 1: return x*x;
182  case 2: return y*y;
183 
184  case 3: return 2.*x*r;
185  case 4: return 2.*x*y;
186  case 5: return 2.*r*y;
187 
188  default: libmesh_error(); return 0;
189  }
190  }
191  case THIRD:
192  {
193  const Real x=p(0);
194  const Real y=p(1);
195  const Real r=1.-x-y;
196  libmesh_assert_less (i, 10);
197 
198  unsigned int shape=i;
199 
200 
201  if((i==3||i==4) && elem->point(0) > elem->point(1)) shape=7-i;
202  if((i==5||i==6) && elem->point(1) > elem->point(2)) shape=11-i;
203  if((i==7||i==8) && elem->point(0) > elem->point(2)) shape=15-i;
204 
205  switch(shape)
206  {
207  case 0: return r*r*r;
208  case 1: return x*x*x;
209  case 2: return y*y*y;
210 
211  case 3: return 3.*x*r*r;
212  case 4: return 3.*x*x*r;
213 
214  case 5: return 3.*y*x*x;
215  case 6: return 3.*y*y*x;
216 
217  case 7: return 3.*y*r*r;
218  case 8: return 3.*y*y*r;
219 
220  case 9: return 6.*x*y*r;
221 
222  default: libmesh_error(); return 0;
223  }
224  }
225  case FOURTH:
226  {
227  const Real x=p(0);
228  const Real y=p(1);
229  const Real r=1-x-y;
230  unsigned int shape=i;
231 
232  libmesh_assert_less (i, 15);
233 
234  if((i==3||i== 5) && elem->point(0) > elem->point(1))shape=8-i;
235  if((i==6||i== 8) && elem->point(1) > elem->point(2))shape=14-i;
236  if((i==9||i==11) && elem->point(0) > elem->point(2))shape=20-i;
237 
238 
239  switch(shape)
240  {
241  // point functions
242  case 0: return r*r*r*r;
243  case 1: return x*x*x*x;
244  case 2: return y*y*y*y;
245 
246  // edge functions
247  case 3: return 4.*x*r*r*r;
248  case 4: return 6.*x*x*r*r;
249  case 5: return 4.*x*x*x*r;
250 
251  case 6: return 4.*y*x*x*x;
252  case 7: return 6.*y*y*x*x;
253  case 8: return 4.*y*y*y*x;
254 
255  case 9: return 4.*y*r*r*r;
256  case 10: return 6.*y*y*r*r;
257  case 11: return 4.*y*y*y*r;
258 
259  // inner functions
260  case 12: return 12.*x*y*r*r;
261  case 13: return 12.*x*x*y*r;
262  case 14: return 12.*x*y*y*r;
263 
264  default: libmesh_error(); return 0;
265  }
266  }
267  case FIFTH:
268  {
269  const Real x=p(0);
270  const Real y=p(1);
271  const Real r=1-x-y;
272  unsigned int shape=i;
273 
274  libmesh_assert_less (i, 21);
275 
276  if((i>= 3&&i<= 6) && elem->point(0) > elem->point(1))shape=9-i;
277  if((i>= 7&&i<=10) && elem->point(1) > elem->point(2))shape=17-i;
278  if((i>=11&&i<=14) && elem->point(0) > elem->point(2))shape=25-i;
279 
280  switch(shape)
281  {
282  //point functions
283  case 0: return pow<5>(r);
284  case 1: return pow<5>(x);
285  case 2: return pow<5>(y);
286 
287  //edge functions
288  case 3: return 5.*x *pow<4>(r);
289  case 4: return 10.*pow<2>(x)*pow<3>(r);
290  case 5: return 10.*pow<3>(x)*pow<2>(r);
291  case 6: return 5.*pow<4>(x)*r;
292 
293  case 7: return 5.*y *pow<4>(x);
294  case 8: return 10.*pow<2>(y)*pow<3>(x);
295  case 9: return 10.*pow<3>(y)*pow<2>(x);
296  case 10: return 5.*pow<4>(y)*x;
297 
298  case 11: return 5.*y *pow<4>(r);
299  case 12: return 10.*pow<2>(y)*pow<3>(r);
300  case 13: return 10.*pow<3>(y)*pow<2>(r);
301  case 14: return 5.*pow<4>(y)*r;
302 
303  //inner functions
304  case 15: return 20.*x*y*pow<3>(r);
305  case 16: return 30.*x*pow<2>(y)*pow<2>(r);
306  case 17: return 30.*pow<2>(x)*y*pow<2>(r);
307  case 18: return 20.*x*pow<3>(y)*r;
308  case 19: return 20.*pow<3>(x)*y*r;
309  case 20: return 30.*pow<2>(x)*pow<2>(y)*r;
310 
311  default: libmesh_error(); return 0;
312  }
313  }
314  case SIXTH:
315  {
316  const Real x=p(0);
317  const Real y=p(1);
318  const Real r=1-x-y;
319  unsigned int shape=i;
320 
321  libmesh_assert_less (i, 28);
322 
323  if((i>= 3&&i<= 7) && elem->point(0) > elem->point(1))shape=10-i;
324  if((i>= 8&&i<=12) && elem->point(1) > elem->point(2))shape=20-i;
325  if((i>=13&&i<=17) && elem->point(0) > elem->point(2))shape=30-i;
326 
327  switch(shape)
328  {
329  //point functions
330  case 0: return pow<6>(r);
331  case 1: return pow<6>(x);
332  case 2: return pow<6>(y);
333 
334  //edge functions
335  case 3: return 6.*x *pow<5>(r);
336  case 4: return 15.*pow<2>(x)*pow<4>(r);
337  case 5: return 20.*pow<3>(x)*pow<3>(r);
338  case 6: return 15.*pow<4>(x)*pow<2>(r);
339  case 7: return 6.*pow<5>(x)*r;
340 
341  case 8: return 6.*y *pow<5>(x);
342  case 9: return 15.*pow<2>(y)*pow<4>(x);
343  case 10: return 20.*pow<3>(y)*pow<3>(x);
344  case 11: return 15.*pow<4>(y)*pow<2>(x);
345  case 12: return 6.*pow<5>(y)*x;
346 
347  case 13: return 6.*y *pow<5>(r);
348  case 14: return 15.*pow<2>(y)*pow<4>(r);
349  case 15: return 20.*pow<3>(y)*pow<3>(r);
350  case 16: return 15.*pow<4>(y)*pow<2>(r);
351  case 17: return 6.*pow<5>(y)*r;
352 
353  //inner functions
354  case 18: return 30.*x*y*pow<4>(r);
355  case 19: return 60.*x*pow<2>(y)*pow<3>(r);
356  case 20: return 60.* pow<2>(x)*y*pow<3>(r);
357  case 21: return 60.*x*pow<3>(y)*pow<2>(r);
358  case 22: return 60.*pow<3>(x)*y*pow<2>(r);
359  case 23: return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r);
360  case 24: return 30.*x*pow<4>(y)*r;
361  case 25: return 60.*pow<2>(x)*pow<3>(y)*r;
362  case 26: return 60.*pow<3>(x)*pow<2>(y)*r;
363  case 27: return 30.*pow<4>(x)*y*r;
364 
365  default: libmesh_error(); return 0;
366  } // switch shape
367  } // case TRI6
368  default:
369  {
370  libMesh::err << "ERROR: element order!" << std::endl;
371  libmesh_error();
372  }
373 
374 
375  } // switch order
376 
377 
378  default:
379  {
380  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
381  libmesh_error();
382  }
383 
384  } // switch type
385 
386 
387  // old code
388 // switch (totalorder)
389 // {
390 
391 // case FIRST:
392 
393 // switch(type)
394 // {
395 
396 // case TRI6:
397 // {
398 // const Real x=p(0);
399 // const Real y=p(1);
400 // const Real r=1.-x-y;
401 
402 // libmesh_assert_less (i, 3);
403 
404 // switch(i)
405 // {
406 // case 0: return r; //f0,0,1
407 // case 1: return x; //f0,1,1
408 // case 2: return y; //f1,0,1
409 // }
410 // }
411 
412 // case QUAD8:
413 // case QUAD9:
414 // {
415 // // Compute quad shape functions as a tensor-product
416 // const Real xi = p(0);
417 // const Real eta = p(1);
418 
419 // libmesh_assert_less (i, 4);
420 
421 // // 0 1 2 3
422 // static const unsigned int i0[] = {0, 1, 1, 0};
423 // static const unsigned int i1[] = {0, 0, 1, 1};
424 
425 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
426 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta));
427 
428 // }
429 // default:
430 // libmesh_error();
431 // }
432 
433 // case SECOND:
434 // switch(type)
435 // {
436 // case TRI6:
437 // {
438 // const Real x=p(0);
439 // const Real y=p(1);
440 // const Real r=1.-x-y;
441 
442 // libmesh_assert_less (i, 6);
443 
444 // switch(i)
445 // {
446 // case 0: return r*r;
447 // case 1: return x*x;
448 // case 2: return y*y;
449 
450 // case 3: return 2.*x*r;
451 // case 4: return 2.*x*y;
452 // case 5: return 2.*r*y;
453 // }
454 // }
455 
456 // // Bernstein shape functions on the 8-noded quadrilateral.
457 // case QUAD8:
458 // {
459 // const Real xi = p(0);
460 // const Real eta = p(1);
461 
462 // libmesh_assert_less (i, 8);
463 
464 // // 0 1 2 3 4 5 6 7 8
465 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
466 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
467 // static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
468 
469 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
470 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)
471 // +scal[i]*
472 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[8], xi)*
473 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[8], eta));
474 
475 // //B_t,i0(i)|xi * B_s,i1(i)|eta + scal(i) * B_t,2|xi * B_t,2|eta
476 // }
477 
478 // // Bernstein shape functions on the 9-noded quadrilateral.
479 // case QUAD9:
480 // {
481 // // Compute quad shape functions as a tensor-product
482 // const Real xi = p(0);
483 // const Real eta = p(1);
484 
485 // libmesh_assert_less (i, 9);
486 
487 // // 0 1 2 3 4 5 6 7 8
488 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
489 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
490 
491 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
492 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta));
493 
494 // }
495 // default:
496 // libmesh_error();
497 // }
498 
499 // case THIRD:
500 // switch(type)
501 // {
502 // case TRI6:
503 // {
504 // const Real x=p(0);
505 // const Real y=p(1);
506 // const Real r=1.-x-y;
507 // libmesh_assert_less (i, 10);
508 
509 // unsigned int shape=i;
510 
511 
512 // if((i==3||i==4) && elem->node(0) > elem->node(1))shape=7-i;
513 // if((i==5||i==6) && elem->node(1) > elem->node(2))shape=11-i;
514 // if((i==7||i==8) && elem->node(0) > elem->node(2))shape=15-i;
515 
516 
517 // switch(shape)
518 // {
519 // case 0: return r*r*r;
520 // case 1: return x*x*x;
521 // case 2: return y*y*y;
522 
523 // case 3: return 3.*x*r*r;
524 // case 4: return 3.*x*x*r;
525 
526 // case 5: return 3.*y*x*x;
527 // case 6: return 3.*y*y*x;
528 
529 // case 7: return 3.*y*r*r;
530 // case 8: return 3.*y*y*r;
531 
532 // case 9: return 6.*x*y*r;
533 // }
534 // }
535 
536 // // Bernstein shape functions on the quadrilateral.
537 // case QUAD8:
538 // case QUAD9:
539 // {
540 // // Compute quad shape functions as a tensor-product
541 // Real xi = p(0);
542 // Real eta = p(1);
543 
544 // libmesh_assert_less (i, 16);
545 
546 // // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
547 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
548 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
549 
550 // unsigned int i0=i0_reg[i];
551 // unsigned int i1=i1_reg[i];
552 
553 // if((i== 4||i== 5) && elem->node(0) > elem->node(1)) i0=5-i0; // 2->3 3->2
554 // if((i== 6||i== 7) && elem->node(1) > elem->node(2)) i1=5-i1;
555 // if((i== 8||i== 9) && elem->node(3) > elem->node(2)) i0=5-i0;
556 // if((i==10||i==11) && elem->node(0) > elem->node(3)) i1=5-i1;
557 
558 // // element dof orientation is needed when used with ifems
559 // // if(i > 11)
560 // // {
561 // // const unsigned int min_node = std::min(elem->node(1),
562 // // std::min(elem->node(2),
563 // // std::min(elem->node(0),
564 // // elem->node(3))));
565 // // if (elem->node(0) == min_node)
566 // // if (elem->node(1) == std::min(elem->node(1), elem->node(3)))
567 // // {
568 // // // Case 1
569 // // xi = xi;
570 // // eta = eta;
571 // // }
572 // // else
573 // // {
574 // // // Case 2
575 // // xi = eta;
576 // // eta = xi;
577 // // }
578 
579 // // else if (elem->node(3) == min_node)
580 // // if (elem->node(0) == std::min(elem->node(0), elem->node(2)))
581 // // {
582 // // // Case 3
583 // // xi = -eta;
584 // // eta = xi;
585 // // }
586 // // else
587 // // {
588 // // // Case 4
589 // // xi = xi;
590 // // eta = -eta;
591 // // }
592 
593 // // else if (elem->node(2) == min_node)
594 // // if (elem->node(3) == std::min(elem->node(3), elem->node(1)))
595 // // {
596 // // // Case 5
597 // // xi = -xi;
598 // // eta = -eta;
599 // // }
600 // // else
601 // // {
602 // // // Case 6
603 // // xi = -eta;
604 // // eta = -xi;
605 // // }
606 
607 // // else if (elem->node(1) == min_node)
608 // // if (elem->node(2) == std::min(elem->node(2), elem->node(0)))
609 // // {
610 // // // Case 7
611 // // xi = eta;
612 // // eta = -xi;
613 // // }
614 // // else
615 // // {
616 // // // Case 8
617 // // xi = -xi;
618 // // eta = eta;
619 // // }
620 // // }
621 
622 
623 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)*
624 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta));
625 
626 // }
627 
628 // default:
629 // libmesh_error();
630 
631 // }
632 
633 // case FOURTH:
634 // switch(type)
635 // {
636 // case TRI6:
637 // {
638 // const Real x=p(0);
639 // const Real y=p(1);
640 // const Real r=1-x-y;
641 // unsigned int shape=i;
642 
643 // libmesh_assert_less (i, 15);
644 
645 // if((i==3||i== 5) && elem->node(0) > elem->node(1))shape=8-i;
646 // if((i==6||i== 8) && elem->node(1) > elem->node(2))shape=14-i;
647 // if((i==9||i==11) && elem->node(0) > elem->node(2))shape=20-i;
648 
649 
650 // switch(shape)
651 // {
652 // // point functions
653 // case 0: return r*r*r*r;
654 // case 1: return x*x*x*x;
655 // case 2: return y*y*y*y;
656 
657 // // edge functions
658 // case 3: return 4.*x*r*r*r;
659 // case 4: return 6.*x*x*r*r;
660 // case 5: return 4.*x*x*x*r;
661 
662 // case 6: return 4.*y*x*x*x;
663 // case 7: return 6.*y*y*x*x;
664 // case 8: return 4.*y*y*y*x;
665 
666 // case 9: return 4.*y*r*r*r;
667 // case 10: return 6.*y*y*r*r;
668 // case 11: return 4.*y*y*y*r;
669 
670 // // inner functions
671 // case 12: return 12.*x*y*r*r;
672 // case 13: return 12.*x*x*y*r;
673 // case 14: return 12.*x*y*y*r;
674 
675 // }
676 // }
677 
678 
679 // // Bernstein shape functions on the quadrilateral.
680 // case QUAD8:
681 // case QUAD9:
682 // {
683 // // Compute quad shape functions as a tensor-product
684 // const Real xi = p(0);
685 // const Real eta = p(1);
686 
687 // libmesh_assert_less (i, 25);
688 
689 // // 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
690 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
691 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
692 
693 // unsigned int i0=i0_reg[i];
694 // unsigned int i1=i1_reg[i];
695 
696 // if((i>= 4&&i<= 6) && elem->node(0) > elem->node(1)) i0=6-i0; // 2->4, 4->2
697 // if((i>= 7&&i<= 9) && elem->node(1) > elem->node(2)) i1=6-i1;
698 // if((i>=10&&i<=12) && elem->node(3) > elem->node(2)) i0=6-i0;
699 // if((i>=13&&i<=15) && elem->node(0) > elem->node(3)) i1=6-i1;
700 
701 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)*
702 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta));
703 
704 // }
705 
706 // default:
707 // libmesh_error();
708 
709 // }
710 
711 // case FIFTH:
712 // switch(type)
713 // {
714 // case TRI6:
715 // {
716 // const Real x=p(0);
717 // const Real y=p(1);
718 // const Real r=1-x-y;
719 // unsigned int shape=i;
720 
721 // libmesh_assert_less (i, 21);
722 
723 // if((i>= 3&&i<= 6) && elem->node(0) > elem->node(1))shape=9-i;
724 // if((i>= 7&&i<=10) && elem->node(1) > elem->node(2))shape=17-i;
725 // if((i>=11&&i<=14) && elem->node(0) > elem->node(2))shape=25-i;
726 
727 // switch(shape)
728 // {
729 // //point functions
730 // case 0: return pow<5>(r);
731 // case 1: return pow<5>(x);
732 // case 2: return pow<5>(y);
733 
734 // //edge functions
735 // case 3: return 5.*x *pow<4>(r);
736 // case 4: return 10.*pow<2>(x)*pow<3>(r);
737 // case 5: return 10.*pow<3>(x)*pow<2>(r);
738 // case 6: return 5.*pow<4>(x)*r;
739 
740 // case 7: return 5.*y *pow<4>(x);
741 // case 8: return 10.*pow<2>(y)*pow<3>(x);
742 // case 9: return 10.*pow<3>(y)*pow<2>(x);
743 // case 10: return 5.*pow<4>(y)*x;
744 
745 // case 11: return 5.*y *pow<4>(r);
746 // case 12: return 10.*pow<2>(y)*pow<3>(r);
747 // case 13: return 10.*pow<3>(y)*pow<2>(r);
748 // case 14: return 5.*pow<4>(y)*r;
749 
750 // //inner functions
751 // case 15: return 20.*x*y*pow<3>(r);
752 // case 16: return 30.*x*pow<2>(y)*pow<2>(r);
753 // case 17: return 30.*pow<2>(x)*y*pow<2>(r);
754 // case 18: return 20.*x*pow<3>(y)*r;
755 // case 19: return 20.*pow<3>(x)*y*r;
756 // case 20: return 30.*pow<2>(x)*pow<2>(y)*r;
757 
758 // }
759 // }
760 
761 
762 // // Bernstein shape functions on the quadrilateral.
763 // case QUAD8:
764 // case QUAD9:
765 // {
766 // // Compute quad shape functions as a tensor-product
767 // const Real xi = p(0);
768 // const Real eta = p(1);
769 
770 // libmesh_assert_less (i, 36);
771 
772 // // 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
773 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
774 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
775 
776 // unsigned int i0=i0_reg[i];
777 // unsigned int i1=i1_reg[i];
778 
779 // if((i>= 4&&i<= 7) && elem->node(0) > elem->node(1)) i0=7-i0;
780 // if((i>= 8&&i<=11) && elem->node(1) > elem->node(2)) i1=7-i1;
781 // if((i>=12&&i<=15) && elem->node(3) > elem->node(2)) i0=7-i0;
782 // if((i>=16&&i<=19) && elem->node(0) > elem->node(3)) i1=7-i1;
783 
784 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)*
785 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta));
786 
787 // }
788 // default:
789 // libmesh_error();
790 // }
791 
792 // case SIXTH:
793 // switch(type)
794 // {
795 // case TRI6:
796 // {
797 // const Real x=p(0);
798 // const Real y=p(1);
799 // const Real r=1-x-y;
800 // unsigned int shape=i;
801 
802 // libmesh_assert_less (i, 28);
803 
804 // if((i>= 3&&i<= 7) && elem->node(0) > elem->node(1))shape=10-i;
805 // if((i>= 8&&i<=12) && elem->node(1) > elem->node(2))shape=20-i;
806 // if((i>=13&&i<=17) && elem->node(0) > elem->node(2))shape=30-i;
807 
808 // switch(shape)
809 // {
810 // //point functions
811 // case 0: return pow<6>(r);
812 // case 1: return pow<6>(x);
813 // case 2: return pow<6>(y);
814 
815 // //edge functions
816 // case 3: return 6.*x *pow<5>(r);
817 // case 4: return 15.*pow<2>(x)*pow<4>(r);
818 // case 5: return 20.*pow<3>(x)*pow<3>(r);
819 // case 6: return 15.*pow<4>(x)*pow<2>(r);
820 // case 7: return 6.*pow<5>(x)*r;
821 
822 // case 8: return 6.*y *pow<5>(x);
823 // case 9: return 15.*pow<2>(y)*pow<4>(x);
824 // case 10: return 20.*pow<3>(y)*pow<3>(x);
825 // case 11: return 15.*pow<4>(y)*pow<2>(x);
826 // case 12: return 6.*pow<5>(y)*x;
827 
828 // case 13: return 6.*y *pow<5>(r);
829 // case 14: return 15.*pow<2>(y)*pow<4>(r);
830 // case 15: return 20.*pow<3>(y)*pow<3>(r);
831 // case 16: return 15.*pow<4>(y)*pow<2>(r);
832 // case 17: return 6.*pow<5>(y)*r;
833 
834 // //inner functions
835 // case 18: return 30.*x*y*pow<4>(r);
836 // case 19: return 60.*x*pow<2>(y)*pow<3>(r);
837 // case 20: return 60.* pow<2>(x)*y*pow<3>(r);
838 // case 21: return 60.*x*pow<3>(y)*pow<2>(r);
839 // case 22: return 60.*pow<3>(x)*y*pow<2>(r);
840 // case 23: return 90.*pow<2>(x)*pow<2>(y)*pow<2>(r);
841 // case 24: return 30.*x*pow<4>(y)*r;
842 // case 25: return 60.*pow<2>(x)*pow<3>(y)*r;
843 // case 26: return 60.*pow<3>(x)*pow<2>(y)*r;
844 // case 27: return 30.*pow<4>(x)*y*r;
845 
846 // } // switch shape
847 // } // case TRI6
848 
849 
850 
851 // // Bernstein shape functions on the quadrilateral.
852 // case QUAD8:
853 // case QUAD9:
854 // {
855 // // Compute quad shape functions as a tensor-product
856 // const Real xi = p(0);
857 // const Real eta = p(1);
858 
859 // libmesh_assert_less (i, 49);
860 
861 // // 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 36 37 38 39 40 41 42 43 44 45 46 47 48
862 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
863 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
864 
865 // unsigned int i0=i0_reg[i];
866 // unsigned int i1=i1_reg[i];
867 
868 // if((i>= 4&&i<= 8) && elem->node(0) > elem->node(1)) i0=8-i0;
869 // if((i>= 9&&i<=13) && elem->node(1) > elem->node(2)) i1=8-i1;
870 // if((i>=14&&i<=18) && elem->node(3) > elem->node(2)) i0=8-i0;
871 // if((i>=19&&i<=23) && elem->node(0) > elem->node(3)) i1=8-i1;
872 
873 // return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0, xi)*
874 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1, eta));
875 
876 // } // case QUAD8/9
877 
878 // default:
879 // libmesh_error();
880 
881 // }
882 // // 7th-order Bernstein.
883 // case SEVENTH:
884 // {
885 // switch (type)
886 // {
887 
888 // // Szabo-Babuska shape functions on the quadrilateral.
889 // case QUAD9:
890 // {
891 // // Compute quad shape functions as a tensor-product
892 // const Real xi = p(0);
893 // const Real eta = p(1);
894 
895 // libmesh_assert_less (i, 64);
896 
897 // // 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
898 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
899 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
900 
901 // Real f=1.;
902 
903 // switch(i)
904 // {
905 // case 5: // edge 0 nodes
906 // case 7:
907 // case 9:
908 // if (elem->node(0) > elem->node(1))f = -1.;
909 // break;
910 // case 11: // edge 1 nodes
911 // case 13:
912 // case 15:
913 // if (elem->node(1) > elem->node(2))f = -1.;
914 // break;
915 // case 17: // edge 2 nodes
916 // case 19:
917 // case 21:
918 // if (elem->node(3) > elem->node(2))f = -1.;
919 // break;
920 // case 23: // edge 3 nodes
921 // case 25:
922 // case 27:
923 // if (elem->node(0) > elem->node(3))f = -1.;
924 // break;
925 // }
926 
927 // return f*(FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
928 // FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta));
929 
930 // } // case QUAD8/QUAD9
931 
932 // default:
933 // libmesh_error();
934 
935 // } // switch type
936 
937 // } // case SEVENTH
938 
939 
940 // // by default throw an error
941 // default:
942 // libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
943 // libmesh_error();
944 // }
945 
946  libmesh_error();
947  return 0.;
948 }
949 
950 
951 
952 template <>
954  const Order,
955  const unsigned int,
956  const unsigned int,
957  const Point&)
958 {
959  libMesh::err << "Bernstein polynomials require the element type\n"
960  << "because edge orientation is needed."
961  << std::endl;
962 
963  libmesh_error();
964  return 0.;
965 }
966 
967 
968 
969 template <>
971  const Order order,
972  const unsigned int i,
973  const unsigned int j,
974  const Point& p)
975 {
976  libmesh_assert(elem);
977 
978  const ElemType type = elem->type();
979 
980  const Order totalorder = static_cast<Order>(order + elem->p_level());
981 
982  switch (type)
983  {
984  // Hierarchic shape functions on the quadrilateral.
985  case QUAD4:
986  case QUAD9:
987  {
988  // Compute quad shape functions as a tensor-product
989  const Real xi = p(0);
990  const Real eta = p(1);
991 
992  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
993 
994  unsigned int i0, i1;
995 
996  // Vertex DoFs
997  if (i == 0)
998  { i0 = 0; i1 = 0; }
999  else if (i == 1)
1000  { i0 = 1; i1 = 0; }
1001  else if (i == 2)
1002  { i0 = 1; i1 = 1; }
1003  else if (i == 3)
1004  { i0 = 0; i1 = 1; }
1005 
1006 
1007  // Edge DoFs
1008  else if (i < totalorder + 3u)
1009  { i0 = i - 2; i1 = 0; }
1010  else if (i < 2u*totalorder + 2)
1011  { i0 = 1; i1 = i - totalorder - 1; }
1012  else if (i < 3u*totalorder + 1)
1013  { i0 = i - 2u*totalorder; i1 = 1; }
1014  else if (i < 4u*totalorder)
1015  { i0 = 0; i1 = i - 3u*totalorder + 1; }
1016  // Interior DoFs
1017  else
1018  {
1019  unsigned int basisnum = i - 4*totalorder;
1020  i0 = square_number_column[basisnum] + 2;
1021  i1 = square_number_row[basisnum] + 2;
1022  }
1023 
1024 
1025  // Flip odd degree of freedom values if necessary
1026  // to keep continuity on sides
1027  if ((i>= 4 && i<= 4+ totalorder-2u) && elem->point(0) > elem->point(1)) i0=totalorder+2-i0; //
1028  else if((i>= 4+ totalorder-1u && i<= 4+2*totalorder-3u) && elem->point(1) > elem->point(2)) i1=totalorder+2-i1;
1029  else if((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->point(3) > elem->point(2)) i0=totalorder+2-i0;
1030  else if((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->point(0) > elem->point(3)) i1=totalorder+2-i1;
1031 
1032  switch (j)
1033  {
1034  // d()/dxi
1035  case 0:
1036  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
1037  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta));
1038 
1039  // d()/deta
1040  case 1:
1041  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)*
1042  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
1043 
1044  }
1045  }
1046 
1047  // Bernstein shape functions on the 8-noded quadrilateral
1048  // is handled separately.
1049  case QUAD8:
1050  {
1051  libmesh_assert_less (totalorder, 3);
1052 
1053  const Real xi = p(0);
1054  const Real eta = p(1);
1055 
1056  libmesh_assert_less (i, 8);
1057 
1058  // 0 1 2 3 4 5 6 7 8
1059  static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
1060  static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
1061  static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
1062  switch (j)
1063  {
1064  // d()/dxi
1065  case 0:
1066  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1067  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)
1068  +scal[i]*
1069  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)*
1070  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[8], eta));
1071 
1072  // d()/deta
1073  case 1:
1074  return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1075  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)
1076  +scal[i]*
1077  FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[8], xi)*
1078  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta));
1079 
1080  default:
1081  libmesh_error();
1082  }
1083  }
1084 
1085  case TRI3:
1086  libmesh_assert_less (totalorder, 2);
1087  case TRI6:
1088  {
1089  // I have been lazy here and am using finite differences
1090  // to compute the derivatives!
1091  const Real eps = 1.e-6;
1092 
1093  switch (j)
1094  {
1095  // d()/dxi
1096  case 0:
1097  {
1098  const Point pp(p(0)+eps, p(1));
1099  const Point pm(p(0)-eps, p(1));
1100 
1101  return (FE<2,BERNSTEIN>::shape(elem, totalorder, i, pp) -
1102  FE<2,BERNSTEIN>::shape(elem, totalorder, i, pm))/2./eps;
1103  }
1104 
1105  // d()/deta
1106  case 1:
1107  {
1108  const Point pp(p(0), p(1)+eps);
1109  const Point pm(p(0), p(1)-eps);
1110 
1111  return (FE<2,BERNSTEIN>::shape(elem, totalorder, i, pp) -
1112  FE<2,BERNSTEIN>::shape(elem, totalorder, i, pm))/2./eps;
1113  }
1114 
1115 
1116  default:
1117  libmesh_error();
1118  }
1119  }
1120 
1121  default:
1122  {
1123  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
1124  libmesh_error();
1125  }
1126 
1127  }
1128 
1129  // old code
1130 // switch (totalorder)
1131 // {
1132 
1133 // // 1st order Bernsteins are aquivalent to Lagrange elements.
1134 // case FIRST:
1135 // {
1136 // switch (type)
1137 // {
1138 // // Bernstein shape functions on the triangle.
1139 // case TRI6:
1140 // {
1141 // // I have been lazy here and am using finite differences
1142 // // to compute the derivatives!
1143 // const Real eps = 1.e-6;
1144 
1145 // libmesh_assert_less (i, 3);
1146 // libmesh_assert_less (j, 2);
1147 
1148 // switch (j)
1149 // {
1150 // // d()/dxi
1151 // case 0:
1152 // {
1153 // const Point pp(p(0)+eps, p(1));
1154 // const Point pm(p(0)-eps, p(1));
1155 
1156 // return (FE<2,BERNSTEIN>::shape(elem, order, i, pp) -
1157 // FE<2,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1158 // }
1159 
1160 // // d()/deta
1161 // case 1:
1162 // {
1163 // const Point pp(p(0), p(1)+eps);
1164 // const Point pm(p(0), p(1)-eps);
1165 
1166 // return (FE<2,BERNSTEIN>::shape(elem, order, i, pp) -
1167 // FE<2,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1168 // }
1169 
1170 
1171 // default:
1172 // libmesh_error();
1173 // }
1174 // }
1175 
1176 // // Bernstein shape functions on the quadrilateral.
1177 // case QUAD8:
1178 // case QUAD9:
1179 // {
1180 // // Compute quad shape functions as a tensor-product
1181 // const Real xi = p(0);
1182 // const Real eta = p(1);
1183 
1184 // libmesh_assert_less (i, 4);
1185 
1186 // // 0 1 2 3
1187 // static const unsigned int i0[] = {0, 1, 1, 0};
1188 // static const unsigned int i1[] = {0, 0, 1, 1};
1189 
1190 // switch (j)
1191 // {
1192 // // d()/dxi
1193 // case 0:
1194 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1195 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta));
1196 
1197 // // d()/deta
1198 // case 1:
1199 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1200 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1201 
1202 // default:
1203 // libmesh_error();
1204 // }
1205 // }
1206 
1207 // default:
1208 // libmesh_error();
1209 // }
1210 // }
1211 
1212 // // second order Bernsteins.
1213 // case SECOND:
1214 // {
1215 // switch (type)
1216 // {
1217 // // Bernstein shape functions on the triangle.
1218 // case TRI6:
1219 // {
1220 // // I have been lazy here and am using finite differences
1221 // // to compute the derivatives!
1222 // const Real eps = 1.e-6;
1223 
1224 // libmesh_assert_less (i, 6);
1225 // libmesh_assert_less (j, 2);
1226 
1227 // switch (j)
1228 // {
1229 // // d()/dxi
1230 // case 0:
1231 // {
1232 // const Point pp(p(0)+eps, p(1));
1233 // const Point pm(p(0)-eps, p(1));
1234 
1235 // return (FE<2,BERNSTEIN>::shape(elem, order, i, pp) -
1236 // FE<2,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1237 // }
1238 
1239 // // d()/deta
1240 // case 1:
1241 // {
1242 // const Point pp(p(0), p(1)+eps);
1243 // const Point pm(p(0), p(1)-eps);
1244 
1245 // return (FE<2,BERNSTEIN>::shape(elem, order, i, pp) -
1246 // FE<2,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
1247 // }
1248 
1249 
1250 // default:
1251 // libmesh_error();
1252 // }
1253 // }
1254 
1255 // // Bernstein shape functions on the 8-noded quadrilateral.
1256 // case QUAD8:
1257 // {
1258 // const Real xi = p(0);
1259 // const Real eta = p(1);
1260 
1261 // libmesh_assert_less (i, 8);
1262 
1263 // // 0 1 2 3 4 5 6 7 8
1264 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
1265 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
1266 // static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
1267 // switch (j)
1268 // {
1269 // // d()/dxi
1270 // case 0:
1271 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1272 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta)
1273 // +scal[i]*
1274 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)*
1275 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[8], eta));
1276 
1277 // // d()/deta
1278 // case 1:
1279 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1280 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)
1281 // +scal[i]*
1282 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[8], xi)*
1283 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta));
1284 
1285 // default:
1286 // libmesh_error();
1287 // }
1288 // }
1289 
1290 
1291 // // Bernstein shape functions on the 9-noded quadrilateral.
1292 // case QUAD9:
1293 // {
1294 // // Compute quad shape functions as a tensor-product
1295 // const Real xi = p(0);
1296 // const Real eta = p(1);
1297 
1298 // libmesh_assert_less (i, 9);
1299 
1300 // // 0 1 2 3 4 5 6 7 8
1301 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
1302 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
1303 
1304 // switch (j)
1305 // {
1306 // // d()/dxi
1307 // case 0:
1308 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1309 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta));
1310 
1311 // // d()/deta
1312 // case 1:
1313 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1314 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1315 
1316 // default:
1317 // libmesh_error();
1318 // }
1319 // }
1320 
1321 // default:
1322 // libmesh_error();
1323 // }
1324 // }
1325 
1326 
1327 // // 3rd-order Bernsteins.
1328 // case THIRD:
1329 // {
1330 // switch (type)
1331 // {
1332 // // Bernstein shape functions on the triangle.
1333 // case TRI6:
1334 // {
1335 // // I have been lazy here and am using finite differences
1336 // // to compute the derivatives!
1337 // const Real eps = 1.e-6;
1338 
1339 // libmesh_assert_less (i, 10);
1340 // libmesh_assert_less (j, 2);
1341 
1342 
1343 // unsigned int shape=i;
1344 
1345 // /**
1346 // * Note, since the derivatives are computed using FE::shape
1347 // * the continuity along neighboring elements (shape_flip)
1348 // * is already considered.
1349 // */
1350 
1351 
1352 // switch (j)
1353 // {
1354 // // d()/dxi
1355 // case 0:
1356 // {
1357 // const Point pp(p(0)+eps, p(1));
1358 // const Point pm(p(0)-eps, p(1));
1359 
1360 // return (FE<2,BERNSTEIN>::shape(elem, totalorder, shape, pp) -
1361 // FE<2,BERNSTEIN>::shape(elem, totalorder, shape, pm))/2./eps;
1362 // }
1363 
1364 // // d()/deta
1365 // case 1:
1366 // {
1367 // const Point pp(p(0), p(1)+eps);
1368 // const Point pm(p(0), p(1)-eps);
1369 
1370 // return (FE<2,BERNSTEIN>::shape(elem, totalorder, shape, pp) -
1371 // FE<2,BERNSTEIN>::shape(elem, totalorder, shape, pm))/2./eps;
1372 // }
1373 
1374 
1375 // default:
1376 // libmesh_error();
1377 // }
1378 // }
1379 
1380 
1381 
1382 // // Bernstein shape functions on the quadrilateral.
1383 // case QUAD8:
1384 // case QUAD9:
1385 // {
1386 // // Compute quad shape functions as a tensor-product
1387 // const Real xi = p(0);
1388 // const Real eta = p(1);
1389 
1390 // libmesh_assert_less (i, 16);
1391 
1392 // // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1393 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3};
1394 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3};
1395 
1396 // unsigned int i0=i0_reg[i];
1397 // unsigned int i1=i1_reg[i];
1398 
1399 // if((i== 4||i== 5) && elem->node(0) > elem->node(1)) i0=5-i0;
1400 // if((i== 6||i== 7) && elem->node(1) > elem->node(2)) i1=5-i1;
1401 // if((i== 8||i== 9) && elem->node(3) > elem->node(2)) i0=5-i0;
1402 // if((i==10||i==11) && elem->node(0) > elem->node(3)) i1=5-i1;
1403 
1404 // switch (j)
1405 // {
1406 // // d()/dxi
1407 // case 0:
1408 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
1409 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta));
1410 
1411 // // d()/deta
1412 // case 1:
1413 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)*
1414 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
1415 
1416 // default:
1417 // libmesh_error();
1418 // }
1419 // }
1420 
1421 // default:
1422 // libmesh_error();
1423 // }
1424 // }
1425 
1426 
1427 
1428 
1429 // // 4th-order Bernsteins.
1430 // case FOURTH:
1431 // {
1432 // switch (type)
1433 // {
1434 // // Bernstein shape functions on the triangle.
1435 // case TRI6:
1436 // {
1437 // // I have been lazy here and am using finite differences
1438 // // to compute the derivatives!
1439 // const Real eps = 1.e-6;
1440 
1441 // libmesh_assert_less (i, 15);
1442 // libmesh_assert_less (j, 2);
1443 
1444 // unsigned int shape=i;
1445 
1446 // switch (j)
1447 // {
1448 // // d()/dxi
1449 // case 0:
1450 // {
1451 // const Point pp(p(0)+eps, p(1));
1452 // const Point pm(p(0)-eps, p(1));
1453 
1454 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) -
1455 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps;
1456 // }
1457 
1458 // // d()/deta
1459 // case 1:
1460 // {
1461 // const Point pp(p(0), p(1)+eps);
1462 // const Point pm(p(0), p(1)-eps);
1463 
1464 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) -
1465 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps;
1466 // }
1467 
1468 
1469 // default:
1470 // libmesh_error();
1471 // }
1472 // }
1473 
1474 
1475 
1476 // // Bernstein shape functions on the quadrilateral.
1477 // case QUAD8:
1478 // case QUAD9:
1479 // {
1480 // // Compute quad shape functions as a tensor-product
1481 // const Real xi = p(0);
1482 // const Real eta = p(1);
1483 
1484 // libmesh_assert_less (i, 25);
1485 
1486 // // 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
1487 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
1488 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
1489 
1490 // unsigned int i0=i0_reg[i];
1491 // unsigned int i1=i1_reg[i];
1492 
1493 // if((i== 4||i== 6) && elem->node(0) > elem->node(1)) i0=6-i0;
1494 // if((i== 7||i== 9) && elem->node(1) > elem->node(2)) i1=6-i1;
1495 // if((i==10||i==12) && elem->node(3) > elem->node(2)) i0=6-i0;
1496 // if((i==13||i==15) && elem->node(0) > elem->node(3)) i1=6-i1;
1497 
1498 
1499 // switch (j)
1500 // {
1501 // // d()/dxi
1502 // case 0:
1503 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
1504 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta));
1505 
1506 // // d()/deta
1507 // case 1:
1508 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)*
1509 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
1510 
1511 // default:
1512 // libmesh_error();
1513 // }
1514 // }
1515 
1516 // default:
1517 // libmesh_error();
1518 // }
1519 // }
1520 
1521 
1522 
1523 
1524 // // 5th-order Bernsteins.
1525 // case FIFTH:
1526 // {
1527 // // Bernstein shape functions on the quadrilateral.
1528 // switch (type)
1529 // {
1530 // // Bernstein shape functions on the triangle.
1531 // case TRI6:
1532 // {
1533 // // I have been lazy here and am using finite differences
1534 // // to compute the derivatives!
1535 // const Real eps = 1.e-6;
1536 
1537 // libmesh_assert_less (i, 21);
1538 // libmesh_assert_less (j, 2);
1539 
1540 // unsigned int shape=i;
1541 
1542 // switch (j)
1543 // {
1544 // // d()/dxi
1545 // case 0:
1546 // {
1547 // const Point pp(p(0)+eps, p(1));
1548 // const Point pm(p(0)-eps, p(1));
1549 
1550 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) -
1551 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps;
1552 // }
1553 
1554 // // d()/deta
1555 // case 1:
1556 // {
1557 // const Point pp(p(0), p(1)+eps);
1558 // const Point pm(p(0), p(1)-eps);
1559 
1560 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) -
1561 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps;
1562 // }
1563 
1564 
1565 // default:
1566 // libmesh_error();
1567 // }
1568 // } // case TRI6
1569 
1570 
1571 
1572 // case QUAD8:
1573 // case QUAD9:
1574 // {
1575 // // Compute quad shape functions as a tensor-product
1576 // const Real xi = p(0);
1577 // const Real eta = p(1);
1578 
1579 // libmesh_assert_less (i, 36);
1580 
1581 // // 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
1582 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
1583 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
1584 
1585 // unsigned int i0=i0_reg[i];
1586 // unsigned int i1=i1_reg[i];
1587 
1588 // if((i>= 4&&i<= 7) && elem->node(0) > elem->node(1)) i0=7-i0;
1589 // if((i>= 8&&i<=11) && elem->node(1) > elem->node(2)) i1=7-i1;
1590 // if((i>=12&&i<=15) && elem->node(3) > elem->node(2)) i0=7-i0;
1591 // if((i>=16&&i<=19) && elem->node(0) > elem->node(3)) i1=7-i1;
1592 
1593 
1594 // switch (j)
1595 // {
1596 // // d()/dxi
1597 // case 0:
1598 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
1599 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta));
1600 
1601 // // d()/deta
1602 // case 1:
1603 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)*
1604 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
1605 
1606 // default:
1607 // libmesh_error();
1608 // }
1609 // } // case QUAD8/9
1610 
1611 // default:
1612 // libmesh_error();
1613 // }
1614 
1615 // }
1616 
1617 
1618 // // 6th-order Bernsteins.
1619 // case SIXTH:
1620 // {
1621 // // Bernstein shape functions on the quadrilateral.
1622 // switch (type)
1623 // {
1624 // // Bernstein shape functions on the triangle.
1625 // case TRI6:
1626 // {
1627 // // I have been lazy here and am using finite differences
1628 // // to compute the derivatives!
1629 // const Real eps = 1.e-6;
1630 
1631 // libmesh_assert_less (i, 28);
1632 // libmesh_assert_less (j, 2);
1633 
1634 // unsigned int shape=i;
1635 
1636 // switch (j)
1637 // {
1638 // // d()/dxi
1639 // case 0:
1640 // {
1641 // const Point pp(p(0)+eps, p(1));
1642 // const Point pm(p(0)-eps, p(1));
1643 
1644 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) -
1645 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps;
1646 // }
1647 
1648 // // d()/deta
1649 // case 1:
1650 // {
1651 // const Point pp(p(0), p(1)+eps);
1652 // const Point pm(p(0), p(1)-eps);
1653 
1654 // return (FE<2,BERNSTEIN>::shape(elem, order, shape, pp) -
1655 // FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps;
1656 // }
1657 
1658 
1659 // default:
1660 // libmesh_error();
1661 // }
1662 // } // case TRI6
1663 
1664 
1665 
1666 // case QUAD8:
1667 // case QUAD9:
1668 // {
1669 // // Compute quad shape functions as a tensor-product
1670 // const Real xi = p(0);
1671 // const Real eta = p(1);
1672 
1673 // libmesh_assert_less (i, 49);
1674 
1675 // // 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 36 37 38 39 40 41 42 43 44 45 46 47 48
1676 // static const unsigned int i0_reg[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
1677 // static const unsigned int i1_reg[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
1678 
1679 // unsigned int i0=i0_reg[i];
1680 // unsigned int i1=i1_reg[i];
1681 
1682 // if((i>= 4&&i<= 8) && elem->node(0) > elem->node(1)) i0=8-i0;
1683 // if((i>= 9&&i<=13) && elem->node(1) > elem->node(2)) i1=8-i1;
1684 // if((i>=14&&i<=18) && elem->node(3) > elem->node(2)) i0=8-i0;
1685 // if((i>=19&&i<=23) && elem->node(0) > elem->node(3)) i1=8-i1;
1686 
1687 
1688 // switch (j)
1689 // {
1690 // // d()/dxi
1691 // case 0:
1692 // return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
1693 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1, eta));
1694 
1695 // // d()/deta
1696 // case 1:
1697 // return (FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0, xi)*
1698 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
1699 
1700 // default:
1701 // libmesh_error();
1702 // }
1703 // } // case QUAD8/9
1704 
1705 // default:
1706 // libmesh_error();
1707 // }
1708 
1709 // // 7th-order Bernstein.
1710 // case SEVENTH:
1711 // {
1712 // // Szabo-Babuska shape functions on the quadrilateral.
1713 // switch (type)
1714 // {
1715 
1716 
1717 // case QUAD9:
1718 // {
1719 // // Compute quad shape functions as a tensor-product
1720 // const Real xi = p(0);
1721 // const Real eta = p(1);
1722 
1723 // libmesh_assert_less (i, 64);
1724 
1725 // // 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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
1726 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
1727 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
1728 
1729 // Real f=1.;
1730 
1731 // switch(i)
1732 // {
1733 // case 5: // edge 0 nodes
1734 // case 7:
1735 // case 9:
1736 // if (elem->node(0) > elem->node(1))f = -1.;
1737 // break;
1738 // case 11: // edge 1 nodes
1739 // case 13:
1740 // case 15:
1741 // if (elem->node(1) > elem->node(2))f = -1.;
1742 // break;
1743 // case 17: // edge 2 nodes
1744 // case 19:
1745 // case 21:
1746 // if (elem->node(3) > elem->node(2))f = -1.;
1747 // break;
1748 // case 23: // edge 3 nodes
1749 // case 25:
1750 // case 27:
1751 // if (elem->node(0) > elem->node(3))f = -1.;
1752 // break;
1753 // }
1754 
1755 
1756 // switch (j)
1757 // {
1758 // // d()/dxi
1759 // case 0:
1760 // return f*(FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
1761 // FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i1[i], eta));
1762 
1763 // // d()/deta
1764 // case 1:
1765 // return f*(FE<1,BERNSTEIN>::shape (EDGE3, totalorder, i0[i], xi)*
1766 // FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
1767 
1768 // default:
1769 // libmesh_error();
1770 // }
1771 // }
1772 
1773 // default:
1774 // libmesh_error();
1775 // }
1776 // }
1777 
1778 // }
1779 // // by default throw an error;call the orientation-independent shape functions
1780 // default:
1781 // libMesh::err << "ERROR: Unsupported polynomial order!" << std::endl;
1782 // libmesh_error();
1783 // }
1784 
1785 
1786  libmesh_error();
1787  return 0.;
1788 }
1789 
1790 
1791 
1792 template <>
1794  const Order,
1795  const unsigned int,
1796  const unsigned int,
1797  const Point&)
1798 {
1799  static bool warning_given = false;
1800 
1801  if (!warning_given)
1802  libMesh::err << "Second derivatives for Bernstein elements "
1803  << "are not yet implemented!"
1804  << std::endl;
1805 
1806  warning_given = true;
1807  return 0.;
1808 }
1809 
1810 
1811 
1812 template <>
1814  const Order,
1815  const unsigned int,
1816  const unsigned int,
1817  const Point&)
1818 {
1819  static bool warning_given = false;
1820 
1821  if (!warning_given)
1822  libMesh::err << "Second derivatives for Bernstein elements "
1823  << "are not yet implemented!"
1824  << std::endl;
1825 
1826  warning_given = true;
1827  return 0.;
1828 }
1829 
1830 } // namespace libMesh
1831 
1832 
1833 #endif // LIBMESH_ENABLE_HIGHER_ORDER_SHAPES

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo