fe_monomial_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 // C++ inlcludes
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();
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 
294  libmesh_error();
295  return 0.;
296 
297 #endif
298 }
299 
300 
301 
302 template <>
304  const Order order,
305  const unsigned int i,
306  const unsigned int j,
307  const Point& p)
308 {
309  libmesh_assert(elem);
310 
311  // by default call the orientation-independent shape functions
312  return FE<2,MONOMIAL>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
313 }
314 
315 
316 
317 template <>
319  const Order libmesh_dbg_var(order),
320  const unsigned int i,
321  const unsigned int j,
322  const Point& p)
323 {
324 #if LIBMESH_DIM > 1
325 
326 
327  libmesh_assert_less_equal (j, 2);
328 
329  libmesh_assert_less (i, (static_cast<unsigned int>(order)+1)*
330  (static_cast<unsigned int>(order)+2)/2);
331 
332  const Real xi = p(0);
333  const Real eta = p(1);
334 
335  // monomials. since they are hierarchic we only need one case block.
336 
337  switch (j)
338  {
339  // d^2()/dxi^2
340  case 0:
341  {
342  switch (i)
343  {
344  // constants
345  case 0:
346  // linears
347  case 1:
348  case 2:
349  return 0.;
350 
351  // quadratics
352  case 3:
353  return 2.;
354 
355  case 4:
356  case 5:
357  return 0.;
358 
359  // cubics
360  case 6:
361  return 6.*xi;
362 
363  case 7:
364  return 2.*eta;
365 
366  case 8:
367  case 9:
368  return 0.;
369 
370  // quartics
371  case 10:
372  return 12.*xi*xi;
373 
374  case 11:
375  return 6.*xi*eta;
376 
377  case 12:
378  return 2.*eta*eta;
379 
380  case 13:
381  case 14:
382  return 0.;
383 
384  default:
385  unsigned int o = 0;
386  for (; i >= (o+1)*(o+2)/2; o++) { }
387  unsigned int ny = i - (o*(o+1)/2);
388  unsigned int nx = o - ny;
389  Real val = nx * (nx - 1);
390  for (unsigned int index=2; index < nx; index++)
391  val *= xi;
392  for (unsigned int index=0; index != ny; index++)
393  val *= eta;
394  return val;
395  }
396  }
397 
398  // d^2()/dxideta
399  case 1:
400  {
401  switch (i)
402  {
403  // constants
404  case 0:
405 
406  // linears
407  case 1:
408  case 2:
409  return 0.;
410 
411  // quadratics
412  case 3:
413  return 0.;
414 
415  case 4:
416  return 1.;
417 
418  case 5:
419  return 0.;
420 
421  // cubics
422  case 6:
423  return 0.;
424  case 7:
425  return 2.*xi;
426 
427  case 8:
428  return 2.*eta;
429 
430  case 9:
431  return 0.;
432 
433  // quartics
434  case 10:
435  return 0.;
436 
437  case 11:
438  return 3.*xi*xi;
439 
440  case 12:
441  return 4.*xi*eta;
442 
443  case 13:
444  return 3.*eta*eta;
445 
446  case 14:
447  return 0.;
448 
449  default:
450  unsigned int o = 0;
451  for (; i >= (o+1)*(o+2)/2; o++) { }
452  unsigned int ny = i - (o*(o+1)/2);
453  unsigned int nx = o - ny;
454  Real val = nx * ny;
455  for (unsigned int index=1; index < nx; index++)
456  val *= xi;
457  for (unsigned int index=1; index < ny; index++)
458  val *= eta;
459  return val;
460  }
461  }
462 
463  // d^2()/deta^2
464  case 2:
465  {
466  switch (i)
467  {
468  // constants
469  case 0:
470 
471  // linears
472  case 1:
473  case 2:
474  return 0.;
475 
476  // quadratics
477  case 3:
478  case 4:
479  return 0.;
480 
481  case 5:
482  return 2.;
483 
484  // cubics
485  case 6:
486  return 0.;
487 
488  case 7:
489  return 0.;
490 
491  case 8:
492  return 2.*xi;
493 
494  case 9:
495  return 6.*eta;
496 
497  // quartics
498  case 10:
499  case 11:
500  return 0.;
501 
502  case 12:
503  return 2.*xi*xi;
504 
505  case 13:
506  return 6.*xi*eta;
507 
508  case 14:
509  return 12.*eta*eta;
510 
511  default:
512  unsigned int o = 0;
513  for (; i >= (o+1)*(o+2)/2; o++) { }
514  unsigned int ny = i - (o*(o+1)/2);
515  unsigned int nx = o - ny;
516  Real val = ny * (ny - 1);
517  for (unsigned int index=0; index != nx; index++)
518  val *= xi;
519  for (unsigned int index=2; index < ny; index++)
520  val *= eta;
521  return val;
522  }
523  }
524  }
525 
526  libmesh_error();
527  return 0.;
528 
529 #endif
530 }
531 
532 
533 
534 template <>
536  const Order order,
537  const unsigned int i,
538  const unsigned int j,
539  const Point& p)
540 {
541  libmesh_assert(elem);
542 
543  // by default call the orientation-independent shape functions
544  return FE<2,MONOMIAL>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
545 }
546 
547 
548 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo