fe_hierarchic.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/elem.h"
22 #include "libmesh/fe.h"
23 #include "libmesh/fe_interface.h"
24 
25 namespace libMesh
26 {
27 
28  // ------------------------------------------------------------
29  // Hierarchic-specific implementations
30 
31  // Anonymous namespace for local helper functions
32  namespace {
33 
34  void hierarchic_nodal_soln(const Elem* elem,
35  const Order order,
36  const std::vector<Number>& elem_soln,
37  std::vector<Number>& nodal_soln,
38  unsigned Dim)
39  {
40  const unsigned int n_nodes = elem->n_nodes();
41 
42  const ElemType elem_type = elem->type();
43 
44  nodal_soln.resize(n_nodes);
45 
46  const Order totalorder = static_cast<Order>(order + elem->p_level());
47 
48  // FEType object to be passed to various FEInterface functions below.
49  FEType fe_type(totalorder, HIERARCHIC);
50 
51  switch (totalorder)
52  {
53  // Constant shape functions
54  case CONSTANT:
55  {
56  libmesh_assert_equal_to (elem_soln.size(), 1);
57 
58  const Number val = elem_soln[0];
59 
60  for (unsigned int n=0; n<n_nodes; n++)
61  nodal_soln[n] = val;
62 
63  return;
64  }
65 
66 
67  // For other orders do interpolation at the nodes
68  // explicitly.
69  default:
70  {
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  }
96  }
97  } // hierarchic_nodal_soln()
98 
99 
100 
101 
102  unsigned int hierarchic_n_dofs(const ElemType t, const Order o)
103  {
104  libmesh_assert_greater (o, 0);
105  switch (t)
106  {
107  case NODEELEM:
108  return 1;
109  case EDGE2:
110  case EDGE3:
111  return (o+1);
112  case QUAD4:
113  libmesh_assert_less (o, 2);
114  case QUAD8:
115  case QUAD9:
116  return ((o+1)*(o+1));
117  case HEX8:
118  libmesh_assert_less (o, 2);
119  case HEX20:
120  libmesh_assert_less (o, 2);
121  case HEX27:
122  return ((o+1)*(o+1)*(o+1));
123  case TRI6:
124  return ((o+1)*(o+2)/2);
125  default:
126  libmesh_error();
127  }
128 
129  libmesh_error();
130  return 0;
131  } // hierarchic_n_dofs()
132 
133 
134 
135 
136  unsigned int hierarchic_n_dofs_at_node(const ElemType t,
137  const Order o,
138  const unsigned int n)
139  {
140  libmesh_assert_greater (o, 0);
141  switch (t)
142  {
143  case NODEELEM:
144  return 1;
145  case EDGE2:
146  case EDGE3:
147  switch (n)
148  {
149  case 0:
150  case 1:
151  return 1;
152  // Internal DoFs are associated with the elem, not its nodes
153  case 2:
154  return 0;
155  default:
156  libmesh_error();
157  }
158  case TRI6:
159  switch (n)
160  {
161  case 0:
162  case 1:
163  case 2:
164  return 1;
165 
166  case 3:
167  case 4:
168  case 5:
169  return (o-1);
170 
171  // Internal DoFs are associated with the elem, not its nodes
172  default:
173  libmesh_error();
174  }
175  case QUAD4:
176  libmesh_assert_less (n, 4);
177  libmesh_assert_less (o, 2);
178  case QUAD8:
179  case QUAD9:
180  switch (n)
181  {
182  case 0:
183  case 1:
184  case 2:
185  case 3:
186  return 1;
187 
188  case 4:
189  case 5:
190  case 6:
191  case 7:
192  return (o-1);
193 
194  // Internal DoFs are associated with the elem, not its nodes
195  case 8:
196  return 0;
197 
198  default:
199  libmesh_error();
200  }
201  case HEX8:
202  libmesh_assert_less (n, 8);
203  libmesh_assert_less (o, 2);
204  case HEX20:
205  libmesh_assert_less (n, 20);
206  libmesh_assert_less (o, 2);
207  case HEX27:
208  switch (n)
209  {
210  case 0:
211  case 1:
212  case 2:
213  case 3:
214  case 4:
215  case 5:
216  case 6:
217  case 7:
218  return 1;
219 
220  case 8:
221  case 9:
222  case 10:
223  case 11:
224  case 12:
225  case 13:
226  case 14:
227  case 15:
228  case 16:
229  case 17:
230  case 18:
231  case 19:
232  return (o-1);
233 
234  case 20:
235  case 21:
236  case 22:
237  case 23:
238  case 24:
239  case 25:
240  return ((o-1)*(o-1));
241 
242  // Internal DoFs are associated with the elem, not its nodes
243  case 26:
244  return 0;
245  default:
246  libmesh_error();
247  }
248  default:
249 #ifdef DEBUG
250  libMesh::err << "ERROR: Bad ElemType = " << t
251  << std::endl;
252 #endif
253  libmesh_error();
254  }
255 
256  libmesh_error();
257 
258  return 0;
259  } // hierarchic_n_dofs_at_node()
260 
261 
262 
263 
264  unsigned int hierarchic_n_dofs_per_elem(const ElemType t,
265  const Order o)
266  {
267  libmesh_assert_greater (o, 0);
268  switch (t)
269  {
270  case NODEELEM:
271  return 0;
272  case EDGE2:
273  case EDGE3:
274  return (o-1);
275  case TRI3:
276  case QUAD4:
277  return 0;
278  case TRI6:
279  return ((o-1)*(o-2)/2);
280  case QUAD8:
281  case QUAD9:
282  return ((o-1)*(o-1));
283  case HEX8:
284  case HEX20:
285  libmesh_assert_less (o, 2);
286  return 0;
287  case HEX27:
288  return ((o-1)*(o-1)*(o-1));
289  default:
290 #ifdef DEBUG
291  libMesh::err << "ERROR: Bad ElemType = " << t
292  << std::endl;
293 #endif
294  libmesh_error();
295  }
296 
297  // Will never get here...
298  libmesh_error();
299  return 0;
300  } // hierarchic_n_dofs_per_elem()
301 
302  } // anonymous namespace
303 
304 
305 
306 
307  // Do full-specialization of nodal_soln() function for every
308  // dimension, instead of explicit instantiation at the end of this
309  // file.
310  // This could be macro-ified so that it fits on one line...
311  template <>
313  const Order order,
314  const std::vector<Number>& elem_soln,
315  std::vector<Number>& nodal_soln)
316  { hierarchic_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/0); }
317 
318  template <>
320  const Order order,
321  const std::vector<Number>& elem_soln,
322  std::vector<Number>& nodal_soln)
323  { hierarchic_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/1); }
324 
325  template <>
327  const Order order,
328  const std::vector<Number>& elem_soln,
329  std::vector<Number>& nodal_soln)
330  { hierarchic_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/2); }
331 
332  template <>
334  const Order order,
335  const std::vector<Number>& elem_soln,
336  std::vector<Number>& nodal_soln)
337  { hierarchic_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/3); }
338 
339 
340  // Full specialization of n_dofs() function for every dimension
341  template <> unsigned int FE<0,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
342  template <> unsigned int FE<1,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
343  template <> unsigned int FE<2,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
344  template <> unsigned int FE<3,HIERARCHIC>::n_dofs(const ElemType t, const Order o) { return hierarchic_n_dofs(t, o); }
345 
346  // Full specialization of n_dofs_at_node() function for every dimension.
347  template <> unsigned int FE<0,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
348  template <> unsigned int FE<1,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
349  template <> unsigned int FE<2,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
350  template <> unsigned int FE<3,HIERARCHIC>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return hierarchic_n_dofs_at_node(t, o, n); }
351 
352  // Full specialization of n_dofs_per_elem() function for every dimension.
353  template <> unsigned int FE<0,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
354  template <> unsigned int FE<1,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
355  template <> unsigned int FE<2,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
356  template <> unsigned int FE<3,HIERARCHIC>::n_dofs_per_elem(const ElemType t, const Order o) { return hierarchic_n_dofs_per_elem(t, o); }
357 
358  // Hierarchic FEMs are C^0 continuous
359  template <> FEContinuity FE<0,HIERARCHIC>::get_continuity() const { return C_ZERO; }
360  template <> FEContinuity FE<1,HIERARCHIC>::get_continuity() const { return C_ZERO; }
361  template <> FEContinuity FE<2,HIERARCHIC>::get_continuity() const { return C_ZERO; }
362  template <> FEContinuity FE<3,HIERARCHIC>::get_continuity() const { return C_ZERO; }
363 
364  // Hierarchic FEMs are hierarchic (duh!)
365  template <> bool FE<0,HIERARCHIC>::is_hierarchic() const { return true; }
366  template <> bool FE<1,HIERARCHIC>::is_hierarchic() const { return true; }
367  template <> bool FE<2,HIERARCHIC>::is_hierarchic() const { return true; }
368  template <> bool FE<3,HIERARCHIC>::is_hierarchic() const { return true; }
369 
370 #ifdef LIBMESH_ENABLE_AMR
371  // compute_constraints() specializations are only needed for 2 and 3D
372  template <>
374  DofMap &dof_map,
375  const unsigned int variable_number,
376  const Elem* elem)
377  { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
378 
379  template <>
381  DofMap &dof_map,
382  const unsigned int variable_number,
383  const Elem* elem)
384  { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
385 #endif // #ifdef LIBMESH_ENABLE_AMR
386 
387  // Hierarchic FEM shapes need reinit
388  template <> bool FE<0,HIERARCHIC>::shapes_need_reinit() const { return true; }
389  template <> bool FE<1,HIERARCHIC>::shapes_need_reinit() const { return true; }
390  template <> bool FE<2,HIERARCHIC>::shapes_need_reinit() const { return true; }
391  template <> bool FE<3,HIERARCHIC>::shapes_need_reinit() const { return true; }
392 
393 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo