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