libMesh
fe_clough.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/elem.h"
22 #include "libmesh/fe.h"
23 #include "libmesh/fe_interface.h"
24 #include "libmesh/string_to_enum.h"
25 
26 
27 namespace libMesh
28 {
29 
30 // ------------------------------------------------------------
31 // Clough-specific implementations
32 
33 // Anonymous namespace for local helper functions
34 namespace {
35 
36 void clough_nodal_soln(const Elem * elem,
37  const Order order,
38  const std::vector<Number> & elem_soln,
39  std::vector<Number> & nodal_soln,
40  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  // FEType object to be passed to various FEInterface functions below.
51  FEType fe_type(totalorder, CLOUGH);
52 
53  switch (totalorder)
54  {
55  // Piecewise cubic shape functions with linear flux edges
56  case SECOND:
57  // Piecewise cubic shape functions
58  case THIRD:
59  {
60 
61  const unsigned int n_sf =
62  // FE<Dim,T>::n_shape_functions(elem_type, totalorder);
63  FEInterface::n_shape_functions(Dim, fe_type, elem_type);
64 
65  std::vector<Point> refspace_nodes;
66  FEBase::get_refspace_nodes(elem_type,refspace_nodes);
67  libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
68 
69  for (unsigned int n=0; n<n_nodes; n++)
70  {
71  libmesh_assert_equal_to (elem_soln.size(), n_sf);
72 
73  // Zero before summation
74  nodal_soln[n] = 0;
75 
76  // u_i = Sum (alpha_i phi_i)
77  for (unsigned int i=0; i<n_sf; i++)
78  nodal_soln[n] += elem_soln[i] *
79  // FE<Dim,T>::shape(elem, order, i, mapped_point);
80  FEInterface::shape(Dim, fe_type, elem, i, refspace_nodes[n]);
81  }
82 
83  return;
84  }
85 
86  default:
87  libmesh_error_msg("ERROR: Invalid total order " << totalorder);
88  }
89 } // clough_nodal_soln()
90 
91 
92 
93 
94 unsigned int clough_n_dofs(const ElemType t, const Order o)
95 {
96  if (t == INVALID_ELEM)
97  return 0;
98 
99  switch (o)
100  {
101  // Piecewise cubic shape functions with linear flux edges
102  case SECOND:
103  {
104  switch (t)
105  {
106  case TRI6:
107  return 9;
108 
109  default:
110  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
111  }
112  }
113  // Piecewise cubic Clough-Tocher element
114  case THIRD:
115  {
116  switch (t)
117  {
118  case TRI6:
119  return 12;
120 
121  default:
122  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
123  }
124  }
125 
126  default:
127  libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for CLOUGH FE family!");
128  }
129 
130  libmesh_error_msg("We'll never get here!");
131  return 0;
132 } // clough_n_dofs()
133 
134 
135 
136 
137 
138 unsigned int clough_n_dofs_at_node(const ElemType t,
139  const Order o,
140  const unsigned int n)
141 {
142  switch (o)
143  {
144  // Piecewise cubic shape functions with linear flux edges
145  case SECOND:
146  {
147  switch (t)
148  {
149  // The 2D Clough-Tocher defined on a 6-noded triangle
150  case TRI6:
151  {
152  switch (n)
153  {
154  case 0:
155  case 1:
156  case 2:
157  return 3;
158 
159  case 3:
160  case 4:
161  case 5:
162  return 0;
163 
164  default:
165  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
166  }
167  }
168 
169  default:
170  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
171  }
172  }
173  // The third-order clough shape functions
174  case THIRD:
175  {
176  switch (t)
177  {
178  // The 2D Clough-Tocher defined on a 6-noded triangle
179  case TRI6:
180  {
181  switch (n)
182  {
183  case 0:
184  case 1:
185  case 2:
186  return 3;
187 
188  case 3:
189  case 4:
190  case 5:
191  return 1;
192 
193  default:
194  libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
195  }
196  }
197 
198  default:
199  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
200  }
201  }
202  default:
203  libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for CLOUGH FE family!");
204  }
205 
206  libmesh_error_msg("We'll never get here!");
207  return 0;
208 } // clough_n_dofs_at_node()
209 
210 
211 
212 
213 unsigned int clough_n_dofs_per_elem(const ElemType t, const Order o)
214 {
215  switch (o)
216  {
217  // Piecewise cubic shape functions with linear flux edges
218  case SECOND:
219  // The third-order Clough-Tocher shape functions
220  case THIRD:
221  {
222  switch (t)
223  {
224  // The 2D clough defined on a 6-noded triangle
225  case TRI6:
226  return 0;
227 
228  default:
229  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
230  }
231  }
232 
233  default:
234  libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for CLOUGH FE family!");
235  }
236 
237  libmesh_error_msg("We'll never get here!");
238  return 0;
239 } // clough_n_dofs_per_elem()
240 
241 } // anonymous
242 
243 
244 
245 
246  // Do full-specialization of nodal_soln() function for every
247  // dimension, instead of explicit instantiation at the end of this
248  // file.
249  // This could be macro-ified so that it fits on one line...
250 template <>
251 void FE<0,CLOUGH>::nodal_soln(const Elem * elem,
252  const Order order,
253  const std::vector<Number> & elem_soln,
254  std::vector<Number> & nodal_soln)
255 { clough_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/0); }
256 
257 template <>
258 void FE<1,CLOUGH>::nodal_soln(const Elem * elem,
259  const Order order,
260  const std::vector<Number> & elem_soln,
261  std::vector<Number> & nodal_soln)
262 { clough_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/1); }
263 
264 template <>
265 void FE<2,CLOUGH>::nodal_soln(const Elem * elem,
266  const Order order,
267  const std::vector<Number> & elem_soln,
268  std::vector<Number> & nodal_soln)
269 { clough_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/2); }
270 
271 template <>
272 void FE<3,CLOUGH>::nodal_soln(const Elem * elem,
273  const Order order,
274  const std::vector<Number> & elem_soln,
275  std::vector<Number> & nodal_soln)
276 { clough_nodal_soln(elem, order, elem_soln, nodal_soln, /*Dim=*/3); }
277 
278 
279 // Full specialization of n_dofs() function for every dimension
280 template <> unsigned int FE<0,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
281 template <> unsigned int FE<1,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
282 template <> unsigned int FE<2,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
283 template <> unsigned int FE<3,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
284 
285 
286 // Full specialization of n_dofs_at_node() function for every dimension.
287 template <> unsigned int FE<0,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
288 template <> unsigned int FE<1,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
289 template <> unsigned int FE<2,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
290 template <> unsigned int FE<3,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
291 
292 // Full specialization of n_dofs_per_elem() function for every dimension.
293 template <> unsigned int FE<0,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
294 template <> unsigned int FE<1,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
295 template <> unsigned int FE<2,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
296 template <> unsigned int FE<3,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
297 
298 // Clough FEMs are C^1 continuous
299 template <> FEContinuity FE<0,CLOUGH>::get_continuity() const { return C_ONE; }
300 template <> FEContinuity FE<1,CLOUGH>::get_continuity() const { return C_ONE; }
301 template <> FEContinuity FE<2,CLOUGH>::get_continuity() const { return C_ONE; }
302 template <> FEContinuity FE<3,CLOUGH>::get_continuity() const { return C_ONE; }
303 
304 // Clough FEMs are not (currently) hierarchic
305 template <> bool FE<0,CLOUGH>::is_hierarchic() const { return false; } // FIXME - this will be changed
306 template <> bool FE<1,CLOUGH>::is_hierarchic() const { return false; } // FIXME - this will be changed
307 template <> bool FE<2,CLOUGH>::is_hierarchic() const { return false; } // FIXME - this will be changed
308 template <> bool FE<3,CLOUGH>::is_hierarchic() const { return false; } // FIXME - this will be changed
309 
310 #ifdef LIBMESH_ENABLE_AMR
311 // compute_constraints() specializations are only needed for 2 and 3D
312 template <>
314  DofMap & dof_map,
315  const unsigned int variable_number,
316  const Elem * elem)
317 { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
318 
319 template <>
321  DofMap & dof_map,
322  const unsigned int variable_number,
323  const Elem * elem)
324 { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
325 #endif // #ifdef LIBMESH_ENABLE_AMR
326 
327 // Clough FEM shapes need reinit
328 template <> bool FE<0,CLOUGH>::shapes_need_reinit() const { return true; }
329 template <> bool FE<1,CLOUGH>::shapes_need_reinit() const { return true; }
330 template <> bool FE<2,CLOUGH>::shapes_need_reinit() const { return true; }
331 template <> bool FE<3,CLOUGH>::shapes_need_reinit() const { return true; }
332 
333 } // 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