libMesh
fe_bernstein.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/libmesh_config.h"
22 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
23 
24 #include "libmesh/fe.h"
25 #include "libmesh/elem.h"
26 #include "libmesh/fe_interface.h"
27 #include "libmesh/string_to_enum.h"
28 
29 namespace libMesh
30 {
31 
32 // ------------------------------------------------------------
33 // Bernstein-specific implementations, Steffen Petersen 2005
34 
35 // Anonymous namespace for local helper functions
36 namespace {
37 
38 void bernstein_nodal_soln(const Elem * elem,
39  const Order order,
40  const std::vector<Number> & elem_soln,
41  std::vector<Number> & nodal_soln,
42  unsigned Dim)
43 {
44  const unsigned int n_nodes = elem->n_nodes();
45 
46  const ElemType elem_type = elem->type();
47 
48  nodal_soln.resize(n_nodes);
49 
50  const Order totalorder = static_cast<Order>(order + elem->p_level());
51 
52  // FEType object to be passed to various FEInterface functions below.
53  FEType fe_type(totalorder, BERNSTEIN);
54 
55  switch (totalorder)
56  {
57  // Constant shape functions
58  case CONSTANT:
59  {
60  libmesh_assert_equal_to (elem_soln.size(), 1);
61 
62  const Number val = elem_soln[0];
63 
64  for (unsigned int n=0; n<n_nodes; n++)
65  nodal_soln[n] = val;
66 
67  return;
68  }
69 
70 
71  // For other bases do interpolation at the nodes
72  // explicitly.
73  case FIRST:
74  case SECOND:
75  case THIRD:
76  case FOURTH:
77  case FIFTH:
78  case SIXTH:
79  {
80 
81  const unsigned int n_sf =
82  // FE<Dim,T>::n_shape_functions(elem_type, totalorder);
83  FEInterface::n_shape_functions(Dim, fe_type, elem_type);
84 
85  std::vector<Point> refspace_nodes;
86  FEBase::get_refspace_nodes(elem_type,refspace_nodes);
87  libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
88 
89  for (unsigned int n=0; n<n_nodes; n++)
90  {
91  libmesh_assert_equal_to (elem_soln.size(), n_sf);
92 
93  // Zero before summation
94  nodal_soln[n] = 0;
95 
96  // u_i = Sum (alpha_i phi_i)
97  for (unsigned int i=0; i<n_sf; i++)
98  nodal_soln[n] += elem_soln[i] *
99  // FE<Dim,T>::shape(elem, order, i, mapped_point);
100  FEInterface::shape(Dim, fe_type, elem, i, refspace_nodes[n]);
101  }
102 
103  return;
104  }
105 
106  default:
107  libmesh_error_msg("ERROR: Invalid total order " << totalorder);
108  }
109 } // bernstein_nodal_soln()
110 
111 
112 
113 unsigned int bernstein_n_dofs(const ElemType t, const Order o)
114 {
115  switch (t)
116  {
117  case NODEELEM:
118  return 1;
119  case EDGE2:
120  case EDGE3:
121  return (o+1);
122  case QUAD4:
123  case QUADSHELL4:
124  libmesh_assert_less (o, 2);
125  libmesh_fallthrough();
126  case QUAD8:
127  case QUADSHELL8:
128  {
129  if (o == 1)
130  return 4;
131  else if (o == 2)
132  return 8;
133  else
134  libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for BERNSTEIN FE family!");
135  }
136  case QUAD9:
137  return ((o+1)*(o+1));
138  case HEX8:
139  libmesh_assert_less (o, 2);
140  libmesh_fallthrough();
141  case HEX20:
142  {
143  if (o == 1)
144  return 8;
145  else if (o == 2)
146  return 20;
147  else
148  libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for BERNSTEIN FE family!");
149  }
150  case HEX27:
151  return ((o+1)*(o+1)*(o+1));
152  case TRI3:
153  case TRISHELL3:
154  libmesh_assert_less (o, 2);
155  libmesh_fallthrough();
156  case TRI6:
157  return ((o+1)*(o+2)/2);
158  case TET4:
159  libmesh_assert_less (o, 2);
160  libmesh_fallthrough();
161  case TET10:
162  {
163  libmesh_assert_less (o, 3);
164  return ((o+1)*(o+2)*(o+3)/6);
165  }
166  case INVALID_ELEM:
167  return 0;
168  default:
169  libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
170  }
171 
172  libmesh_error_msg("We'll never get here!");
173  return 0;
174 } // bernstein_n_dofs()
175 
176 
177 
178 
179 unsigned int bernstein_n_dofs_at_node(const ElemType t,
180  const Order o,
181  const unsigned int n)
182 {
183  switch (t)
184  {
185  case NODEELEM:
186  return 1;
187 
188  case EDGE2:
189  case EDGE3:
190  switch (n)
191  {
192  case 0:
193  case 1:
194  return 1;
195  case 2:
196  libmesh_assert (o>1);
197  return (o-1);
198  default:
199  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for EDGE2/3!");
200  }
201  case TRI3:
202  libmesh_assert_less (n, 3);
203  libmesh_assert_less (o, 2);
204  libmesh_fallthrough();
205  case TRI6:
206  switch (n)
207  {
208  case 0:
209  case 1:
210  case 2:
211  return 1;
212 
213  case 3:
214  case 4:
215  case 5:
216  return (o-1);
217  // Internal DoFs are associated with the elem, not its nodes
218  default:
219  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
220  }
221  case QUAD4:
222  libmesh_assert_less (n, 4);
223  libmesh_assert_less (o, 2);
224  libmesh_fallthrough();
225  case QUAD8:
226  libmesh_assert_less (n, 8);
227  libmesh_assert_less (o, 3);
228  libmesh_fallthrough();
229  case QUAD9:
230  {
231  switch (n)
232  {
233  case 0:
234  case 1:
235  case 2:
236  case 3:
237  return 1;
238 
239  case 4:
240  case 5:
241  case 6:
242  case 7:
243  return (o-1);
244 
245  // Internal DoFs are associated with the elem, not its nodes
246  case 8:
247  return 0;
248 
249  default:
250  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for QUAD9!");
251  }
252  }
253  case HEX8:
254  libmesh_assert_less (n, 8);
255  libmesh_assert_less (o, 2);
256  libmesh_fallthrough();
257  case HEX20:
258  libmesh_assert_less (n, 20);
259  libmesh_assert_less (o, 3);
260  libmesh_fallthrough();
261  case HEX27:
262  switch (n)
263  {
264  case 0:
265  case 1:
266  case 2:
267  case 3:
268  case 4:
269  case 5:
270  case 6:
271  case 7:
272  return 1;
273 
274  case 8:
275  case 9:
276  case 10:
277  case 11:
278  case 12:
279  case 13:
280  case 14:
281  case 15:
282  case 16:
283  case 17:
284  case 18:
285  case 19:
286  return (o-1);
287 
288  case 20:
289  case 21:
290  case 22:
291  case 23:
292  case 24:
293  case 25:
294  return ((o-1)*(o-1));
295 
296  // Internal DoFs are associated with the elem, not its nodes
297  case 26:
298  return 0;
299 
300  default:
301  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for HEX27!");
302  }
303  case TET4:
304  libmesh_assert_less (n, 4);
305  libmesh_assert_less (o, 2);
306  libmesh_fallthrough();
307  case TET10:
308  libmesh_assert_less (o, 3);
309  libmesh_assert_less (n, 10);
310  switch (n)
311  {
312  case 0:
313  case 1:
314  case 2:
315  case 3:
316  return 1;
317 
318  case 4:
319  case 5:
320  case 6:
321  case 7:
322  case 8:
323  case 9:
324  return (o-1);
325 
326  default:
327  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TET10!");
328  }
329  case INVALID_ELEM:
330  return 0;
331  default:
332  libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
333  }
334 
335  libmesh_error_msg("We'll never get here!");
336  return 0;
337 } // bernstein_n_dofs_at_node()
338 
339 
340 
341 
342 unsigned int bernstein_n_dofs_per_elem(const ElemType t, const Order o)
343 {
344  switch (t)
345  {
346  case NODEELEM:
347  case EDGE2:
348  case EDGE3:
349  return 0;
350  case TRI3:
351  case TRISHELL3:
352  case QUAD4:
353  case QUADSHELL4:
354  return 0;
355  case TRI6:
356  return ((o-1)*(o-2)/2);
357  case QUAD8:
358  case QUADSHELL8:
359  if (o <= 2)
360  return 0;
361  libmesh_fallthrough();
362  case QUAD9:
363  return ((o-1)*(o-1));
364  case HEX8:
365  libmesh_assert_less (o, 2);
366  libmesh_fallthrough();
367  case HEX20:
368  libmesh_assert_less (o, 3);
369  return 0;
370  case HEX27:
371  return ((o-1)*(o-1)*(o-1));
372  case TET4:
373  libmesh_assert_less (o, 2);
374  libmesh_fallthrough();
375  case TET10:
376  libmesh_assert_less (o, 3);
377  return 0;
378  case INVALID_ELEM:
379  return 0;
380  default:
381  libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for BERNSTEIN FE family!");
382  }
383 
384  libmesh_error_msg("We'll never get here!");
385  return 0;
386 } // bernstein_n_dofs_per_elem
387 
388 } // anonymous namespace
389 
390 
391 
392 
393  // Do full-specialization of nodal_soln() function for every
394  // dimension, instead of explicit instantiation at the end of this
395  // file.
396  // This could be macro-ified so that it fits on one line...
397 template <>
399  const Order order,
400  const std::vector<Number> & elem_soln,
401  std::vector<Number> & nodal_soln)
402 { bernstein_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/0); }
403 
404 template <>
406  const Order order,
407  const std::vector<Number> & elem_soln,
408  std::vector<Number> & nodal_soln)
409 { bernstein_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/1); }
410 
411 template <>
413  const Order order,
414  const std::vector<Number> & elem_soln,
415  std::vector<Number> & nodal_soln)
416 { bernstein_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/2); }
417 
418 template <>
420  const Order order,
421  const std::vector<Number> & elem_soln,
422  std::vector<Number> & nodal_soln)
423 { bernstein_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/3); }
424 
425 
426 // Full specialization of n_dofs() function for every dimension
427 template <> unsigned int FE<0,BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return bernstein_n_dofs(t, o); }
428 template <> unsigned int FE<1,BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return bernstein_n_dofs(t, o); }
429 template <> unsigned int FE<2,BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return bernstein_n_dofs(t, o); }
430 template <> unsigned int FE<3,BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return bernstein_n_dofs(t, o); }
431 
432 // Full specialization of n_dofs_at_node() function for every dimension.
433 template <> unsigned int FE<0,BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return bernstein_n_dofs_at_node(t, o, n); }
434 template <> unsigned int FE<1,BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return bernstein_n_dofs_at_node(t, o, n); }
435 template <> unsigned int FE<2,BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return bernstein_n_dofs_at_node(t, o, n); }
436 template <> unsigned int FE<3,BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return bernstein_n_dofs_at_node(t, o, n); }
437 
438 // Full specialization of n_dofs_per_elem() function for every dimension.
439 template <> unsigned int FE<0,BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return bernstein_n_dofs_per_elem(t, o); }
440 template <> unsigned int FE<1,BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return bernstein_n_dofs_per_elem(t, o); }
441 template <> unsigned int FE<2,BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return bernstein_n_dofs_per_elem(t, o); }
442 template <> unsigned int FE<3,BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return bernstein_n_dofs_per_elem(t, o); }
443 
444 // Bernstein FEMs are C^0 continuous
445 template <> FEContinuity FE<0,BERNSTEIN>::get_continuity() const { return C_ZERO; }
446 template <> FEContinuity FE<1,BERNSTEIN>::get_continuity() const { return C_ZERO; }
447 template <> FEContinuity FE<2,BERNSTEIN>::get_continuity() const { return C_ZERO; }
448 template <> FEContinuity FE<3,BERNSTEIN>::get_continuity() const { return C_ZERO; }
449 
450 // Bernstein FEMs are not hierarchic
451 template <> bool FE<0,BERNSTEIN>::is_hierarchic() const { return false; }
452 template <> bool FE<1,BERNSTEIN>::is_hierarchic() const { return false; }
453 template <> bool FE<2,BERNSTEIN>::is_hierarchic() const { return false; }
454 template <> bool FE<3,BERNSTEIN>::is_hierarchic() const { return false; }
455 
456 #ifdef LIBMESH_ENABLE_AMR
457 // compute_constraints() specializations are only needed for 2 and 3D
458 template <>
460  DofMap & dof_map,
461  const unsigned int variable_number,
462  const Elem * elem)
463 { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
464 
465 template <>
467  DofMap & dof_map,
468  const unsigned int variable_number,
469  const Elem * elem)
470 { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
471 #endif // #ifdef LIBMESH_ENABLE_AMR
472 
473 // Bernstein shapes need reinit only for approximation orders >= 3,
474 // but we might reach that with p refinement
475 template <> bool FE<0,BERNSTEIN>::shapes_need_reinit() const { return true; }
476 template <> bool FE<1,BERNSTEIN>::shapes_need_reinit() const { return true; }
477 template <> bool FE<2,BERNSTEIN>::shapes_need_reinit() const { return true; }
478 template <> bool FE<3,BERNSTEIN>::shapes_need_reinit() const { return true; }
479 
480 } // namespace libMesh
481 
482 #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
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.
libmesh_assert(j)
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