fe_lagrange.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 
27 namespace libMesh
28 {
29 
30  // ------------------------------------------------------------
31  // Lagrange-specific implementations
32 
33 
34  // Anonymous namespace for local helper functions
35  namespace {
36  void lagrange_nodal_soln(const Elem* elem,
37  const Order order,
38  const std::vector<Number>& elem_soln,
39  std::vector<Number>& nodal_soln)
40  {
41  const unsigned int n_nodes = elem->n_nodes();
42  const ElemType type = elem->type();
43 
44  const Order totalorder = static_cast<Order>(order+elem->p_level());
45 
46  nodal_soln.resize(n_nodes);
47 
48 
49 
50  switch (totalorder)
51  {
52  // linear Lagrange shape functions
53  case FIRST:
54  {
55  switch (type)
56  {
57  case EDGE3:
58  {
59  libmesh_assert_equal_to (elem_soln.size(), 2);
60  libmesh_assert_equal_to (nodal_soln.size(), 3);
61 
62  nodal_soln[0] = elem_soln[0];
63  nodal_soln[1] = elem_soln[1];
64  nodal_soln[2] = .5*(elem_soln[0] + elem_soln[1]);
65 
66  return;
67  }
68 
69  case EDGE4:
70  {
71  libmesh_assert_equal_to (elem_soln.size(), 2);
72  libmesh_assert_equal_to (nodal_soln.size(), 4);
73 
74  nodal_soln[0] = elem_soln[0];
75  nodal_soln[1] = elem_soln[1];
76  nodal_soln[2] = (2.*elem_soln[0] + elem_soln[1])/3.;
77  nodal_soln[3] = (elem_soln[0] + 2.*elem_soln[1])/3.;
78 
79  return;
80  }
81 
82 
83  case TRI6:
84  {
85  libmesh_assert_equal_to (elem_soln.size(), 3);
86  libmesh_assert_equal_to (nodal_soln.size(), 6);
87 
88  nodal_soln[0] = elem_soln[0];
89  nodal_soln[1] = elem_soln[1];
90  nodal_soln[2] = elem_soln[2];
91  nodal_soln[3] = .5*(elem_soln[0] + elem_soln[1]);
92  nodal_soln[4] = .5*(elem_soln[1] + elem_soln[2]);
93  nodal_soln[5] = .5*(elem_soln[2] + elem_soln[0]);
94 
95  return;
96  }
97 
98 
99  case QUAD8:
100  case QUAD9:
101  {
102  libmesh_assert_equal_to (elem_soln.size(), 4);
103 
104  if (type == QUAD8)
105  libmesh_assert_equal_to (nodal_soln.size(), 8);
106  else
107  libmesh_assert_equal_to (nodal_soln.size(), 9);
108 
109 
110  nodal_soln[0] = elem_soln[0];
111  nodal_soln[1] = elem_soln[1];
112  nodal_soln[2] = elem_soln[2];
113  nodal_soln[3] = elem_soln[3];
114  nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
115  nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
116  nodal_soln[6] = .5*(elem_soln[2] + elem_soln[3]);
117  nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
118 
119  if (type == QUAD9)
120  nodal_soln[8] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
121 
122  return;
123  }
124 
125 
126  case TET10:
127  {
128  libmesh_assert_equal_to (elem_soln.size(), 4);
129  libmesh_assert_equal_to (nodal_soln.size(), 10);
130 
131  nodal_soln[0] = elem_soln[0];
132  nodal_soln[1] = elem_soln[1];
133  nodal_soln[2] = elem_soln[2];
134  nodal_soln[3] = elem_soln[3];
135  nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
136  nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
137  nodal_soln[6] = .5*(elem_soln[2] + elem_soln[0]);
138  nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
139  nodal_soln[8] = .5*(elem_soln[3] + elem_soln[1]);
140  nodal_soln[9] = .5*(elem_soln[3] + elem_soln[2]);
141 
142  return;
143  }
144 
145 
146  case HEX20:
147  case HEX27:
148  {
149  libmesh_assert_equal_to (elem_soln.size(), 8);
150 
151  if (type == HEX20)
152  libmesh_assert_equal_to (nodal_soln.size(), 20);
153  else
154  libmesh_assert_equal_to (nodal_soln.size(), 27);
155 
156  nodal_soln[0] = elem_soln[0];
157  nodal_soln[1] = elem_soln[1];
158  nodal_soln[2] = elem_soln[2];
159  nodal_soln[3] = elem_soln[3];
160  nodal_soln[4] = elem_soln[4];
161  nodal_soln[5] = elem_soln[5];
162  nodal_soln[6] = elem_soln[6];
163  nodal_soln[7] = elem_soln[7];
164  nodal_soln[8] = .5*(elem_soln[0] + elem_soln[1]);
165  nodal_soln[9] = .5*(elem_soln[1] + elem_soln[2]);
166  nodal_soln[10] = .5*(elem_soln[2] + elem_soln[3]);
167  nodal_soln[11] = .5*(elem_soln[3] + elem_soln[0]);
168  nodal_soln[12] = .5*(elem_soln[0] + elem_soln[4]);
169  nodal_soln[13] = .5*(elem_soln[1] + elem_soln[5]);
170  nodal_soln[14] = .5*(elem_soln[2] + elem_soln[6]);
171  nodal_soln[15] = .5*(elem_soln[3] + elem_soln[7]);
172  nodal_soln[16] = .5*(elem_soln[4] + elem_soln[5]);
173  nodal_soln[17] = .5*(elem_soln[5] + elem_soln[6]);
174  nodal_soln[18] = .5*(elem_soln[6] + elem_soln[7]);
175  nodal_soln[19] = .5*(elem_soln[4] + elem_soln[7]);
176 
177  if (type == HEX27)
178  {
179  nodal_soln[20] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
180  nodal_soln[21] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[5]);
181  nodal_soln[22] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[6]);
182  nodal_soln[23] = .25*(elem_soln[2] + elem_soln[3] + elem_soln[6] + elem_soln[7]);
183  nodal_soln[24] = .25*(elem_soln[3] + elem_soln[0] + elem_soln[7] + elem_soln[4]);
184  nodal_soln[25] = .25*(elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
185 
186  nodal_soln[26] = .125*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3] +
187  elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
188  }
189 
190  return;
191  }
192 
193 
194  case PRISM15:
195  case PRISM18:
196  {
197  libmesh_assert_equal_to (elem_soln.size(), 6);
198 
199  if (type == PRISM15)
200  libmesh_assert_equal_to (nodal_soln.size(), 15);
201  else
202  libmesh_assert_equal_to (nodal_soln.size(), 18);
203 
204  nodal_soln[0] = elem_soln[0];
205  nodal_soln[1] = elem_soln[1];
206  nodal_soln[2] = elem_soln[2];
207  nodal_soln[3] = elem_soln[3];
208  nodal_soln[4] = elem_soln[4];
209  nodal_soln[5] = elem_soln[5];
210  nodal_soln[6] = .5*(elem_soln[0] + elem_soln[1]);
211  nodal_soln[7] = .5*(elem_soln[1] + elem_soln[2]);
212  nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]);
213  nodal_soln[9] = .5*(elem_soln[0] + elem_soln[3]);
214  nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
215  nodal_soln[11] = .5*(elem_soln[2] + elem_soln[5]);
216  nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
217  nodal_soln[13] = .5*(elem_soln[4] + elem_soln[5]);
218  nodal_soln[14] = .5*(elem_soln[3] + elem_soln[5]);
219 
220  if (type == PRISM18)
221  {
222  nodal_soln[15] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[3]);
223  nodal_soln[16] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[4]);
224  nodal_soln[17] = .25*(elem_soln[2] + elem_soln[0] + elem_soln[3] + elem_soln[5]);
225  }
226 
227  return;
228  }
229 
230  case PYRAMID14:
231  {
232  libmesh_assert_equal_to (elem_soln.size(), 5);
233  libmesh_assert_equal_to (nodal_soln.size(), 14);
234 
235  nodal_soln[0] = elem_soln[0];
236  nodal_soln[1] = elem_soln[1];
237  nodal_soln[2] = elem_soln[2];
238  nodal_soln[3] = elem_soln[3];
239  nodal_soln[4] = elem_soln[4];
240  nodal_soln[5] = .5*(elem_soln[0] + elem_soln[1]);
241  nodal_soln[6] = .5*(elem_soln[1] + elem_soln[2]);
242  nodal_soln[7] = .5*(elem_soln[2] + elem_soln[3]);
243  nodal_soln[8] = .5*(elem_soln[3] + elem_soln[0]);
244  nodal_soln[9] = .5*(elem_soln[0] + elem_soln[4]);
245  nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
246  nodal_soln[11] = .5*(elem_soln[2] + elem_soln[4]);
247  nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
248  nodal_soln[13] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
249 
250  return;
251  }
252 
253  default:
254  {
255  // By default the element solution _is_ nodal,
256  // so just copy it.
257  nodal_soln = elem_soln;
258 
259  return;
260  }
261  }
262  }
263 
264  case SECOND:
265  {
266  switch (type)
267  {
268  case EDGE4:
269  {
270  libmesh_assert_equal_to (elem_soln.size(), 3);
271  libmesh_assert_equal_to (nodal_soln.size(), 4);
272 
273  // Project quadratic solution onto cubic element nodes
274  nodal_soln[0] = elem_soln[0];
275  nodal_soln[1] = elem_soln[1];
276  nodal_soln[2] = (2.*elem_soln[0] - elem_soln[1] +
277  8.*elem_soln[2])/9.;
278  nodal_soln[3] = (-elem_soln[0] + 2.*elem_soln[1] +
279  8.*elem_soln[2])/9.;
280  return;
281  }
282 
283  default:
284  {
285  // By default the element solution _is_ nodal,
286  // so just copy it.
287  nodal_soln = elem_soln;
288 
289  return;
290  }
291  }
292  }
293 
294 
295 
296 
297  default:
298  {
299  // By default the element solution _is_ nodal,
300  // so just copy it.
301  nodal_soln = elem_soln;
302 
303  return;
304  }
305  }
306  }
307 
308 
309 
310  unsigned int lagrange_n_dofs(const ElemType t, const Order o)
311  {
312  switch (o)
313  {
314 
315  // linear Lagrange shape functions
316  case FIRST:
317  {
318  switch (t)
319  {
320  case NODEELEM:
321  return 1;
322 
323  case EDGE2:
324  case EDGE3:
325  case EDGE4:
326  return 2;
327 
328  case TRI3:
329  case TRI6:
330  return 3;
331 
332  case QUAD4:
333  case QUAD8:
334  case QUAD9:
335  return 4;
336 
337  case TET4:
338  case TET10:
339  return 4;
340 
341  case HEX8:
342  case HEX20:
343  case HEX27:
344  return 8;
345 
346  case PRISM6:
347  case PRISM15:
348  case PRISM18:
349  return 6;
350 
351  case PYRAMID5:
352  case PYRAMID14:
353  return 5;
354 
355  default:
356  {
357 #ifdef DEBUG
358  libMesh::err << "ERROR: Bad ElemType = " << t
359  << " for FIRST order approximation!"
360  << std::endl;
361 #endif
362  libmesh_error();
363  }
364  }
365  }
366 
367 
368  // quadratic Lagrange shape functions
369  case SECOND:
370  {
371  switch (t)
372  {
373  case NODEELEM:
374  return 1;
375 
376  case EDGE3:
377  return 3;
378 
379  case TRI6:
380  return 6;
381 
382  case QUAD8:
383  return 8;
384 
385  case QUAD9:
386  return 9;
387 
388  case TET10:
389  return 10;
390 
391  case HEX20:
392  return 20;
393 
394  case HEX27:
395  return 27;
396 
397  case PRISM15:
398  return 15;
399 
400  case PRISM18:
401  return 18;
402 
403  case PYRAMID14:
404  return 14;
405 
406  default:
407  {
408 #ifdef DEBUG
409  libMesh::err << "ERROR: Bad ElemType = " << t
410  << " for SECOND order approximation!"
411  << std::endl;
412 #endif
413  libmesh_error();
414  }
415  }
416  }
417 
418  case THIRD:
419  {
420  switch (t)
421  {
422  case NODEELEM:
423  return 1;
424 
425  case EDGE4:
426  return 4;
427 
428  default:
429  {
430 #ifdef DEBUG
431  libMesh::err << "ERROR: Bad ElemType = " << t
432  << " for THIRD order approximation!"
433  << std::endl;
434 #endif
435  libmesh_error();
436  }
437  }
438  }
439 
440  default:
441  libmesh_error();
442  }
443 
444  libmesh_error();
445  return 0;
446  }
447 
448 
449 
450 
451  unsigned int lagrange_n_dofs_at_node(const ElemType t,
452  const Order o,
453  const unsigned int n)
454  {
455  switch (o)
456  {
457 
458  // linear Lagrange shape functions
459  case FIRST:
460  {
461  switch (t)
462  {
463  case NODEELEM:
464  return 1;
465 
466  case EDGE2:
467  case EDGE3:
468  case EDGE4:
469  {
470  switch (n)
471  {
472  case 0:
473  case 1:
474  return 1;
475 
476  default:
477  return 0;
478  }
479  }
480 
481  case TRI3:
482  case TRI6:
483  {
484  switch (n)
485  {
486  case 0:
487  case 1:
488  case 2:
489  return 1;
490 
491  default:
492  return 0;
493  }
494  }
495 
496  case QUAD4:
497  case QUAD8:
498  case QUAD9:
499  {
500  switch (n)
501  {
502  case 0:
503  case 1:
504  case 2:
505  case 3:
506  return 1;
507 
508  default:
509  return 0;
510  }
511  }
512 
513 
514  case TET4:
515  case TET10:
516  {
517  switch (n)
518  {
519  case 0:
520  case 1:
521  case 2:
522  case 3:
523  return 1;
524 
525  default:
526  return 0;
527  }
528  }
529 
530  case HEX8:
531  case HEX20:
532  case HEX27:
533  {
534  switch (n)
535  {
536  case 0:
537  case 1:
538  case 2:
539  case 3:
540  case 4:
541  case 5:
542  case 6:
543  case 7:
544  return 1;
545 
546  default:
547  return 0;
548  }
549  }
550 
551  case PRISM6:
552  case PRISM15:
553  case PRISM18:
554  {
555  switch (n)
556  {
557  case 0:
558  case 1:
559  case 2:
560  case 3:
561  case 4:
562  case 5:
563  return 1;
564 
565  default:
566  return 0;
567  }
568  }
569 
570  case PYRAMID5:
571  case PYRAMID14:
572  {
573  switch (n)
574  {
575  case 0:
576  case 1:
577  case 2:
578  case 3:
579  case 4:
580  return 1;
581 
582  default:
583  return 0;
584  }
585  }
586 
587  default:
588  {
589 #ifdef DEBUG
590  libMesh::err << "ERROR: Bad ElemType = " << t
591  << " for FIRST order approximation!"
592  << std::endl;
593 #endif
594  libmesh_error();
595  }
596  }
597  }
598 
599  // quadratic Lagrange shape functions
600  case SECOND:
601  {
602  switch (t)
603  {
604  // quadratic lagrange has one dof at each node
605  case NODEELEM:
606  case EDGE3:
607  case TRI6:
608  case QUAD8:
609  case QUAD9:
610  case TET10:
611  case HEX20:
612  case HEX27:
613  case PRISM15:
614  case PRISM18:
615  case PYRAMID14:
616  return 1;
617 
618  default:
619  {
620 #ifdef DEBUG
621  libMesh::err << "ERROR: Bad ElemType = " << t
622  << " for SECOND order approximation!"
623  << std::endl;
624 #endif
625  libmesh_error();
626  }
627  }
628  }
629 
630  case THIRD:
631  {
632  switch (t)
633  {
634  case NODEELEM:
635  case EDGE4:
636  return 1;
637 
638  default:
639  {
640 #ifdef DEBUG
641  libMesh::err << "ERROR: Bad ElemType = " << t
642  << " for THIRD order approximation!"
643  << std::endl;
644 #endif
645  libmesh_error();
646  }
647  }
648  }
649 
650 
651 
652  default:
653  libmesh_error();
654  }
655 
656  libmesh_error();
657  return 0;
658 
659  }
660 
661 
662 
663 #ifdef LIBMESH_ENABLE_AMR
664  void lagrange_compute_constraints (DofConstraints &constraints,
665  DofMap &dof_map,
666  const unsigned int variable_number,
667  const Elem* elem,
668  const unsigned Dim)
669  {
670  // Only constrain elements in 2,3D.
671  if (Dim == 1)
672  return;
673 
674  libmesh_assert(elem);
675 
676  // Only constrain active and ancestor elements
677  if (elem->subactive())
678  return;
679 
680  FEType fe_type = dof_map.variable_type(variable_number);
681  fe_type.order = static_cast<Order>(fe_type.order + elem->p_level());
682 
683  std::vector<dof_id_type> my_dof_indices, parent_dof_indices;
684 
685  // Look at the element faces. Check to see if we need to
686  // build constraints.
687  for (unsigned int s=0; s<elem->n_sides(); s++)
688  if (elem->neighbor(s) != NULL)
689  if (elem->neighbor(s)->level() < elem->level()) // constrain dofs shared between
690  { // this element and ones coarser
691  // than this element.
692  // Get pointers to the elements of interest and its parent.
693  const Elem* parent = elem->parent();
694 
695  // This can't happen... Only level-0 elements have NULL
696  // parents, and no level-0 elements can be at a higher
697  // level than their neighbors!
698  libmesh_assert(parent);
699 
700  const AutoPtr<Elem> my_side (elem->build_side(s));
701  const AutoPtr<Elem> parent_side (parent->build_side(s));
702 
703  // This function gets called element-by-element, so there
704  // will be a lot of memory allocation going on. We can
705  // at least minimize this for the case of the dof indices
706  // by efficiently preallocating the requisite storage.
707  my_dof_indices.reserve (my_side->n_nodes());
708  parent_dof_indices.reserve (parent_side->n_nodes());
709 
710  dof_map.dof_indices (my_side.get(), my_dof_indices,
711  variable_number);
712  dof_map.dof_indices (parent_side.get(), parent_dof_indices,
713  variable_number);
714 
715  for (unsigned int my_dof=0;
716  my_dof<FEInterface::n_dofs(Dim-1, fe_type, my_side->type());
717  my_dof++)
718  {
719  libmesh_assert_less (my_dof, my_side->n_nodes());
720 
721  // My global dof index.
722  const dof_id_type my_dof_g = my_dof_indices[my_dof];
723 
724  // Hunt for "constraining against myself" cases before
725  // we bother creating a constraint row
726  bool self_constraint = false;
727  for (unsigned int their_dof=0;
728  their_dof<FEInterface::n_dofs(Dim-1, fe_type, parent_side->type());
729  their_dof++)
730  {
731  libmesh_assert_less (their_dof, parent_side->n_nodes());
732 
733  // Their global dof index.
734  const dof_id_type their_dof_g =
735  parent_dof_indices[their_dof];
736 
737  if (their_dof_g == my_dof_g)
738  {
739  self_constraint = true;
740  break;
741  }
742  }
743 
744  if (self_constraint)
745  continue;
746 
747  DofConstraintRow* constraint_row;
748 
749  // we may be running constraint methods concurretly
750  // on multiple threads, so we need a lock to
751  // ensure that this constraint is "ours"
752  {
753  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
754 
755  if (dof_map.is_constrained_dof(my_dof_g))
756  continue;
757 
758  constraint_row = &(constraints[my_dof_g]);
759  libmesh_assert(constraint_row->empty());
760  }
761 
762  // The support point of the DOF
763  const Point& support_point = my_side->point(my_dof);
764 
765  // Figure out where my node lies on their reference element.
766  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
767  parent_side.get(),
768  support_point);
769 
770  // Compute the parent's side shape function values.
771  for (unsigned int their_dof=0;
772  their_dof<FEInterface::n_dofs(Dim-1, fe_type, parent_side->type());
773  their_dof++)
774  {
775  libmesh_assert_less (their_dof, parent_side->n_nodes());
776 
777  // Their global dof index.
778  const dof_id_type their_dof_g =
779  parent_dof_indices[their_dof];
780 
781  const Real their_dof_value = FEInterface::shape(Dim-1,
782  fe_type,
783  parent_side->type(),
784  their_dof,
785  mapped_point);
786 
787  // Only add non-zero and non-identity values
788  // for Lagrange basis functions.
789  if ((std::abs(their_dof_value) > 1.e-5) &&
790  (std::abs(their_dof_value) < .999))
791  {
792  constraint_row->insert(std::make_pair (their_dof_g,
793  their_dof_value));
794  }
795 #ifdef DEBUG
796  // Protect for the case u_i = 0.999 u_j,
797  // in which case i better equal j.
798  else if (their_dof_value >= .999)
799  libmesh_assert_equal_to (my_dof_g, their_dof_g);
800 #endif
801  }
802  }
803  }
804  } // lagrange_compute_constrants()
805 #endif // #ifdef LIBMESH_ENABLE_AMR
806 
807  } // anonymous namespace
808 
809 
810  // Do full-specialization for every dimension, instead
811  // of explicit instantiation at the end of this file.
812  // This could be macro-ified so that it fits on one line...
813  template <>
815  const Order order,
816  const std::vector<Number>& elem_soln,
817  std::vector<Number>& nodal_soln)
818  { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
819 
820  template <>
822  const Order order,
823  const std::vector<Number>& elem_soln,
824  std::vector<Number>& nodal_soln)
825  { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
826 
827  template <>
829  const Order order,
830  const std::vector<Number>& elem_soln,
831  std::vector<Number>& nodal_soln)
832  { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
833 
834  template <>
836  const Order order,
837  const std::vector<Number>& elem_soln,
838  std::vector<Number>& nodal_soln)
839  { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
840 
841 
842  // Do full-specialization for every dimension, instead
843  // of explicit instantiation at the end of this function.
844  // This could be macro-ified.
845  template <> unsigned int FE<0,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
846  template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
847  template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
848  template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
849 
850 
851  // Do full-specialization for every dimension, instead
852  // of explicit instantiation at the end of this function.
853  template <> unsigned int FE<0,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
854  template <> unsigned int FE<1,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
855  template <> unsigned int FE<2,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
856  template <> unsigned int FE<3,LAGRANGE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return lagrange_n_dofs_at_node(t, o, n); }
857 
858 
859  // Lagrange elements have no dofs per element
860  // (just at the nodes)
861  template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
862  template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
863  template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
864  template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
865 
866  // Lagrange FEMs are always C^0 continuous
867  template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; }
868  template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; }
869  template <> FEContinuity FE<2,LAGRANGE>::get_continuity() const { return C_ZERO; }
870  template <> FEContinuity FE<3,LAGRANGE>::get_continuity() const { return C_ZERO; }
871 
872  // Lagrange FEMs are not hierarchic
873  template <> bool FE<0,LAGRANGE>::is_hierarchic() const { return false; }
874  template <> bool FE<1,LAGRANGE>::is_hierarchic() const { return false; }
875  template <> bool FE<2,LAGRANGE>::is_hierarchic() const { return false; }
876  template <> bool FE<3,LAGRANGE>::is_hierarchic() const { return false; }
877 
878  // Lagrange FEM shapes do not need reinit (is this always true?)
879  template <> bool FE<0,LAGRANGE>::shapes_need_reinit() const { return false; }
880  template <> bool FE<1,LAGRANGE>::shapes_need_reinit() const { return false; }
881  template <> bool FE<2,LAGRANGE>::shapes_need_reinit() const { return false; }
882  template <> bool FE<3,LAGRANGE>::shapes_need_reinit() const { return false; }
883 
884  // Methods for computing Lagrange constraints. Note: we pass the
885  // dimension as the last argument to the anonymous helper function.
886  // Also note: we only need instantiations of this function for
887  // Dim==2 and 3.
888 #ifdef LIBMESH_ENABLE_AMR
889  template <>
891  DofMap &dof_map,
892  const unsigned int variable_number,
893  const Elem* elem)
894  { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
895 
896  template <>
898  DofMap &dof_map,
899  const unsigned int variable_number,
900  const Elem* elem)
901  { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
902 #endif // LIBMESH_ENABLE_AMR
903 
904 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo