fe_hierarchic_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++ includes
20 
21 // Local includes
22 #include "libmesh/fe.h"
23 #include "libmesh/elem.h"
24 #include "libmesh/number_lookups.h"
25 #include "libmesh/utility.h"
26 
27 
28 namespace libMesh
29 {
30 
31 template <>
33  const Order,
34  const unsigned int,
35  const Point&)
36 {
37  libMesh::err << "Hierarchic polynomials require the element type\n"
38  << "because edge orientation is needed."
39  << std::endl;
40 
41  libmesh_error();
42  return 0.;
43 }
44 
45 
46 
47 template <>
49  const Order order,
50  const unsigned int i,
51  const Point& p)
52 {
53  libmesh_assert(elem);
54 
55  const Order totalorder = static_cast<Order>(order+elem->p_level());
56  libmesh_assert_greater (totalorder, 0);
57 
58  switch (elem->type())
59  {
60  case TRI3:
61  case TRI6:
62  {
63  const Real zeta1 = p(0);
64  const Real zeta2 = p(1);
65  const Real zeta0 = 1. - zeta1 - zeta2;
66 
67  libmesh_assert_less (i, (totalorder+1u)*(totalorder+2u)/2);
68  libmesh_assert (elem->type() == TRI6 || totalorder < 2);
69 
70  // Vertex DoFs
71  if (i == 0)
72  return zeta0;
73  else if (i == 1)
74  return zeta1;
75  else if (i == 2)
76  return zeta2;
77  // Edge DoFs
78  else if (i < totalorder + 2u)
79  {
80  // Avoid returning NaN on vertices!
81  if (zeta0 + zeta1 == 0.)
82  return 0.;
83 
84  const unsigned int basisorder = i - 1;
85  // Get factors to account for edge-flipping
86  Real f0 = 1;
87  if (basisorder%2 && (elem->point(0) > elem->point(1)))
88  f0 = -1.;
89 
90  Real edgeval = (zeta1 - zeta0) / (zeta1 + zeta0);
91  Real crossfunc = zeta0 + zeta1;
92  for (unsigned int n=1; n != basisorder; ++n)
93  crossfunc *= (zeta0 + zeta1);
94 
95  return f0 * crossfunc *
96  FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
97  basisorder, edgeval);
98  }
99  else if (i < 2u*totalorder + 1)
100  {
101  // Avoid returning NaN on vertices!
102  if (zeta1 + zeta2 == 0.)
103  return 0.;
104 
105  const unsigned int basisorder = i - totalorder;
106  // Get factors to account for edge-flipping
107  Real f1 = 1;
108  if (basisorder%2 && (elem->point(1) > elem->point(2)))
109  f1 = -1.;
110 
111  Real edgeval = (zeta2 - zeta1) / (zeta2 + zeta1);
112  Real crossfunc = zeta2 + zeta1;
113  for (unsigned int n=1; n != basisorder; ++n)
114  crossfunc *= (zeta2 + zeta1);
115 
116  return f1 * crossfunc *
117  FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
118  basisorder, edgeval);
119  }
120  else if (i < 3u*totalorder)
121  {
122  // Avoid returning NaN on vertices!
123  if (zeta0 + zeta2 == 0.)
124  return 0.;
125 
126  const unsigned int basisorder = i - (2u*totalorder) + 1;
127  // Get factors to account for edge-flipping
128  Real f2 = 1;
129  if (basisorder%2 && (elem->point(2) > elem->point(0)))
130  f2 = -1.;
131 
132  Real edgeval = (zeta0 - zeta2) / (zeta0 + zeta2);
133  Real crossfunc = zeta0 + zeta2;
134  for (unsigned int n=1; n != basisorder; ++n)
135  crossfunc *= (zeta0 + zeta2);
136 
137  return f2 * crossfunc *
138  FE<1,HIERARCHIC>::shape(EDGE3, totalorder,
139  basisorder, edgeval);
140  }
141  // Interior DoFs
142  else
143  {
144  const unsigned int basisnum = i - (3u*totalorder);
145  unsigned int exp0 = triangular_number_column[basisnum] + 1;
146  unsigned int exp1 = triangular_number_row[basisnum] + 1 -
147  triangular_number_column[basisnum];
148 
149  Real returnval = 1;
150  for (unsigned int n = 0; n != exp0; ++n)
151  returnval *= zeta0;
152  for (unsigned int n = 0; n != exp1; ++n)
153  returnval *= zeta1;
154  returnval *= zeta2;
155  return returnval;
156  }
157  }
158 
159  // Hierarchic shape functions on the quadrilateral.
160  case QUAD4:
161  libmesh_assert_less (totalorder, 2);
162  case QUAD8:
163  case QUAD9:
164  {
165  // Compute quad shape functions as a tensor-product
166  const Real xi = p(0);
167  const Real eta = p(1);
168 
169  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
170 
171 // Example i, i0, i1 values for totalorder = 5:
172 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
173 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 3, 2, 4, 4, 4, 3, 2, 5, 5, 5, 5, 4, 3, 2};
174 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 3, 3, 2, 3, 4, 4, 4, 2, 3, 4, 5, 5, 5, 5};
175 
176  unsigned int i0, i1;
177 
178  // Vertex DoFs
179  if (i == 0)
180  { i0 = 0; i1 = 0; }
181  else if (i == 1)
182  { i0 = 1; i1 = 0; }
183  else if (i == 2)
184  { i0 = 1; i1 = 1; }
185  else if (i == 3)
186  { i0 = 0; i1 = 1; }
187  // Edge DoFs
188  else if (i < totalorder + 3u)
189  { i0 = i - 2; i1 = 0; }
190  else if (i < 2u*totalorder + 2)
191  { i0 = 1; i1 = i - totalorder - 1; }
192  else if (i < 3u*totalorder + 1)
193  { i0 = i - 2u*totalorder; i1 = 1; }
194  else if (i < 4u*totalorder)
195  { i0 = 0; i1 = i - 3u*totalorder + 1; }
196  // Interior DoFs
197  else
198  {
199  unsigned int basisnum = i - 4*totalorder;
200  i0 = square_number_column[basisnum] + 2;
201  i1 = square_number_row[basisnum] + 2;
202  }
203 
204  // Flip odd degree of freedom values if necessary
205  // to keep continuity on sides
206  Real f = 1.;
207 
208  if ((i0%2) && (i0 > 2) && (i1 == 0))
209  f = (elem->point(0) > elem->point(1))?-1.:1.;
210  else if ((i0%2) && (i0>2) && (i1 == 1))
211  f = (elem->point(3) > elem->point(2))?-1.:1.;
212  else if ((i0 == 0) && (i1%2) && (i1>2))
213  f = (elem->point(0) > elem->point(3))?-1.:1.;
214  else if ((i0 == 1) && (i1%2) && (i1>2))
215  f = (elem->point(1) > elem->point(2))?-1.:1.;
216 
217  return f*(FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i0, xi)*
218  FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i1, eta));
219  }
220 
221  default:
222  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
223  libmesh_error();
224  }
225 
226  return 0.;
227 }
228 
229 
230 
231 template <>
233  const Order,
234  const unsigned int,
235  const unsigned int,
236  const Point&)
237 {
238  libMesh::err << "Hierarchic polynomials require the element type\n"
239  << "because edge orientation is needed."
240  << std::endl;
241 
242  libmesh_error();
243  return 0.;
244 }
245 
246 
247 
248 template <>
250  const Order order,
251  const unsigned int i,
252  const unsigned int j,
253  const Point& p)
254 {
255  libmesh_assert(elem);
256 
257  const ElemType type = elem->type();
258 
259  const Order totalorder = static_cast<Order>(order+elem->p_level());
260 
261  libmesh_assert_greater (totalorder, 0);
262 
263  switch (type)
264  {
265  // 1st & 2nd-order Hierarchics.
266  case TRI3:
267  case TRI6:
268  {
269  const Real eps = 1.e-6;
270 
271  libmesh_assert_less (j, 2);
272 
273  switch (j)
274  {
275  // d()/dxi
276  case 0:
277  {
278  const Point pp(p(0)+eps, p(1));
279  const Point pm(p(0)-eps, p(1));
280 
281  return (FE<2,HIERARCHIC>::shape(elem, order, i, pp) -
282  FE<2,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
283  }
284 
285  // d()/deta
286  case 1:
287  {
288  const Point pp(p(0), p(1)+eps);
289  const Point pm(p(0), p(1)-eps);
290 
291  return (FE<2,HIERARCHIC>::shape(elem, order, i, pp) -
292  FE<2,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
293  }
294 
295 
296  default:
297  libmesh_error();
298  }
299  }
300 
301  case QUAD4:
302  libmesh_assert_less (totalorder, 2);
303  case QUAD8:
304  case QUAD9:
305  {
306  // Compute quad shape functions as a tensor-product
307  const Real xi = p(0);
308  const Real eta = p(1);
309 
310  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
311 
312 // Example i, i0, i1 values for totalorder = 5:
313 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
314 // static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
315 // static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
316 
317  unsigned int i0, i1;
318 
319  // Vertex DoFs
320  if (i == 0)
321  { i0 = 0; i1 = 0; }
322  else if (i == 1)
323  { i0 = 1; i1 = 0; }
324  else if (i == 2)
325  { i0 = 1; i1 = 1; }
326  else if (i == 3)
327  { i0 = 0; i1 = 1; }
328  // Edge DoFs
329  else if (i < totalorder + 3u)
330  { i0 = i - 2; i1 = 0; }
331  else if (i < 2u*totalorder + 2)
332  { i0 = 1; i1 = i - totalorder - 1; }
333  else if (i < 3u*totalorder + 1u)
334  { i0 = i - 2u*totalorder; i1 = 1; }
335  else if (i < 4u*totalorder)
336  { i0 = 0; i1 = i - 3u*totalorder + 1; }
337  // Interior DoFs
338  else
339  {
340  unsigned int basisnum = i - 4*totalorder;
341  i0 = square_number_column[basisnum] + 2;
342  i1 = square_number_row[basisnum] + 2;
343  }
344 
345  // Flip odd degree of freedom values if necessary
346  // to keep continuity on sides
347  Real f = 1.;
348 
349  if ((i0%2) && (i0 > 2) && (i1 == 0))
350  f = (elem->point(0) > elem->point(1))?-1.:1.;
351  else if ((i0%2) && (i0>2) && (i1 == 1))
352  f = (elem->point(3) > elem->point(2))?-1.:1.;
353  else if ((i0 == 0) && (i1%2) && (i1>2))
354  f = (elem->point(0) > elem->point(3))?-1.:1.;
355  else if ((i0 == 1) && (i1%2) && (i1>2))
356  f = (elem->point(1) > elem->point(2))?-1.:1.;
357 
358  switch (j)
359  {
360  // d()/dxi
361  case 0:
362  return f*(FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
363  FE<1,HIERARCHIC>::shape (EDGE3, totalorder, i1, eta));
364 
365  // d()/deta
366  case 1:
367  return f*(FE<1,HIERARCHIC>::shape (EDGE3, totalorder, i0, xi)*
368  FE<1,HIERARCHIC>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
369 
370  default:
371  libmesh_error();
372  }
373 
374  }
375 
376  default:
377  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
378  libmesh_error();
379  }
380 
381  return 0.;
382 }
383 
384 
385 
386 template <>
388  const Order,
389  const unsigned int,
390  const unsigned int,
391  const Point&)
392 {
393  libMesh::err << "Hierarchic polynomials require the element type\n"
394  << "because edge orientation is needed."
395  << std::endl;
396 
397  libmesh_error();
398  return 0.;
399 }
400 
401 
402 
403 template <>
405  const Order order,
406  const unsigned int i,
407  const unsigned int j,
408  const Point& p)
409 {
410  libmesh_assert(elem);
411 
412  // I have been lazy here and am using finite differences
413  // to compute the derivatives!
414  const Real eps = 1.e-6;
415  Point pp, pm;
416  unsigned int prevj = libMesh::invalid_uint;
417 
418  switch (j)
419  {
420  // d^2()/dxi^2
421  case 0:
422  {
423  pp = Point(p(0)+eps, p(1));
424  pm = Point(p(0)-eps, p(1));
425  prevj = 0;
426  break;
427  }
428 
429  // d^2()/dxideta
430  case 1:
431  {
432  pp = Point(p(0), p(1)+eps);
433  pm = Point(p(0), p(1)-eps);
434  prevj = 0;
435  break;
436  }
437 
438  // d^2()/deta^2
439  case 2:
440  {
441  pp = Point(p(0), p(1)+eps);
442  pm = Point(p(0), p(1)-eps);
443  prevj = 1;
444  break;
445  }
446  default:
447  libmesh_error();
448  }
449  return (FE<2,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp) -
450  FE<2,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pm)
451  )/2./eps;
452 }
453 
454 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo