libMesh
fe_nedelec_one.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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  UniquePtr<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  libmesh_error_msg("We'll never get here!");
188  return 0;
189 }
190 
191 
192 
193 
194 unsigned int nedelec_one_n_dofs_at_node(const ElemType t,
195  const Order o,
196  const unsigned int n)
197 {
198  switch (o)
199  {
200  case FIRST:
201  {
202  switch (t)
203  {
204  case TRI6:
205  {
206  switch (n)
207  {
208  case 0:
209  case 1:
210  case 2:
211  return 0;
212  case 3:
213  case 4:
214  case 5:
215  return 1;
216 
217  default:
218  libmesh_error_msg("ERROR: Invalid node ID " << n);
219  }
220  }
221  case QUAD8:
222  {
223  switch (n)
224  {
225  case 0:
226  case 1:
227  case 2:
228  case 3:
229  return 0;
230  case 4:
231  case 5:
232  case 6:
233  case 7:
234  return 1;
235 
236  default:
237  libmesh_error_msg("ERROR: Invalid node ID " << n);
238  }
239  }
240  case QUAD9:
241  {
242  switch (n)
243  {
244  case 0:
245  case 1:
246  case 2:
247  case 3:
248  case 8:
249  return 0;
250  case 4:
251  case 5:
252  case 6:
253  case 7:
254  return 1;
255 
256  default:
257  libmesh_error_msg("ERROR: Invalid node ID " << n);
258  }
259  }
260  case TET10:
261  {
262  switch (n)
263  {
264  case 0:
265  case 1:
266  case 2:
267  case 3:
268  return 0;
269  case 4:
270  case 5:
271  case 6:
272  case 7:
273  case 8:
274  case 9:
275  return 1;
276 
277  default:
278  libmesh_error_msg("ERROR: Invalid node ID " << n);
279  }
280  }
281 
282  case HEX20:
283  {
284  switch (n)
285  {
286  case 0:
287  case 1:
288  case 2:
289  case 3:
290  case 4:
291  case 5:
292  case 6:
293  case 7:
294  return 0;
295  case 8:
296  case 9:
297  case 10:
298  case 11:
299  case 12:
300  case 13:
301  case 14:
302  case 15:
303  case 16:
304  case 17:
305  case 18:
306  case 19:
307  return 1;
308 
309  default:
310  libmesh_error_msg("ERROR: Invalid node ID " << n);
311  }
312  }
313  case HEX27:
314  {
315  switch (n)
316  {
317  case 0:
318  case 1:
319  case 2:
320  case 3:
321  case 4:
322  case 5:
323  case 6:
324  case 7:
325  case 20:
326  case 21:
327  case 22:
328  case 23:
329  case 24:
330  case 25:
331  case 26:
332  return 0;
333  case 8:
334  case 9:
335  case 10:
336  case 11:
337  case 12:
338  case 13:
339  case 14:
340  case 15:
341  case 16:
342  case 17:
343  case 18:
344  case 19:
345  return 1;
346 
347  default:
348  libmesh_error_msg("ERROR: Invalid node ID " << n);
349  }
350  }
351 
352  case INVALID_ELEM:
353  return 0;
354 
355  default:
356  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
357  }
358  }
359 
360  default:
361  libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for NEDELEC_ONE FE family!");
362  }
363 
364  libmesh_error_msg("We'll never get here!");
365  return 0;
366 }
367 
368 
369 
370 #ifdef LIBMESH_ENABLE_AMR
371 void nedelec_one_compute_constraints (DofConstraints & /*constraints*/,
372  DofMap & /*dof_map*/,
373  const unsigned int /*variable_number*/,
374  const Elem * libmesh_dbg_var(elem),
375  const unsigned Dim)
376 {
377  // Only constrain elements in 2,3D.
378  if (Dim == 1)
379  return;
380 
381  libmesh_assert(elem);
382 
383  libmesh_not_implemented();
384 
385  /*
386  // Only constrain active and ancestor elements
387  if (elem->subactive())
388  return;
389 
390  FEType fe_type = dof_map.variable_type(variable_number);
391  fe_type.order = static_cast<Order>(fe_type.order + elem->p_level());
392 
393  std::vector<unsigned int> my_dof_indices, parent_dof_indices;
394 
395  // Look at the element faces. Check to see if we need to
396  // build constraints.
397  for (unsigned int s=0; s<elem->n_sides(); s++)
398  if (elem->neighbor(s) != libmesh_nullptr)
399  if (elem->neighbor(s)->level() < elem->level()) // constrain dofs shared between
400  { // this element and ones coarser
401  // than this element.
402  // Get pointers to the elements of interest and its parent.
403  const Elem * parent = elem->parent();
404 
405  // This can't happen... Only level-0 elements have NULL
406  // parents, and no level-0 elements can be at a higher
407  // level than their neighbors!
408  libmesh_assert(parent);
409 
410  const UniquePtr<const Elem> my_side (elem->build_side_ptr(s));
411  const UniquePtr<const Elem> parent_side (parent->build_side_ptr(s));
412 
413  // This function gets called element-by-element, so there
414  // will be a lot of memory allocation going on. We can
415  // at least minimize this for the case of the dof indices
416  // by efficiently preallocating the requisite storage.
417  my_dof_indices.reserve (my_side->n_nodes());
418  parent_dof_indices.reserve (parent_side->n_nodes());
419 
420  dof_map.dof_indices (my_side.get(), my_dof_indices,
421  variable_number);
422  dof_map.dof_indices (parent_side.get(), parent_dof_indices,
423  variable_number);
424 
425  for (unsigned int my_dof=0;
426  my_dof<FEInterface::n_dofs(Dim-1, fe_type, my_side->type());
427  my_dof++)
428  {
429  libmesh_assert_less (my_dof, my_side->n_nodes());
430 
431  // My global dof index.
432  const unsigned int my_dof_g = my_dof_indices[my_dof];
433 
434  // The support point of the DOF
435  const Point & support_point = my_side->point(my_dof);
436 
437  // Figure out where my node lies on their reference element.
438  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
439  parent_side.get(),
440  support_point);
441 
442  // Compute the parent's side shape function values.
443  for (unsigned int their_dof=0;
444  their_dof<FEInterface::n_dofs(Dim-1, fe_type, parent_side->type());
445  their_dof++)
446  {
447  libmesh_assert_less (their_dof, parent_side->n_nodes());
448 
449  // Their global dof index.
450  const unsigned int their_dof_g =
451  parent_dof_indices[their_dof];
452 
453  const Real their_dof_value = FEInterface::shape(Dim-1,
454  fe_type,
455  parent_side->type(),
456  their_dof,
457  mapped_point);
458 
459  // Only add non-zero and non-identity values
460  // for Lagrange basis functions.
461  if ((std::abs(their_dof_value) > 1.e-5) &&
462  (std::abs(their_dof_value) < .999))
463  {
464  // since we may be running this method concurrently
465  // on multiple threads we need to acquire a lock
466  // before modifying the shared constraint_row object.
467  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
468 
469  // A reference to the constraint row.
470  DofConstraintRow & constraint_row = constraints[my_dof_g].first;
471 
472  constraint_row.insert(std::make_pair (their_dof_g,
473  their_dof_value));
474  }
475  #ifdef DEBUG
476  // Protect for the case u_i = 0.999 u_j,
477  // in which case i better equal j.
478  else if (their_dof_value >= .999)
479  libmesh_assert_equal_to (my_dof_g, their_dof_g);
480  #endif
481  }
482  }
483  }
484  */
485 } // nedelec_one_compute_constraints()
486 #endif // #ifdef LIBMESH_ENABLE_AMR
487 
488 } // anonymous namespace
489 
490 #define NEDELEC_LOW_D_ERROR_MESSAGE \
491  libmesh_error_msg("ERROR: This method makes no sense for low-D elements!");
492 
493 
494 // Do full-specialization for every dimension, instead
495 // of explicit instantiation at the end of this file.
496 template <>
498  const Order,
499  const std::vector<Number> &,
500  std::vector<Number> &)
501 { NEDELEC_LOW_D_ERROR_MESSAGE }
502 
503 template <>
505  const Order,
506  const std::vector<Number> &,
507  std::vector<Number> &)
508 { NEDELEC_LOW_D_ERROR_MESSAGE }
509 
510 template <>
512  const Order order,
513  const std::vector<Number> & elem_soln,
514  std::vector<Number> & nodal_soln)
515 { nedelec_one_nodal_soln(elem, order, elem_soln, 2 /*dim*/, nodal_soln); }
516 
517 template <>
519  const Order order,
520  const std::vector<Number> & elem_soln,
521  std::vector<Number> & nodal_soln)
522 { nedelec_one_nodal_soln(elem, order, elem_soln, 3 /*dim*/, nodal_soln); }
523 
524 
525 // Do full-specialization for every dimension, instead
526 // of explicit instantiation at the end of this function.
527 // This could be macro-ified.
528 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
529 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
530 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); }
531 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); }
532 
533 
534 // Do full-specialization for every dimension, instead
535 // of explicit instantiation at the end of this function.
536 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
537 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE }
538 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); }
539 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); }
540 
541 
542 // Nedelec first type elements have no dofs per element
543 // FIXME: Only for first order!
544 template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
545 template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE }
546 template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
547 template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
548 
549 // Nedelec first type FEMs are always tangentially continuous
550 template <> FEContinuity FE<0,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE }
551 template <> FEContinuity FE<1,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE }
554 
555 // Nedelec first type FEMs are not hierarchic
556 template <> bool FE<0,NEDELEC_ONE>::is_hierarchic() const { NEDELEC_LOW_D_ERROR_MESSAGE }
557 template <> bool FE<1,NEDELEC_ONE>::is_hierarchic() const { NEDELEC_LOW_D_ERROR_MESSAGE }
558 template <> bool FE<2,NEDELEC_ONE>::is_hierarchic() const { return false; }
559 template <> bool FE<3,NEDELEC_ONE>::is_hierarchic() const { return false; }
560 
561 // Nedelec first type FEM shapes always need to be reinit'ed (because of orientation dependence)
562 template <> bool FE<0,NEDELEC_ONE>::shapes_need_reinit() const { NEDELEC_LOW_D_ERROR_MESSAGE }
563 template <> bool FE<1,NEDELEC_ONE>::shapes_need_reinit() const { NEDELEC_LOW_D_ERROR_MESSAGE }
564 template <> bool FE<2,NEDELEC_ONE>::shapes_need_reinit() const { return true; }
565 template <> bool FE<3,NEDELEC_ONE>::shapes_need_reinit() const { return true; }
566 
567 #ifdef LIBMESH_ENABLE_AMR
568 template <>
570  DofMap &,
571  const unsigned int,
572  const Elem *)
573 { NEDELEC_LOW_D_ERROR_MESSAGE }
574 
575 template <>
577  DofMap &,
578  const unsigned int,
579  const Elem *)
580 { NEDELEC_LOW_D_ERROR_MESSAGE }
581 
582 template <>
584  DofMap & dof_map,
585  const unsigned int variable_number,
586  const Elem * elem)
587 { nedelec_one_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
588 
589 template <>
591  DofMap & dof_map,
592  const unsigned int variable_number,
593  const Elem * elem)
594 { nedelec_one_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
595 #endif // LIBMESH_ENABLE_AMR
596 
597 // Specialize useless shape function methods
598 template <>
599 RealGradient FE<0,NEDELEC_ONE>::shape(const ElemType, const Order,const unsigned int,const Point &)
600 { NEDELEC_LOW_D_ERROR_MESSAGE }
601 template <>
602 RealGradient FE<0,NEDELEC_ONE>::shape(const Elem *,const Order,const unsigned int,const Point &)
603 { NEDELEC_LOW_D_ERROR_MESSAGE }
604 template <>
606  const unsigned int,const Point &)
607 { NEDELEC_LOW_D_ERROR_MESSAGE }
608 template <>
609 RealGradient FE<0,NEDELEC_ONE>::shape_deriv(const Elem *,const Order,const unsigned int,
610  const unsigned int,const Point &)
611 { NEDELEC_LOW_D_ERROR_MESSAGE }
612 template <>
614  const unsigned int,const Point &)
615 { NEDELEC_LOW_D_ERROR_MESSAGE }
616 template <>
618  const unsigned int,const Point &)
619 { NEDELEC_LOW_D_ERROR_MESSAGE }
620 
621 template <>
622 RealGradient FE<1,NEDELEC_ONE>::shape(const ElemType, const Order,const unsigned int,const Point &)
623 { NEDELEC_LOW_D_ERROR_MESSAGE }
624 template <>
625 RealGradient FE<1,NEDELEC_ONE>::shape(const Elem *,const Order,const unsigned int,const Point &)
626 { NEDELEC_LOW_D_ERROR_MESSAGE }
627 template <>
629  const unsigned int,const Point &)
630 { NEDELEC_LOW_D_ERROR_MESSAGE }
631 template <>
632 RealGradient FE<1,NEDELEC_ONE>::shape_deriv(const Elem *,const Order,const unsigned int,
633  const unsigned int,const Point &)
634 { NEDELEC_LOW_D_ERROR_MESSAGE }
635 template <>
637  const unsigned int,const Point &)
638 { NEDELEC_LOW_D_ERROR_MESSAGE }
639 template <>
641  const unsigned int,const Point &)
642 { NEDELEC_LOW_D_ERROR_MESSAGE }
643 
644 } // 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)
virtual bool shapes_need_reinit() const libmesh_override
unsigned int dim
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
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)
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
const dof_id_type n_nodes
Definition: tecplot_io.C:67
virtual FEContinuity get_continuity() const libmesh_override
static unsigned int n_shape_functions(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:385
static void get_refspace_nodes(const ElemType t, std::vector< Point > &nodes)
Definition: fe_abstract.C:259
virtual bool is_hierarchic() const libmesh_override
PetscErrorCode Vec Mat libmesh_dbg_var(j)
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
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...
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
Build the nodal soln from the element soln.
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
The constraint matrix storage format.
Definition: dof_map.h:96
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.