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