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

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

Hosted By:
SourceForge.net Logo