1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 {
lagrange_nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)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
lagrange_n_dofs(const ElemType t,const Order o)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
lagrange_n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)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
lagrange_compute_constraints(DofConstraints & constraints,DofMap & dof_map,const unsigned int variable_number,const Elem * elem,const unsigned Dim)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
685 // Pull objects out of the loop to reduce heap operations
686 std::vector<dof_id_type> my_dof_indices, parent_dof_indices;
687 std::unique_ptr<const Elem> my_side, parent_side;
688
689 // Look at the element faces. Check to see if we need to
690 // build constraints.
691 for (auto s : elem->side_index_range())
692 if (elem->neighbor_ptr(s) != nullptr &&
693 elem->neighbor_ptr(s) != remote_elem)
694 if (elem->neighbor_ptr(s)->level() < elem->level()) // constrain dofs shared between
695 { // this element and ones coarser
696 // than this element.
697 // Get pointers to the elements of interest and its parent.
698 const Elem * parent = elem->parent();
699
700 // This can't happen... Only level-0 elements have nullptr
701 // parents, and no level-0 elements can be at a higher
702 // level than their neighbors!
703 libmesh_assert(parent);
704
705 elem->build_side_ptr(my_side, s);
706 parent->build_side_ptr(parent_side, s);
707
708 // This function gets called element-by-element, so there
709 // will be a lot of memory allocation going on. We can
710 // at least minimize this for the case of the dof indices
711 // by efficiently preallocating the requisite storage.
712 my_dof_indices.reserve (my_side->n_nodes());
713 parent_dof_indices.reserve (parent_side->n_nodes());
714
715 dof_map.dof_indices (my_side.get(), my_dof_indices,
716 variable_number);
717 dof_map.dof_indices (parent_side.get(), parent_dof_indices,
718 variable_number);
719
720 const unsigned int n_side_dofs =
721 FEInterface::n_dofs(fe_type, my_side.get());
722 const unsigned int n_parent_side_dofs =
723 FEInterface::n_dofs(fe_type, parent_side.get());
724 for (unsigned int my_dof=0; my_dof != n_side_dofs; my_dof++)
725 {
726 libmesh_assert_less (my_dof, my_side->n_nodes());
727
728 // My global dof index.
729 const dof_id_type my_dof_g = my_dof_indices[my_dof];
730
731 // Hunt for "constraining against myself" cases before
732 // we bother creating a constraint row
733 bool self_constraint = false;
734 for (unsigned int their_dof=0;
735 their_dof != n_parent_side_dofs; their_dof++)
736 {
737 libmesh_assert_less (their_dof, parent_side->n_nodes());
738
739 // Their global dof index.
740 const dof_id_type their_dof_g =
741 parent_dof_indices[their_dof];
742
743 if (their_dof_g == my_dof_g)
744 {
745 self_constraint = true;
746 break;
747 }
748 }
749
750 if (self_constraint)
751 continue;
752
753 DofConstraintRow * constraint_row;
754
755 // we may be running constraint methods concurrently
756 // on multiple threads, so we need a lock to
757 // ensure that this constraint is "ours"
758 {
759 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
760
761 if (dof_map.is_constrained_dof(my_dof_g))
762 continue;
763
764 constraint_row = &(constraints[my_dof_g]);
765 libmesh_assert(constraint_row->empty());
766 }
767
768 // The support point of the DOF
769 const Point & support_point = my_side->point(my_dof);
770
771 // Figure out where my node lies on their reference element.
772 const Point mapped_point = FEMap::inverse_map(Dim-1,
773 parent_side.get(),
774 support_point);
775
776 // Compute the parent's side shape function values.
777 for (unsigned int their_dof=0;
778 their_dof != n_parent_side_dofs; their_dof++)
779 {
780 libmesh_assert_less (their_dof, parent_side->n_nodes());
781
782 // Their global dof index.
783 const dof_id_type their_dof_g =
784 parent_dof_indices[their_dof];
785
786 const Real their_dof_value = FEInterface::shape(Dim-1,
787 fe_type,
788 parent_side.get(),
789 their_dof,
790 mapped_point);
791
792 // Only add non-zero and non-identity values
793 // for Lagrange basis functions.
794 if ((std::abs(their_dof_value) > 1.e-5) &&
795 (std::abs(their_dof_value) < .999))
796 {
797 constraint_row->emplace(their_dof_g, their_dof_value);
798 }
799 #ifdef DEBUG
800 // Protect for the case u_i = 0.999 u_j,
801 // in which case i better equal j.
802 else if (their_dof_value >= .999)
803 libmesh_assert_equal_to (my_dof_g, their_dof_g);
804 #endif
805 }
806 }
807 }
808 } // lagrange_compute_constraints()
809 #endif // #ifdef LIBMESH_ENABLE_AMR
810
811 } // anonymous namespace
812
813
814 // Do full-specialization for every dimension, instead
815 // of explicit instantiation at the end of this file.
816 // This could be macro-ified so that it fits on one line...
817 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)818 void FE<0,LAGRANGE>::nodal_soln(const Elem * elem,
819 const Order order,
820 const std::vector<Number> & elem_soln,
821 std::vector<Number> & nodal_soln)
822 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
823
824 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)825 void FE<1,LAGRANGE>::nodal_soln(const Elem * elem,
826 const Order order,
827 const std::vector<Number> & elem_soln,
828 std::vector<Number> & nodal_soln)
829 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
830
831 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)832 void FE<2,LAGRANGE>::nodal_soln(const Elem * elem,
833 const Order order,
834 const std::vector<Number> & elem_soln,
835 std::vector<Number> & nodal_soln)
836 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
837
838 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)839 void FE<3,LAGRANGE>::nodal_soln(const Elem * elem,
840 const Order order,
841 const std::vector<Number> & elem_soln,
842 std::vector<Number> & nodal_soln)
843 { lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
844
845
846 // Do full-specialization for every dimension, instead
847 // of explicit instantiation at the end of this function.
848 // This could be macro-ified.
n_dofs(const ElemType t,const Order o)849 template <> unsigned int FE<0,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
n_dofs(const ElemType t,const Order o)850 template <> unsigned int FE<1,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
n_dofs(const ElemType t,const Order o)851 template <> unsigned int FE<2,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
n_dofs(const ElemType t,const Order o)852 template <> unsigned int FE<3,LAGRANGE>::n_dofs(const ElemType t, const Order o) { return lagrange_n_dofs(t, o); }
853
854
855 // Do full-specialization for every dimension, instead
856 // of explicit instantiation at the end of this function.
n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)857 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); }
n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)858 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); }
n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)859 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); }
n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)860 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); }
861
862
863 // Lagrange elements have no dofs per element
864 // (just at the nodes)
n_dofs_per_elem(const ElemType,const Order)865 template <> unsigned int FE<0,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
n_dofs_per_elem(const ElemType,const Order)866 template <> unsigned int FE<1,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
n_dofs_per_elem(const ElemType,const Order)867 template <> unsigned int FE<2,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
n_dofs_per_elem(const ElemType,const Order)868 template <> unsigned int FE<3,LAGRANGE>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
869
870 // Lagrange FEMs are always C^0 continuous
get_continuity()871 template <> FEContinuity FE<0,LAGRANGE>::get_continuity() const { return C_ZERO; }
get_continuity()872 template <> FEContinuity FE<1,LAGRANGE>::get_continuity() const { return C_ZERO; }
get_continuity()873 template <> FEContinuity FE<2,LAGRANGE>::get_continuity() const { return C_ZERO; }
get_continuity()874 template <> FEContinuity FE<3,LAGRANGE>::get_continuity() const { return C_ZERO; }
875
876 // Lagrange FEMs are not hierarchic
is_hierarchic()877 template <> bool FE<0,LAGRANGE>::is_hierarchic() const { return false; }
is_hierarchic()878 template <> bool FE<1,LAGRANGE>::is_hierarchic() const { return false; }
is_hierarchic()879 template <> bool FE<2,LAGRANGE>::is_hierarchic() const { return false; }
is_hierarchic()880 template <> bool FE<3,LAGRANGE>::is_hierarchic() const { return false; }
881
882 // Lagrange FEM shapes do not need reinit (is this always true?)
shapes_need_reinit()883 template <> bool FE<0,LAGRANGE>::shapes_need_reinit() const { return false; }
shapes_need_reinit()884 template <> bool FE<1,LAGRANGE>::shapes_need_reinit() const { return false; }
shapes_need_reinit()885 template <> bool FE<2,LAGRANGE>::shapes_need_reinit() const { return false; }
shapes_need_reinit()886 template <> bool FE<3,LAGRANGE>::shapes_need_reinit() const { return false; }
887
888 // Methods for computing Lagrange constraints. Note: we pass the
889 // dimension as the last argument to the anonymous helper function.
890 // Also note: we only need instantiations of this function for
891 // Dim==2 and 3.
892 #ifdef LIBMESH_ENABLE_AMR
893 template <>
compute_constraints(DofConstraints & constraints,DofMap & dof_map,const unsigned int variable_number,const Elem * elem)894 void FE<2,LAGRANGE>::compute_constraints (DofConstraints & constraints,
895 DofMap & dof_map,
896 const unsigned int variable_number,
897 const Elem * elem)
898 { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/2); }
899
900 template <>
compute_constraints(DofConstraints & constraints,DofMap & dof_map,const unsigned int variable_number,const Elem * elem)901 void FE<3,LAGRANGE>::compute_constraints (DofConstraints & constraints,
902 DofMap & dof_map,
903 const unsigned int variable_number,
904 const Elem * elem)
905 { lagrange_compute_constraints(constraints, dof_map, variable_number, elem, /*Dim=*/3); }
906 #endif // LIBMESH_ENABLE_AMR
907
908 } // namespace libMesh
909