libMesh
fe_hierarchic.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 // Hierarchic-specific implementations
31 
32 // Anonymous namespace for local helper functions
33 namespace {
34 
35 void hierarchic_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, HIERARCHIC);
51 
52  switch (totalorder)
53  {
54  // Constant shape functions
55  case CONSTANT:
56  {
57  libmesh_assert_equal_to (elem_soln.size(), 1);
58 
59  const Number val = elem_soln[0];
60 
61  for (unsigned int n=0; n<n_nodes; n++)
62  nodal_soln[n] = val;
63 
64  return;
65  }
66 
67 
68  // For other orders do interpolation at the nodes
69  // explicitly.
70  default:
71  {
72 
73  const unsigned int n_sf =
74  // FE<Dim,T>::n_shape_functions(elem_type, totalorder);
75  FEInterface::n_shape_functions(Dim, fe_type, elem_type);
76 
77  std::vector<Point> refspace_nodes;
78  FEBase::get_refspace_nodes(elem_type,refspace_nodes);
79  libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
80 
81  for (unsigned int n=0; n<n_nodes; n++)
82  {
83  libmesh_assert_equal_to (elem_soln.size(), n_sf);
84 
85  // Zero before summation
86  nodal_soln[n] = 0;
87 
88  // u_i = Sum (alpha_i phi_i)
89  for (unsigned int i=0; i<n_sf; i++)
90  nodal_soln[n] += elem_soln[i] *
91  // FE<Dim,T>::shape(elem, order, i, mapped_point);
92  FEInterface::shape(Dim, fe_type, elem, i, refspace_nodes[n]);
93  }
94 
95  return;
96  }
97  }
98 } // hierarchic_nodal_soln()
99 
100 
101 
102 
103 unsigned int hierarchic_n_dofs(const ElemType t, const Order o)
104 {
105  libmesh_assert_greater (o, 0);
106  switch (t)
107  {
108  case NODEELEM:
109  return 1;
110  case EDGE2:
111  case EDGE3:
112  return (o+1);
113  case QUAD4:
114  case QUADSHELL4:
115  libmesh_assert_less (o, 2);
116  libmesh_fallthrough();
117  case QUAD8:
118  case QUADSHELL8:
119  case QUAD9:
120  return ((o+1)*(o+1));
121  case HEX8:
122  libmesh_assert_less (o, 2);
123  libmesh_fallthrough();
124  case HEX20:
125  libmesh_assert_less (o, 2);
126  libmesh_fallthrough();
127  case HEX27:
128  return ((o+1)*(o+1)*(o+1));
129  case TRI3:
130  libmesh_assert_less (o, 2);
131  libmesh_fallthrough();
132  case TRI6:
133  return ((o+1)*(o+2)/2);
134  case INVALID_ELEM:
135  return 0;
136  default:
137  libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for HIERARCHIC FE family!");
138  }
139 
140  libmesh_error_msg("We'll never get here!");
141  return 0;
142 } // hierarchic_n_dofs()
143 
144 
145 
146 
147 unsigned int hierarchic_n_dofs_at_node(const ElemType t,
148  const Order o,
149  const unsigned int n)
150 {
151  libmesh_assert_greater (o, 0);
152  switch (t)
153  {
154  case NODEELEM:
155  return 1;
156  case EDGE2:
157  case EDGE3:
158  switch (n)
159  {
160  case 0:
161  case 1:
162  return 1;
163  // Internal DoFs are associated with the elem, not its nodes
164  case 2:
165  return 0;
166  default:
167  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
168  }
169  case TRI3:
170  libmesh_assert_less (n, 3);
171  libmesh_assert_less (o, 2);
172  libmesh_fallthrough();
173  case TRI6:
174  switch (n)
175  {
176  case 0:
177  case 1:
178  case 2:
179  return 1;
180 
181  case 3:
182  case 4:
183  case 5:
184  return (o-1);
185 
186  // Internal DoFs are associated with the elem, not its nodes
187  default:
188  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
189  }
190  case QUAD4:
191  case QUADSHELL4:
192  libmesh_assert_less (n, 4);
193  libmesh_assert_less (o, 2);
194  libmesh_fallthrough();
195  case QUAD8:
196  case QUADSHELL8:
197  case QUAD9:
198  switch (n)
199  {
200  case 0:
201  case 1:
202  case 2:
203  case 3:
204  return 1;
205 
206  case 4:
207  case 5:
208  case 6:
209  case 7:
210  return (o-1);
211 
212  // Internal DoFs are associated with the elem, not its nodes
213  case 8:
214  return 0;
215 
216  default:
217  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD4/8/9!");
218  }
219  case HEX8:
220  libmesh_assert_less (n, 8);
221  libmesh_assert_less (o, 2);
222  libmesh_fallthrough();
223  case HEX20:
224  libmesh_assert_less (n, 20);
225  libmesh_assert_less (o, 2);
226  libmesh_fallthrough();
227  case HEX27:
228  switch (n)
229  {
230  case 0:
231  case 1:
232  case 2:
233  case 3:
234  case 4:
235  case 5:
236  case 6:
237  case 7:
238  return 1;
239 
240  case 8:
241  case 9:
242  case 10:
243  case 11:
244  case 12:
245  case 13:
246  case 14:
247  case 15:
248  case 16:
249  case 17:
250  case 18:
251  case 19:
252  return (o-1);
253 
254  case 20:
255  case 21:
256  case 22:
257  case 23:
258  case 24:
259  case 25:
260  return ((o-1)*(o-1));
261 
262  // Internal DoFs are associated with the elem, not its nodes
263  case 26:
264  return 0;
265  default:
266  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for HEX8/20/27!");
267  }
268 
269  case INVALID_ELEM:
270  return 0;
271 
272  default:
273  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
274  }
275 
276  libmesh_error_msg("We'll never get here!");
277  return 0;
278 } // hierarchic_n_dofs_at_node()
279 
280 
281 
282 
283 unsigned int hierarchic_n_dofs_per_elem(const ElemType t,
284  const Order o)
285 {
286  libmesh_assert_greater (o, 0);
287  switch (t)
288  {
289  case NODEELEM:
290  return 0;
291  case EDGE2:
292  case EDGE3:
293  return (o-1);
294  case TRI3:
295  case TRISHELL3:
296  case QUAD4:
297  case QUADSHELL4:
298  return 0;
299  case TRI6:
300  return ((o-1)*(o-2)/2);
301  case QUAD8:
302  case QUADSHELL8:
303  case QUAD9:
304  return ((o-1)*(o-1));
305  case HEX8:
306  case HEX20:
307  libmesh_assert_less (o, 2);
308  return 0;
309  case HEX27:
310  return ((o-1)*(o-1)*(o-1));
311  case INVALID_ELEM:
312  return 0;
313  default:
314  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
315  }
316 
317  libmesh_error_msg("We'll never get here!");
318  return 0;
319 } // hierarchic_n_dofs_per_elem()
320 
321 } // anonymous namespace
322 
323 
324 
325 
326  // Do full-specialization of nodal_soln() function for every
327  // dimension, instead of explicit instantiation at the end of this
328  // file.
329  // This could be macro-ified so that it fits on one line...
330 template <>
332  const Order order,
333  const std::vector<Number> & elem_soln,
334  std::vector<Number> & nodal_soln)
335 { hierarchic_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/0); }
336 
337 template <>
339  const Order order,
340  const std::vector<Number> & elem_soln,
341  std::vector<Number> & nodal_soln)
342 { hierarchic_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/1); }
343 
344 template <>
346  const Order order,
347  const std::vector<Number> & elem_soln,
348  std::vector<Number> & nodal_soln)
349 { hierarchic_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/2); }
350 
351 template <>
353  const Order order,
354  const std::vector<Number> & elem_soln,
355  std::vector<Number> & nodal_soln)
356 { hierarchic_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/3); }
357 
358 
359 // Full specialization of n_dofs() function for every dimension
360 template <> unsigned int FE<0,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
361 template <> unsigned int FE<1,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
362 template <> unsigned int FE<2,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
363 template <> unsigned int FE<3,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
364 
365 // Full specialization of n_dofs_at_node() function for every dimension.
366 template <> unsigned int FE<0,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
367 template <> unsigned int FE<1,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
368 template <> unsigned int FE<2,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
369 template <> unsigned int FE<3,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
370 
371 // Full specialization of n_dofs_per_elem() function for every dimension.
372 template <> unsigned int FE<0,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
373 template <> unsigned int FE<1,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
374 template <> unsigned int FE<2,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
375 template <> unsigned int FE<3,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
376 
377 // Hierarchic FEMs are C^0 continuous
378 template <> FEContinuity FE<0,HIERARCHIC>::get_continuity() const { return C_ZERO; }
379 template <> FEContinuity FE<1,HIERARCHIC>::get_continuity() const { return C_ZERO; }
380 template <> FEContinuity FE<2,HIERARCHIC>::get_continuity() const { return C_ZERO; }
381 template <> FEContinuity FE<3,HIERARCHIC>::get_continuity() const { return C_ZERO; }
382 
383 // Hierarchic FEMs are hierarchic (duh!)
384 template <> bool FE<0,HIERARCHIC>::is_hierarchic() const { return true; }
385 template <> bool FE<1,HIERARCHIC>::is_hierarchic() const { return true; }
386 template <> bool FE<2,HIERARCHIC>::is_hierarchic() const { return true; }
387 template <> bool FE<3,HIERARCHIC>::is_hierarchic() const { return true; }
388 
389 #ifdef LIBMESH_ENABLE_AMR
390 // compute_constraints() specializations are only needed for 2 and 3D
391 template <>
393  DofMap & dof_map,
394  const unsigned int variable_number,
395  const Elem * elem)
396 { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
397 
398 template <>
400  DofMap & dof_map,
401  const unsigned int variable_number,
402  const Elem * elem)
403 { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
404 #endif // #ifdef LIBMESH_ENABLE_AMR
405 
406 // Hierarchic FEM shapes need reinit
407 template <> bool FE<0,HIERARCHIC>::shapes_need_reinit() const { return true; }
408 template <> bool FE<1,HIERARCHIC>::shapes_need_reinit() const { return true; }
409 template <> bool FE<2,HIERARCHIC>::shapes_need_reinit() const { return true; }
410 template <> bool FE<3,HIERARCHIC>::shapes_need_reinit() const { return true; }
411 
412 } // 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