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