fe_nedelec_one.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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/fe.h"
23 #include "libmesh/fe_interface.h"
24 #include "libmesh/elem.h"
25 #include "libmesh/threads.h"
26 #include "libmesh/tensor_value.h"
27 #include "libmesh/string_to_enum.h"
28 
29 namespace libMesh
30 {
31 
32 // ------------------------------------------------------------
33 // Nedelec first kind specific implementations
34 
35 
36 // Anonymous namespace for local helper functions
37 namespace {
38 void nedelec_one_nodal_soln(const Elem * elem,
39  const Order order,
40  const std::vector<Number> & elem_soln,
41  const int dim,
42  std::vector<Number> & nodal_soln)
43 {
44  const unsigned int n_nodes = elem->n_nodes();
45  const ElemType elem_type = elem->type();
46 
47  const Order totalorder = static_cast<Order>(order+elem->p_level());
48 
49  nodal_soln.resize(n_nodes*dim);
50 
51  FEType fe_type(totalorder, NEDELEC_ONE);
52 
53  switch (totalorder)
54  {
55  case FIRST:
56  {
57  switch (elem_type)
58  {
59  case TRI6:
60  {
61  libmesh_assert_equal_to (elem_soln.size(), 3);
62  libmesh_assert_equal_to (nodal_soln.size(), 6*2);
63  break;
64  }
65  case QUAD8:
66  case QUAD9:
67  {
68  libmesh_assert_equal_to (elem_soln.size(), 4);
69 
70  if (elem_type == QUAD8)
71  libmesh_assert_equal_to (nodal_soln.size(), 8*2);
72  else
73  libmesh_assert_equal_to (nodal_soln.size(), 9*2);
74  break;
75  }
76  case TET10:
77  {
78  libmesh_assert_equal_to (elem_soln.size(), 6);
79  libmesh_assert_equal_to (nodal_soln.size(), 10*3);
80 
81  libmesh_not_implemented();
82 
83  break;
84  }
85 
86 
87  case HEX20:
88  case HEX27:
89  {
90  libmesh_assert_equal_to (elem_soln.size(), 12);
91 
92  if (elem_type == HEX20)
93  libmesh_assert_equal_to (nodal_soln.size(), 20*3);
94  else
95  libmesh_assert_equal_to (nodal_soln.size(), 27*3);
96 
97  break;
98  }
99 
100  default:
101  libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(elem_type) << " selected for NEDELEC_ONE FE family!");
102 
103  } // switch(elem_type)
104 
105  const unsigned int n_sf =
106  FEInterface::n_shape_functions(dim, fe_type, elem_type);
107 
108  std::vector<Point> refspace_nodes;
109  FEVectorBase::get_refspace_nodes(elem_type,refspace_nodes);
110  libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
111 
112 
113  // Need to create new fe object so the shape function as the FETransformation
114  // applied to it.
115  std::unique_ptr<FEVectorBase> vis_fe = FEVectorBase::build(dim,fe_type);
116 
117  const std::vector<std::vector<RealGradient>> & vis_phi = vis_fe->get_phi();
118 
119  vis_fe->reinit(elem,&refspace_nodes);
120 
121  for (unsigned int n = 0; n < n_nodes; n++)
122  {
123  libmesh_assert_equal_to (elem_soln.size(), n_sf);
124 
125  // Zero before summation
126  for (int d = 0; d < dim; d++)
127  {
128  nodal_soln[dim*n+d] = 0;
129  }
130 
131  // u = Sum (u_i phi_i)
132  for (unsigned int i=0; i<n_sf; i++)
133  {
134  for (int d = 0; d < dim; d++)
135  {
136  nodal_soln[dim*n+d] += elem_soln[i]*(vis_phi[i][n](d));
137  }
138  }
139  }
140 
141  return;
142  } // case FIRST
143 
144  default:
145  libmesh_error_msg("ERROR: Invalid total order " << Utility::enum_to_string(totalorder) << " selected for NEDELEC_ONE FE family!");
146 
147  }//switch (totalorder)
148 
149  return;
150 } // nedelec_one_nodal_soln
151 
152 
153 unsigned int nedelec_one_n_dofs(const ElemType t, const Order o)
154 {
155  switch (o)
156  {
157  case FIRST:
158  {
159  switch (t)
160  {
161  case TRI6:
162  return 3;
163 
164  case QUAD8:
165  case QUAD9:
166  return 4;
167 
168  case TET10:
169  return 6;
170 
171  case HEX20:
172  case HEX27:
173  return 12;
174 
175  case INVALID_ELEM:
176  return 0;
177 
178  default:
179  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
180  }
181  }
182 
183  default:
184  libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for NEDELEC_ONE FE family!");
185  }
186 }
187 
188 
189 
190 
191 unsigned int nedelec_one_n_dofs_at_node(const ElemType t,
192  const Order o,
193  const unsigned int n)
194 {
195  switch (o)
196  {
197  case FIRST:
198  {
199  switch (t)
200  {
201  case TRI6:
202  {
203  switch (n)
204  {
205  case 0:
206  case 1:
207  case 2:
208  return 0;
209  case 3:
210  case 4:
211  case 5:
212  return 1;
213 
214  default:
215  libmesh_error_msg("ERROR: Invalid node ID " << n);
216  }
217  }
218  case QUAD8:
219  {
220  switch (n)
221  {
222  case 0:
223  case 1:
224  case 2:
225  case 3:
226  return 0;
227  case 4:
228  case 5:
229  case 6:
230  case 7:
231  return 1;
232 
233  default:
234  libmesh_error_msg("ERROR: Invalid node ID " << n);
235  }
236  }
237  case QUAD9:
238  {
239  switch (n)
240  {
241  case 0:
242  case 1:
243  case 2:
244  case 3:
245  case 8:
246  return 0;
247  case 4:
248  case 5:
249  case 6:
250  case 7:
251  return 1;
252 
253  default:
254  libmesh_error_msg("ERROR: Invalid node ID " << n);
255  }
256  }
257  case TET10:
258  {
259  switch (n)
260  {
261  case 0:
262  case 1:
263  case 2:
264  case 3:
265  return 0;
266  case 4:
267  case 5:
268  case 6:
269  case 7:
270  case 8:
271  case 9:
272  return 1;
273 
274  default:
275  libmesh_error_msg("ERROR: Invalid node ID " << n);
276  }
277  }
278 
279  case HEX20:
280  {
281  switch (n)
282  {
283  case 0:
284  case 1:
285  case 2:
286  case 3:
287  case 4:
288  case 5:
289  case 6:
290  case 7:
291  return 0;
292  case 8:
293  case 9:
294  case 10:
295  case 11:
296  case 12:
297  case 13:
298  case 14:
299  case 15:
300  case 16:
301  case 17:
302  case 18:
303  case 19:
304  return 1;
305 
306  default:
307  libmesh_error_msg("ERROR: Invalid node ID " << n);
308  }
309  }
310  case HEX27:
311  {
312  switch (n)
313  {
314  case 0:
315  case 1:
316  case 2:
317  case 3:
318  case 4:
319  case 5:
320  case 6:
321  case 7:
322  case 20:
323  case 21:
324  case 22:
325  case 23:
326  case 24:
327  case 25:
328  case 26:
329  return 0;
330  case 8:
331  case 9:
332  case 10:
333  case 11:
334  case 12:
335  case 13:
336  case 14:
337  case 15:
338  case 16:
339  case 17:
340  case 18:
341  case 19:
342  return 1;
343 
344  default:
345  libmesh_error_msg("ERROR: Invalid node ID " << n);
346  }
347  }
348 
349  case INVALID_ELEM:
350  return 0;
351 
352  default:
353  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
354  }
355  }
356 
357  default:
358  libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for NEDELEC_ONE FE family!");
359  }
360 }
361 
362 
363 
364 #ifdef LIBMESH_ENABLE_AMR
365 void nedelec_one_compute_constraints (DofConstraints & /*constraints*/,
366  DofMap & /*dof_map*/,
367  const unsigned int /*variable_number*/,
368  const Elem * libmesh_dbg_var(elem),
369  const unsigned Dim)
370 {
371  // Only constrain elements in 2,3D.
372  if (Dim == 1)
373  return;
374 
375  libmesh_assert(elem);
376 
377  libmesh_not_implemented();
378 
379  /*
380  // Only constrain active and ancestor elements
381  if (elem->subactive())
382  return;
383 
384  FEType fe_type = dof_map.variable_type(variable_number);
385  fe_type.order = static_cast<Order>(fe_type.order + elem->p_level());
386 
387  std::vector<unsigned int> my_dof_indices, parent_dof_indices;
388 
389  // Look at the element faces. Check to see if we need to
390  // build constraints.
391  for (unsigned int s=0; s<elem->n_sides(); s++)
392  if (elem->neighbor(s) != nullptr)
393  if (elem->neighbor(s)->level() < elem->level()) // constrain dofs shared between
394  { // this element and ones coarser
395  // than this element.
396  // Get pointers to the elements of interest and its parent.
397  const Elem * parent = elem->parent();
398 
399  // This can't happen... Only level-0 elements have nullptr
400  // parents, and no level-0 elements can be at a higher
401  // level than their neighbors!
402  libmesh_assert(parent);
403 
404  const std::unique_ptr<const Elem> my_side (elem->build_side_ptr(s));
405  const std::unique_ptr<const Elem> parent_side (parent->build_side_ptr(s));
406 
407  // This function gets called element-by-element, so there
408  // will be a lot of memory allocation going on. We can
409  // at least minimize this for the case of the dof indices
410  // by efficiently preallocating the requisite storage.
411  my_dof_indices.reserve (my_side->n_nodes());
412  parent_dof_indices.reserve (parent_side->n_nodes());
413 
414  dof_map.dof_indices (my_side.get(), my_dof_indices,
415  variable_number);
416  dof_map.dof_indices (parent_side.get(), parent_dof_indices,
417  variable_number);
418 
419  for (unsigned int my_dof=0;
420  my_dof<FEInterface::n_dofs(Dim-1, fe_type, my_side->type());
421  my_dof++)
422  {
423  libmesh_assert_less (my_dof, my_side->n_nodes());
424 
425  // My global dof index.
426  const unsigned int my_dof_g = my_dof_indices[my_dof];
427 
428  // The support point of the DOF
429  const Point & support_point = my_side->point(my_dof);
430 
431  // Figure out where my node lies on their reference element.
432  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
433  parent_side.get(),
434  support_point);
435 
436  // Compute the parent's side shape function values.
437  for (unsigned int their_dof=0;
438  their_dof<FEInterface::n_dofs(Dim-1, fe_type, parent_side->type());
439  their_dof++)
440  {
441  libmesh_assert_less (their_dof, parent_side->n_nodes());
442 
443  // Their global dof index.
444  const unsigned int their_dof_g =
445  parent_dof_indices[their_dof];
446 
447  const Real their_dof_value = FEInterface::shape(Dim-1,
448  fe_type,
449  parent_side->type(),
450  their_dof,
451  mapped_point);
452 
453  // Only add non-zero and non-identity values
454  // for Lagrange basis functions.
455  if ((std::abs(their_dof_value) > 1.e-5) &&
456  (std::abs(their_dof_value) < .999))
457  {
458  // since we may be running this method concurrently
459  // on multiple threads we need to acquire a lock
460  // before modifying the shared constraint_row object.
461  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
462 
463  // A reference to the constraint row.
464  DofConstraintRow & constraint_row = constraints[my_dof_g].first;
465 
466  constraint_row.insert(std::make_pair (their_dof_g,
467  their_dof_value));
468  }
469  #ifdef DEBUG
470  // Protect for the case u_i = 0.999 u_j,
471  // in which case i better equal j.
472  else if (their_dof_value >= .999)
473  libmesh_assert_equal_to (my_dof_g, their_dof_g);
474  #endif
475  }
476  }
477  }
478  */
479 } // nedelec_one_compute_constraints()
480 #endif // #ifdef LIBMESH_ENABLE_AMR
481 
482 } // anonymous namespace
483 
484 #define NEDELEC_LOW_D_ERROR_MESSAGE \
485  libmesh_error_msg("ERROR: This method makes no sense for low-D elements!");
486 
487 
488 // Do full-specialization for every dimension, instead
489 // of explicit instantiation at the end of this file.
490 template <>
492  const Order,
493  const std::vector<Number> &,
494  std::vector<Number> &)
495 { NEDELEC_LOW_D_ERROR_MESSAGE }
496 
497 template <>
499  const Order,
500  const std::vector<Number> &,
501  std::vector<Number> &)
502 { NEDELEC_LOW_D_ERROR_MESSAGE }
503 
504 template <>
506  const Order order,
507  const std::vector<Number> & elem_soln,
508  std::vector<Number> & nodal_soln)
509 { nedelec_one_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln); }
510 
511 template <>
513  const Order order,
514  const std::vector<Number> & elem_soln,
515  std::vector<Number> & nodal_soln)
516 { nedelec_one_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln); }
517 
518 
519 // Do full-specialization for every dimension, instead
520 // of explicit instantiation at the end of this function.
521 // This could be macro-ified.
522 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
523 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
524 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); }
525 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); }
526 
527 
528 // Do full-specialization for every dimension, instead
529 // of explicit instantiation at the end of this function.
530 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
531 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
532 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); }
533 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); }
534 
535 
536 // Nedelec first type elements have no dofs per element
537 // FIXME: Only for first order!
538 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
539 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
540 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
541 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
542 
543 // Nedelec first type FEMs are always tangentially continuous
544 template <> FEContinuity FE<0,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE }
545 template <> FEContinuity FE<1,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE }
548 
549 // Nedelec first type FEMs are not hierarchic
550 template <> bool FE<0,NEDELEC_ONE>::is_hierarchic() const { NEDELEC_LOW_D_ERROR_MESSAGE }
551 template <> bool FE<1,NEDELEC_ONE>::is_hierarchic() const { NEDELEC_LOW_D_ERROR_MESSAGE }
552 template <> bool FE<2,NEDELEC_ONE>::is_hierarchic() const { return false; }
553 template <> bool FE<3,NEDELEC_ONE>::is_hierarchic() const { return false; }
554 
555 // Nedelec first type FEM shapes always need to be reinit'ed (because of orientation dependence)
556 template <> bool FE<0,NEDELEC_ONE>::shapes_need_reinit() const { NEDELEC_LOW_D_ERROR_MESSAGE }
557 template <> bool FE<1,NEDELEC_ONE>::shapes_need_reinit() const { NEDELEC_LOW_D_ERROR_MESSAGE }
558 template <> bool FE<2,NEDELEC_ONE>::shapes_need_reinit() const { return true; }
559 template <> bool FE<3,NEDELEC_ONE>::shapes_need_reinit() const { return true; }
560 
561 #ifdef LIBMESH_ENABLE_AMR
562 template <>
564  DofMap &,
565  const unsigned int,
566  const Elem *)
567 { NEDELEC_LOW_D_ERROR_MESSAGE }
568 
569 template <>
571  DofMap &,
572  const unsigned int,
573  const Elem *)
574 { NEDELEC_LOW_D_ERROR_MESSAGE }
575 
576 template <>
578  DofMap & dof_map,
579  const unsigned int variable_number,
580  const Elem * elem)
581 { nedelec_one_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
582 
583 template <>
585  DofMap & dof_map,
586  const unsigned int variable_number,
587  const Elem * elem)
588 { nedelec_one_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
589 #endif // LIBMESH_ENABLE_AMR
590 
591 // Specialize useless shape function methods
592 template <>
593 RealGradient FE<0,NEDELEC_ONE>::shape(const ElemType, const Order,const unsigned int,const Point &)
594 { NEDELEC_LOW_D_ERROR_MESSAGE }
595 template <>
596 RealGradient FE<0,NEDELEC_ONE>::shape(const Elem *,const Order,const unsigned int,const Point &)
597 { NEDELEC_LOW_D_ERROR_MESSAGE }
598 template <>
600  const unsigned int,const Point &)
601 { NEDELEC_LOW_D_ERROR_MESSAGE }
602 template <>
603 RealGradient FE<0,NEDELEC_ONE>::shape_deriv(const Elem *,const Order,const unsigned int,
604  const unsigned int,const Point &)
605 { NEDELEC_LOW_D_ERROR_MESSAGE }
606 template <>
608  const unsigned int,const Point &)
609 { NEDELEC_LOW_D_ERROR_MESSAGE }
610 template <>
612  const unsigned int,const Point &)
613 { NEDELEC_LOW_D_ERROR_MESSAGE }
614 
615 template <>
616 RealGradient FE<1,NEDELEC_ONE>::shape(const ElemType, const Order,const unsigned int,const Point &)
617 { NEDELEC_LOW_D_ERROR_MESSAGE }
618 template <>
619 RealGradient FE<1,NEDELEC_ONE>::shape(const Elem *,const Order,const unsigned int,const Point &)
620 { NEDELEC_LOW_D_ERROR_MESSAGE }
621 template <>
623  const unsigned int,const Point &)
624 { NEDELEC_LOW_D_ERROR_MESSAGE }
625 template <>
626 RealGradient FE<1,NEDELEC_ONE>::shape_deriv(const Elem *,const Order,const unsigned int,
627  const unsigned int,const Point &)
628 { NEDELEC_LOW_D_ERROR_MESSAGE }
629 template <>
631  const unsigned int,const Point &)
632 { NEDELEC_LOW_D_ERROR_MESSAGE }
633 template <>
635  const unsigned int,const Point &)
636 { NEDELEC_LOW_D_ERROR_MESSAGE }
637 
638 } // namespace libMesh
static unsigned int n_dofs(const ElemType t, const Order o)
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
The base class for all geometric element types.
Definition: elem.h:100
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)
virtual bool shapes_need_reinit() const override
virtual bool is_hierarchic() const override
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
const dof_id_type n_nodes
Definition: tecplot_io.C:68
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:428
static void get_refspace_nodes(const ElemType t, std::vector< Point > &nodes)
Definition: fe_abstract.C:283
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
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)
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
A geometric point in (x,y,z) space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)