fe_hermite_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 #include "libmesh/number_lookups.h"
25 
26 namespace
27 {
28 using namespace libMesh;
29 
30 // Compute the static coefficients for an element
31 void hermite_compute_coefs(const Elem* elem, std::vector<std::vector<Real> > & dxdxi
32 #ifdef DEBUG
33  , std::vector<Real> & dxdeta, std::vector<Real> & dydxi
34 #endif
35  )
36 {
37  const Order mapping_order (elem->default_order());
38  const ElemType mapping_elem_type (elem->type());
39  const int n_mapping_shape_functions =
40  FE<2,LAGRANGE>::n_shape_functions(mapping_elem_type,
41  mapping_order);
42 
43  std::vector<Point> dofpt;
44  dofpt.push_back(Point(-1,-1));
45  dofpt.push_back(Point(1,1));
46 
47  for (int p = 0; p != 2; ++p)
48  {
49  dxdxi[0][p] = 0;
50  dxdxi[1][p] = 0;
51 #ifdef DEBUG
52  dxdeta[p] = 0;
53  dydxi[p] = 0;
54 #endif
55  for (int i = 0; i != n_mapping_shape_functions; ++i)
56  {
58  (mapping_elem_type, mapping_order, i, 0, dofpt[p]);
59  const Real ddeta = FE<2,LAGRANGE>::shape_deriv
60  (mapping_elem_type, mapping_order, i, 1, dofpt[p]);
61 
62  dxdxi[0][p] += elem->point(i)(0) * ddxi;
63  dxdxi[1][p] += elem->point(i)(1) * ddeta;
64 // dxdeta and dydxi should be 0!
65 #ifdef DEBUG
66  dxdeta[p] += elem->point(i)(0) * ddeta;
67  dydxi[p] += elem->point(i)(1) * ddxi;
68 #endif
69  }
70  // No singular elements!
71  libmesh_assert(dxdxi[0][p]);
72  libmesh_assert(dxdxi[1][p]);
73  // No non-rectilinear or non-axis-aligned elements!
74 #ifdef DEBUG
75  libmesh_assert_less (std::abs(dxdeta[p]), 1e-9);
76  libmesh_assert_less (std::abs(dydxi[p]), 1e-9);
77 #endif
78  }
79 }
80 
81 
82 
83 Real hermite_bases_2D
84  (std::vector<unsigned int> &bases1D,
85  const std::vector<std::vector<Real> > &dxdxi,
86  const Order &o,
87  unsigned int i)
88 {
89  bases1D.clear();
90  bases1D.resize(2,0);
91  Real coef = 1.0;
92 
93  unsigned int e = o-3;
94 
95  // Nodes
96  if (i < 16)
97  {
98  switch (i / 4)
99  {
100  case 0:
101  break;
102  case 1:
103  bases1D[0] = 1;
104  break;
105  case 2:
106  bases1D[0] = 1;
107  bases1D[1] = 1;
108  break;
109  case 3:
110  bases1D[1] = 1;
111  break;
112  }
113 
114  unsigned int basisnum = i%4;
115  switch (basisnum)
116  {
117  case 0: // DoF = value at node
118  coef = 1.0;
119  break;
120  case 1: // DoF = x derivative at node
121  coef = dxdxi[0][bases1D[0]];
122  bases1D[0] += 2; break;
123  case 2: // DoF = y derivative at node
124  coef = dxdxi[1][bases1D[1]];
125  bases1D[1] += 2; break;
126  case 3: // DoF = xy derivative at node
127  coef = dxdxi[0][bases1D[0]] * dxdxi[1][bases1D[1]];
128  bases1D[0] += 2; bases1D[1] += 2; break;
129  default:
130  libmesh_error();
131  }
132  }
133  // Edges
134  else if (i < 16 + 4*2*e)
135  {
136  unsigned int basisnum = (i - 16) % (2*e);
137  switch ((i - 16) / (2*e))
138  {
139  case 0:
140  bases1D[0] = basisnum/2 + 4;
141  bases1D[1] = basisnum%2*2;
142  if (basisnum%2)
143  coef = dxdxi[1][0];
144  break;
145  case 1:
146  bases1D[0] = basisnum%2*2 + 1;
147  bases1D[1] = basisnum/2 + 4;
148  if (basisnum%2)
149  coef = dxdxi[0][1];
150  break;
151  case 2:
152  bases1D[0] = basisnum/2 + 4;
153  bases1D[1] = basisnum%2*2 + 1;
154  if (basisnum%2)
155  coef = dxdxi[1][1];
156  break;
157  case 3:
158  bases1D[0] = basisnum%2*2;
159  bases1D[1] = basisnum/2 + 4;
160  if (basisnum%2)
161  coef = dxdxi[0][0];
162  break;
163  default:
164  libmesh_error();
165  }
166  }
167  // Interior
168  else
169  {
170  unsigned int basisnum = i - 16 - 8*e;
171  bases1D[0] = square_number_row[basisnum]+4;
172  bases1D[1] = square_number_column[basisnum]+4;
173  }
174 
175  // No singular elements
176  libmesh_assert(coef);
177  return coef;
178 }
179 
180 } // end anonymous namespace
181 
182 
183 namespace libMesh
184 {
185 
186 
187 template <>
189  const Order,
190  const unsigned int,
191  const Point&)
192 {
193  libMesh::err << "Hermite elements require the real element\n"
194  << "to construct gradient-based degrees of freedom."
195  << std::endl;
196 
197  libmesh_error();
198  return 0.;
199 }
200 
201 
202 
203 template <>
205  const Order order,
206  const unsigned int i,
207  const Point& p)
208 {
209  libmesh_assert(elem);
210 
211  std::vector<std::vector<Real> > dxdxi(2, std::vector<Real>(2, 0));
212 
213 #ifdef DEBUG
214  std::vector<Real> dxdeta(2), dydxi(2);
215 #endif
216 
217  hermite_compute_coefs(elem,dxdxi
218 #ifdef DEBUG
219  ,dxdeta,dydxi
220 #endif
221  );
222 
223  const ElemType type = elem->type();
224 
225  const Order totalorder = static_cast<Order>(order + elem->p_level());
226 
227  switch (type)
228  {
229  case QUAD4:
230  libmesh_assert_less (totalorder, 4);
231  case QUAD8:
232  case QUAD9:
233  {
234  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
235 
236  std::vector<unsigned int> bases1D;
237 
238  Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
239 
240  return coef * FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
241  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
242  }
243  default:
244  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
245  libmesh_error();
246  }
247 
248  libmesh_error();
249  return 0.;
250 }
251 
252 
253 
254 template <>
256  const Order,
257  const unsigned int,
258  const unsigned int,
259  const Point&)
260 {
261  libMesh::err << "Hermite elements require the real element\n"
262  << "to construct gradient-based degrees of freedom."
263  << std::endl;
264 
265  libmesh_error();
266  return 0.;
267 }
268 
269 
270 
271 template <>
273  const Order order,
274  const unsigned int i,
275  const unsigned int j,
276  const Point& p)
277 {
278  libmesh_assert(elem);
279  libmesh_assert (j == 0 || j == 1);
280 
281  std::vector<std::vector<Real> > dxdxi(2, std::vector<Real>(2, 0));
282 
283 #ifdef DEBUG
284  std::vector<Real> dxdeta(2), dydxi(2);
285 #endif
286 
287  hermite_compute_coefs(elem,dxdxi
288 #ifdef DEBUG
289  ,dxdeta,dydxi
290 #endif
291  );
292 
293  const ElemType type = elem->type();
294 
295  const Order totalorder = static_cast<Order>(order + elem->p_level());
296 
297  switch (type)
298  {
299  case QUAD4:
300  libmesh_assert_less (totalorder, 4);
301  case QUAD8:
302  case QUAD9:
303  {
304  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
305 
306  std::vector<unsigned int> bases1D;
307 
308  Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
309 
310  switch (j)
311  {
312  case 0:
313  return coef *
314  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
315  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
316  case 1:
317  return coef *
318  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
319  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1));
320  default:
321  libmesh_error();
322  }
323  }
324  default:
325  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
326  libmesh_error();
327  }
328 
329  libmesh_error();
330  return 0.;
331 }
332 
333 
334 
335 template <>
337  const Order order,
338  const unsigned int i,
339  const unsigned int j,
340  const Point& p)
341 {
342  libmesh_assert(elem);
343  libmesh_assert (j == 0 || j == 1 || j == 2);
344 
345  std::vector<std::vector<Real> > dxdxi(2, std::vector<Real>(2, 0));
346 
347 #ifdef DEBUG
348  std::vector<Real> dxdeta(2), dydxi(2);
349 #endif
350 
351  hermite_compute_coefs(elem,dxdxi
352 #ifdef DEBUG
353  ,dxdeta,dydxi
354 #endif
355  );
356 
357  const ElemType type = elem->type();
358 
359  const Order totalorder = static_cast<Order>(order + elem->p_level());
360 
361  switch (type)
362  {
363  case QUAD4:
364  libmesh_assert_less (totalorder, 4);
365  case QUAD8:
366  case QUAD9:
367  {
368  libmesh_assert_less (i, (totalorder+1u)*(totalorder+1u));
369 
370  std::vector<unsigned int> bases1D;
371 
372  Real coef = hermite_bases_2D(bases1D, dxdxi, totalorder, i);
373 
374  switch (j)
375  {
376  case 0:
377  return coef *
379  FEHermite<1>::hermite_raw_shape(bases1D[1],p(1));
380  case 1:
381  return coef *
382  FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) *
383  FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1));
384  case 2:
385  return coef *
386  FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
388  default:
389  libmesh_error();
390  }
391  }
392  default:
393  libMesh::err << "ERROR: Unsupported element type!" << std::endl;
394  libmesh_error();
395  }
396 
397  libmesh_error();
398  return 0.;
399 }
400 
401 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo