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