libMesh
fe_monomial_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 // C++ includes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 
31 template <>
33  const Order libmesh_dbg_var(order),
34  const unsigned int i,
35  const Point & p)
36 {
37 #if LIBMESH_DIM > 1
38 
39  libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
40  (static_cast<unsigned int>(order)+2)/2);
41 
42  const Real xi = p(0);
43  const Real eta = p(1);
44 
45  switch (i)
46  {
47  // constant
48  case 0:
49  return 1.;
50 
51  // linear
52  case 1:
53  return xi;
54 
55  case 2:
56  return eta;
57 
58  // quadratics
59  case 3:
60  return xi*xi;
61 
62  case 4:
63  return xi*eta;
64 
65  case 5:
66  return eta*eta;
67 
68  // cubics
69  case 6:
70  return xi*xi*xi;
71 
72  case 7:
73  return xi*xi*eta;
74 
75  case 8:
76  return xi*eta*eta;
77 
78  case 9:
79  return eta*eta*eta;
80 
81  // quartics
82  case 10:
83  return xi*xi*xi*xi;
84 
85  case 11:
86  return xi*xi*xi*eta;
87 
88  case 12:
89  return xi*xi*eta*eta;
90 
91  case 13:
92  return xi*eta*eta*eta;
93 
94  case 14:
95  return eta*eta*eta*eta;
96 
97  default:
98  unsigned int o = 0;
99  for (; i >= (o+1)*(o+2)/2; o++) { }
100  unsigned int ny = i - (o*(o+1)/2);
101  unsigned int nx = o - ny;
102  Real val = 1.;
103  for (unsigned int index=0; index != nx; index++)
104  val *= xi;
105  for (unsigned int index=0; index != ny; index++)
106  val *= eta;
107  return val;
108  }
109 
110  libmesh_error_msg("We'll never get here!");
111  return 0.;
112 
113 #endif
114 }
115 
116 
117 
118 template <>
120  const Order order,
121  const unsigned int i,
122  const Point & p)
123 {
124  libmesh_assert(elem);
125 
126  // by default call the orientation-independent shape functions
127  return FE<2,MONOMIAL>::shape(elem->type(), static_cast<Order>(order + elem->p_level()), i, p);
128 }
129 
130 
131 
132 template <>
134  const Order libmesh_dbg_var(order),
135  const unsigned int i,
136  const unsigned int j,
137  const Point & p)
138 {
139 #if LIBMESH_DIM > 1
140 
141 
142  libmesh_assert_less (j, 2);
143 
144  libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
145  (static_cast<unsigned int>(order)+2)/2);
146 
147  const Real xi = p(0);
148  const Real eta = p(1);
149 
150  // monomials. since they are hierarchic we only need one case block.
151 
152  switch (j)
153  {
154  // d()/dxi
155  case 0:
156  {
157  switch (i)
158  {
159  // constants
160  case 0:
161  return 0.;
162 
163  // linears
164  case 1:
165  return 1.;
166 
167  case 2:
168  return 0.;
169 
170  // quadratics
171  case 3:
172  return 2.*xi;
173 
174  case 4:
175  return eta;
176 
177  case 5:
178  return 0.;
179 
180  // cubics
181  case 6:
182  return 3.*xi*xi;
183 
184  case 7:
185  return 2.*xi*eta;
186 
187  case 8:
188  return eta*eta;
189 
190  case 9:
191  return 0.;
192 
193  // quartics
194  case 10:
195  return 4.*xi*xi*xi;
196 
197  case 11:
198  return 3.*xi*xi*eta;
199 
200  case 12:
201  return 2.*xi*eta*eta;
202 
203  case 13:
204  return eta*eta*eta;
205 
206  case 14:
207  return 0.;
208 
209  default:
210  unsigned int o = 0;
211  for (; i >= (o+1)*(o+2)/2; o++) { }
212  unsigned int ny = i - (o*(o+1)/2);
213  unsigned int nx = o - ny;
214  Real val = nx;
215  for (unsigned int index=1; index < nx; index++)
216  val *= xi;
217  for (unsigned int index=0; index != ny; index++)
218  val *= eta;
219  return val;
220  }
221  }
222 
223 
224  // d()/deta
225  case 1:
226  {
227  switch (i)
228  {
229  // constants
230  case 0:
231  return 0.;
232 
233  // linears
234  case 1:
235  return 0.;
236 
237  case 2:
238  return 1.;
239 
240  // quadratics
241  case 3:
242  return 0.;
243 
244  case 4:
245  return xi;
246 
247  case 5:
248  return 2.*eta;
249 
250  // cubics
251  case 6:
252  return 0.;
253 
254  case 7:
255  return xi*xi;
256 
257  case 8:
258  return 2.*xi*eta;
259 
260  case 9:
261  return 3.*eta*eta;
262 
263  // quartics
264  case 10:
265  return 0.;
266 
267  case 11:
268  return xi*xi*xi;
269 
270  case 12:
271  return 2.*xi*xi*eta;
272 
273  case 13:
274  return 3.*xi*eta*eta;
275 
276  case 14:
277  return 4.*eta*eta*eta;
278 
279  default:
280  unsigned int o = 0;
281  for (; i >= (o+1)*(o+2)/2; o++) { }
282  unsigned int ny = i - (o*(o+1)/2);
283  unsigned int nx = o - ny;
284  Real val = ny;
285  for (unsigned int index=0; index != nx; index++)
286  val *= xi;
287  for (unsigned int index=1; index < ny; index++)
288  val *= eta;
289  return val;
290  }
291  }
292 
293  default:
294  libmesh_error_msg("Invalid shape function derivative j = " << j);
295  }
296 
297  libmesh_error_msg("We'll never get here!");
298  return 0.;
299 
300 #endif
301 }
302 
303 
304 
305 template <>
307  const Order order,
308  const unsigned int i,
309  const unsigned int j,
310  const Point & p)
311 {
312  libmesh_assert(elem);
313 
314  // by default call the orientation-independent shape functions
315  return FE<2,MONOMIAL>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
316 }
317 
318 
319 
320 template <>
322  const Order libmesh_dbg_var(order),
323  const unsigned int i,
324  const unsigned int j,
325  const Point & p)
326 {
327 #if LIBMESH_DIM > 1
328 
329 
330  libmesh_assert_less_equal (j, 2);
331 
332  libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
333  (static_cast<unsigned int>(order)+2)/2);
334 
335  const Real xi = p(0);
336  const Real eta = p(1);
337 
338  // monomials. since they are hierarchic we only need one case block.
339 
340  switch (j)
341  {
342  // d^2()/dxi^2
343  case 0:
344  {
345  switch (i)
346  {
347  // constants
348  case 0:
349  // linears
350  case 1:
351  case 2:
352  return 0.;
353 
354  // quadratics
355  case 3:
356  return 2.;
357 
358  case 4:
359  case 5:
360  return 0.;
361 
362  // cubics
363  case 6:
364  return 6.*xi;
365 
366  case 7:
367  return 2.*eta;
368 
369  case 8:
370  case 9:
371  return 0.;
372 
373  // quartics
374  case 10:
375  return 12.*xi*xi;
376 
377  case 11:
378  return 6.*xi*eta;
379 
380  case 12:
381  return 2.*eta*eta;
382 
383  case 13:
384  case 14:
385  return 0.;
386 
387  default:
388  unsigned int o = 0;
389  for (; i >= (o+1)*(o+2)/2; o++) { }
390  unsigned int ny = i - (o*(o+1)/2);
391  unsigned int nx = o - ny;
392  Real val = nx * (nx - 1);
393  for (unsigned int index=2; index < nx; index++)
394  val *= xi;
395  for (unsigned int index=0; index != ny; index++)
396  val *= eta;
397  return val;
398  }
399  }
400 
401  // d^2()/dxideta
402  case 1:
403  {
404  switch (i)
405  {
406  // constants
407  case 0:
408 
409  // linears
410  case 1:
411  case 2:
412  return 0.;
413 
414  // quadratics
415  case 3:
416  return 0.;
417 
418  case 4:
419  return 1.;
420 
421  case 5:
422  return 0.;
423 
424  // cubics
425  case 6:
426  return 0.;
427  case 7:
428  return 2.*xi;
429 
430  case 8:
431  return 2.*eta;
432 
433  case 9:
434  return 0.;
435 
436  // quartics
437  case 10:
438  return 0.;
439 
440  case 11:
441  return 3.*xi*xi;
442 
443  case 12:
444  return 4.*xi*eta;
445 
446  case 13:
447  return 3.*eta*eta;
448 
449  case 14:
450  return 0.;
451 
452  default:
453  unsigned int o = 0;
454  for (; i >= (o+1)*(o+2)/2; o++) { }
455  unsigned int ny = i - (o*(o+1)/2);
456  unsigned int nx = o - ny;
457  Real val = nx * ny;
458  for (unsigned int index=1; index < nx; index++)
459  val *= xi;
460  for (unsigned int index=1; index < ny; index++)
461  val *= eta;
462  return val;
463  }
464  }
465 
466  // d^2()/deta^2
467  case 2:
468  {
469  switch (i)
470  {
471  // constants
472  case 0:
473 
474  // linears
475  case 1:
476  case 2:
477  return 0.;
478 
479  // quadratics
480  case 3:
481  case 4:
482  return 0.;
483 
484  case 5:
485  return 2.;
486 
487  // cubics
488  case 6:
489  return 0.;
490 
491  case 7:
492  return 0.;
493 
494  case 8:
495  return 2.*xi;
496 
497  case 9:
498  return 6.*eta;
499 
500  // quartics
501  case 10:
502  case 11:
503  return 0.;
504 
505  case 12:
506  return 2.*xi*xi;
507 
508  case 13:
509  return 6.*xi*eta;
510 
511  case 14:
512  return 12.*eta*eta;
513 
514  default:
515  unsigned int o = 0;
516  for (; i >= (o+1)*(o+2)/2; o++) { }
517  unsigned int ny = i - (o*(o+1)/2);
518  unsigned int nx = o - ny;
519  Real val = ny * (ny - 1);
520  for (unsigned int index=0; index != nx; index++)
521  val *= xi;
522  for (unsigned int index=2; index < ny; index++)
523  val *= eta;
524  return val;
525  }
526  }
527 
528  default:
529  libmesh_error_msg("Invalid shape function derivative j = " << j);
530  }
531 
532  libmesh_error_msg("We'll never get here!");
533  return 0.;
534 
535 #endif
536 }
537 
538 
539 
540 template <>
542  const Order order,
543  const unsigned int i,
544  const unsigned int j,
545  const Point & p)
546 {
547  libmesh_assert(elem);
548 
549  // by default call the orientation-independent shape functions
550  return FE<2,MONOMIAL>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
551 }
552 
553 
554 } // namespace libMesh
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)
PetscErrorCode Vec Mat libmesh_dbg_var(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
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)