fe_lagrange.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/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 
466 
467 
468 
469 unsigned int lagrange_n_dofs_at_node(const ElemType t,
470  const Order o,
471  const unsigned int n)
472 {
473  switch (o)
474  {
475 
476  // linear Lagrange shape functions
477  case FIRST:
478  {
479  switch (t)
480  {
481  case NODEELEM:
482  return 1;
483 
484  case EDGE2:
485  case EDGE3:
486  case EDGE4:
487  {
488  switch (n)
489  {
490  case 0:
491  case 1:
492  return 1;
493 
494  default:
495  return 0;
496  }
497  }
498 
499  case TRI3:
500  case TRISHELL3:
501  case TRI3SUBDIVISION:
502  case TRI6:
503  {
504  switch (n)
505  {
506  case 0:
507  case 1:
508  case 2:
509  return 1;
510 
511  default:
512  return 0;
513  }
514  }
515 
516  case QUAD4:
517  case QUADSHELL4:
518  case QUAD8:
519  case QUADSHELL8:
520  case QUAD9:
521  {
522  switch (n)
523  {
524  case 0:
525  case 1:
526  case 2:
527  case 3:
528  return 1;
529 
530  default:
531  return 0;
532  }
533  }
534 
535 
536  case TET4:
537  case TET10:
538  {
539  switch (n)
540  {
541  case 0:
542  case 1:
543  case 2:
544  case 3:
545  return 1;
546 
547  default:
548  return 0;
549  }
550  }
551 
552  case HEX8:
553  case HEX20:
554  case HEX27:
555  {
556  switch (n)
557  {
558  case 0:
559  case 1:
560  case 2:
561  case 3:
562  case 4:
563  case 5:
564  case 6:
565  case 7:
566  return 1;
567 
568  default:
569  return 0;
570  }
571  }
572 
573  case PRISM6:
574  case PRISM15:
575  case PRISM18:
576  {
577  switch (n)
578  {
579  case 0:
580  case 1:
581  case 2:
582  case 3:
583  case 4:
584  case 5:
585  return 1;
586 
587  default:
588  return 0;
589  }
590  }
591 
592  case PYRAMID5:
593  case PYRAMID13:
594  case PYRAMID14:
595  {
596  switch (n)
597  {
598  case 0:
599  case 1:
600  case 2:
601  case 3:
602  case 4:
603  return 1;
604 
605  default:
606  return 0;
607  }
608  }
609 
610  case INVALID_ELEM:
611  return 0;
612 
613  default:
614  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
615  }
616  }
617 
618  // quadratic Lagrange shape functions
619  case SECOND:
620  {
621  switch (t)
622  {
623  // quadratic lagrange has one dof at each node
624  case NODEELEM:
625  case EDGE3:
626  case TRI6:
627  case QUAD8:
628  case QUADSHELL8:
629  case QUAD9:
630  case TET10:
631  case HEX20:
632  case HEX27:
633  case PRISM15:
634  case PRISM18:
635  case PYRAMID13:
636  case PYRAMID14:
637  return 1;
638 
639  case INVALID_ELEM:
640  return 0;
641 
642  default:
643  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
644  }
645  }
646 
647  case THIRD:
648  {
649  switch (t)
650  {
651  case NODEELEM:
652  case EDGE4:
653  return 1;
654 
655  case INVALID_ELEM:
656  return 0;
657 
658  default:
659  libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
660  }
661  }
662 
663  default:
664  libmesh_error_msg("Unsupported order: " << o );
665  }
666 }
667 
668 
669 
670 #ifdef LIBMESH_ENABLE_AMR
671 void lagrange_compute_constraints (DofConstraints & constraints,
672  DofMap & dof_map,
673  const unsigned int variable_number,
674  const Elem * elem,
675  const unsigned Dim)
676 {
677  // Only constrain elements in 2,3D.
678  if (Dim == 1)
679  return;
680 
681  libmesh_assert(elem);
682 
683  // Only constrain active and ancestor elements
684  if (elem->subactive())
685  return;
686 
687  FEType fe_type = dof_map.variable_type(variable_number);
688  fe_type.order = static_cast<Order>(fe_type.order + elem->p_level());
689 
690  // Pull objects out of the loop to reduce heap operations
691  std::vector<dof_id_type> my_dof_indices, parent_dof_indices;
692  std::unique_ptr<const Elem> my_side, parent_side;
693 
694  // Look at the element faces. Check to see if we need to
695  // build constraints.
696  for (auto s : elem->side_index_range())
697  if (elem->neighbor_ptr(s) != nullptr &&
698  elem->neighbor_ptr(s) != remote_elem)
699  if (elem->neighbor_ptr(s)->level() < elem->level()) // constrain dofs shared between
700  { // this element and ones coarser
701  // than this element.
702  // Get pointers to the elements of interest and its parent.
703  const Elem * parent = elem->parent();
704 
705  // This can't happen... Only level-0 elements have nullptr
706  // parents, and no level-0 elements can be at a higher
707  // level than their neighbors!
708  libmesh_assert(parent);
709 
710  elem->build_side_ptr(my_side, s);
711  parent->build_side_ptr(parent_side, s);
712 
713  // This function gets called element-by-element, so there
714  // will be a lot of memory allocation going on. We can
715  // at least minimize this for the case of the dof indices
716  // by efficiently preallocating the requisite storage.
717  my_dof_indices.reserve (my_side->n_nodes());
718  parent_dof_indices.reserve (parent_side->n_nodes());
719 
720  dof_map.dof_indices (my_side.get(), my_dof_indices,
721  variable_number);
722  dof_map.dof_indices (parent_side.get(), parent_dof_indices,
723  variable_number);
724 
725  const unsigned int n_side_dofs =
726  FEInterface::n_dofs(Dim-1, fe_type, my_side->type());
727  const unsigned int n_parent_side_dofs =
728  FEInterface::n_dofs(Dim-1, fe_type, parent_side->type());
729  for (unsigned int my_dof=0; my_dof != n_side_dofs; my_dof++)
730  {
731  libmesh_assert_less (my_dof, my_side->n_nodes());
732 
733  // My global dof index.
734  const dof_id_type my_dof_g = my_dof_indices[my_dof];
735 
736  // Hunt for "constraining against myself" cases before
737  // we bother creating a constraint row
738  bool self_constraint = false;
739  for (unsigned int their_dof=0;
740  their_dof != n_parent_side_dofs; their_dof++)
741  {
742  libmesh_assert_less (their_dof, parent_side->n_nodes());
743 
744  // Their global dof index.
745  const dof_id_type their_dof_g =
746  parent_dof_indices[their_dof];
747 
748  if (their_dof_g == my_dof_g)
749  {
750  self_constraint = true;
751  break;
752  }
753  }
754 
755  if (self_constraint)
756  continue;
757 
758  DofConstraintRow * constraint_row;
759 
760  // we may be running constraint methods concurrently
761  // on multiple threads, so we need a lock to
762  // ensure that this constraint is "ours"
763  {
764  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
765 
766  if (dof_map.is_constrained_dof(my_dof_g))
767  continue;
768 
769  constraint_row = &(constraints[my_dof_g]);
770  libmesh_assert(constraint_row->empty());
771  }
772 
773  // The support point of the DOF
774  const Point & support_point = my_side->point(my_dof);
775 
776  // Figure out where my node lies on their reference element.
777  const Point mapped_point = FEInterface::inverse_map(Dim-1, fe_type,
778  parent_side.get(),
779  support_point);
780 
781  // Compute the parent's side shape function values.
782  for (unsigned int their_dof=0;
783  their_dof != n_parent_side_dofs; their_dof++)
784  {
785  libmesh_assert_less (their_dof, parent_side->n_nodes());
786 
787  // Their global dof index.
788  const dof_id_type their_dof_g =
789  parent_dof_indices[their_dof];
790 
791  const Real their_dof_value = FEInterface::shape(Dim-1,
792  fe_type,
793  parent_side->type(),
794  their_dof,
795  mapped_point);
796 
797  // Only add non-zero and non-identity values
798  // for Lagrange basis functions.
799  if ((std::abs(their_dof_value) > 1.e-5) &&
800  (std::abs(their_dof_value) < .999))
801  {
802  constraint_row->insert(std::make_pair (their_dof_g,
803  their_dof_value));
804  }
805 #ifdef DEBUG
806  // Protect for the case u_i = 0.999 u_j,
807  // in which case i better equal j.
808  else if (their_dof_value >= .999)
809  libmesh_assert_equal_to (my_dof_g, their_dof_g);
810 #endif
811  }
812  }
813  }
814 } // lagrange_compute_constraints()
815 #endif // #ifdef LIBMESH_ENABLE_AMR
816 
817 } // anonymous namespace
818 
819 
820  // Do full-specialization for every dimension, instead
821  // of explicit instantiation at the end of this file.
822  // This could be macro-ified so that it fits on one line...
823 template <>
824 void FE<0,LAGRANGE>::nodal_soln(const Elem * elem,
825  const Order order,
826  const std::vector<Number> & elem_soln,
827  std::vector<Number> & nodal_soln)
828 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
829 
830 template <>
831 void FE<1,LAGRANGE>::nodal_soln(const Elem * elem,
832  const Order order,
833  const std::vector<Number> & elem_soln,
834  std::vector<Number> & nodal_soln)
835 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
836 
837 template <>
838 void FE<2,LAGRANGE>::nodal_soln(const Elem * elem,
839  const Order order,
840  const std::vector<Number> & elem_soln,
841  std::vector<Number> & nodal_soln)
842 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
843 
844 template <>
845 void FE<3,LAGRANGE>::nodal_soln(const Elem * elem,
846  const Order order,
847  const std::vector<Number> & elem_soln,
848  std::vector<Number> & nodal_soln)
849 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
850 
851 
852 // Do full-specialization for every dimension, instead
853 // of explicit instantiation at the end of this function.
854 // This could be macro-ified.
855 template <> unsigned int FE<0,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
856 template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
857 template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
858 template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
859 
860 
861 // Do full-specialization for every dimension, instead
862 // of explicit instantiation at the end of this function.
863 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); }
864 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); }
865 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); }
866 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); }
867 
868 
869 // Lagrange elements have no dofs per element
870 // (just at the nodes)
871 template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
872 template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
873 template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
874 template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
875 
876 // Lagrange FEMs are always C^0 continuous
877 template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; }
878 template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; }
879 template <> FEContinuity FE<2,LAGRANGE>::get_continuity() const { return C_ZERO; }
880 template <> FEContinuity FE<3,LAGRANGE>::get_continuity() const { return C_ZERO; }
881 
882 // Lagrange FEMs are not hierarchic
883 template <> bool FE<0,LAGRANGE>::is_hierarchic() const { return false; }
884 template <> bool FE<1,LAGRANGE>::is_hierarchic() const { return false; }
885 template <> bool FE<2,LAGRANGE>::is_hierarchic() const { return false; }
886 template <> bool FE<3,LAGRANGE>::is_hierarchic() const { return false; }
887 
888 // Lagrange FEM shapes do not need reinit (is this always true?)
889 template <> bool FE<0,LAGRANGE>::shapes_need_reinit() const { return false; }
890 template <> bool FE<1,LAGRANGE>::shapes_need_reinit() const { return false; }
891 template <> bool FE<2,LAGRANGE>::shapes_need_reinit() const { return false; }
892 template <> bool FE<3,LAGRANGE>::shapes_need_reinit() const { return false; }
893 
894 // Methods for computing Lagrange constraints. Note: we pass the
895 // dimension as the last argument to the anonymous helper function.
896 // Also note: we only need instantiations of this function for
897 // Dim==2 and 3.
898 #ifdef LIBMESH_ENABLE_AMR
899 template <>
900 void FE<2,LAGRANGE>::compute_constraints (DofConstraints & constraints,
901  DofMap & dof_map,
902  const unsigned int variable_number,
903  const Elem * elem)
904 { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
905 
906 template <>
907 void FE<3,LAGRANGE>::compute_constraints (DofConstraints & constraints,
908  DofMap & dof_map,
909  const unsigned int variable_number,
910  const Elem * elem)
911 { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
912 #endif // LIBMESH_ENABLE_AMR
913 
914 } // namespace libMesh
libMesh::FE::n_dofs
static unsigned int n_dofs(const ElemType t, const Order o)
std::abs
double abs(double a)
Definition: libmesh_augment_std_namespace.h:59
libMesh::FEInterface::n_dofs
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:454
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
libMesh::FE::n_dofs_at_node
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
remote_elem.h
fe_interface.h
libMesh::FE::shapes_need_reinit
virtual bool shapes_need_reinit() const override
libMesh::FE::is_hierarchic
virtual bool is_hierarchic() const override
libMesh::DofMap
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
n_nodes
const dof_id_type n_nodes
Definition: tecplot_io.C:68
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
Definition: threads.C:29
libMesh::FEInterface::shape
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:657
libMesh::FEInterface::inverse_map
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:590
libMesh::FE::n_dofs_per_elem
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
libMesh::FE::get_continuity
virtual FEContinuity get_continuity() const override
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::FE::compute_constraints
static void compute_constraints(DofConstraints &constraints, DofMap &dof_map, const unsigned int variable_number, const Elem *elem)
dof_map.h
string_to_enum.h
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::FE::nodal_soln
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
libMesh::DofConstraintRow
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
Definition: dof_map.h:97
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:64
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57
generated by doxygen