libMesh
Loading...
Searching...
No Matches
fe_nedelec_one.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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/dof_map.h"
22#include "libmesh/elem.h"
23#include "libmesh/enum_to_string.h"
24#include "libmesh/fe.h"
25#include "libmesh/fe_interface.h"
26#include "libmesh/fe_macro.h"
27#include "libmesh/tensor_value.h"
28
29
30namespace libMesh
31{
32
33
38
39
40// Anonymous namespace for local helper functions
41namespace {
42void nedelec_one_nodal_soln(const Elem * elem,
43 const Order order,
44 const std::vector<Number> & elem_soln,
45 const int dim,
46 const int vdim,
47 std::vector<Number> & nodal_soln,
48 const bool add_p_level)
49{
50 const unsigned int n_nodes = elem->n_nodes();
51 const ElemType elem_type = elem->type();
52
53 const Order totalorder = order + add_p_level*elem->p_level();
54
55 nodal_soln.resize(n_nodes*vdim);
56
57 FEType p_refined_fe_type(totalorder, NEDELEC_ONE);
58
59 if (elem_type != TRI6 && elem_type != TRI7 &&
60 elem_type != QUAD8 && elem_type != QUAD9 &&
61 elem_type != TET10 && elem_type != TET14 &&
62 elem_type != HEX20 && elem_type != HEX27)
63 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(elem_type) << " selected for NEDELEC_ONE FE family!");
64
65 const unsigned int n_sf = FEInterface::n_shape_functions(p_refined_fe_type, elem, false);
66
67 std::vector<Point> refspace_nodes;
68 FEVectorBase::get_refspace_nodes(elem_type,refspace_nodes);
69 libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
70 libmesh_assert_equal_to (elem_soln.size(), n_sf);
71
72 // Need to create new fe object so the shape function has the FETransformation
73 // applied to it.
74 std::unique_ptr<FEVectorBase> vis_fe = FEVectorBase::build(dim, p_refined_fe_type);
75
76 const std::vector<std::vector<RealGradient>> & vis_phi = vis_fe->get_phi();
77
78 vis_fe->reinit(elem,&refspace_nodes);
79
80 // Zero before summation
81 std::fill(nodal_soln.begin(), nodal_soln.end(), 0);
82
83 for (unsigned int n = 0; n < n_nodes; n++)
84 // u = Sum (u_i phi_i)
85 for (unsigned int i=0; i<n_sf; i++)
86 for (int d = 0; d < vdim; d++)
87 nodal_soln[vdim*n+d] += elem_soln[i]*(vis_phi[i][n](d));
88
89 return;
90} // nedelec_one_nodal_soln
91
92
93
94unsigned int nedelec_one_n_dofs(const ElemType t, const Order o)
95{
96 libmesh_assert_greater (o, 0);
97 switch (t)
98 {
99 case TRI6:
100 case TRI7:
101 return o*(o+2);
102 case QUAD8:
103 case QUAD9:
104 return 2*o*(o+1);
105 case TET10:
106 libmesh_assert_less (o, 2);
107 libmesh_fallthrough();
108 case TET14:
109 return o*(o+2)*(o+3)/2;
110 case HEX20:
111 libmesh_assert_less (o, 2);
112 libmesh_fallthrough();
113 case HEX27:
114 return 3*o*(o+1)*(o+1);
115 case INVALID_ELEM:
116 return 0;
117 default:
118 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for NEDELEC_ONE FE family!");
119 }
120}
121
122
123
124unsigned int nedelec_one_n_dofs(const Elem * e, const Order o)
125{
127 return nedelec_one_n_dofs(e->type(), o);
128}
129
130
131
132unsigned int nedelec_one_n_dofs_at_node(const ElemType t,
133 const Order o,
134 const unsigned int n)
135{
136 libmesh_assert_greater (o, 0);
137 switch (t)
138 {
139 case TRI6:
140 case TRI7:
141 {
142 switch (n)
143 {
144 case 0:
145 case 1:
146 case 2:
147 return 0;
148 case 3:
149 case 4:
150 case 5:
151 return o;
152 case 6:
153 libmesh_assert_equal_to(t, TRI7);
154 return 0;
155 default:
156 libmesh_error_msg("ERROR: Invalid node ID " << n);
157 }
158 }
159 case QUAD8:
160 case QUAD9:
161 {
162 switch (n)
163 {
164 case 0:
165 case 1:
166 case 2:
167 case 3:
168 return 0;
169 case 4:
170 case 5:
171 case 6:
172 case 7:
173 return o;
174 case 8:
175 libmesh_assert_equal_to(t, QUAD9);
176 return 0;
177 default:
178 libmesh_error_msg("ERROR: Invalid node ID " << n);
179 }
180 }
181 case TET10:
182 libmesh_assert_less (o, 2);
183 libmesh_fallthrough();
184 case TET14:
185 {
186 switch (n)
187 {
188 case 0:
189 case 1:
190 case 2:
191 case 3:
192 return 0;
193 case 4:
194 case 5:
195 case 6:
196 case 7:
197 case 8:
198 case 9:
199 return o;
200 case 10:
201 case 11:
202 case 12:
203 case 13:
204 libmesh_assert_equal_to(t, TET14);
205 return o*(o-1);
206 default:
207 libmesh_error_msg("ERROR: Invalid node ID " << n);
208 }
209 }
210 case HEX20:
211 libmesh_assert_less (o, 2);
212 libmesh_fallthrough();
213 case HEX27:
214 {
215 switch (n)
216 {
217 case 0:
218 case 1:
219 case 2:
220 case 3:
221 case 4:
222 case 5:
223 case 6:
224 case 7:
225 return 0;
226 case 8:
227 case 9:
228 case 10:
229 case 11:
230 case 12:
231 case 13:
232 case 14:
233 case 15:
234 case 16:
235 case 17:
236 case 18:
237 case 19:
238 return o;
239 case 20:
240 case 21:
241 case 22:
242 case 23:
243 case 24:
244 case 25:
245 libmesh_assert_equal_to(t, HEX27);
246 return 2*o*(o-1);
247 case 26:
248 libmesh_assert_equal_to(t, HEX27);
249 return 0;
250 default:
251 libmesh_error_msg("ERROR: Invalid node ID " << n);
252 }
253 }
254 case INVALID_ELEM:
255 return 0;
256 default:
257 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for NEDELEC_ONE FE family!");
258 }
259}
260
261
262
263unsigned int nedelec_one_n_dofs_at_node(const Elem & e,
264 const Order o,
265 const unsigned int n)
266{
267 return nedelec_one_n_dofs_at_node(e.type(), o, n);
268}
269
270
271
272unsigned int nedelec_one_n_dofs_per_elem(const ElemType t,
273 const Order o)
274{
275 libmesh_assert_greater (o, 0);
276 switch (t)
277 {
278 case TRI6:
279 case TRI7:
280 return o*(o-1);
281 case QUAD8:
282 case QUAD9:
283 return 2*o*(o-1);
284 case TET10:
285 libmesh_assert_less (o, 2);
286 libmesh_fallthrough();
287 case TET14:
288 return o*(o-1)*(o-2)/2;
289 case HEX20:
290 libmesh_assert_less (o, 2);
291 libmesh_fallthrough();
292 case HEX27:
293 return 3*o*(o-1)*(o-1);
294 case INVALID_ELEM:
295 return 0;
296 default:
297 libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(t) << " selected for NEDELEC_ONE FE family!");
298 }
299}
300
301
302
303unsigned int nedelec_one_n_dofs_per_elem(const Elem & e,
304 const Order o)
305{
306 return nedelec_one_n_dofs_per_elem(e.type(), o);
307}
308
309
310
311#ifdef LIBMESH_ENABLE_AMR
312void nedelec_one_compute_constraints (DofConstraints & /*constraints*/,
313 DofMap & /*dof_map*/,
314 const unsigned int /*variable_number*/,
315 const Elem * libmesh_dbg_var(elem),
316 const unsigned Dim)
317{
318 // Only constrain elements in 2,3D.
319 if (Dim == 1)
320 return;
321
322 libmesh_assert(elem);
323
324 libmesh_not_implemented();
325} // nedelec_one_compute_constraints()
326#endif // #ifdef LIBMESH_ENABLE_AMR
327
328} // anonymous namespace
329
330#define NEDELEC_LOW_D_ERROR_MESSAGE \
331 libmesh_error_msg("ERROR: This method makes no sense for low-D elements!");
332
333
334// Do full-specialization for every dimension, instead
335// of explicit instantiation at the end of this file.
336template <>
338 const Order,
339 const std::vector<Number> &,
340 std::vector<Number> &,
341 bool,
342 const unsigned)
343{ NEDELEC_LOW_D_ERROR_MESSAGE }
344
345template <>
347 const Order,
348 const std::vector<Number> &,
349 std::vector<Number> &,
350 bool,
351 const unsigned)
352{ NEDELEC_LOW_D_ERROR_MESSAGE }
353
354template <>
356 const Order order,
357 const std::vector<Number> & elem_soln,
358 std::vector<Number> & nodal_soln,
359 const bool add_p_level,
360 const unsigned vdim)
361{ nedelec_one_nodal_soln(elem, order, elem_soln, 2 /*dim*/, vdim, nodal_soln, add_p_level); }
362
363template <>
365 const Order order,
366 const std::vector<Number> & elem_soln,
367 std::vector<Number> & nodal_soln,
368 const bool add_p_level,
369 const unsigned)
370{ nedelec_one_nodal_soln(elem, order, elem_soln, 3 /*dim*/, 3 /*vdim*/, nodal_soln, add_p_level); }
371
373
374
375// Do full-specialization for every dimension, instead
376// of explicit instantiation at the end of this function.
377template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
378template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
379template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); }
380template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); }
381
382template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs(const Elem *, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
383template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs(const Elem *, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
384template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs(const Elem * e, const Order o) { return nedelec_one_n_dofs(e, o); }
385template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs(const Elem * e, const Order o) { return nedelec_one_n_dofs(e, o); }
386
387template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
388template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
389template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(t, o, n); }
390template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(t, o, n); }
391
392template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
393template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_at_node(const Elem &, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
394template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(e, o, n); }
395template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(e, o, n); }
396
397template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
398template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
399template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_per_elem(const ElemType t, const Order o) { return nedelec_one_n_dofs_per_elem(t, o); }
400template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_per_elem(const ElemType t, const Order o) { return nedelec_one_n_dofs_per_elem(t, o); }
401
402template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_per_elem(const Elem &, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
403template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_per_elem(const Elem &, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
404template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_per_elem(const Elem & e, const Order o) { return nedelec_one_n_dofs_per_elem(e, o); }
405template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_per_elem(const Elem & e, const Order o) { return nedelec_one_n_dofs_per_elem(e, o); }
406
407// Nedelec first type FEMs are always tangentially continuous
408template <> FEContinuity FE<0,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE }
409template <> FEContinuity FE<1,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE }
412
413// Nedelec first type FEMs are not hierarchic
414template <> bool FE<0,NEDELEC_ONE>::is_hierarchic() const { NEDELEC_LOW_D_ERROR_MESSAGE }
415template <> bool FE<1,NEDELEC_ONE>::is_hierarchic() const { NEDELEC_LOW_D_ERROR_MESSAGE }
416template <> bool FE<2,NEDELEC_ONE>::is_hierarchic() const { return false; }
417template <> bool FE<3,NEDELEC_ONE>::is_hierarchic() const { return false; }
418
419// Nedelec first type FEM shapes always need to be reinit'ed (because of orientation dependence)
420template <> bool FE<0,NEDELEC_ONE>::shapes_need_reinit() const { NEDELEC_LOW_D_ERROR_MESSAGE }
421template <> bool FE<1,NEDELEC_ONE>::shapes_need_reinit() const { NEDELEC_LOW_D_ERROR_MESSAGE }
422template <> bool FE<2,NEDELEC_ONE>::shapes_need_reinit() const { return true; }
423template <> bool FE<3,NEDELEC_ONE>::shapes_need_reinit() const { return true; }
424
425#ifdef LIBMESH_ENABLE_AMR
426template <>
428 DofMap &,
429 const unsigned int,
430 const Elem *)
431{ NEDELEC_LOW_D_ERROR_MESSAGE }
432
433template <>
435 DofMap &,
436 const unsigned int,
437 const Elem *)
438{ NEDELEC_LOW_D_ERROR_MESSAGE }
439
440template <>
442 DofMap & dof_map,
443 const unsigned int variable_number,
444 const Elem * elem)
445{ nedelec_one_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
446
447template <>
449 DofMap & dof_map,
450 const unsigned int variable_number,
451 const Elem * elem)
452{ nedelec_one_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
453#endif // LIBMESH_ENABLE_AMR
454
455// Specialize useless shape function methods
456template <>
457RealGradient FE<0,NEDELEC_ONE>::shape(const ElemType, const Order,const unsigned int,const Point &)
458{ NEDELEC_LOW_D_ERROR_MESSAGE }
459template <>
460RealGradient FE<0,NEDELEC_ONE>::shape(const Elem *,const Order,const unsigned int,const Point &,const bool)
461{ NEDELEC_LOW_D_ERROR_MESSAGE }
462template <>
464 const unsigned int,const Point &)
465{ NEDELEC_LOW_D_ERROR_MESSAGE }
466template <>
467RealGradient FE<0,NEDELEC_ONE>::shape_deriv(const Elem *,const Order,const unsigned int,
468 const unsigned int,const Point &,const bool)
469{ NEDELEC_LOW_D_ERROR_MESSAGE }
470
471#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
472template <>
474 const unsigned int,const Point &)
475{ NEDELEC_LOW_D_ERROR_MESSAGE }
476template <>
478 const unsigned int,const Point &,const bool)
479{ NEDELEC_LOW_D_ERROR_MESSAGE }
480
481#endif
482
483template <>
484RealGradient FE<1,NEDELEC_ONE>::shape(const ElemType, const Order,const unsigned int,const Point &)
485{ NEDELEC_LOW_D_ERROR_MESSAGE }
486template <>
487RealGradient FE<1,NEDELEC_ONE>::shape(const Elem *,const Order,const unsigned int,const Point &,const bool)
488{ NEDELEC_LOW_D_ERROR_MESSAGE }
489template <>
491 const unsigned int,const Point &)
492{ NEDELEC_LOW_D_ERROR_MESSAGE }
493template <>
494RealGradient FE<1,NEDELEC_ONE>::shape_deriv(const Elem *,const Order,const unsigned int,
495 const unsigned int,const Point &,const bool)
496{ NEDELEC_LOW_D_ERROR_MESSAGE }
497
498#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
499template <>
501 const unsigned int,const Point &)
502{ NEDELEC_LOW_D_ERROR_MESSAGE }
503template <>
505 const unsigned int,const Point &,const bool)
506{ NEDELEC_LOW_D_ERROR_MESSAGE }
507
508#endif
509
510} // namespace libMesh
unsigned int dim
The constraint matrix storage format.
Definition dof_map.h:112
This class handles the numbering of degrees of freedom on a mesh.
Definition dof_map.h:181
This is the base class from which all geometric element types are derived.
Definition elem.h:96
virtual unsigned int n_nodes() const =0
virtual ElemType type() const =0
unsigned int p_level() const
Definition elem.h:3125
static void get_refspace_nodes(const ElemType t, std::vector< Point > &nodes)
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual bool is_hierarchic() const override
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln, bool add_p_level=true, const unsigned vdim=1)
Build the nodal soln from the element soln.
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
virtual FEContinuity get_continuity() const override
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 unsigned int n_dofs(const ElemType t, const Order o)
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
virtual bool shapes_need_reinit() const override
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
std::string enum_to_string(const T e)
The libMesh namespace provides an interface to certain functionality in the library.
LIBMESH_DEFAULT_VECTORIZED_FE(template<>Real FE< 0, BERNSTEIN)
ElemType
Defines an enum for geometric element types.
libmesh_assert(ctx)
LIBMESH_FE_SIDE_NODAL_SOLN(HIERARCHIC_VEC)
const dof_id_type n_nodes
Definition tecplot_io.C:67