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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo