libMesh
fe_xyz_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 
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 
93  // quadratics
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  libmesh_error_msg("We'll never get here!");
145  return 0.;
146 
147 #endif
148 }
149 
150 
151 
152 template <>
154  const Order,
155  const unsigned int,
156  const unsigned int,
157  const Point &)
158 {
159  libmesh_error_msg("XYZ polynomials require the element \nbecause the centroid is needed.");
160  return 0.;
161 }
162 
163 
164 
165 template <>
167  const Order libmesh_dbg_var(order),
168  const unsigned int i,
169  const unsigned int j,
170  const Point & point_in)
171 {
172 #if LIBMESH_DIM > 1
173 
174 
175  libmesh_assert_less (j, 2);
176  libmesh_assert(elem);
177 
178  Point centroid = elem->centroid();
179  Point max_distance = Point(0.,0.,0.);
180  for (unsigned int p = 0; p < elem->n_nodes(); p++)
181  for (unsigned int d = 0; d < 2; d++)
182  {
183  const Real distance = std::abs(centroid(d) - elem->point(p)(d));
184  max_distance(d) = std::max(distance, max_distance(d));
185  }
186 
187  const Real x = point_in(0);
188  const Real y = point_in(1);
189  const Real xc = centroid(0);
190  const Real yc = centroid(1);
191  const Real distx = max_distance(0);
192  const Real disty = max_distance(1);
193  const Real dx = (x - xc)/distx;
194  const Real dy = (y - yc)/disty;
195 
196 #ifndef NDEBUG
197  // totalorder is only used in the assertion below, so
198  // we avoid declaring it when asserts are not active.
199  const unsigned int totalorder = order + elem->p_level();
200 #endif
201  libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
202 
203  // monomials. since they are hierarchic we only need one case block.
204 
205  switch (j)
206  {
207  // d()/dx
208  case 0:
209  {
210  switch (i)
211  {
212  // constants
213  case 0:
214  return 0.;
215 
216  // linears
217  case 1:
218  return 1./distx;
219 
220  case 2:
221  return 0.;
222 
223  // quadratics
224  case 3:
225  return 2.*dx/distx;
226 
227  case 4:
228  return dy/distx;
229 
230  case 5:
231  return 0.;
232 
233  // cubics
234  case 6:
235  return 3.*dx*dx/distx;
236 
237  case 7:
238  return 2.*dx*dy/distx;
239 
240  case 8:
241  return dy*dy/distx;
242 
243  case 9:
244  return 0.;
245 
246  // quartics
247  case 10:
248  return 4.*dx*dx*dx/distx;
249 
250  case 11:
251  return 3.*dx*dx*dy/distx;
252 
253  case 12:
254  return 2.*dx*dy*dy/distx;
255 
256  case 13:
257  return dy*dy*dy/distx;
258 
259  case 14:
260  return 0.;
261 
262  default:
263  unsigned int o = 0;
264  for (; i >= (o+1)*(o+2)/2; o++) { }
265  unsigned int i2 = i - (o*(o+1)/2);
266  Real val = o - i2;
267  for (unsigned int index=i2+1; index < o; index++)
268  val *= dx;
269  for (unsigned int index=0; index != i2; index++)
270  val *= dy;
271  return val/distx;
272  }
273  }
274 
275 
276  // d()/dy
277  case 1:
278  {
279  switch (i)
280  {
281  // constants
282  case 0:
283  return 0.;
284 
285  // linears
286  case 1:
287  return 0.;
288 
289  case 2:
290  return 1./disty;
291 
292  // quadratics
293  case 3:
294  return 0.;
295 
296  case 4:
297  return dx/disty;
298 
299  case 5:
300  return 2.*dy/disty;
301 
302  // cubics
303  case 6:
304  return 0.;
305 
306  case 7:
307  return dx*dx/disty;
308 
309  case 8:
310  return 2.*dx*dy/disty;
311 
312  case 9:
313  return 3.*dy*dy/disty;
314 
315  // quartics
316  case 10:
317  return 0.;
318 
319  case 11:
320  return dx*dx*dx/disty;
321 
322  case 12:
323  return 2.*dx*dx*dy/disty;
324 
325  case 13:
326  return 3.*dx*dy*dy/disty;
327 
328  case 14:
329  return 4.*dy*dy*dy/disty;
330 
331  default:
332  unsigned int o = 0;
333  for (; i >= (o+1)*(o+2)/2; o++) { }
334  unsigned int i2 = i - (o*(o+1)/2);
335  Real val = i2;
336  for (unsigned int index=i2; index != o; index++)
337  val *= dx;
338  for (unsigned int index=1; index <= i2; index++)
339  val *= dy;
340  return val/disty;
341  }
342  }
343 
344 
345  default:
346  libmesh_error_msg("Invalid j = " << j);
347  }
348 
349  libmesh_error_msg("We'll never get here!");
350  return 0.;
351 
352 #endif
353 }
354 
355 
356 
357 template <>
359  const Order,
360  const unsigned int,
361  const unsigned int,
362  const Point &)
363 {
364  libmesh_error_msg("XYZ polynomials require the element \nbecause the centroid is needed.");
365  return 0.;
366 }
367 
368 
369 
370 template <>
372  const Order libmesh_dbg_var(order),
373  const unsigned int i,
374  const unsigned int j,
375  const Point & point_in)
376 {
377 #if LIBMESH_DIM > 1
378 
379  libmesh_assert_less_equal (j, 2);
380  libmesh_assert(elem);
381 
382  Point centroid = elem->centroid();
383  Point max_distance = Point(0.,0.,0.);
384  for (unsigned int p = 0; p < elem->n_nodes(); p++)
385  for (unsigned int d = 0; d < 2; d++)
386  {
387  const Real distance = std::abs(centroid(d) - elem->point(p)(d));
388  max_distance(d) = std::max(distance, max_distance(d));
389  }
390 
391  const Real x = point_in(0);
392  const Real y = point_in(1);
393  const Real xc = centroid(0);
394  const Real yc = centroid(1);
395  const Real distx = max_distance(0);
396  const Real disty = max_distance(1);
397  const Real dx = (x - xc)/distx;
398  const Real dy = (y - yc)/disty;
399  const Real dist2x = pow(distx,2.);
400  const Real dist2y = pow(disty,2.);
401  const Real distxy = distx * disty;
402 
403 #ifndef NDEBUG
404  // totalorder is only used in the assertion below, so
405  // we avoid declaring it when asserts are not active.
406  const unsigned int totalorder = order + elem->p_level();
407 #endif
408  libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
409 
410  // monomials. since they are hierarchic we only need one case block.
411 
412  switch (j)
413  {
414  // d^2()/dx^2
415  case 0:
416  {
417  switch (i)
418  {
419  // constants
420  case 0:
421  // linears
422  case 1:
423  case 2:
424  return 0.;
425 
426  // quadratics
427  case 3:
428  return 2./dist2x;
429 
430  case 4:
431  case 5:
432  return 0.;
433 
434  // cubics
435  case 6:
436  return 6.*dx/dist2x;
437 
438  case 7:
439  return 2.*dy/dist2x;
440 
441  case 8:
442  case 9:
443  return 0.;
444 
445  // quartics
446  case 10:
447  return 12.*dx*dx/dist2x;
448 
449  case 11:
450  return 6.*dx*dy/dist2x;
451 
452  case 12:
453  return 2.*dy*dy/dist2x;
454 
455  case 13:
456  case 14:
457  return 0.;
458 
459  default:
460  unsigned int o = 0;
461  for (; i >= (o+1)*(o+2)/2; o++) { }
462  unsigned int i2 = i - (o*(o+1)/2);
463  Real val = (o - i2) * (o - i2 - 1);
464  for (unsigned int index=i2+2; index < o; index++)
465  val *= dx;
466  for (unsigned int index=0; index != i2; index++)
467  val *= dy;
468  return val/dist2x;
469  }
470  }
471 
472  // d^2()/dxdy
473  case 1:
474  {
475  switch (i)
476  {
477  // constants
478  case 0:
479 
480  // linears
481  case 1:
482  case 2:
483  return 0.;
484 
485  // quadratics
486  case 3:
487  return 0.;
488 
489  case 4:
490  return 1./distxy;
491 
492  case 5:
493  return 0.;
494 
495  // cubics
496  case 6:
497  return 0.;
498  case 7:
499  return 2.*dx/distxy;
500 
501  case 8:
502  return 2.*dy/distxy;
503 
504  case 9:
505  return 0.;
506 
507  // quartics
508  case 10:
509  return 0.;
510 
511  case 11:
512  return 3.*dx*dx/distxy;
513 
514  case 12:
515  return 4.*dx*dy/distxy;
516 
517  case 13:
518  return 3.*dy*dy/distxy;
519 
520  case 14:
521  return 0.;
522 
523  default:
524  unsigned int o = 0;
525  for (; i >= (o+1)*(o+2)/2; o++) { }
526  unsigned int i2 = i - (o*(o+1)/2);
527  Real val = (o - i2) * i2;
528  for (unsigned int index=i2+1; index < o; index++)
529  val *= dx;
530  for (unsigned int index=1; index < i2; index++)
531  val *= dy;
532  return val/distxy;
533  }
534  }
535 
536  // d^2()/dy^2
537  case 2:
538  {
539  switch (i)
540  {
541  // constants
542  case 0:
543 
544  // linears
545  case 1:
546  case 2:
547  return 0.;
548 
549  // quadratics
550  case 3:
551  case 4:
552  return 0.;
553 
554  case 5:
555  return 2./dist2y;
556 
557  // cubics
558  case 6:
559  return 0.;
560 
561  case 7:
562  return 0.;
563 
564  case 8:
565  return 2.*dx/dist2y;
566 
567  case 9:
568  return 6.*dy/dist2y;
569 
570  // quartics
571  case 10:
572  case 11:
573  return 0.;
574 
575  case 12:
576  return 2.*dx*dx/dist2y;
577 
578  case 13:
579  return 6.*dx*dy/dist2y;
580 
581  case 14:
582  return 12.*dy*dy/dist2y;
583 
584  default:
585  unsigned int o = 0;
586  for (; i >= (o+1)*(o+2)/2; o++) { }
587  unsigned int i2 = i - (o*(o+1)/2);
588  Real val = i2 * (i2 - 1);
589  for (unsigned int index=i2; index != o; index++)
590  val *= dx;
591  for (unsigned int index=2; index < i2; index++)
592  val *= dy;
593  return val/dist2y;
594  }
595  }
596 
597  default:
598  libmesh_error_msg("Invalid shape function derivative j = " << j);
599  }
600 
601  libmesh_error_msg("We'll never get here!");
602  return 0.;
603 
604 #endif
605 }
606 
607 } // namespace libMesh
double abs(double a)
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
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.
long double max(long double a, double b)
Real distance(const Point &p)
libmesh_assert(j)
virtual unsigned int n_nodes() const =0
PetscErrorCode Vec x
PetscErrorCode Vec Mat libmesh_dbg_var(j)
double pow(double a, int b)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
virtual Point centroid() const
Definition: elem.C:446
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)