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/fe.h"
22 #include "libmesh/elem.h"
23 #include "libmesh/fe_interface.h"
24 #include "libmesh/enum_to_string.h"
25
26 namespace libMesh
27 {
28
29 // Anonymous namespace for local helper functions
30 namespace {
31
monomial_nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)32 void monomial_nodal_soln(const Elem * elem,
33 const Order order,
34 const std::vector<Number> & elem_soln,
35 std::vector<Number> & nodal_soln)
36 {
37 const unsigned int n_nodes = elem->n_nodes();
38
39 const ElemType elem_type = elem->type();
40
41 nodal_soln.resize(n_nodes);
42
43 const Order totalorder = static_cast<Order>(order+elem->p_level());
44
45 switch (totalorder)
46 {
47 // Constant shape functions
48 case CONSTANT:
49 {
50 libmesh_assert_equal_to (elem_soln.size(), 1);
51
52 const Number val = elem_soln[0];
53
54 for (unsigned int n=0; n<n_nodes; n++)
55 nodal_soln[n] = val;
56
57 return;
58 }
59
60
61 // For other orders, do interpolation at the nodes
62 // explicitly.
63 default:
64 {
65 // FEType object to be passed to various FEInterface functions below.
66 FEType fe_type(order, MONOMIAL);
67
68 const unsigned int n_sf =
69 FEInterface::n_shape_functions(fe_type, elem);
70
71 std::vector<Point> refspace_nodes;
72 FEBase::get_refspace_nodes(elem_type,refspace_nodes);
73 libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
74
75 for (unsigned int n=0; n<n_nodes; n++)
76 {
77 libmesh_assert_equal_to (elem_soln.size(), n_sf);
78
79 // Zero before summation
80 nodal_soln[n] = 0;
81
82 // u_i = Sum (alpha_i phi_i)
83 for (unsigned int i=0; i<n_sf; i++)
84 nodal_soln[n] += elem_soln[i] *
85 FEInterface::shape(fe_type, elem, i, refspace_nodes[n]);
86 }
87
88 return;
89 } // default
90 } // switch
91 } // monomial_nodal_soln()
92
93
94
95
monomial_n_dofs(const ElemType t,const Order o)96 unsigned int monomial_n_dofs(const ElemType t, const Order o)
97 {
98 switch (o)
99 {
100
101 // constant shape functions
102 // no matter what shape there is only one DOF.
103 case CONSTANT:
104 return (t != INVALID_ELEM) ? 1 : 0;
105
106
107 // Discontinuous linear shape functions
108 // expressed in the monomials.
109 case FIRST:
110 {
111 switch (t)
112 {
113 case NODEELEM:
114 return 1;
115
116 case EDGE2:
117 case EDGE3:
118 case EDGE4:
119 return 2;
120
121 case TRI3:
122 case TRISHELL3:
123 case TRI6:
124 case QUAD4:
125 case QUADSHELL4:
126 case QUAD8:
127 case QUADSHELL8:
128 case QUAD9:
129 return 3;
130
131 case TET4:
132 case TET10:
133 case HEX8:
134 case HEX20:
135 case HEX27:
136 case PRISM6:
137 case PRISM15:
138 case PRISM18:
139 case PYRAMID5:
140 case PYRAMID13:
141 case PYRAMID14:
142 return 4;
143
144 case INVALID_ELEM:
145 return 0;
146
147 default:
148 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
149 }
150 }
151
152
153 // Discontinuous quadratic shape functions
154 // expressed in the monomials.
155 case SECOND:
156 {
157 switch (t)
158 {
159 case NODEELEM:
160 return 1;
161
162 case EDGE2:
163 case EDGE3:
164 case EDGE4:
165 return 3;
166
167 case TRI3:
168 case TRISHELL3:
169 case TRI6:
170 case QUAD4:
171 case QUADSHELL4:
172 case QUAD8:
173 case QUADSHELL8:
174 case QUAD9:
175 return 6;
176
177 case TET4:
178 case TET10:
179 case HEX8:
180 case HEX20:
181 case HEX27:
182 case PRISM6:
183 case PRISM15:
184 case PRISM18:
185 case PYRAMID5:
186 case PYRAMID13:
187 case PYRAMID14:
188 return 10;
189
190 case INVALID_ELEM:
191 return 0;
192
193 default:
194 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
195 }
196 }
197
198
199 // Discontinuous cubic shape functions
200 // expressed in the monomials.
201 case THIRD:
202 {
203 switch (t)
204 {
205 case NODEELEM:
206 return 1;
207
208 case EDGE2:
209 case EDGE3:
210 case EDGE4:
211 return 4;
212
213 case TRI3:
214 case TRISHELL3:
215 case TRI6:
216 case QUAD4:
217 case QUADSHELL4:
218 case QUAD8:
219 case QUADSHELL8:
220 case QUAD9:
221 return 10;
222
223 case TET4:
224 case TET10:
225 case HEX8:
226 case HEX20:
227 case HEX27:
228 case PRISM6:
229 case PRISM15:
230 case PRISM18:
231 case PYRAMID5:
232 case PYRAMID13:
233 case PYRAMID14:
234 return 20;
235
236 case INVALID_ELEM:
237 return 0;
238
239 default:
240 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
241 }
242 }
243
244
245 // Discontinuous quartic shape functions
246 // expressed in the monomials.
247 case FOURTH:
248 {
249 switch (t)
250 {
251 case NODEELEM:
252 return 1;
253
254 case EDGE2:
255 case EDGE3:
256 return 5;
257
258 case TRI3:
259 case TRISHELL3:
260 case TRI6:
261 case QUAD4:
262 case QUADSHELL4:
263 case QUAD8:
264 case QUADSHELL8:
265 case QUAD9:
266 return 15;
267
268 case TET4:
269 case TET10:
270 case HEX8:
271 case HEX20:
272 case HEX27:
273 case PRISM6:
274 case PRISM15:
275 case PRISM18:
276 case PYRAMID5:
277 case PYRAMID13:
278 case PYRAMID14:
279 return 35;
280
281 case INVALID_ELEM:
282 return 0;
283
284 default:
285 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
286 }
287 }
288
289
290 default:
291 {
292 const unsigned int order = static_cast<unsigned int>(o);
293 switch (t)
294 {
295 case NODEELEM:
296 return 1;
297
298 case EDGE2:
299 case EDGE3:
300 return (order+1);
301
302 case TRI3:
303 case TRISHELL3:
304 case TRI6:
305 case QUAD4:
306 case QUADSHELL4:
307 case QUAD8:
308 case QUADSHELL8:
309 case QUAD9:
310 return (order+1)*(order+2)/2;
311
312 case TET4:
313 case TET10:
314 case HEX8:
315 case HEX20:
316 case HEX27:
317 case PRISM6:
318 case PRISM15:
319 case PRISM18:
320 case PYRAMID5:
321 case PYRAMID13:
322 case PYRAMID14:
323 return (order+1)*(order+2)*(order+3)/6;
324
325 case INVALID_ELEM:
326 return 0;
327
328 default:
329 libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
330 }
331 }
332 }
333 } // monomial_n_dofs()
334
335
336 } // anonymous namespace
337
338
339
340
341
342 // Do full-specialization for every dimension, instead
343 // of explicit instantiation at the end of this file.
344 // This could be macro-ified so that it fits on one line...
345 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)346 void FE<0,MONOMIAL>::nodal_soln(const Elem * elem,
347 const Order order,
348 const std::vector<Number> & elem_soln,
349 std::vector<Number> & nodal_soln)
350 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln); }
351
352 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)353 void FE<1,MONOMIAL>::nodal_soln(const Elem * elem,
354 const Order order,
355 const std::vector<Number> & elem_soln,
356 std::vector<Number> & nodal_soln)
357 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln); }
358
359 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)360 void FE<2,MONOMIAL>::nodal_soln(const Elem * elem,
361 const Order order,
362 const std::vector<Number> & elem_soln,
363 std::vector<Number> & nodal_soln)
364 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln); }
365
366 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)367 void FE<3,MONOMIAL>::nodal_soln(const Elem * elem,
368 const Order order,
369 const std::vector<Number> & elem_soln,
370 std::vector<Number> & nodal_soln)
371 { monomial_nodal_soln(elem, order, elem_soln, nodal_soln); }
372
373
374 // Full specialization of n_dofs() function for every dimension
n_dofs(const ElemType t,const Order o)375 template <> unsigned int FE<0,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
n_dofs(const ElemType t,const Order o)376 template <> unsigned int FE<1,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
n_dofs(const ElemType t,const Order o)377 template <> unsigned int FE<2,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
n_dofs(const ElemType t,const Order o)378 template <> unsigned int FE<3,MONOMIAL>::n_dofs(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
379
380 // Full specialization of n_dofs_at_node() function for every dimension.
381 // Monomials have no dofs at nodes, only element dofs.
n_dofs_at_node(const ElemType,const Order,const unsigned int)382 template <> unsigned int FE<0,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
n_dofs_at_node(const ElemType,const Order,const unsigned int)383 template <> unsigned int FE<1,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
n_dofs_at_node(const ElemType,const Order,const unsigned int)384 template <> unsigned int FE<2,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
n_dofs_at_node(const ElemType,const Order,const unsigned int)385 template <> unsigned int FE<3,MONOMIAL>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
386
387 // Full specialization of n_dofs_per_elem() function for every dimension.
n_dofs_per_elem(const ElemType t,const Order o)388 template <> unsigned int FE<0,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
n_dofs_per_elem(const ElemType t,const Order o)389 template <> unsigned int FE<1,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
n_dofs_per_elem(const ElemType t,const Order o)390 template <> unsigned int FE<2,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
n_dofs_per_elem(const ElemType t,const Order o)391 template <> unsigned int FE<3,MONOMIAL>::n_dofs_per_elem(const ElemType t, const Order o) { return monomial_n_dofs(t, o); }
392
393
394 // Full specialization of get_continuity() function for every dimension.
get_continuity()395 template <> FEContinuity FE<0,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
get_continuity()396 template <> FEContinuity FE<1,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
get_continuity()397 template <> FEContinuity FE<2,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
get_continuity()398 template <> FEContinuity FE<3,MONOMIAL>::get_continuity() const { return DISCONTINUOUS; }
399
400 // Full specialization of is_hierarchic() function for every dimension.
401 // The monomials are hierarchic!
is_hierarchic()402 template <> bool FE<0,MONOMIAL>::is_hierarchic() const { return true; }
is_hierarchic()403 template <> bool FE<1,MONOMIAL>::is_hierarchic() const { return true; }
is_hierarchic()404 template <> bool FE<2,MONOMIAL>::is_hierarchic() const { return true; }
is_hierarchic()405 template <> bool FE<3,MONOMIAL>::is_hierarchic() const { return true; }
406
407 #ifdef LIBMESH_ENABLE_AMR
408
409 // Full specialization of compute_constraints() function for 2D and
410 // 3D only. There are no constraints for discontinuous elements, so
411 // we do nothing.
compute_constraints(DofConstraints &,DofMap &,const unsigned int,const Elem *)412 template <> void FE<2,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
compute_constraints(DofConstraints &,DofMap &,const unsigned int,const Elem *)413 template <> void FE<3,MONOMIAL>::compute_constraints (DofConstraints &, DofMap &, const unsigned int, const Elem *) {}
414
415 #endif // #ifdef LIBMESH_ENABLE_AMR
416
417 // Full specialization of shapes_need_reinit() function for every dimension.
shapes_need_reinit()418 template <> bool FE<0,MONOMIAL>::shapes_need_reinit() const { return false; }
shapes_need_reinit()419 template <> bool FE<1,MONOMIAL>::shapes_need_reinit() const { return false; }
shapes_need_reinit()420 template <> bool FE<2,MONOMIAL>::shapes_need_reinit() const { return false; }
shapes_need_reinit()421 template <> bool FE<3,MONOMIAL>::shapes_need_reinit() const { return false; }
422
423 } // namespace libMesh
424