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