fe_monomial.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
25 namespace libMesh
26 {
27 
28  // ------------------------------------------------------------
29  // Monomials-specific implementations
30 
31 
32  // Anonymous namespace for local helper functions
33  namespace {
34 
35  void monomial_nodal_soln(const Elem* elem,
36  const Order order,
37  const std::vector<Number>& elem_soln,
38  std::vector<Number>& nodal_soln,
39  const 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  switch (totalorder)
50  {
51  // Constant shape functions
52  case CONSTANT:
53  {
54  libmesh_assert_equal_to (elem_soln.size(), 1);
55 
56  const Number val = elem_soln[0];
57 
58  for (unsigned int n=0; n<n_nodes; n++)
59  nodal_soln[n] = val;
60 
61  return;
62  }
63 
64 
65  // For other orders, do interpolation at the nodes
66  // explicitly.
67  default:
68  {
69  // FEType object to be passed to various FEInterface functions below.
70  FEType fe_type(totalorder, MONOMIAL);
71 
72  const unsigned int n_sf =
73  // FE<Dim,T>::n_shape_functions(elem_type, totalorder);
74  FEInterface::n_shape_functions(Dim, fe_type, elem_type);
75 
76  std::vector<Point> refspace_nodes;
77  FEBase::get_refspace_nodes(elem_type,refspace_nodes);
78  libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
79 
80  for (unsigned int n=0; n<n_nodes; n++)
81  {
82  libmesh_assert_equal_to (elem_soln.size(), n_sf);
83 
84  // Zero before summation
85  nodal_soln[n] = 0;
86 
87  // u_i = Sum (alpha_i phi_i)
88  for (unsigned int i=0; i<n_sf; i++)
89  nodal_soln[n] += elem_soln[i] *
90  // FE<Dim,T>::shape(elem, order, i, mapped_point);
91  FEInterface::shape(Dim, fe_type, elem, i, refspace_nodes[n]);
92  }
93 
94  return;
95  } // default
96  } // switch
97  } // monomial_nodal_soln()
98 
99 
100 
101 
102  unsigned int monomial_n_dofs(const ElemType t, const Order o)
103  {
104  switch (o)
105  {
106 
107  // constant shape functions
108  // no matter what shape there is only one DOF.
109  case CONSTANT:
110  return 1;
111 
112 
113  // Discontinuous linear shape functions
114  // expressed in the monomials.
115  case FIRST:
116  {
117  switch (t)
118  {
119  case NODEELEM:
120  return 1;
121 
122  case EDGE2:
123  case EDGE3:
124  case EDGE4:
125  return 2;
126 
127  case TRI3:
128  case TRI6:
129  case QUAD4:
130  case QUAD8:
131  case QUAD9:
132  return 3;
133 
134  case TET4:
135  case TET10:
136  case HEX8:
137  case HEX20:
138  case HEX27:
139  case PRISM6:
140  case PRISM15:
141  case PRISM18:
142  case PYRAMID5:
143  case PYRAMID14:
144  return 4;
145 
146  default:
147  {
148 #ifdef DEBUG
149  libMesh::err << "ERROR: Bad ElemType = " << t
150  << " for " << o << "th order approximation!"
151  << std::endl;
152 #endif
153  libmesh_error();
154  }
155  }
156  }
157 
158 
159  // Discontinuous quadratic shape functions
160  // expressed in the monomials.
161  case SECOND:
162  {
163  switch (t)
164  {
165  case NODEELEM:
166  return 1;
167 
168  case EDGE2:
169  case EDGE3:
170  case EDGE4:
171  return 3;
172 
173  case TRI3:
174  case TRI6:
175  case QUAD4:
176  case QUAD8:
177  case QUAD9:
178  return 6;
179 
180  case TET4:
181  case TET10:
182  case HEX8:
183  case HEX20:
184  case HEX27:
185  case PRISM6:
186  case PRISM15:
187  case PRISM18:
188  case PYRAMID5:
189  case PYRAMID14:
190  return 10;
191 
192  default:
193  {
194 #ifdef DEBUG
195  libMesh::err << "ERROR: Bad ElemType = " << t
196  << " for " << o << "th order approximation!"
197  << std::endl;
198 #endif
199  libmesh_error();
200  }
201  }
202  }
203 
204 
205  // Discontinuous cubic shape functions
206  // expressed in the monomials.
207  case THIRD:
208  {
209  switch (t)
210  {
211  case NODEELEM:
212  return 1;
213 
214  case EDGE2:
215  case EDGE3:
216  case EDGE4:
217  return 4;
218 
219  case TRI3:
220  case TRI6:
221  case QUAD4:
222  case QUAD8:
223  case QUAD9:
224  return 10;
225 
226  case TET4:
227  case TET10:
228  case HEX8:
229  case HEX20:
230  case HEX27:
231  case PRISM6:
232  case PRISM15:
233  case PRISM18:
234  case PYRAMID5:
235  case PYRAMID14:
236  return 20;
237 
238  default:
239  {
240 #ifdef DEBUG
241  libMesh::err << "ERROR: Bad ElemType = " << t
242  << " for " << o << "th order approximation!"
243  << std::endl;
244 #endif
245  libmesh_error();
246  }
247  }
248  }
249 
250 
251  // Discontinuous quartic shape functions
252  // expressed in the monomials.
253  case FOURTH:
254  {
255  switch (t)
256  {
257  case NODEELEM:
258  return 1;
259 
260  case EDGE2:
261  case EDGE3:
262  return 5;
263 
264  case TRI3:
265  case TRI6:
266  case QUAD4:
267  case QUAD8:
268  case QUAD9:
269  return 15;
270 
271  case TET4:
272  case TET10:
273  case HEX8:
274  case HEX20:
275  case HEX27:
276  case PRISM6:
277  case PRISM15:
278  case PRISM18:
279  case PYRAMID5:
280  case PYRAMID14:
281  return 35;
282 
283  default:
284  {
285 #ifdef DEBUG
286  libMesh::err << "ERROR: Bad ElemType = " << t
287  << " for " << o << "th order approximation!"
288  << std::endl;
289 #endif
290  libmesh_error();
291  }
292  }
293  }
294 
295 
296  default:
297  {
298  const unsigned int order = static_cast<unsigned int>(o);
299  switch (t)
300  {
301  case NODEELEM:
302  return 1;
303 
304  case EDGE2:
305  case EDGE3:
306  return (order+1);
307 
308  case TRI3:
309  case TRI6:
310  case QUAD4:
311  case QUAD8:
312  case QUAD9:
313  return (order+1)*(order+2)/2;
314 
315  case TET4:
316  case TET10:
317  case HEX8:
318  case HEX20:
319  case HEX27:
320  case PRISM6:
321  case PRISM15:
322  case PRISM18:
323  case PYRAMID5:
324  case PYRAMID14:
325  return (order+1)*(order+2)*(order+3)/6;
326 
327  default:
328  {
329 #ifdef DEBUG
330  libMesh::err << "ERROR: Bad ElemType = " << t
331  << " for " << o << "th order approximation!"
332  << std::endl;
333 #endif
334  libmesh_error();
335  }
336  }
337  }
338  }
339 
340  libmesh_error();
341 
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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo