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