fe_hermite.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 
20 // Local includes
21 #include "libmesh/elem.h"
22 #include "libmesh/fe.h"
23 #include "libmesh/fe_interface.h"
24 
25 namespace libMesh
26 {
27 
28  // ------------------------------------------------------------
29  // Hermite-specific implementations
30 
31  // Anonymous namespace for local helper functions
32  namespace {
33 
34  void hermite_nodal_soln(const Elem* elem,
35  const Order order,
36  const std::vector<Number>& elem_soln,
37  std::vector<Number>& nodal_soln,
38  unsigned Dim)
39  {
40  const unsigned int n_nodes = elem->n_nodes();
41 
42  const ElemType elem_type = elem->type();
43 
44  nodal_soln.resize(n_nodes);
45 
46  const Order totalorder = static_cast<Order>(order + elem->p_level());
47 
48  // FEType object to be passed to various FEInterface functions below.
49  FEType fe_type(totalorder, HERMITE);
50 
51  const unsigned int n_sf =
52  // FE<Dim,T>::n_shape_functions(elem_type, totalorder);
53  FEInterface::n_shape_functions(Dim, fe_type, elem_type);
54 
55  std::vector<Point> refspace_nodes;
56  FEBase::get_refspace_nodes(elem_type,refspace_nodes);
57  libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
58 
59  for (unsigned int n=0; n<n_nodes; n++)
60  {
61  libmesh_assert_equal_to (elem_soln.size(), n_sf);
62 
63  // Zero before summation
64  nodal_soln[n] = 0;
65 
66  // u_i = Sum (alpha_i phi_i)
67  for (unsigned int i=0; i<n_sf; i++)
68  nodal_soln[n] += elem_soln[i] *
69  // FE<Dim,T>::shape(elem, order, i, mapped_point);
70  FEInterface::shape(Dim, fe_type, elem, i, refspace_nodes[n]);
71  }
72  } // hermite_nodal_soln()
73 
74 
75 
76  unsigned int hermite_n_dofs(const ElemType t, const Order o)
77  {
79  // Piecewise (bi/tri)cubic C1 Hermite splines
80  switch (t)
81  {
82  case NODEELEM:
83  return 1;
84  case EDGE2:
85  libmesh_assert_less (o, 4);
86  case EDGE3:
87  return (o+1);
88 
89  case QUAD4:
90  case QUAD8:
91  libmesh_assert_less (o, 4);
92  case QUAD9:
93  return ((o+1)*(o+1));
94 
95  case HEX8:
96  case HEX20:
97  libmesh_assert_less (o, 4);
98  case HEX27:
99  return ((o+1)*(o+1)*(o+1));
100 
101  default:
102  {
103 #ifdef DEBUG
104  libMesh::err << "ERROR: Bad ElemType = " << t
105  << " for " << o << "th order approximation!"
106  << std::endl;
107 #endif
108  libmesh_error();
109  }
110  }
111 
112  libmesh_error();
113  return 0;
114  } // hermite_n_dofs()
115 
116 
117 
118 
119  unsigned int hermite_n_dofs_at_node(const ElemType t,
120  const Order o,
121  const unsigned int n)
122  {
123  libmesh_assert_greater (o, 2);
124  // Piecewise (bi/tri)cubic C1 Hermite splines
125  switch (t)
126  {
127  case NODEELEM:
128  return 1;
129  case EDGE2:
130  case EDGE3:
131  {
132  switch (n)
133  {
134  case 0:
135  case 1:
136  return 2;
137  case 2:
138  // Interior DoFs are carried on Elems
139  // return (o-3);
140  return 0;
141 
142  default:
143  libmesh_error();
144  }
145  }
146 
147  case QUAD4:
148  libmesh_assert_less (o, 4);
149  case QUAD8:
150  case QUAD9:
151  {
152  switch (n)
153  {
154  // Vertices
155  case 0:
156  case 1:
157  case 2:
158  case 3:
159  return 4;
160  // Edges
161  case 4:
162  case 5:
163  case 6:
164  case 7:
165  return (2*(o-3));
166  case 8:
167  // Interior DoFs are carried on Elems
168  // return ((o-3)*(o-3));
169  return 0;
170 
171  default:
172  libmesh_error();
173  }
174  }
175 
176  case HEX8:
177  case HEX20:
178  libmesh_assert_less (o, 4);
179  case HEX27:
180  {
181  switch (n)
182  {
183  // Vertices
184  case 0:
185  case 1:
186  case 2:
187  case 3:
188  case 4:
189  case 5:
190  case 6:
191  case 7:
192  return 8;
193  // Edges
194  case 8:
195  case 9:
196  case 10:
197  case 11:
198  case 12:
199  case 13:
200  case 14:
201  case 15:
202  case 16:
203  case 17:
204  case 18:
205  case 19:
206  return (4*(o-3));
207  // Faces
208  case 20:
209  case 21:
210  case 22:
211  case 23:
212  case 24:
213  case 25:
214  return (2*(o-3)*(o-3));
215  case 26:
216  // Interior DoFs are carried on Elems
217  // return ((o-3)*(o-3)*(o-3));
218  return 0;
219 
220  default:
221  libmesh_error();
222  }
223  }
224 
225  default:
226  {
227 #ifdef DEBUG
228  libMesh::err << "ERROR: Bad ElemType = " << t
229  << " for " << o << "th order approximation!"
230  << std::endl;
231 #endif
232  libmesh_error();
233  }
234 
235  }
236 
237  libmesh_error();
238 
239  return 0;
240  } // hermite_n_dofs_at_node()
241 
242 
243 
244  unsigned int hermite_n_dofs_per_elem(const ElemType t,
245  const Order o)
246  {
247  libmesh_assert_greater (o, 2);
248 
249  switch (t)
250  {
251  case NODEELEM:
252  return 0;
253  case EDGE2:
254  case EDGE3:
255  return (o-3);
256  case QUAD4:
257  libmesh_assert_less (o, 4);
258  case QUAD8:
259  case QUAD9:
260  return ((o-3)*(o-3));
261  case HEX8:
262  libmesh_assert_less (o, 4);
263  case HEX20:
264  case HEX27:
265  return ((o-3)*(o-3)*(o-3));
266 
267  default:
268  {
269 #ifdef DEBUG
270  libMesh::err << "ERROR: Bad ElemType = " << t
271  << " for " << o << "th order approximation!"
272  << std::endl;
273 #endif
274  libmesh_error();
275  }
276  }
277 
278  // Will never get here...
279  libmesh_error();
280  return 0;
281  } // hermite_n_dofs_per_elem()
282 
283 
284  } // anonymous namespace
285 
286 
287 
288 
289  // Do full-specialization of nodal_soln() function for every
290  // dimension, instead of explicit instantiation at the end of this
291  // file.
292  // This could be macro-ified so that it fits on one line...
293  template <>
294  void FE<0,HERMITE>::nodal_soln(const Elem* elem,
295  const Order order,
296  const std::vector<Number>& elem_soln,
297  std::vector<Number>& nodal_soln)
298  { hermite_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/0); }
299 
300  template <>
301  void FE<1,HERMITE>::nodal_soln(const Elem* elem,
302  const Order order,
303  const std::vector<Number>& elem_soln,
304  std::vector<Number>& nodal_soln)
305  { hermite_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/1); }
306 
307  template <>
308  void FE<2,HERMITE>::nodal_soln(const Elem* elem,
309  const Order order,
310  const std::vector<Number>& elem_soln,
311  std::vector<Number>& nodal_soln)
312  { hermite_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/2); }
313 
314  template <>
315  void FE<3,HERMITE>::nodal_soln(const Elem* elem,
316  const Order order,
317  const std::vector<Number>& elem_soln,
318  std::vector<Number>& nodal_soln)
319  { hermite_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/3); }
320 
321 
322 
323 
324  // Do full-specialization for every dimension, instead
325  // of explicit instantiation at the end of this function.
326  // This could be macro-ified.
327  template <> unsigned int FE<0,HERMITE>::n_dofs(const ElemType t, const Order o) { return hermite_n_dofs(t, o); }
328  template <> unsigned int FE<1,HERMITE>::n_dofs(const ElemType t, const Order o) { return hermite_n_dofs(t, o); }
329  template <> unsigned int FE<2,HERMITE>::n_dofs(const ElemType t, const Order o) { return hermite_n_dofs(t, o); }
330  template <> unsigned int FE<3,HERMITE>::n_dofs(const ElemType t, const Order o) { return hermite_n_dofs(t, o); }
331 
332  // Do full-specialization for every dimension, instead
333  // of explicit instantiation at the end of this function.
334  template <> unsigned int FE<0,HERMITE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hermite_n_dofs_at_node(t, o, n); }
335  template <> unsigned int FE<1,HERMITE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hermite_n_dofs_at_node(t, o, n); }
336  template <> unsigned int FE<2,HERMITE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hermite_n_dofs_at_node(t, o, n); }
337  template <> unsigned int FE<3,HERMITE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hermite_n_dofs_at_node(t, o, n); }
338 
339  // Full specialization of n_dofs_per_elem() function for every dimension.
340  template <> unsigned int FE<0,HERMITE>::n_dofs_per_elem(const ElemType t, const Order o) { return hermite_n_dofs_per_elem(t, o); }
341  template <> unsigned int FE<1,HERMITE>::n_dofs_per_elem(const ElemType t, const Order o) { return hermite_n_dofs_per_elem(t, o); }
342  template <> unsigned int FE<2,HERMITE>::n_dofs_per_elem(const ElemType t, const Order o) { return hermite_n_dofs_per_elem(t, o); }
343  template <> unsigned int FE<3,HERMITE>::n_dofs_per_elem(const ElemType t, const Order o) { return hermite_n_dofs_per_elem(t, o); }
344 
345  // Hermite FEMs are C^1 continuous
346  template <> FEContinuity FE<0,HERMITE>::get_continuity() const { return C_ONE; }
347  template <> FEContinuity FE<1,HERMITE>::get_continuity() const { return C_ONE; }
348  template <> FEContinuity FE<2,HERMITE>::get_continuity() const { return C_ONE; }
349  template <> FEContinuity FE<3,HERMITE>::get_continuity() const { return C_ONE; }
350 
351  // Hermite FEMs are hierarchic
352  template <> bool FE<0,HERMITE>::is_hierarchic() const { return true; }
353  template <> bool FE<1,HERMITE>::is_hierarchic() const { return true; }
354  template <> bool FE<2,HERMITE>::is_hierarchic() const { return true; }
355  template <> bool FE<3,HERMITE>::is_hierarchic() const { return true; }
356 
357 
358 #ifdef LIBMESH_ENABLE_AMR
359  // compute_constraints() specializations are only needed for 2 and 3D
360  template <>
362  DofMap &dof_map,
363  const unsigned int variable_number,
364  const Elem* elem)
365  { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
366 
367  template <>
369  DofMap &dof_map,
370  const unsigned int variable_number,
371  const Elem* elem)
372  { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
373 #endif // #ifdef LIBMESH_ENABLE_AMR
374 
375  // Hermite FEM shapes need reinit
376  template <> bool FE<0,HERMITE>::shapes_need_reinit() const { return true; }
377  template <> bool FE<1,HERMITE>::shapes_need_reinit() const { return true; }
378  template <> bool FE<2,HERMITE>::shapes_need_reinit() const { return true; }
379  template <> bool FE<3,HERMITE>::shapes_need_reinit() const { return true; }
380 
381 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo