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