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