libMesh
fe_monomial.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/fe.h"
22 #include "libmesh/elem.h"
23 #include "libmesh/fe_interface.h"
24 #include "libmesh/string_to_enum.h"
25 
26 namespace libMesh
27 {
28 
29 // ------------------------------------------------------------
30 // Monomials-specific implementations
31 
32 
33 // Anonymous namespace for local helper functions
34 namespace {
35 
36 void monomial_nodal_soln(const Elem * elem,
37  const Order order,
38  const std::vector<Number> & elem_soln,
39  std::vector<Number> & nodal_soln,
40  const unsigned Dim)
41 {
42  const unsigned int n_nodes = elem->n_nodes();
43 
44  const ElemType elem_type = elem->type();
45 
46  nodal_soln.resize(n_nodes);
47 
48  const Order totalorder = static_cast<Order>(order+elem->p_level());
49 
50  switch (totalorder)
51  {
52  // Constant shape functions
53  case CONSTANT:
54  {
55  libmesh_assert_equal_to (elem_soln.size(), 1);
56 
57  const Number val = elem_soln[0];
58 
59  for (unsigned int n=0; n<n_nodes; n++)
60  nodal_soln[n] = val;
61 
62  return;
63  }
64 
65 
66  // For other orders, do interpolation at the nodes
67  // explicitly.
68  default:
69  {
70  // FEType object to be passed to various FEInterface functions below.
71  FEType fe_type(totalorder, MONOMIAL);
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  } // default
97  } // switch
98 } // monomial_nodal_soln()
99 
100 
101 
102 
103 unsigned int monomial_n_dofs(const ElemType t, const Order o)
104 {
105  switch (o)
106  {
107 
108  // constant shape functions
109  // no matter what shape there is only one DOF.
110  case CONSTANT:
111  return (t != INVALID_ELEM) ? 1 : 0;
112 
113 
114  // Discontinuous linear shape functions
115  // expressed in the monomials.
116  case FIRST:
117  {
118  switch (t)
119  {
120  case NODEELEM:
121  return 1;
122 
123  case EDGE2:
124  case EDGE3:
125  case EDGE4:
126  return 2;
127 
128  case TRI3:
129  case TRISHELL3:
130  case TRI6:
131  case QUAD4:
132  case QUADSHELL4:
133  case QUAD8:
134  case QUADSHELL8:
135  case QUAD9:
136  return 3;
137 
138  case TET4:
139  case TET10:
140  case HEX8:
141  case HEX20:
142  case HEX27:
143  case PRISM6:
144  case PRISM15:
145  case PRISM18:
146  case PYRAMID5:
147  case PYRAMID13:
148  case PYRAMID14:
149  return 4;
150 
151  case INVALID_ELEM:
152  return 0;
153 
154  default:
155  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
156  }
157  }
158 
159 
160  // Discontinuous quadratic shape functions
161  // expressed in the monomials.
162  case SECOND:
163  {
164  switch (t)
165  {
166  case NODEELEM:
167  return 1;
168 
169  case EDGE2:
170  case EDGE3:
171  case EDGE4:
172  return 3;
173 
174  case TRI3:
175  case TRISHELL3:
176  case TRI6:
177  case QUAD4:
178  case QUADSHELL4:
179  case QUAD8:
180  case QUADSHELL8:
181  case QUAD9:
182  return 6;
183 
184  case TET4:
185  case TET10:
186  case HEX8:
187  case HEX20:
188  case HEX27:
189  case PRISM6:
190  case PRISM15:
191  case PRISM18:
192  case PYRAMID5:
193  case PYRAMID13:
194  case PYRAMID14:
195  return 10;
196 
197  case INVALID_ELEM:
198  return 0;
199 
200  default:
201  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
202  }
203  }
204 
205 
206  // Discontinuous cubic shape functions
207  // expressed in the monomials.
208  case THIRD:
209  {
210  switch (t)
211  {
212  case NODEELEM:
213  return 1;
214 
215  case EDGE2:
216  case EDGE3:
217  case EDGE4:
218  return 4;
219 
220  case TRI3:
221  case TRISHELL3:
222  case TRI6:
223  case QUAD4:
224  case QUADSHELL4:
225  case QUAD8:
226  case QUADSHELL8:
227  case QUAD9:
228  return 10;
229 
230  case TET4:
231  case TET10:
232  case HEX8:
233  case HEX20:
234  case HEX27:
235  case PRISM6:
236  case PRISM15:
237  case PRISM18:
238  case PYRAMID5:
239  case PYRAMID13:
240  case PYRAMID14:
241  return 20;
242 
243  case INVALID_ELEM:
244  return 0;
245 
246  default:
247  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
248  }
249  }
250 
251 
252  // Discontinuous quartic shape functions
253  // expressed in the monomials.
254  case FOURTH:
255  {
256  switch (t)
257  {
258  case NODEELEM:
259  return 1;
260 
261  case EDGE2:
262  case EDGE3:
263  return 5;
264 
265  case TRI3:
266  case TRISHELL3:
267  case TRI6:
268  case QUAD4:
269  case QUADSHELL4:
270  case QUAD8:
271  case QUADSHELL8:
272  case QUAD9:
273  return 15;
274 
275  case TET4:
276  case TET10:
277  case HEX8:
278  case HEX20:
279  case HEX27:
280  case PRISM6:
281  case PRISM15:
282  case PRISM18:
283  case PYRAMID5:
284  case PYRAMID13:
285  case PYRAMID14:
286  return 35;
287 
288  case INVALID_ELEM:
289  return 0;
290 
291  default:
292  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
293  }
294  }
295 
296 
297  default:
298  {
299  const unsigned int order = static_cast<unsigned int>(o);
300  switch (t)
301  {
302  case NODEELEM:
303  return 1;
304 
305  case EDGE2:
306  case EDGE3:
307  return (order+1);
308 
309  case TRI3:
310  case TRISHELL3:
311  case TRI6:
312  case QUAD4:
313  case QUADSHELL4:
314  case QUAD8:
315  case QUADSHELL8:
316  case QUAD9:
317  return (order+1)*(order+2)/2;
318 
319  case TET4:
320  case TET10:
321  case HEX8:
322  case HEX20:
323  case HEX27:
324  case PRISM6:
325  case PRISM15:
326  case PRISM18:
327  case PYRAMID5:
328  case PYRAMID13:
329  case PYRAMID14:
330  return (order+1)*(order+2)*(order+3)/6;
331 
332  case INVALID_ELEM:
333  return 0;
334 
335  default:
336  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
337  }
338  }
339  }
340 
341  libmesh_error_msg("We'll never get here!");
342  return 0;
343 } // monomial_n_dofs()
344 
345 
346 } // anonymous namespace
347 
348 
349 
350 
351 
352  // Do full-specialization for every dimension, instead
353  // of explicit instantiation at the end of this file.
354  // This could be macro-ified so that it fits on one line...
355 template <>
357  const Order order,
358  const std::vector<Number> & elem_soln,
359  std::vector<Number> & nodal_soln)
360 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/0); }
361 
362 template <>
364  const Order order,
365  const std::vector<Number> & elem_soln,
366  std::vector<Number> & nodal_soln)
367 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/1); }
368 
369 template <>
371  const Order order,
372  const std::vector<Number> & elem_soln,
373  std::vector<Number> & nodal_soln)
374 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/2); }
375 
376 template <>
378  const Order order,
379  const std::vector<Number> & elem_soln,
380  std::vector<Number> & nodal_soln)
381 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/3); }
382 
383 
384 // Full specialization of n_dofs() function for every dimension
385 template <> unsigned int FE<0,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
386 template <> unsigned int FE<1,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
387 template <> unsigned int FE<2,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
388 template <> unsigned int FE<3,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
389 
390 // Full specialization of n_dofs_at_node() function for every dimension.
391 // Monomials have no dofs at nodes, only element dofs.
392 template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
393 template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
394 template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
395 template <> unsigned int FE<3,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
396 
397 // Full specialization of n_dofs_per_elem() function for every dimension.
398 template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
399 template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
400 template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
401 template <> unsigned int FE<3,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
402 
403 
404 // Full specialization of get_continuity() function for every dimension.
409 
410 // Full specialization of is_hierarchic() function for every dimension.
411 // The monomials are hierarchic!
412 template <> bool FE<0,MONOMIAL>::is_hierarchic() const { return true; }
413 template <> bool FE<1,MONOMIAL>::is_hierarchic() const { return true; }
414 template <> bool FE<2,MONOMIAL>::is_hierarchic() const { return true; }
415 template <> bool FE<3,MONOMIAL>::is_hierarchic() const { return true; }
416 
417 #ifdef LIBMESH_ENABLE_AMR
418 
419 // Full specialization of compute_constraints() function for 2D and
420 // 3D only. There are no constraints for discontinuous elements, so
421 // we do nothing.
422 template <> void FE<2,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
423 template <> void FE<3,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
424 
425 #endif // #ifdef LIBMESH_ENABLE_AMR
426 
427 // Full specialization of shapes_need_reinit() function for every dimension.
428 template <> bool FE<0,MONOMIAL>::shapes_need_reinit() const { return false; }
429 template <> bool FE<1,MONOMIAL>::shapes_need_reinit() const { return false; }
430 template <> bool FE<2,MONOMIAL>::shapes_need_reinit() const { return false; }
431 template <> bool FE<3,MONOMIAL>::shapes_need_reinit() const { return false; }
432 
433 } // 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