fe_monomial.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 } // monomial_n_dofs()
341 
342 
343 } // anonymous namespace
344 
345 
346 
347 
348 
349  // Do full-specialization for every dimension, instead
350  // of explicit instantiation at the end of this file.
351  // This could be macro-ified so that it fits on one line...
352 template <>
354  const Order order,
355  const std::vector<Number> & elem_soln,
356  std::vector<Number> & nodal_soln)
357 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/0); }
358 
359 template <>
361  const Order order,
362  const std::vector<Number> & elem_soln,
363  std::vector<Number> & nodal_soln)
364 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/1); }
365 
366 template <>
368  const Order order,
369  const std::vector<Number> & elem_soln,
370  std::vector<Number> & nodal_soln)
371 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/2); }
372 
373 template <>
375  const Order order,
376  const std::vector<Number> & elem_soln,
377  std::vector<Number> & nodal_soln)
378 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/3); }
379 
380 
381 // Full specialization of n_dofs() function for every dimension
382 template <> unsigned int FE<0,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
383 template <> unsigned int FE<1,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
384 template <> unsigned int FE<2,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
385 template <> unsigned int FE<3,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
386 
387 // Full specialization of n_dofs_at_node() function for every dimension.
388 // Monomials have no dofs at nodes, only element dofs.
389 template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
390 template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
391 template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
392 template <> unsigned int FE<3,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
393 
394 // Full specialization of n_dofs_per_elem() function for every dimension.
395 template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
396 template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
397 template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
398 template <> unsigned int FE<3,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
399 
400 
401 // Full specialization of get_continuity() function for every dimension.
406 
407 // Full specialization of is_hierarchic() function for every dimension.
408 // The monomials are hierarchic!
409 template <> bool FE<0,MONOMIAL>::is_hierarchic() const { return true; }
410 template <> bool FE<1,MONOMIAL>::is_hierarchic() const { return true; }
411 template <> bool FE<2,MONOMIAL>::is_hierarchic() const { return true; }
412 template <> bool FE<3,MONOMIAL>::is_hierarchic() const { return true; }
413 
414 #ifdef LIBMESH_ENABLE_AMR
415 
416 // Full specialization of compute_constraints() function for 2D and
417 // 3D only. There are no constraints for discontinuous elements, so
418 // we do nothing.
419 template <> void FE<2,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
420 template <> void FE<3,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
421 
422 #endif // #ifdef LIBMESH_ENABLE_AMR
423 
424 // Full specialization of shapes_need_reinit() function for every dimension.
425 template <> bool FE<0,MONOMIAL>::shapes_need_reinit() const { return false; }
426 template <> bool FE<1,MONOMIAL>::shapes_need_reinit() const { return false; }
427 template <> bool FE<2,MONOMIAL>::shapes_need_reinit() const { return false; }
428 template <> bool FE<3,MONOMIAL>::shapes_need_reinit() const { return false; }
429 
430 } // namespace libMesh
static unsigned int n_dofs(const ElemType t, const Order o)
The base class for all geometric element types.
Definition: elem.h:100
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
virtual bool shapes_need_reinit() const override
virtual bool is_hierarchic() const override
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
const dof_id_type n_nodes
Definition: tecplot_io.C:68
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:428
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:657
static void get_refspace_nodes(const ElemType t, std::vector< Point > &nodes)
Definition: fe_abstract.C:283
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
virtual FEContinuity get_continuity() const override
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)
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)