fe_xyz_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 
26 
27 
28 // Anonymous namespace for persistant variables.
29 // This allows us to determine when the centroid needs
30 // to be recalculated.
31 namespace
32 {
33  using namespace libMesh;
34 
35  static dof_id_type old_elem_id = DofObject::invalid_id;
36  static Point centroid;
37  static Point max_distance;
38 }
39 
40 
41 namespace libMesh
42 {
43 
44 
45 template <>
47  const Order,
48  const unsigned int,
49  const Point&)
50 {
51  libMesh::err << "XYZ polynomials require the element\n"
52  << "because the centroid is needed."
53  << std::endl;
54 
55  libmesh_error();
56  return 0.;
57 }
58 
59 
60 
61 template <>
63  const Order libmesh_dbg_var(order),
64  const unsigned int i,
65  const Point& point_in)
66 {
67 #if LIBMESH_DIM > 1
68 
69  libmesh_assert(elem);
70 
71  // Only recompute the centroid if the element
72  // has changed from the last one we computed.
73  // This avoids repeated centroid calculations
74  // when called in succession with the same element.
75  if (elem->id() != old_elem_id)
76  {
77  centroid = elem->centroid();
78  old_elem_id = elem->id();
79  max_distance = Point(0.,0.,0.);
80  for (unsigned int p = 0; p < elem->n_nodes(); p++)
81  for (unsigned int d = 0; d < 2; d++)
82  {
83  const Real distance = std::abs(centroid(d) - elem->point(p)(d));
84  max_distance(d) = std::max(distance, max_distance(d));
85  }
86  }
87 
88  // Using static globals for old_elem_id, etc. will fail
89  // horribly with more than one thread.
90  libmesh_assert_equal_to (libMesh::n_threads(), 1);
91 
92  const Real x = point_in(0);
93  const Real y = point_in(1);
94  const Real xc = centroid(0);
95  const Real yc = centroid(1);
96  const Real distx = max_distance(0);
97  const Real disty = max_distance(1);
98  const Real dx = (x - xc)/distx;
99  const Real dy = (y - yc)/disty;
100 
101 #ifndef NDEBUG
102  // totalorder is only used in the assertion below, so
103  // we avoid declaring it when asserts are not active.
104  const unsigned int totalorder = order + elem->p_level();
105 #endif
106  libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
107 
108 
109  // monomials. since they are hierarchic we only need one case block.
110  switch (i)
111  {
112  // constant
113  case 0:
114  return 1.;
115 
116  // linear
117  case 1:
118  return dx;
119 
120  case 2:
121  return dy;
122 
123  // quadratics
124  case 3:
125  return dx*dx;
126 
127  case 4:
128  return dx*dy;
129 
130  case 5:
131  return dy*dy;
132 
133  // cubics
134  case 6:
135  return dx*dx*dx;
136 
137  case 7:
138  return dx*dx*dy;
139 
140  case 8:
141  return dx*dy*dy;
142 
143  case 9:
144  return dy*dy*dy;
145 
146  // quartics
147  case 10:
148  return dx*dx*dx*dx;
149 
150  case 11:
151  return dx*dx*dx*dy;
152 
153  case 12:
154  return dx*dx*dy*dy;
155 
156  case 13:
157  return dx*dy*dy*dy;
158 
159  case 14:
160  return dy*dy*dy*dy;
161 
162  default:
163  unsigned int o = 0;
164  for (; i >= (o+1)*(o+2)/2; o++) { }
165  unsigned int i2 = i - (o*(o+1)/2);
166  Real val = 1.;
167  for (unsigned int index=i2; index != o; index++)
168  val *= dx;
169  for (unsigned int index=0; index != i2; index++)
170  val *= dy;
171  return val;
172  }
173 
174  libmesh_error();
175  return 0.;
176 
177 #endif
178 }
179 
180 
181 
182 template <>
184  const Order,
185  const unsigned int,
186  const unsigned int,
187  const Point&)
188 {
189  libMesh::err << "XYZ polynomials require the element\n"
190  << "because the centroid is needed."
191  << std::endl;
192 
193  libmesh_error();
194  return 0.;
195 }
196 
197 
198 
199 template <>
201  const Order libmesh_dbg_var(order),
202  const unsigned int i,
203  const unsigned int j,
204  const Point& point_in)
205 {
206 #if LIBMESH_DIM > 1
207 
208 
209  libmesh_assert_less (j, 2);
210  libmesh_assert(elem);
211 
212  // Only recompute the centroid if the element
213  // has changed from the last one we computed.
214  // This avoids repeated centroid calculations
215  // when called in succession with the same element.
216  if (elem->id() != old_elem_id)
217  {
218  centroid = elem->centroid();
219  old_elem_id = elem->id();
220  max_distance = Point(0.,0.,0.);
221  for (unsigned int p = 0; p < elem->n_nodes(); p++)
222  for (unsigned int d = 0; d < 2; d++)
223  {
224  const Real distance = std::abs(centroid(d) - elem->point(p)(d));
225  max_distance(d) = std::max(distance, max_distance(d));
226  }
227  }
228 
229  // Using static globals for old_elem_id, etc. will fail
230  // horribly with more than one thread.
231  libmesh_assert_equal_to (libMesh::n_threads(), 1);
232 
233  const Real x = point_in(0);
234  const Real y = point_in(1);
235  const Real xc = centroid(0);
236  const Real yc = centroid(1);
237  const Real distx = max_distance(0);
238  const Real disty = max_distance(1);
239  const Real dx = (x - xc)/distx;
240  const Real dy = (y - yc)/disty;
241 
242 #ifndef NDEBUG
243  // totalorder is only used in the assertion below, so
244  // we avoid declaring it when asserts are not active.
245  const unsigned int totalorder = order + elem->p_level();
246 #endif
247  libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
248 
249  // monomials. since they are hierarchic we only need one case block.
250 
251 switch (j)
252  {
253  // d()/dx
254  case 0:
255  {
256  switch (i)
257  {
258  // constants
259  case 0:
260  return 0.;
261 
262  // linears
263  case 1:
264  return 1./distx;
265 
266  case 2:
267  return 0.;
268 
269  // quadratics
270  case 3:
271  return 2.*dx/distx;
272 
273  case 4:
274  return dy/distx;
275 
276  case 5:
277  return 0.;
278 
279  // cubics
280  case 6:
281  return 3.*dx*dx/distx;
282 
283  case 7:
284  return 2.*dx*dy/distx;
285 
286  case 8:
287  return dy*dy/distx;
288 
289  case 9:
290  return 0.;
291 
292  // quartics
293  case 10:
294  return 4.*dx*dx*dx/distx;
295 
296  case 11:
297  return 3.*dx*dx*dy/distx;
298 
299  case 12:
300  return 2.*dx*dy*dy/distx;
301 
302  case 13:
303  return dy*dy*dy/distx;
304 
305  case 14:
306  return 0.;
307 
308  default:
309  unsigned int o = 0;
310  for (; i >= (o+1)*(o+2)/2; o++) { }
311  unsigned int i2 = i - (o*(o+1)/2);
312  Real val = o - i2;
313  for (unsigned int index=i2+1; index < o; index++)
314  val *= dx;
315  for (unsigned int index=0; index != i2; index++)
316  val *= dy;
317  return val/distx;
318  }
319  }
320 
321 
322  // d()/dy
323  case 1:
324  {
325  switch (i)
326  {
327  // constants
328  case 0:
329  return 0.;
330 
331  // linears
332  case 1:
333  return 0.;
334 
335  case 2:
336  return 1./disty;
337 
338  // quadratics
339  case 3:
340  return 0.;
341 
342  case 4:
343  return dx/disty;
344 
345  case 5:
346  return 2.*dy/disty;
347 
348  // cubics
349  case 6:
350  return 0.;
351 
352  case 7:
353  return dx*dx/disty;
354 
355  case 8:
356  return 2.*dx*dy/disty;
357 
358  case 9:
359  return 3.*dy*dy/disty;
360 
361  // quartics
362  case 10:
363  return 0.;
364 
365  case 11:
366  return dx*dx*dx/disty;
367 
368  case 12:
369  return 2.*dx*dx*dy/disty;
370 
371  case 13:
372  return 3.*dx*dy*dy/disty;
373 
374  case 14:
375  return 4.*dy*dy*dy/disty;
376 
377  default:
378  unsigned int o = 0;
379  for (; i >= (o+1)*(o+2)/2; o++) { }
380  unsigned int i2 = i - (o*(o+1)/2);
381  Real val = i2;
382  for (unsigned int index=i2; index != o; index++)
383  val *= dx;
384  for (unsigned int index=1; index <= i2; index++)
385  val *= dy;
386  return val/disty;
387  }
388  }
389 
390 
391  default:
392  libmesh_error();
393  }
394 
395  libmesh_error();
396  return 0.;
397 
398 #endif
399 }
400 
401 
402 
403 template <>
405  const Order,
406  const unsigned int,
407  const unsigned int,
408  const Point&)
409 {
410  libMesh::err << "XYZ polynomials require the element\n"
411  << "because the centroid is needed."
412  << std::endl;
413 
414  libmesh_error();
415  return 0.;
416 }
417 
418 
419 
420 template <>
422  const Order libmesh_dbg_var(order),
423  const unsigned int i,
424  const unsigned int j,
425  const Point& point_in)
426 {
427 #if LIBMESH_DIM > 1
428 
429  libmesh_assert_less_equal (j, 2);
430  libmesh_assert(elem);
431 
432  // Only recompute the centroid if the element
433  // has changed from the last one we computed.
434  // This avoids repeated centroid calculations
435  // when called in succession with the same element.
436  if (elem->id() != old_elem_id)
437  {
438  centroid = elem->centroid();
439  old_elem_id = elem->id();
440  max_distance = Point(0.,0.,0.);
441  for (unsigned int p = 0; p < elem->n_nodes(); p++)
442  for (unsigned int d = 0; d < 2; d++)
443  {
444  const Real distance = std::abs(centroid(d) - elem->point(p)(d));
445  max_distance(d) = std::max(distance, max_distance(d));
446  }
447  }
448 
449  // Using static globals for old_elem_id, etc. will fail
450  // horribly with more than one thread.
451  libmesh_assert_equal_to (libMesh::n_threads(), 1);
452 
453  const Real x = point_in(0);
454  const Real y = point_in(1);
455  const Real xc = centroid(0);
456  const Real yc = centroid(1);
457  const Real distx = max_distance(0);
458  const Real disty = max_distance(1);
459  const Real dx = (x - xc)/distx;
460  const Real dy = (y - yc)/disty;
461  const Real dist2x = pow(distx,2.);
462  const Real dist2y = pow(disty,2.);
463  const Real distxy = distx * disty;
464 
465 #ifndef NDEBUG
466  // totalorder is only used in the assertion below, so
467  // we avoid declaring it when asserts are not active.
468  const unsigned int totalorder = order + elem->p_level();
469 #endif
470  libmesh_assert_less (i, (totalorder+1)*(totalorder+2)/2);
471 
472  // monomials. since they are hierarchic we only need one case block.
473 
474  switch (j)
475  {
476  // d^2()/dx^2
477  case 0:
478  {
479  switch (i)
480  {
481  // constants
482  case 0:
483  // linears
484  case 1:
485  case 2:
486  return 0.;
487 
488  // quadratics
489  case 3:
490  return 2./dist2x;
491 
492  case 4:
493  case 5:
494  return 0.;
495 
496  // cubics
497  case 6:
498  return 6.*dx/dist2x;
499 
500  case 7:
501  return 2.*dy/dist2x;
502 
503  case 8:
504  case 9:
505  return 0.;
506 
507  // quartics
508  case 10:
509  return 12.*dx*dx/dist2x;
510 
511  case 11:
512  return 6.*dx*dy/dist2x;
513 
514  case 12:
515  return 2.*dy*dy/dist2x;
516 
517  case 13:
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) * (o - i2 - 1);
526  for (unsigned int index=i2+2; index < o; index++)
527  val *= dx;
528  for (unsigned int index=0; index != i2; index++)
529  val *= dy;
530  return val/dist2x;
531  }
532  }
533 
534  // d^2()/dxdy
535  case 1:
536  {
537  switch (i)
538  {
539  // constants
540  case 0:
541 
542  // linears
543  case 1:
544  case 2:
545  return 0.;
546 
547  // quadratics
548  case 3:
549  return 0.;
550 
551  case 4:
552  return 1./distxy;
553 
554  case 5:
555  return 0.;
556 
557  // cubics
558  case 6:
559  return 0.;
560  case 7:
561  return 2.*dx/distxy;
562 
563  case 8:
564  return 2.*dy/distxy;
565 
566  case 9:
567  return 0.;
568 
569  // quartics
570  case 10:
571  return 0.;
572 
573  case 11:
574  return 3.*dx*dx/distxy;
575 
576  case 12:
577  return 4.*dx*dy/distxy;
578 
579  case 13:
580  return 3.*dy*dy/distxy;
581 
582  case 14:
583  return 0.;
584 
585  default:
586  unsigned int o = 0;
587  for (; i >= (o+1)*(o+2)/2; o++) { }
588  unsigned int i2 = i - (o*(o+1)/2);
589  Real val = (o - i2) * i2;
590  for (unsigned int index=i2+1; index < o; index++)
591  val *= dx;
592  for (unsigned int index=1; index < i2; index++)
593  val *= dy;
594  return val/distxy;
595  }
596  }
597 
598  // d^2()/dy^2
599  case 2:
600  {
601  switch (i)
602  {
603  // constants
604  case 0:
605 
606  // linears
607  case 1:
608  case 2:
609  return 0.;
610 
611  // quadratics
612  case 3:
613  case 4:
614  return 0.;
615 
616  case 5:
617  return 2./dist2y;
618 
619  // cubics
620  case 6:
621  return 0.;
622 
623  case 7:
624  return 0.;
625 
626  case 8:
627  return 2.*dx/dist2y;
628 
629  case 9:
630  return 6.*dy/dist2y;
631 
632  // quartics
633  case 10:
634  case 11:
635  return 0.;
636 
637  case 12:
638  return 2.*dx*dx/dist2y;
639 
640  case 13:
641  return 6.*dx*dy/dist2y;
642 
643  case 14:
644  return 12.*dy*dy/dist2y;
645 
646  default:
647  unsigned int o = 0;
648  for (; i >= (o+1)*(o+2)/2; o++) { }
649  unsigned int i2 = i - (o*(o+1)/2);
650  Real val = i2 * (i2 - 1);
651  for (unsigned int index=i2; index != o; index++)
652  val *= dx;
653  for (unsigned int index=2; index < i2; index++)
654  val *= dy;
655  return val/dist2y;
656  }
657  }
658  }
659 
660  libmesh_error();
661  return 0.;
662 
663 #endif
664 }
665 
666 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo