libMesh
Loading...
Searching...
No Matches
fe_base.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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/fe.h"
22#include "libmesh/inf_fe.h"
23#include "libmesh/libmesh_logging.h"
24
25// For projection code:
26#include "libmesh/boundary_info.h"
27#include "libmesh/mesh_base.h"
28#include "libmesh/dense_matrix.h"
29#include "libmesh/dense_vector.h"
30#include "libmesh/dof_map.h"
31#include "libmesh/elem.h"
32#include "libmesh/fe_interface.h"
33#include "libmesh/int_range.h"
34#include "libmesh/numeric_vector.h"
35#include "libmesh/periodic_boundary_base.h"
36#include "libmesh/periodic_boundaries.h"
37#include "libmesh/quadrature.h"
38#include "libmesh/quadrature_gauss.h"
39#include "libmesh/tensor_value.h"
40#include "libmesh/threads.h"
41#include "libmesh/fe_type.h"
42#include "libmesh/enum_to_string.h"
43
44// C++ Includes
45#include <memory>
46
47// Anonymous namespace, for a helper function for periodic boundary
48// constraint calculations
49namespace
50{
51using namespace libMesh;
52
53#ifdef LIBMESH_ENABLE_PERIODIC
54
55// Find the "primary" element around a boundary point:
56const Elem * primary_boundary_point_neighbor(const Elem * elem,
57 const Point & p,
58 const BoundaryInfo & boundary_info,
59 const std::set<boundary_id_type> & boundary_ids)
60{
61 // If we don't find a better alternative, the user will have
62 // provided the primary element
63 const Elem * primary = elem;
64
65 // Container to catch boundary IDs passed back by BoundaryInfo.
66 std::vector<boundary_id_type> bc_ids;
67
68 // Pull object out of the loop to reduce heap operations
69 std::unique_ptr<const Elem> periodic_side;
70
71 std::set<const Elem *> point_neighbors;
72 elem->find_point_neighbors(p, point_neighbors);
73 for (const auto & pt_neighbor : point_neighbors)
74 {
75 // If this point neighbor isn't at least
76 // as coarse as the current primary elem, or if it is at
77 // the same level but has a lower id, then
78 // we won't defer to it.
79 if ((pt_neighbor->level() > primary->level()) ||
80 (pt_neighbor->level() == primary->level() &&
81 pt_neighbor->id() < primary->id()))
82 continue;
83
84 // Otherwise, we will defer to the point neighbor, but only if
85 // one of its sides is on a relevant boundary and that side
86 // contains this vertex
87 bool vertex_on_periodic_side = false;
88 for (auto ns : pt_neighbor->side_index_range())
89 {
90 boundary_info.boundary_ids (pt_neighbor, ns, bc_ids);
91
92 bool on_relevant_boundary = false;
93 for (const auto & id : boundary_ids)
94 if (std::find(bc_ids.begin(), bc_ids.end(), id) != bc_ids.end())
95 on_relevant_boundary = true;
96
97 if (!on_relevant_boundary)
98 continue;
99
100 pt_neighbor->build_side_ptr(periodic_side, ns);
101 if (!periodic_side->contains_point(p))
102 continue;
103
104 vertex_on_periodic_side = true;
105 break;
106 }
107
108 if (vertex_on_periodic_side)
109 primary = pt_neighbor;
110 }
111
112 return primary;
113}
114
115// Find the "primary" element around a boundary edge:
116const Elem * primary_boundary_edge_neighbor(const Elem * elem,
117 const Point & p1,
118 const Point & p2,
119 const BoundaryInfo & boundary_info,
120 const std::set<boundary_id_type> & boundary_ids)
121{
122 // If we don't find a better alternative, the user will have
123 // provided the primary element
124 const Elem * primary = elem;
125
126 std::set<const Elem *> edge_neighbors;
127 elem->find_edge_neighbors(p1, p2, edge_neighbors);
128
129 // Container to catch boundary IDs handed back by BoundaryInfo
130 std::vector<boundary_id_type> bc_ids;
131
132 // Pull object out of the loop to reduce heap operations
133 std::unique_ptr<const Elem> periodic_side;
134
135 for (const auto & e_neighbor : edge_neighbors)
136 {
137 // If this edge neighbor isn't at least
138 // as coarse as the current primary elem, or if it is at
139 // the same level but has a lower id, then
140 // we won't defer to it.
141 if ((e_neighbor->level() > primary->level()) ||
142 (e_neighbor->level() == primary->level() &&
143 e_neighbor->id() < primary->id()))
144 continue;
145
146 // Otherwise, we will defer to the edge neighbor, but only if
147 // one of its sides is on this periodic boundary and that
148 // side contains this edge
149 bool vertex_on_periodic_side = false;
150 for (auto ns : e_neighbor->side_index_range())
151 {
152 boundary_info.boundary_ids (e_neighbor, ns, bc_ids);
153
154 bool on_relevant_boundary = false;
155 for (const auto & id : boundary_ids)
156 if (std::find(bc_ids.begin(), bc_ids.end(), id) != bc_ids.end())
157 on_relevant_boundary = true;
158
159 if (!on_relevant_boundary)
160 continue;
161
162 e_neighbor->build_side_ptr(periodic_side, ns);
163 if (!(periodic_side->contains_point(p1) &&
164 periodic_side->contains_point(p2)))
165 continue;
166
167 vertex_on_periodic_side = true;
168 break;
169 }
170
171 if (vertex_on_periodic_side)
172 primary = e_neighbor;
173 }
174
175 return primary;
176}
177
178#endif // LIBMESH_ENABLE_PERIODIC
179
180}
181
182namespace libMesh
183{
184
185
186
187// ------------------------------------------------------------
188// FEBase class members
189template <>
190std::unique_ptr<FEGenericBase<Real>>
191FEGenericBase<Real>::build (const unsigned int dim,
192 const FEType & fet)
193{
194 switch (dim)
195 {
196 // 0D
197 case 0:
198 {
199 switch (fet.family)
200 {
201 case CLOUGH:
202 return std::make_unique<FE<0,CLOUGH>>(fet);
203
204 case HERMITE:
205 return std::make_unique<FE<0,HERMITE>>(fet);
206
207 case LAGRANGE:
208 return std::make_unique<FE<0,LAGRANGE>>(fet);
209
210 case L2_LAGRANGE:
211 return std::make_unique<FE<0,L2_LAGRANGE>>(fet);
212
213 case L2_LAGRANGE_GLL:
214 return std::make_unique<FE<0,L2_LAGRANGE_GLL>>(fet);
215
216 case HIERARCHIC:
217 return std::make_unique<FE<0,HIERARCHIC>>(fet);
218
219 case L2_HIERARCHIC:
220 return std::make_unique<FE<0,L2_HIERARCHIC>>(fet);
221
222 case SIDE_HIERARCHIC:
223 return std::make_unique<FE<0,SIDE_HIERARCHIC>>(fet);
224
225 case MONOMIAL:
226 return std::make_unique<FE<0,MONOMIAL>>(fet);
227
228#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
229 case SZABAB:
230 return std::make_unique<FE<0,SZABAB>>(fet);
231
232 case BERNSTEIN:
233 return std::make_unique<FE<0,BERNSTEIN>>(fet);
234
235 case RATIONAL_BERNSTEIN:
236 return std::make_unique<FE<0,RATIONAL_BERNSTEIN>>(fet);
237#endif
238
239 case XYZ:
240 return std::make_unique<FEXYZ<0>>(fet);
241
242 case SCALAR:
243 return std::make_unique<FEScalar<0>>(fet);
244
245 default:
246 libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
247 }
248 }
249 // 1D
250 case 1:
251 {
252 switch (fet.family)
253 {
254 case CLOUGH:
255 return std::make_unique<FE<1,CLOUGH>>(fet);
256
257 case HERMITE:
258 return std::make_unique<FE<1,HERMITE>>(fet);
259
260 case LAGRANGE:
261 return std::make_unique<FE<1,LAGRANGE>>(fet);
262
263 case L2_LAGRANGE:
264 return std::make_unique<FE<1,L2_LAGRANGE>>(fet);
265
266 case L2_LAGRANGE_GLL:
267 return std::make_unique<FE<1,L2_LAGRANGE_GLL>>(fet);
268
269 case HIERARCHIC:
270 return std::make_unique<FE<1,HIERARCHIC>>(fet);
271
272 case L2_HIERARCHIC:
273 return std::make_unique<FE<1,L2_HIERARCHIC>>(fet);
274
275 case SIDE_HIERARCHIC:
276 return std::make_unique<FE<1,SIDE_HIERARCHIC>>(fet);
277
278 case MONOMIAL:
279 return std::make_unique<FE<1,MONOMIAL>>(fet);
280
281#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
282 case SZABAB:
283 return std::make_unique<FE<1,SZABAB>>(fet);
284
285 case BERNSTEIN:
286 return std::make_unique<FE<1,BERNSTEIN>>(fet);
287
288 case RATIONAL_BERNSTEIN:
289 return std::make_unique<FE<1,RATIONAL_BERNSTEIN>>(fet);
290#endif
291
292 case XYZ:
293 return std::make_unique<FEXYZ<1>>(fet);
294
295 case SCALAR:
296 return std::make_unique<FEScalar<1>>(fet);
297
298 default:
299 libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
300 }
301 }
302
303
304 // 2D
305 case 2:
306 {
307 switch (fet.family)
308 {
309 case CLOUGH:
310 return std::make_unique<FE<2,CLOUGH>>(fet);
311
312 case HERMITE:
313 return std::make_unique<FE<2,HERMITE>>(fet);
314
315 case LAGRANGE:
316 return std::make_unique<FE<2,LAGRANGE>>(fet);
317
318 case L2_LAGRANGE:
319 return std::make_unique<FE<2,L2_LAGRANGE>>(fet);
320
321 case L2_LAGRANGE_GLL:
322 return std::make_unique<FE<2,L2_LAGRANGE_GLL>>(fet);
323
324 case HIERARCHIC:
325 return std::make_unique<FE<2,HIERARCHIC>>(fet);
326
327 case L2_HIERARCHIC:
328 return std::make_unique<FE<2,L2_HIERARCHIC>>(fet);
329
330 case SIDE_HIERARCHIC:
331 return std::make_unique<FE<2,SIDE_HIERARCHIC>>(fet);
332
333 case MONOMIAL:
334 return std::make_unique<FE<2,MONOMIAL>>(fet);
335
336#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
337 case SZABAB:
338 return std::make_unique<FE<2,SZABAB>>(fet);
339
340 case BERNSTEIN:
341 return std::make_unique<FE<2,BERNSTEIN>>(fet);
342
343 case RATIONAL_BERNSTEIN:
344 return std::make_unique<FE<2,RATIONAL_BERNSTEIN>>(fet);
345#endif
346
347 case XYZ:
348 return std::make_unique<FEXYZ<2>>(fet);
349
350 case SCALAR:
351 return std::make_unique<FEScalar<2>>(fet);
352
353 case SUBDIVISION:
354 return std::make_unique<FESubdivision>(fet);
355
356 default:
357 libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
358 }
359 }
360
361
362 // 3D
363 case 3:
364 {
365 switch (fet.family)
366 {
367 case CLOUGH:
368 libmesh_error_msg("ERROR: Clough-Tocher elements currently only support 1D and 2D");
369
370 case HERMITE:
371 return std::make_unique<FE<3,HERMITE>>(fet);
372
373 case LAGRANGE:
374 return std::make_unique<FE<3,LAGRANGE>>(fet);
375
376 case L2_LAGRANGE:
377 return std::make_unique<FE<3,L2_LAGRANGE>>(fet);
378
379 case L2_LAGRANGE_GLL:
380 return std::make_unique<FE<3,L2_LAGRANGE_GLL>>(fet);
381
382 case HIERARCHIC:
383 return std::make_unique<FE<3,HIERARCHIC>>(fet);
384
385 case L2_HIERARCHIC:
386 return std::make_unique<FE<3,L2_HIERARCHIC>>(fet);
387
388 case SIDE_HIERARCHIC:
389 return std::make_unique<FE<3,SIDE_HIERARCHIC>>(fet);
390
391 case MONOMIAL:
392 return std::make_unique<FE<3,MONOMIAL>>(fet);
393
394#ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
395 case SZABAB:
396 return std::make_unique<FE<3,SZABAB>>(fet);
397
398 case BERNSTEIN:
399 return std::make_unique<FE<3,BERNSTEIN>>(fet);
400
401 case RATIONAL_BERNSTEIN:
402 return std::make_unique<FE<3,RATIONAL_BERNSTEIN>>(fet);
403#endif
404
405 case XYZ:
406 return std::make_unique<FEXYZ<3>>(fet);
407
408 case SCALAR:
409 return std::make_unique<FEScalar<3>>(fet);
410
411 default:
412 libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
413 }
414 }
415
416 default:
417 libmesh_error_msg("Invalid dimension dim = " << dim);
418 }
419}
420
421
422
423template <>
424std::unique_ptr<FEGenericBase<RealGradient>>
425FEGenericBase<RealGradient>::build (const unsigned int dim,
426 const FEType & fet)
427{
428 switch (dim)
429 {
430 // 0D
431 case 0:
432 {
433 switch (fet.family)
434 {
435 case HIERARCHIC_VEC:
436 return std::make_unique<FEHierarchicVec<0>>(fet);
437
438 case L2_HIERARCHIC_VEC:
439 return std::make_unique<FEL2HierarchicVec<0>>(fet);
440
441 case LAGRANGE_VEC:
442 return std::make_unique<FELagrangeVec<0>>(fet);
443
444 case L2_LAGRANGE_VEC:
445 return std::make_unique<FEL2LagrangeVec<0>>(fet);
446
447 case MONOMIAL_VEC:
448 return std::make_unique<FEMonomialVec<0>>(fet);
449
450 default:
451 libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
452 }
453 }
454 case 1:
455 {
456 switch (fet.family)
457 {
458 case HIERARCHIC_VEC:
459 return std::make_unique<FEHierarchicVec<1>>(fet);
460
461 case L2_HIERARCHIC_VEC:
462 return std::make_unique<FEL2HierarchicVec<1>>(fet);
463
464 case LAGRANGE_VEC:
465 return std::make_unique<FELagrangeVec<1>>(fet);
466
467 case L2_LAGRANGE_VEC:
468 return std::make_unique<FEL2LagrangeVec<1>>(fet);
469
470 case MONOMIAL_VEC:
471 return std::make_unique<FEMonomialVec<1>>(fet);
472
473 default:
474 libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
475 }
476 }
477 case 2:
478 {
479 switch (fet.family)
480 {
481 case HIERARCHIC_VEC:
482 return std::make_unique<FEHierarchicVec<2>>(fet);
483
484 case L2_HIERARCHIC_VEC:
485 return std::make_unique<FEL2HierarchicVec<2>>(fet);
486
487 case LAGRANGE_VEC:
488 return std::make_unique<FELagrangeVec<2>>(fet);
489
490 case L2_LAGRANGE_VEC:
491 return std::make_unique<FEL2LagrangeVec<2>>(fet);
492
493 case MONOMIAL_VEC:
494 return std::make_unique<FEMonomialVec<2>>(fet);
495
496 case NEDELEC_ONE:
497 return std::make_unique<FENedelecOne<2>>(fet);
498
499 case RAVIART_THOMAS:
500 return std::make_unique<FERaviartThomas<2>>(fet);
501
502 case L2_RAVIART_THOMAS:
503 return std::make_unique<FEL2RaviartThomas<2>>(fet);
504
505 default:
506 libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
507 }
508 }
509 case 3:
510 {
511 switch (fet.family)
512 {
513 case HIERARCHIC_VEC:
514 return std::make_unique<FEHierarchicVec<3>>(fet);
515
516 case L2_HIERARCHIC_VEC:
517 return std::make_unique<FEL2HierarchicVec<3>>(fet);
518
519 case LAGRANGE_VEC:
520 return std::make_unique<FELagrangeVec<3>>(fet);
521
522 case L2_LAGRANGE_VEC:
523 return std::make_unique<FEL2LagrangeVec<3>>(fet);
524
525 case MONOMIAL_VEC:
526 return std::make_unique<FEMonomialVec<3>>(fet);
527
528 case NEDELEC_ONE:
529 return std::make_unique<FENedelecOne<3>>(fet);
530
531 case RAVIART_THOMAS:
532 return std::make_unique<FERaviartThomas<3>>(fet);
533
534 case L2_RAVIART_THOMAS:
535 return std::make_unique<FEL2RaviartThomas<3>>(fet);
536
537 default:
538 libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
539 }
540 }
541
542 default:
543 libmesh_error_msg("Invalid dimension dim = " << dim);
544 } // switch(dim)
545}
546
547
548
549
550
551
552
553#ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
554
555
556template <>
557std::unique_ptr<FEGenericBase<Real>>
558FEGenericBase<Real>::build_InfFE (const unsigned int dim,
559 const FEType & fet)
560{
561 switch (dim)
562 {
563
564 // 1D
565 case 1:
566 {
567 switch (fet.radial_family)
568 {
569 case INFINITE_MAP:
570 libmesh_error_msg("ERROR: Can't build an infinite element with FEFamily = " << Utility::enum_to_string(fet.radial_family));
571
572 case JACOBI_20_00:
573 {
574 switch (fet.inf_map)
575 {
576 case CARTESIAN:
577 return std::make_unique<InfFE<1,JACOBI_20_00,CARTESIAN>>(fet);
578
579 default:
580 libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
581 }
582 }
583
584 case JACOBI_30_00:
585 {
586 switch (fet.inf_map)
587 {
588 case CARTESIAN:
589 return std::make_unique<InfFE<1,JACOBI_30_00,CARTESIAN>>(fet);
590
591 default:
592 libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
593 }
594 }
595
596 case LEGENDRE:
597 {
598 switch (fet.inf_map)
599 {
600 case CARTESIAN:
601 return std::make_unique<InfFE<1,LEGENDRE,CARTESIAN>>(fet);
602
603 default:
604 libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
605 }
606 }
607
608 case LAGRANGE:
609 {
610 switch (fet.inf_map)
611 {
612 case CARTESIAN:
613 return std::make_unique<InfFE<1,LAGRANGE,CARTESIAN>>(fet);
614
615 default:
616 libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
617 }
618 }
619
620 default:
621 libmesh_error_msg("ERROR: Bad FEType.radial_family= " << Utility::enum_to_string(fet.radial_family));
622 }
623 }
624
625
626
627
628 // 2D
629 case 2:
630 {
631 switch (fet.radial_family)
632 {
633 case INFINITE_MAP:
634 libmesh_error_msg("ERROR: Can't build an infinite element with FEFamily = " << Utility::enum_to_string(fet.radial_family));
635
636 case JACOBI_20_00:
637 {
638 switch (fet.inf_map)
639 {
640 case CARTESIAN:
641 return std::make_unique<InfFE<2,JACOBI_20_00,CARTESIAN>>(fet);
642
643 default:
644 libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
645 }
646 }
647
648 case JACOBI_30_00:
649 {
650 switch (fet.inf_map)
651 {
652 case CARTESIAN:
653 return std::make_unique<InfFE<2,JACOBI_30_00,CARTESIAN>>(fet);
654
655 default:
656 libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
657 }
658 }
659
660 case LEGENDRE:
661 {
662 switch (fet.inf_map)
663 {
664 case CARTESIAN:
665 return std::make_unique<InfFE<2,LEGENDRE,CARTESIAN>>(fet);
666
667 default:
668 libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
669 }
670 }
671
672 case LAGRANGE:
673 {
674 switch (fet.inf_map)
675 {
676 case CARTESIAN:
677 return std::make_unique<InfFE<2,LAGRANGE,CARTESIAN>>(fet);
678
679 default:
680 libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
681 }
682 }
683
684 default:
685 libmesh_error_msg("ERROR: Bad FEType.radial_family= " << Utility::enum_to_string(fet.radial_family));
686 }
687 }
688
689
690
691
692 // 3D
693 case 3:
694 {
695 switch (fet.radial_family)
696 {
697 case INFINITE_MAP:
698 libmesh_error_msg("ERROR: Don't build an infinite element with FEFamily = " << Utility::enum_to_string(fet.radial_family));
699
700 case JACOBI_20_00:
701 {
702 switch (fet.inf_map)
703 {
704 case CARTESIAN:
705 return std::make_unique<InfFE<3,JACOBI_20_00,CARTESIAN>>(fet);
706
707 default:
708 libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
709 }
710 }
711
712 case JACOBI_30_00:
713 {
714 switch (fet.inf_map)
715 {
716 case CARTESIAN:
717 return std::make_unique<InfFE<3,JACOBI_30_00,CARTESIAN>>(fet);
718
719 default:
720 libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
721 }
722 }
723
724 case LEGENDRE:
725 {
726 switch (fet.inf_map)
727 {
728 case CARTESIAN:
729 return std::make_unique<InfFE<3,LEGENDRE,CARTESIAN>>(fet);
730
731 default:
732 libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
733 }
734 }
735
736 case LAGRANGE:
737 {
738 switch (fet.inf_map)
739 {
740 case CARTESIAN:
741 return std::make_unique<InfFE<3,LAGRANGE,CARTESIAN>>(fet);
742
743 default:
744 libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
745 }
746 }
747
748 default:
749 libmesh_error_msg("ERROR: Bad FEType.radial_family= " << Utility::enum_to_string(fet.radial_family));
750 }
751 }
752
753 default:
754 libmesh_error_msg("Invalid dimension dim = " << dim);
755 }
756}
757
758
759
760template <>
761std::unique_ptr<FEGenericBase<RealGradient>>
762FEGenericBase<RealGradient>::build_InfFE (const unsigned int,
763 const FEType & )
764{
765 // No vector types defined... YET.
766 libmesh_not_implemented();
767 return std::unique_ptr<FEVectorBase>();
768}
769
770#endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
771
772
773template <typename OutputType>
774void FEGenericBase<OutputType>::compute_shape_functions (const Elem * elem,
775 const std::vector<Point> & qp)
776{
777 //-------------------------------------------------------------------------
778 // Compute the shape function values (and derivatives)
779 // at the Quadrature points. Note that the actual values
780 // have already been computed via init_shape_functions
781
782 // Start logging the shape function computation
783 LOG_ONE_ENTITY_SCOPE("compute_shape_functions()", "FE");
784
785 this->determine_calculations();
786
787 if (calculate_phi)
788 this->_fe_trans->map_phi(this->dim, elem, qp, (*this), this->phi, this->_add_p_level_in_reinit);
789
790 if (calculate_dphi)
791 this->_fe_trans->map_dphi(this->dim, elem, qp, (*this), this->dphi,
792 this->dphidx, this->dphidy, this->dphidz);
793
794#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
795 if (calculate_d2phi)
796 this->_fe_trans->map_d2phi(this->dim, qp, (*this), this->d2phi,
797 this->d2phidx2, this->d2phidxdy, this->d2phidxdz,
798 this->d2phidy2, this->d2phidydz, this->d2phidz2);
799#endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
800
801 // Only compute curl for vector-valued elements
802 if (calculate_curl_phi && TypesEqual<OutputType,RealGradient>::value)
803 this->_fe_trans->map_curl(this->dim, elem, qp, (*this), this->curl_phi);
804
805 // Only compute div for vector-valued elements
806 if (calculate_div_phi && TypesEqual<OutputType,RealGradient>::value)
807 this->_fe_trans->map_div(this->dim, elem, qp, (*this), this->div_phi);
808}
809
810
811// Here, we rely on the input \p phi_vals for accurate integration of the mass matrix.
812// This is because in contact, we often have customized qrule for mortar segments
813// due to deformation of the element, and the size of \p phi_vals for the secondary
814// element changes accordingly.
815template <>
816void FEGenericBase<Real>::compute_dual_shape_coeffs (const std::vector<Real> & JxW, const std::vector<std::vector<OutputShape>> & phi_vals)
817{
818 // Start logging the dual coeff computation
819 LOG_SCOPE("compute_dual_shape_coeffs()", "FE");
820
821 const unsigned int sz=phi_vals.size();
822 libmesh_error_msg_if(!sz, "ERROR: cannot compute dual shape coefficients with empty phi values");
823
824 //compute dual basis coefficient (dual_coeff)
825 dual_coeff.resize(sz, sz);
826 DenseMatrix<Real> A(sz, sz), D(sz, sz);
827
828 for (const auto i : index_range(phi_vals))
829 for (const auto qp : index_range(phi_vals[i]))
830 {
831 D(i,i) += JxW[qp]*phi_vals[i][qp];
832 for (const auto j : index_range(phi_vals))
833 A(i,j) += JxW[qp]*phi_vals[i][qp]*phi_vals[j][qp];
834 }
835
836 // The integral positivity condition on the dual basis (Popp, Wohlmuth, Gee and Wall, SIAM J. Sci.
837 // Comput. 34(4):B421-B446, 2012, Eq. (4.2)) fails on QUAD8 and TRI6: D(k,k) = \int N_k is -1/3 at a
838 // QUAD8 corner and exactly 0 at a TRI6 vertex (their Eqs. (4.3) and (4.4)), the latter leaving those
839 // dual shape functions identically zero. Biorthogonalize instead against the locally quadratic
840 // transformed basis Ntilde = T N of their Sec. 4.4.1, in which each vertex absorbs a fraction alpha
841 // of its adjacent mid-edge shapes. That replaces D by T^-1 diag(T d), which for this T is the sparse
842 // update below. alpha = 1/5 is recommended there, and makes the weights strictly positive (QUAD8
843 // 1/5 and 4/5, TRI6 1/15 and 1/10) while preserving the partition of unity. T = I elsewhere.
844 // SECOND order on TRI7 is the TRI6 basis, with no degree of freedom at the centroid node, and
845 // L2_LAGRANGE shares the LAGRANGE shape functions.
846 if (_elem && (get_family() == LAGRANGE || get_family() == L2_LAGRANGE) &&
847 (((_elem->type() == TRI6 || _elem->type() == TRI7) && sz == 6) ||
848 (_elem->type() == QUAD8 && sz == 8)))
849 {
850 const Real alpha = Real(1)/5;
851 // Mid-edge nodes are the trailing indices, and only vertex entries are written, so D(m,m) here
852 // is never a value an earlier iteration modified.
853 for (const auto m : make_range(_elem->n_vertices(), sz))
854 for (const auto v : make_range(_elem->n_second_order_adjacent_vertices(m)))
855 {
856 const auto vertex = _elem->second_order_adjacent_vertex(m, v);
857 D(vertex, vertex) += alpha*D(m,m);
858 D(vertex, m) -= alpha*D(m,m);
859 }
860 }
861
862 // dual_coeff = A^-1*D
863 for (const auto j : index_range(phi_vals))
864 {
865 DenseVector<Real> Dcol(sz), coeffcol(sz);
866 for (const auto i : index_range(phi_vals))
867 Dcol(i) = D(i, j);
868 A.cholesky_solve(Dcol, coeffcol);
869
870 for (const auto row : index_range(phi_vals))
871 dual_coeff(row, j)=coeffcol(row);
872 }
873}
874
875template <>
876void FEGenericBase<Real>::compute_dual_shape_functions ()
877{
878 // Start logging the shape function computation
879 LOG_SCOPE("compute_dual_shape_functions()", "FE");
880
881 // The dual coeffs matrix should have the same size as phi
882 libmesh_assert(dual_coeff.m() == phi.size());
883 libmesh_assert(dual_coeff.n() == phi.size());
884
885 // initialize dual basis
886 for (const auto j : index_range(phi))
887 for (const auto qp : index_range(phi[j]))
888 {
889 dual_phi[j][qp] = 0;
890 if (calculate_dphi)
891 dual_dphi[j][qp] = 0;
892#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
893 if (calculate_d2phi)
894 dual_d2phi[j][qp] = 0;
895#endif
896 }
897
898 // compute dual basis
899 for (const auto j : index_range(phi))
900 for (const auto i : index_range(phi))
901 for (const auto qp : index_range(phi[j]))
902 {
903 dual_phi[j][qp] += dual_coeff(i, j) * phi[i][qp];
904 if (calculate_dphi)
905 dual_dphi[j][qp] += dual_coeff(i, j) * dphi[i][qp];
906#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
907 if (calculate_d2phi)
908 dual_d2phi[j][qp] += dual_coeff(i, j) * d2phi[i][qp];
909#endif
910 }
911}
912
913template <typename OutputType>
914void FEGenericBase<OutputType>::print_phi(std::ostream & os) const
915{
916 for (auto i : index_range(phi))
917 for (auto j : index_range(phi[i]))
918 os << " phi[" << i << "][" << j << "]=" << phi[i][j] << std::endl;
919}
920
921template <typename OutputType>
922void FEGenericBase<OutputType>::print_dual_phi(std::ostream & os) const
923{
924 for (auto i : index_range(dual_phi))
925 for (auto j : index_range(dual_phi[i]))
926 os << " dual_phi[" << i << "][" << j << "]=" << dual_phi[i][j] << std::endl;
927}
928
929
930
931
932template <typename OutputType>
933void FEGenericBase<OutputType>::print_dphi(std::ostream & os) const
934{
935 for (auto i : index_range(dphi))
936 for (auto j : index_range(dphi[i]))
937 os << " dphi[" << i << "][" << j << "]=" << dphi[i][j];
938}
939
940template <typename OutputType>
941void FEGenericBase<OutputType>::print_dual_dphi(std::ostream & os) const
942{
943 for (auto i : index_range(dphi))
944 for (auto j : index_range(dphi[i]))
945 os << " dual_dphi[" << i << "][" << j << "]=" << dual_dphi[i][j];
946}
947
948
949
950template <typename OutputType>
951void FEGenericBase<OutputType>::determine_calculations()
952{
953 this->calculations_started = true;
954
955 // If the user did not explicitly pre-request something (or nothing)
956 // to be computed, then we throw an error here.
957 bool requested_ok =
958 this->calculate_nothing || this->calculate_phi || this->calculate_dphi ||
959 this->calculate_dphiref || this->calculate_curl_phi || this->calculate_div_phi ||
960 this->calculate_map;
961
962#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
963 requested_ok = requested_ok || this->calculate_d2phi;
964#endif
965
966 libmesh_error_msg_if(
967 !requested_ok,
968 "You must call one or more of the FE accessors "
969 "(e.g. get_phi(), get_dphi(), get_nothing()) "
970 "_before_ calling reinit()!");
971
972 // Request whichever terms are necessary from the FEMap
973 if (this->calculate_phi)
974 this->_fe_trans->init_map_phi(*this);
975
976 if (this->calculate_dphiref)
977 this->_fe_trans->init_map_dphi(*this);
978
979#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
980 if (this->calculate_d2phi)
981 this->_fe_trans->init_map_d2phi(*this);
982#endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
983}
984
985
986
987#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
988
989
990template <typename OutputType>
991void FEGenericBase<OutputType>::print_d2phi(std::ostream & os) const
992{
993 for (auto i : index_range(dphi))
994 for (auto j : index_range(dphi[i]))
995 os << " d2phi[" << i << "][" << j << "]=" << d2phi[i][j];
996}
997
998template <typename OutputType>
999void FEGenericBase<OutputType>::print_dual_d2phi(std::ostream & os) const
1000{
1001 for (auto i : index_range(dual_d2phi))
1002 for (auto j : index_range(dual_d2phi[i]))
1003 os << " dual_d2phi[" << i << "][" << j << "]=" << dual_d2phi[i][j];
1004}
1005
1006#endif
1007
1008
1009
1010#ifdef LIBMESH_ENABLE_AMR
1011
1012template <typename OutputType>
1013void
1014FEGenericBase<OutputType>::coarsened_dof_values(const NumericVector<Number> & old_vector,
1015 const DofMap & dof_map,
1016 const Elem * elem,
1017 DenseVector<Number> & Ue,
1018 const unsigned int var,
1019 const bool use_old_dof_indices)
1020{
1021 // Side/edge local DOF indices
1022 std::vector<unsigned int> new_side_dofs, old_side_dofs;
1023
1024 // FIXME: what about 2D shells in 3D space?
1025 unsigned int dim = elem->dim();
1026
1027 // Cache n_children(); it's a virtual call but it's const.
1028 const unsigned int n_children = elem->n_children();
1029
1030 // We use local FE objects for now
1031 // FIXME: we should use more, external objects instead for efficiency
1032 const FEType & base_fe_type = dof_map.variable_type(var);
1033 std::unique_ptr<FEGenericBase<OutputShape>> fe
1034 (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1035 std::unique_ptr<FEGenericBase<OutputShape>> fe_coarse
1036 (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1037
1038 std::unique_ptr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
1039 std::unique_ptr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
1040 std::unique_ptr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
1041 std::vector<Point> coarse_qpoints;
1042
1043 // The values of the shape functions at the quadrature
1044 // points
1045 const std::vector<std::vector<OutputShape>> & phi_values =
1046 fe->get_phi();
1047 const std::vector<std::vector<OutputShape>> & phi_coarse =
1048 fe_coarse->get_phi();
1049
1050 // The gradients of the shape functions at the quadrature
1051 // points on the child element.
1052 const std::vector<std::vector<OutputGradient>> * dphi_values =
1053 nullptr;
1054 const std::vector<std::vector<OutputGradient>> * dphi_coarse =
1055 nullptr;
1056
1057 const FEContinuity cont = fe->get_continuity();
1058
1059 if (cont == C_ONE)
1060 {
1061 const std::vector<std::vector<OutputGradient>> &
1062 ref_dphi_values = fe->get_dphi();
1063 dphi_values = &ref_dphi_values;
1064 const std::vector<std::vector<OutputGradient>> &
1065 ref_dphi_coarse = fe_coarse->get_dphi();
1066 dphi_coarse = &ref_dphi_coarse;
1067 }
1068
1069 // The Jacobian * quadrature weight at the quadrature points
1070 const std::vector<Real> & JxW =
1071 fe->get_JxW();
1072
1073 // The XYZ locations of the quadrature points on the
1074 // child element
1075 const std::vector<Point> & xyz_values =
1076 fe->get_xyz();
1077
1078 // Number of nodes on parent element
1079 const unsigned int n_nodes = elem->n_nodes();
1080
1081 // Number of dofs on parent element
1082 const unsigned int new_n_dofs =
1083 FEInterface::n_dofs(base_fe_type, elem->max_descendant_p_level(), elem);
1084
1085 // Fixed vs. free DoFs on edge/face projections
1086 std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
1087 std::vector<int> free_dof(new_n_dofs, 0);
1088
1089 DenseMatrix<Real> Ke;
1090 DenseVector<Number> Fe;
1091 Ue.resize(new_n_dofs); Ue.zero();
1092
1093
1094 // When coarsening, in general, we need a series of
1095 // projections to ensure a unique and continuous
1096 // solution. We start by interpolating nodes, then
1097 // hold those fixed and project edges, then
1098 // hold those fixed and project faces, then
1099 // hold those fixed and project interiors
1100
1101 // Copy node values first
1102 {
1103 std::vector<dof_id_type> node_dof_indices;
1104 if (use_old_dof_indices)
1105 dof_map.old_dof_indices (elem, node_dof_indices, var);
1106 else
1107 dof_map.dof_indices (elem, node_dof_indices, var);
1108
1109 unsigned int current_dof = 0;
1110 for (unsigned int n=0; n!= n_nodes; ++n)
1111 {
1112 // FIXME: this should go through the DofMap,
1113 // not duplicate dof_indices code badly!
1114 const unsigned int my_nc =
1115 FEInterface::n_dofs_at_node (base_fe_type, elem->max_descendant_p_level(), elem, n);
1116 if (!elem->is_vertex(n))
1117 {
1118 current_dof += my_nc;
1119 continue;
1120 }
1121
1122 // We're assuming here that child n shares vertex n,
1123 // which is wrong on non-simplices right now
1124 // ... but this code isn't necessary except on elements
1125 // where p refinement creates more vertex dofs; we have
1126 // no such elements yet.
1127 int extra_order = 0;
1128 // if (elem->child_ptr(n)->p_level() < elem->p_level())
1129 // extra_order = elem->child_ptr(n)->p_level();
1130 const unsigned int nc =
1131 FEInterface::n_dofs_at_node (base_fe_type, extra_order, elem, n);
1132 for (unsigned int i=0; i!= nc; ++i)
1133 {
1134 Ue(current_dof) =
1135 old_vector(node_dof_indices[current_dof]);
1136 dof_is_fixed[current_dof] = true;
1137 current_dof++;
1138 }
1139 }
1140 }
1141
1142 FEType fe_type = base_fe_type, temp_fe_type;
1143 fe_type.order = fe_type.order + elem->max_descendant_p_level();
1144
1145 // In 3D, project any edge values next
1146 if (dim > 2 && cont != DISCONTINUOUS)
1147 for (auto e : elem->edge_index_range())
1148 {
1149 FEInterface::dofs_on_edge(elem, dim, fe_type,
1150 e, new_side_dofs);
1151
1152 const unsigned int n_new_side_dofs =
1153 cast_int<unsigned int>(new_side_dofs.size());
1154
1155 // Some edge dofs are on nodes and already
1156 // fixed, others are free to calculate
1157 unsigned int free_dofs = 0;
1158 for (unsigned int i=0; i != n_new_side_dofs; ++i)
1159 if (!dof_is_fixed[new_side_dofs[i]])
1160 free_dof[free_dofs++] = i;
1161 Ke.resize (free_dofs, free_dofs); Ke.zero();
1162 Fe.resize (free_dofs); Fe.zero();
1163 // The new edge coefficients
1164 DenseVector<Number> Uedge(free_dofs);
1165
1166 // Add projection terms from each child sharing
1167 // this edge
1168 for (unsigned int c=0; c != n_children; ++c)
1169 {
1170 if (!elem->is_child_on_edge(c,e))
1171 continue;
1172 const Elem * child = elem->child_ptr(c);
1173
1174 std::vector<dof_id_type> child_dof_indices;
1175 if (use_old_dof_indices)
1176 dof_map.old_dof_indices (child,
1177 child_dof_indices, var);
1178 else
1179 dof_map.dof_indices (child,
1180 child_dof_indices, var);
1181 const unsigned int child_n_dofs =
1182 cast_int<unsigned int>
1183 (child_dof_indices.size());
1184
1185 temp_fe_type = base_fe_type;
1186 temp_fe_type.order = temp_fe_type.order + child->p_level();
1187
1188 FEInterface::dofs_on_edge(child, dim,
1189 temp_fe_type, e, old_side_dofs);
1190
1191 // Initialize both child and parent FE data
1192 // on the child's edge
1193 fe->attach_quadrature_rule (qedgerule.get());
1194 fe->edge_reinit (child, e);
1195 const unsigned int n_qp = qedgerule->n_points();
1196
1197 FEMap::inverse_map (dim, elem, xyz_values,
1198 coarse_qpoints);
1199
1200 fe_coarse->reinit(elem, &coarse_qpoints);
1201
1202 // Loop over the quadrature points
1203 for (unsigned int qp=0; qp<n_qp; qp++)
1204 {
1205 // solution value at the quadrature point
1206 OutputNumber fineval = libMesh::zero;
1207 // solution grad at the quadrature point
1208 OutputNumberGradient finegrad;
1209
1210 // Sum the solution values * the DOF
1211 // values at the quadrature point to
1212 // get the solution value and gradient.
1213 for (unsigned int i=0; i<child_n_dofs;
1214 i++)
1215 {
1216 fineval +=
1217 (old_vector(child_dof_indices[i])*
1218 phi_values[i][qp]);
1219 if (cont == C_ONE)
1220 finegrad += (*dphi_values)[i][qp] *
1221 old_vector(child_dof_indices[i]);
1222 }
1223
1224 // Form edge projection matrix
1225 for (unsigned int sidei=0, freei=0; sidei != n_new_side_dofs; ++sidei)
1226 {
1227 unsigned int i = new_side_dofs[sidei];
1228 // fixed DoFs aren't test functions
1229 if (dof_is_fixed[i])
1230 continue;
1231 for (unsigned int sidej=0, freej=0; sidej != n_new_side_dofs; ++sidej)
1232 {
1233 unsigned int j =
1234 new_side_dofs[sidej];
1235 if (dof_is_fixed[j])
1236 Fe(freei) -=
1237 TensorTools::inner_product(phi_coarse[i][qp],
1238 phi_coarse[j][qp]) *
1239 JxW[qp] * Ue(j);
1240 else
1241 Ke(freei,freej) +=
1242 TensorTools::inner_product(phi_coarse[i][qp],
1243 phi_coarse[j][qp]) *
1244 JxW[qp];
1245 if (cont == C_ONE)
1246 {
1247 if (dof_is_fixed[j])
1248 Fe(freei) -=
1249 TensorTools::inner_product((*dphi_coarse)[i][qp],
1250 (*dphi_coarse)[j][qp]) *
1251 JxW[qp] * Ue(j);
1252 else
1253 Ke(freei,freej) +=
1254 TensorTools::inner_product((*dphi_coarse)[i][qp],
1255 (*dphi_coarse)[j][qp]) *
1256 JxW[qp];
1257 }
1258 if (!dof_is_fixed[j])
1259 freej++;
1260 }
1261 Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp],
1262 fineval) * JxW[qp];
1263 if (cont == C_ONE)
1264 Fe(freei) +=
1265 TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1266 freei++;
1267 }
1268 }
1269 }
1270 Ke.cholesky_solve(Fe, Uedge);
1271
1272 // Transfer new edge solutions to element
1273 for (unsigned int i=0; i != free_dofs; ++i)
1274 {
1275 Number & ui = Ue(new_side_dofs[free_dof[i]]);
1276 libmesh_assert(std::abs(ui) < TOLERANCE ||
1277 std::abs(ui - Uedge(i)) < TOLERANCE);
1278 ui = Uedge(i);
1279 dof_is_fixed[new_side_dofs[free_dof[i]]] = true;
1280 }
1281 }
1282
1283 // Project any side values (edges in 2D, faces in 3D)
1284 if (dim > 1 && cont != DISCONTINUOUS)
1285 for (auto s : elem->side_index_range())
1286 {
1287 FEInterface::dofs_on_side(elem, dim, fe_type,
1288 s, new_side_dofs);
1289
1290 const unsigned int n_new_side_dofs =
1291 cast_int<unsigned int>(new_side_dofs.size());
1292
1293 // Some side dofs are on nodes/edges and already
1294 // fixed, others are free to calculate
1295 unsigned int free_dofs = 0;
1296 for (unsigned int i=0; i != n_new_side_dofs; ++i)
1297 if (!dof_is_fixed[new_side_dofs[i]])
1298 free_dof[free_dofs++] = i;
1299 Ke.resize (free_dofs, free_dofs); Ke.zero();
1300 Fe.resize (free_dofs); Fe.zero();
1301 // The new side coefficients
1302 DenseVector<Number> Uside(free_dofs);
1303
1304 // Add projection terms from each child sharing
1305 // this side
1306 for (unsigned int c=0; c != n_children; ++c)
1307 {
1308 if (!elem->is_child_on_side(c,s))
1309 continue;
1310 const Elem * child = elem->child_ptr(c);
1311
1312 std::vector<dof_id_type> child_dof_indices;
1313 if (use_old_dof_indices)
1314 dof_map.old_dof_indices (child,
1315 child_dof_indices, var);
1316 else
1317 dof_map.dof_indices (child,
1318 child_dof_indices, var);
1319 const unsigned int child_n_dofs =
1320 cast_int<unsigned int>
1321 (child_dof_indices.size());
1322
1323 temp_fe_type = base_fe_type;
1324 temp_fe_type.order = temp_fe_type.order + child->p_level();
1325
1326 FEInterface::dofs_on_side(child, dim,
1327 temp_fe_type, s, old_side_dofs);
1328
1329 // Initialize both child and parent FE data
1330 // on the child's side
1331 fe->attach_quadrature_rule (qsiderule.get());
1332 fe->reinit (child, s);
1333 const unsigned int n_qp = qsiderule->n_points();
1334
1335 FEMap::inverse_map (dim, elem, xyz_values,
1336 coarse_qpoints);
1337
1338 fe_coarse->reinit(elem, &coarse_qpoints);
1339
1340 // Loop over the quadrature points
1341 for (unsigned int qp=0; qp<n_qp; qp++)
1342 {
1343 // solution value at the quadrature point
1344 OutputNumber fineval = libMesh::zero;
1345 // solution grad at the quadrature point
1346 OutputNumberGradient finegrad;
1347
1348 // Sum the solution values * the DOF
1349 // values at the quadrature point to
1350 // get the solution value and gradient.
1351 for (unsigned int i=0; i<child_n_dofs;
1352 i++)
1353 {
1354 fineval +=
1355 old_vector(child_dof_indices[i]) *
1356 phi_values[i][qp];
1357 if (cont == C_ONE)
1358 finegrad += (*dphi_values)[i][qp] *
1359 old_vector(child_dof_indices[i]);
1360 }
1361
1362 // Form side projection matrix
1363 for (unsigned int sidei=0, freei=0; sidei != n_new_side_dofs; ++sidei)
1364 {
1365 unsigned int i = new_side_dofs[sidei];
1366 // fixed DoFs aren't test functions
1367 if (dof_is_fixed[i])
1368 continue;
1369 for (unsigned int sidej=0, freej=0; sidej != n_new_side_dofs; ++sidej)
1370 {
1371 unsigned int j =
1372 new_side_dofs[sidej];
1373 if (dof_is_fixed[j])
1374 Fe(freei) -=
1375 TensorTools::inner_product(phi_coarse[i][qp],
1376 phi_coarse[j][qp]) *
1377 JxW[qp] * Ue(j);
1378 else
1379 Ke(freei,freej) +=
1380 TensorTools::inner_product(phi_coarse[i][qp],
1381 phi_coarse[j][qp]) *
1382 JxW[qp];
1383 if (cont == C_ONE)
1384 {
1385 if (dof_is_fixed[j])
1386 Fe(freei) -=
1387 TensorTools::inner_product((*dphi_coarse)[i][qp],
1388 (*dphi_coarse)[j][qp]) *
1389 JxW[qp] * Ue(j);
1390 else
1391 Ke(freei,freej) +=
1392 TensorTools::inner_product((*dphi_coarse)[i][qp],
1393 (*dphi_coarse)[j][qp]) *
1394 JxW[qp];
1395 }
1396 if (!dof_is_fixed[j])
1397 freej++;
1398 }
1399 Fe(freei) += TensorTools::inner_product(fineval, phi_coarse[i][qp]) * JxW[qp];
1400 if (cont == C_ONE)
1401 Fe(freei) +=
1402 TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1403 freei++;
1404 }
1405 }
1406 }
1407 Ke.cholesky_solve(Fe, Uside);
1408
1409 // Transfer new side solutions to element
1410 for (unsigned int i=0; i != free_dofs; ++i)
1411 {
1412 Number & ui = Ue(new_side_dofs[free_dof[i]]);
1413 libmesh_assert(std::abs(ui) < TOLERANCE ||
1414 std::abs(ui - Uside(i)) < TOLERANCE);
1415 ui = Uside(i);
1416 dof_is_fixed[new_side_dofs[free_dof[i]]] = true;
1417 }
1418 }
1419
1420 // Project the interior values, finally
1421
1422 // Some interior dofs are on nodes/edges/sides and
1423 // already fixed, others are free to calculate
1424 unsigned int free_dofs = 0;
1425 for (unsigned int i=0; i != new_n_dofs; ++i)
1426 if (!dof_is_fixed[i])
1427 free_dof[free_dofs++] = i;
1428 Ke.resize (free_dofs, free_dofs); Ke.zero();
1429 Fe.resize (free_dofs); Fe.zero();
1430 // The new interior coefficients
1431 DenseVector<Number> Uint(free_dofs);
1432
1433 // Add projection terms from each child
1434 for (auto & child : elem->child_ref_range())
1435 {
1436 std::vector<dof_id_type> child_dof_indices;
1437 if (use_old_dof_indices)
1438 dof_map.old_dof_indices (&child,
1439 child_dof_indices, var);
1440 else
1441 dof_map.dof_indices (&child,
1442 child_dof_indices, var);
1443 const unsigned int child_n_dofs =
1444 cast_int<unsigned int>
1445 (child_dof_indices.size());
1446
1447 // Initialize both child and parent FE data
1448 // on the child's quadrature points
1449 fe->attach_quadrature_rule (qrule.get());
1450 fe->reinit (&child);
1451 const unsigned int n_qp = qrule->n_points();
1452
1453 FEMap::inverse_map (dim, elem, xyz_values, coarse_qpoints);
1454
1455 fe_coarse->reinit(elem, &coarse_qpoints);
1456
1457 // Loop over the quadrature points
1458 for (unsigned int qp=0; qp<n_qp; qp++)
1459 {
1460 // solution value at the quadrature point
1461 OutputNumber fineval = libMesh::zero;
1462 // solution grad at the quadrature point
1463 OutputNumberGradient finegrad;
1464
1465 // Sum the solution values * the DOF
1466 // values at the quadrature point to
1467 // get the solution value and gradient.
1468 for (unsigned int i=0; i<child_n_dofs; i++)
1469 {
1470 fineval +=
1471 (old_vector(child_dof_indices[i]) *
1472 phi_values[i][qp]);
1473 if (cont == C_ONE)
1474 finegrad += (*dphi_values)[i][qp] *
1475 old_vector(child_dof_indices[i]);
1476 }
1477
1478 // Form interior projection matrix
1479 for (unsigned int i=0, freei=0;
1480 i != new_n_dofs; ++i)
1481 {
1482 // fixed DoFs aren't test functions
1483 if (dof_is_fixed[i])
1484 continue;
1485 for (unsigned int j=0, freej=0; j !=
1486 new_n_dofs; ++j)
1487 {
1488 if (dof_is_fixed[j])
1489 Fe(freei) -=
1490 TensorTools::inner_product(phi_coarse[i][qp],
1491 phi_coarse[j][qp]) *
1492 JxW[qp] * Ue(j);
1493 else
1494 Ke(freei,freej) +=
1495 TensorTools::inner_product(phi_coarse[i][qp],
1496 phi_coarse[j][qp]) *
1497 JxW[qp];
1498 if (cont == C_ONE)
1499 {
1500 if (dof_is_fixed[j])
1501 Fe(freei) -=
1502 TensorTools::inner_product((*dphi_coarse)[i][qp],
1503 (*dphi_coarse)[j][qp]) *
1504 JxW[qp] * Ue(j);
1505 else
1506 Ke(freei,freej) +=
1507 TensorTools::inner_product((*dphi_coarse)[i][qp],
1508 (*dphi_coarse)[j][qp]) *
1509 JxW[qp];
1510 }
1511 if (!dof_is_fixed[j])
1512 freej++;
1513 }
1514 Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp], fineval) *
1515 JxW[qp];
1516 if (cont == C_ONE)
1517 Fe(freei) += TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1518 freei++;
1519 }
1520 }
1521 }
1522 Ke.cholesky_solve(Fe, Uint);
1523
1524 // Transfer new interior solutions to element
1525 for (unsigned int i=0; i != free_dofs; ++i)
1526 {
1527 Number & ui = Ue(free_dof[i]);
1528 libmesh_assert(std::abs(ui) < TOLERANCE ||
1529 std::abs(ui - Uint(i)) < TOLERANCE);
1530 ui = Uint(i);
1531 // We should be fixing all dofs by now; no need to keep track of
1532 // that unless we're debugging
1533#ifndef NDEBUG
1534 dof_is_fixed[free_dof[i]] = true;
1535#endif
1536 }
1537
1538#ifndef NDEBUG
1539 // Make sure every DoF got reached!
1540 for (unsigned int i=0; i != new_n_dofs; ++i)
1541 libmesh_assert(dof_is_fixed[i]);
1542#endif
1543}
1544
1545
1546
1547template <typename OutputType>
1548void
1549FEGenericBase<OutputType>::coarsened_dof_values(const NumericVector<Number> & old_vector,
1550 const DofMap & dof_map,
1551 const Elem * elem,
1552 DenseVector<Number> & Ue,
1553 const bool use_old_dof_indices)
1554{
1555 Ue.resize(0);
1556
1557 for (auto v : make_range(dof_map.n_variables()))
1558 {
1559 DenseVector<Number> Usub;
1560
1561 coarsened_dof_values(old_vector, dof_map, elem, Usub,
1562 v, use_old_dof_indices);
1563
1564 Ue.append (Usub);
1565 }
1566}
1567
1568
1569
1570template <typename OutputType>
1571void
1572FEGenericBase<OutputType>::compute_proj_constraints (DofConstraints & constraints,
1573 DofMap & dof_map,
1574 const unsigned int variable_number,
1575 const Elem * elem)
1576{
1577 libmesh_assert(elem);
1578
1579 const unsigned int Dim = elem->dim();
1580
1581 // Only constrain elements in 2,3D.
1582 if (Dim == 1)
1583 return;
1584
1585 // Only constrain active elements with this method
1586 if (!elem->active())
1587 return;
1588
1589 const Variable & var = dof_map.variable(variable_number);
1590 const FEType & base_fe_type = var.type();
1591 const bool add_p_level = base_fe_type.p_refinement;
1592
1593 // Construct FE objects for this element and its neighbors.
1594 std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1595 (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1596 my_fe->add_p_level_in_reinit(add_p_level);
1597 const FEContinuity cont = my_fe->get_continuity();
1598
1599 // We don't need to constrain discontinuous elements
1600 if (cont == DISCONTINUOUS)
1601 return;
1602 libmesh_assert (cont == C_ZERO || cont == C_ONE ||
1603 cont == SIDE_DISCONTINUOUS);
1604
1605 // this would require some generalisation:
1606 // - e.g. the 'my_fe'-object needs generalisation
1607 // - due to lack of one-to-one correspondence of DOFs and nodes,
1608 // this doesn't work easily.
1609 if (elem->infinite())
1610 libmesh_not_implemented();
1611
1612 std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1613 (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1614 neigh_fe->add_p_level_in_reinit(add_p_level);
1615
1616 QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1617 my_fe->attach_quadrature_rule (&my_qface);
1618 std::vector<Point> neigh_qface;
1619
1620 const std::vector<Real> & JxW = my_fe->get_JxW();
1621 const std::vector<Point> & q_point = my_fe->get_xyz();
1622 const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1623 const std::vector<std::vector<OutputShape>> & neigh_phi =
1624 neigh_fe->get_phi();
1625 const std::vector<Point> * face_normals = nullptr;
1626 const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1627 const std::vector<std::vector<OutputGradient>> * neigh_dphi = nullptr;
1628
1629 std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1630 std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1631
1632 if (cont == C_ONE)
1633 {
1634 const std::vector<Point> & ref_face_normals =
1635 my_fe->get_normals();
1636 face_normals = &ref_face_normals;
1637 const std::vector<std::vector<OutputGradient>> & ref_dphi =
1638 my_fe->get_dphi();
1639 dphi = &ref_dphi;
1640 const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1641 neigh_fe->get_dphi();
1642 neigh_dphi = &ref_neigh_dphi;
1643 }
1644
1645 DenseMatrix<Real> Ke;
1646 DenseVector<Real> Fe;
1647 std::vector<DenseVector<Real>> Ue;
1648
1649 // Look at the element faces. Check to see if we need to
1650 // build constraints.
1651 for (auto s : elem->side_index_range())
1652 {
1653 // Get pointers to the element's neighbor.
1654 const Elem * neigh = elem->neighbor_ptr(s);
1655
1656 if (!neigh)
1657 continue;
1658
1659 if (!var.active_on_subdomain(neigh->subdomain_id()))
1660 continue;
1661
1662 // h refinement constraints:
1663 // constrain dofs shared between
1664 // this element and ones coarser
1665 // than this element.
1666 if (neigh->level() < elem->level())
1667 {
1668 unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
1669 libmesh_assert_less (s_neigh, neigh->n_neighbors());
1670
1671 // Find the minimum p level; we build the h constraint
1672 // matrix with this and then constrain away all higher p
1673 // DoFs.
1674 libmesh_assert(neigh->active());
1675 const unsigned int min_p_level = add_p_level *
1676 std::min(elem->p_level(), neigh->p_level());
1677 // we may need to make the FE objects reinit with the
1678 // minimum shared p_level
1679 const unsigned int old_elem_level = add_p_level * elem->p_level();
1680 if (old_elem_level != min_p_level)
1681 my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + min_p_level - old_elem_level);
1682 const unsigned int old_neigh_level = add_p_level * neigh->p_level();
1683 if (old_neigh_level != min_p_level)
1684 neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + min_p_level - old_neigh_level);
1685
1686 my_fe->reinit(elem, s);
1687
1688 // This function gets called element-by-element, so there
1689 // will be a lot of memory allocation going on. We can
1690 // at least minimize this for the case of the dof indices
1691 // by efficiently preallocating the requisite storage.
1692 // n_nodes is not necessarily n_dofs, but it is better
1693 // than nothing!
1694 my_dof_indices.reserve (elem->n_nodes());
1695 neigh_dof_indices.reserve (neigh->n_nodes());
1696
1697 dof_map.dof_indices (elem, my_dof_indices,
1698 variable_number,
1699 min_p_level);
1700 dof_map.dof_indices (neigh, neigh_dof_indices,
1701 variable_number,
1702 min_p_level);
1703
1704 const unsigned int n_qp = my_qface.n_points();
1705
1706 FEMap::inverse_map (Dim, neigh, q_point, neigh_qface);
1707
1708 neigh_fe->reinit(neigh, &neigh_qface);
1709
1710 // We're only concerned with DOFs whose values (and/or first
1711 // derivatives for C1 elements) are supported on side nodes
1712 FEType elem_fe_type = base_fe_type;
1713 if (old_elem_level != min_p_level)
1714 elem_fe_type.order = base_fe_type.order.get_order() + min_p_level - old_elem_level;
1715 FEType neigh_fe_type = base_fe_type;
1716 if (old_neigh_level != min_p_level)
1717 neigh_fe_type.order = base_fe_type.order.get_order() + min_p_level - old_neigh_level;
1718 FEInterface::dofs_on_side(elem, Dim, elem_fe_type, s, my_side_dofs);
1719 FEInterface::dofs_on_side(neigh, Dim, neigh_fe_type, s_neigh, neigh_side_dofs);
1720
1721 const unsigned int n_side_dofs =
1722 cast_int<unsigned int>(my_side_dofs.size());
1723 libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
1724
1725#ifndef NDEBUG
1726 for (auto i : my_side_dofs)
1727 libmesh_assert_less(i, my_dof_indices.size());
1728 for (auto i : neigh_side_dofs)
1729 libmesh_assert_less(i, neigh_dof_indices.size());
1730#endif
1731
1732 Ke.resize (n_side_dofs, n_side_dofs);
1733 Ue.resize(n_side_dofs);
1734
1735 // Form the projection matrix, (inner product of fine basis
1736 // functions against fine test functions)
1737 for (unsigned int is = 0; is != n_side_dofs; ++is)
1738 {
1739 const unsigned int i = my_side_dofs[is];
1740 for (unsigned int js = 0; js != n_side_dofs; ++js)
1741 {
1742 const unsigned int j = my_side_dofs[js];
1743 for (unsigned int qp = 0; qp != n_qp; ++qp)
1744 {
1745 Ke(is,js) += JxW[qp] * TensorTools::inner_product(phi[i][qp], phi[j][qp]);
1746 if (cont == C_ONE)
1747 Ke(is,js) += JxW[qp] *
1748 TensorTools::inner_product((*dphi)[i][qp] *
1749 (*face_normals)[qp],
1750 (*dphi)[j][qp] *
1751 (*face_normals)[qp]);
1752 }
1753 }
1754 }
1755
1756 // Form the right hand sides, (inner product of coarse basis
1757 // functions against fine test functions)
1758 for (unsigned int is = 0; is != n_side_dofs; ++is)
1759 {
1760 const unsigned int i = neigh_side_dofs[is];
1761 Fe.resize (n_side_dofs);
1762 for (unsigned int js = 0; js != n_side_dofs; ++js)
1763 {
1764 const unsigned int j = my_side_dofs[js];
1765 for (unsigned int qp = 0; qp != n_qp; ++qp)
1766 {
1767 Fe(js) += JxW[qp] *
1768 TensorTools::inner_product(neigh_phi[i][qp],
1769 phi[j][qp]);
1770 if (cont == C_ONE)
1771 Fe(js) += JxW[qp] *
1772 TensorTools::inner_product((*neigh_dphi)[i][qp] *
1773 (*face_normals)[qp],
1774 (*dphi)[j][qp] *
1775 (*face_normals)[qp]);
1776 }
1777 }
1778 Ke.cholesky_solve(Fe, Ue[is]);
1779 }
1780
1781 for (unsigned int js = 0; js != n_side_dofs; ++js)
1782 {
1783 const unsigned int j = my_side_dofs[js];
1784 const dof_id_type my_dof_g = my_dof_indices[j];
1785 libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
1786
1787 // Hunt for "constraining against myself" cases before
1788 // we bother creating a constraint row
1789 bool self_constraint = false;
1790 for (unsigned int is = 0; is != n_side_dofs; ++is)
1791 {
1792 const unsigned int i = neigh_side_dofs[is];
1793 const dof_id_type their_dof_g = neigh_dof_indices[i];
1794 libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1795
1796 if (their_dof_g == my_dof_g)
1797 {
1798#ifndef NDEBUG
1799 const Real their_dof_value = Ue[is](js);
1800 libmesh_assert_less (std::abs(their_dof_value-1.),
1801 10*TOLERANCE);
1802
1803 for (unsigned int k = 0; k != n_side_dofs; ++k)
1804 libmesh_assert(k == is ||
1805 std::abs(Ue[k](js)) <
1806 10*TOLERANCE);
1807#endif
1808
1809 self_constraint = true;
1810 break;
1811 }
1812 }
1813
1814 if (self_constraint)
1815 continue;
1816
1817 DofConstraintRow * constraint_row;
1818
1819 // we may be running constraint methods concurrently
1820 // on multiple threads, so we need a lock to
1821 // ensure that this constraint is "ours"
1822 {
1823 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1824
1825 if (dof_map.is_constrained_dof(my_dof_g))
1826 continue;
1827
1828 constraint_row = &(constraints[my_dof_g]);
1829 libmesh_assert(constraint_row->empty());
1830 }
1831
1832 for (unsigned int is = 0; is != n_side_dofs; ++is)
1833 {
1834 const unsigned int i = neigh_side_dofs[is];
1835 const dof_id_type their_dof_g = neigh_dof_indices[i];
1836 libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1837 libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
1838
1839 const Real their_dof_value = Ue[is](js);
1840
1841 if (std::abs(their_dof_value) < 10*TOLERANCE)
1842 continue;
1843
1844 constraint_row->emplace(their_dof_g, their_dof_value);
1845 }
1846 }
1847
1848 my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + old_elem_level - min_p_level);
1849 neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + old_neigh_level - min_p_level);
1850 }
1851
1852 if (add_p_level)
1853 {
1854 // p refinement constraints:
1855 // constrain dofs shared between
1856 // active elements and neighbors with
1857 // lower polynomial degrees
1858 const unsigned int min_p_level =
1859 neigh->min_p_level_by_neighbor(elem, elem->p_level());
1860 if (min_p_level < elem->p_level())
1861 {
1862 // Adaptive p refinement of non-hierarchic bases will
1863 // require more coding
1864 libmesh_assert(my_fe->is_hierarchic());
1865 dof_map.constrain_p_dofs(variable_number, elem,
1866 s, min_p_level);
1867 }
1868 }
1869 }
1870}
1871
1872#endif // #ifdef LIBMESH_ENABLE_AMR
1873
1874
1875
1876#ifdef LIBMESH_ENABLE_PERIODIC
1877template <typename OutputType>
1878void
1879FEGenericBase<OutputType>::
1880compute_periodic_constraints (DofConstraints & constraints,
1881 DofMap & dof_map,
1882 const PeriodicBoundaries & boundaries,
1883 const MeshBase & mesh,
1884 const PointLocatorBase * point_locator,
1885 const unsigned int variable_number,
1886 const Elem * elem)
1887{
1888 // Only bother if we truly have periodic boundaries
1889 if (boundaries.empty())
1890 return;
1891
1892 libmesh_assert(elem);
1893
1894 // Only constrain active elements with this method
1895 if (!elem->active())
1896 return;
1897
1898 if (elem->infinite())
1899 libmesh_not_implemented();
1900
1901 const unsigned int Dim = elem->dim();
1902
1903 // We need sys_number and variable_number for DofObject methods
1904 // later
1905 const unsigned int sys_number = dof_map.sys_number();
1906
1907 const FEType & base_fe_type = dof_map.variable_type(variable_number);
1908
1909 // Construct FE objects for this element and its pseudo-neighbors.
1910 std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1911 (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1912 const FEContinuity cont = my_fe->get_continuity();
1913
1914 // We don't need to constrain discontinuous elements
1915 if (cont == DISCONTINUOUS)
1916 return;
1917 libmesh_assert (cont == C_ZERO || cont == C_ONE);
1918
1919 // We'll use element size to generate relative tolerances later
1920 const Real primary_hmin = elem->hmin();
1921
1922 std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1923 (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1924
1925 QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1926 my_fe->attach_quadrature_rule (&my_qface);
1927 std::vector<Point> neigh_qface;
1928
1929 const std::vector<Real> & JxW = my_fe->get_JxW();
1930 const std::vector<Point> & q_point = my_fe->get_xyz();
1931 const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1932 const std::vector<std::vector<OutputShape>> & neigh_phi =
1933 neigh_fe->get_phi();
1934 const std::vector<Point> * face_normals = nullptr;
1935 const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1936 const std::vector<std::vector<OutputGradient>> * neigh_dphi = nullptr;
1937 std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1938 std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1939
1940 if (cont != C_ZERO)
1941 {
1942 const std::vector<Point> & ref_face_normals =
1943 my_fe->get_normals();
1944 face_normals = &ref_face_normals;
1945 const std::vector<std::vector<OutputGradient>> & ref_dphi =
1946 my_fe->get_dphi();
1947 dphi = &ref_dphi;
1948 const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1949 neigh_fe->get_dphi();
1950 neigh_dphi = &ref_neigh_dphi;
1951 }
1952
1953 DenseMatrix<Real> Ke;
1954 DenseVector<Real> Fe;
1955 std::vector<DenseVector<Real>> Ue;
1956
1957 // Container to catch the boundary ids that BoundaryInfo hands us.
1958 std::vector<boundary_id_type> bc_ids;
1959
1960 // Look at the element faces. Check to see if we need to
1961 // build constraints.
1962 const unsigned short int max_ns = elem->n_sides();
1963 for (unsigned short int s = 0; s != max_ns; ++s)
1964 {
1965 if (elem->neighbor_ptr(s))
1966 continue;
1967
1968 mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
1969
1970 for (const auto & boundary_id : bc_ids)
1971 {
1972 const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
1973 if (!periodic || !periodic->is_my_variable(variable_number))
1974 continue;
1975
1976 libmesh_assert(point_locator);
1977
1978 // Get pointers to the element's neighbor.
1979 unsigned int s_neigh;
1980 const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s, &s_neigh);
1981
1982 libmesh_error_msg_if(neigh == nullptr,
1983 "PeriodicBoundaries point locator object returned nullptr!");
1984
1985 // periodic (and possibly h refinement) constraints:
1986 // constrain dofs shared between
1987 // this element and ones as coarse
1988 // as or coarser than this element.
1989 if (neigh->level() <= elem->level())
1990 {
1991#ifdef LIBMESH_ENABLE_AMR
1992 // Find the minimum p level; we build the h constraint
1993 // matrix with this and then constrain away all higher p
1994 // DoFs.
1995 libmesh_assert(neigh->active());
1996 const unsigned int min_p_level =
1997 std::min(elem->p_level(), neigh->p_level());
1998
1999 // we may need to make the FE objects reinit with the
2000 // minimum shared p_level
2001 // FIXME - I hate using const_cast<> and avoiding
2002 // accessor functions; there's got to be a
2003 // better way to do this!
2004 const unsigned int old_elem_level = elem->p_level();
2005 if (old_elem_level != min_p_level)
2006 (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
2007 const unsigned int old_neigh_level = neigh->p_level();
2008 if (old_neigh_level != min_p_level)
2009 (const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
2010#endif // #ifdef LIBMESH_ENABLE_AMR
2011
2012 // We can do a projection with a single integration,
2013 // due to the assumption of nested finite element
2014 // subspaces.
2015 // FIXME: it might be more efficient to do nodes,
2016 // then edges, then side, to reduce the size of the
2017 // Cholesky factorization(s)
2018 my_fe->reinit(elem, s);
2019
2020 dof_map.dof_indices (elem, my_dof_indices,
2021 variable_number);
2022 dof_map.dof_indices (neigh, neigh_dof_indices,
2023 variable_number);
2024
2025 // We use neigh_dof_indices_all_variables in the case that the
2026 // periodic boundary condition involves mappings between multiple
2027 // variables.
2028 std::vector<std::vector<dof_id_type>> neigh_dof_indices_all_variables;
2029 if(periodic->has_transformation_matrix())
2030 {
2031 const std::set<unsigned int> & variables = periodic->get_variables();
2032 neigh_dof_indices_all_variables.resize(variables.size());
2033 unsigned int index = 0;
2034 for(unsigned int var : variables)
2035 {
2036 dof_map.dof_indices (neigh, neigh_dof_indices_all_variables[index],
2037 var);
2038 index++;
2039 }
2040 }
2041
2042 const unsigned int n_qp = my_qface.n_points();
2043
2044 // Translate the quadrature points over to the
2045 // neighbor's boundary
2046 std::vector<Point> neigh_point(q_point.size());
2047 for (auto i : index_range(neigh_point))
2048 neigh_point[i] = periodic->get_corresponding_pos(q_point[i]);
2049
2050 FEMap::inverse_map (Dim, neigh, neigh_point,
2051 neigh_qface);
2052
2053 neigh_fe->reinit(neigh, &neigh_qface);
2054
2055 // We're only concerned with DOFs whose values (and/or first
2056 // derivatives for C1 elements) are supported on side nodes
2057 FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
2058 FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
2059
2060 // We're done with functions that examine Elem::p_level(),
2061 // so let's unhack those levels
2062#ifdef LIBMESH_ENABLE_AMR
2063 if (elem->p_level() != old_elem_level)
2064 (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
2065 if (neigh->p_level() != old_neigh_level)
2066 (const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
2067#endif // #ifdef LIBMESH_ENABLE_AMR
2068
2069 const unsigned int n_side_dofs =
2070 cast_int<unsigned int>
2071 (my_side_dofs.size());
2072 libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
2073
2074 Ke.resize (n_side_dofs, n_side_dofs);
2075 Ue.resize(n_side_dofs);
2076
2077 // Form the projection matrix, (inner product of fine basis
2078 // functions against fine test functions)
2079 for (unsigned int is = 0; is != n_side_dofs; ++is)
2080 {
2081 const unsigned int i = my_side_dofs[is];
2082 for (unsigned int js = 0; js != n_side_dofs; ++js)
2083 {
2084 const unsigned int j = my_side_dofs[js];
2085 for (unsigned int qp = 0; qp != n_qp; ++qp)
2086 {
2087 Ke(is,js) += JxW[qp] *
2088 TensorTools::inner_product(phi[i][qp],
2089 phi[j][qp]);
2090 if (cont != C_ZERO)
2091 Ke(is,js) += JxW[qp] *
2092 TensorTools::inner_product((*dphi)[i][qp] *
2093 (*face_normals)[qp],
2094 (*dphi)[j][qp] *
2095 (*face_normals)[qp]);
2096 }
2097 }
2098 }
2099
2100 // Form the right hand sides, (inner product of coarse basis
2101 // functions against fine test functions)
2102 for (unsigned int is = 0; is != n_side_dofs; ++is)
2103 {
2104 const unsigned int i = neigh_side_dofs[is];
2105 Fe.resize (n_side_dofs);
2106 for (unsigned int js = 0; js != n_side_dofs; ++js)
2107 {
2108 const unsigned int j = my_side_dofs[js];
2109 for (unsigned int qp = 0; qp != n_qp; ++qp)
2110 {
2111 Fe(js) += JxW[qp] *
2112 TensorTools::inner_product(neigh_phi[i][qp],
2113 phi[j][qp]);
2114 if (cont != C_ZERO)
2115 Fe(js) += JxW[qp] *
2116 TensorTools::inner_product((*neigh_dphi)[i][qp] *
2117 (*face_normals)[qp],
2118 (*dphi)[j][qp] *
2119 (*face_normals)[qp]);
2120 }
2121 }
2122 Ke.cholesky_solve(Fe, Ue[is]);
2123 }
2124
2125 // Make sure we're not adding recursive constraints
2126 // due to the redundancy in the way we add periodic
2127 // boundary constraints
2128 //
2129 // In order for this to work while threaded or on
2130 // distributed meshes, we need a rigorous way to
2131 // avoid recursive constraints. Here it is:
2132 //
2133 // For vertex DoFs, if there is a "prior" element
2134 // (i.e. a coarser element or an equally refined
2135 // element with a lower id) on this boundary which
2136 // contains the vertex point, then we will avoid
2137 // generating constraints; the prior element (or
2138 // something prior to it) may do so. If we are the
2139 // most prior (or "primary") element on this
2140 // boundary sharing this point, then we look at the
2141 // boundary periodic to us, we find the primary
2142 // element there, and if that primary is coarser or
2143 // equal-but-lower-id, then our vertex dofs are
2144 // constrained in terms of that element.
2145 //
2146 // For edge DoFs, if there is a coarser element
2147 // on this boundary sharing this edge, then we will
2148 // avoid generating constraints (we will be
2149 // constrained indirectly via AMR constraints
2150 // connecting us to the coarser element's DoFs). If
2151 // we are the coarsest element sharing this edge,
2152 // then we generate constraints if and only if we
2153 // are finer than the coarsest element on the
2154 // boundary periodic to us sharing the corresponding
2155 // periodic edge, or if we are at equal level but
2156 // our edge nodes have higher ids than the periodic
2157 // edge nodes (sorted from highest to lowest, then
2158 // compared lexicographically)
2159 //
2160 // For face DoFs, we generate constraints if we are
2161 // finer than our periodic neighbor, or if we are at
2162 // equal level but our element id is higher than its
2163 // element id.
2164 //
2165 // If the primary neighbor is also the current elem
2166 // (a 1-element-thick mesh) then we choose which
2167 // vertex dofs to constrain via lexicographic
2168 // ordering on point locations
2169
2170 // FIXME: This code doesn't yet properly handle
2171 // cases where multiple different periodic BCs
2172 // intersect.
2173 std::set<dof_id_type> my_constrained_dofs;
2174
2175 // Container to catch boundary IDs handed back by BoundaryInfo.
2176 std::vector<boundary_id_type> new_bc_ids;
2177
2178 for (auto n : elem->node_index_range())
2179 {
2180 if (!elem->is_node_on_side(n,s))
2181 continue;
2182
2183 const Node & my_node = elem->node_ref(n);
2184
2185 if (elem->is_vertex(n))
2186 {
2187 // Find all boundary ids that include this
2188 // point and have periodic boundary
2189 // conditions for this variable
2190 std::set<boundary_id_type> point_bcids;
2191
2192 for (unsigned int new_s = 0;
2193 new_s != max_ns; ++new_s)
2194 {
2195 if (!elem->is_node_on_side(n,new_s))
2196 continue;
2197
2198 mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
2199
2200 for (const auto & new_boundary_id : new_bc_ids)
2201 {
2202 const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2203 if (new_periodic && new_periodic->is_my_variable(variable_number))
2204 point_bcids.insert(new_boundary_id);
2205 }
2206 }
2207
2208 // See if this vertex has point neighbors to
2209 // defer to
2210 if (primary_boundary_point_neighbor
2211 (elem, my_node, mesh.get_boundary_info(), point_bcids)
2212 != elem)
2213 continue;
2214
2215 // Find the complementary boundary id set
2216 std::set<boundary_id_type> point_pairedids;
2217 for (const auto & new_boundary_id : point_bcids)
2218 {
2219 const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2220 point_pairedids.insert(new_periodic->pairedboundary);
2221 }
2222
2223 // What do we want to constrain against?
2224 const Elem * primary_elem = nullptr;
2225 const Elem * main_neigh = nullptr;
2226 Point main_pt = my_node,
2227 primary_pt = my_node;
2228
2229 for (const auto & new_boundary_id : point_bcids)
2230 {
2231 // Find the corresponding periodic point and
2232 // its primary neighbor
2233 const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2234
2235 const Point neigh_pt =
2236 new_periodic->get_corresponding_pos(my_node);
2237
2238 // If the point is getting constrained
2239 // to itself by this PBC then we don't
2240 // generate any constraints
2241 if (neigh_pt.absolute_fuzzy_equals
2242 (my_node, primary_hmin*TOLERANCE))
2243 continue;
2244
2245 // Otherwise we'll have a constraint in
2246 // one direction or another
2247 if (!primary_elem)
2248 primary_elem = elem;
2249
2250 const Elem * primary_neigh =
2251 primary_boundary_point_neighbor(neigh, neigh_pt,
2252 mesh.get_boundary_info(),
2253 point_pairedids);
2254
2255 libmesh_assert(primary_neigh);
2256
2257 if (new_boundary_id == boundary_id)
2258 {
2259 main_neigh = primary_neigh;
2260 main_pt = neigh_pt;
2261 }
2262
2263 // Finer elements will get constrained in
2264 // terms of coarser neighbors, not the
2265 // other way around
2266 if ((primary_neigh->level() > primary_elem->level()) ||
2267
2268 // For equal-level elements, the one with
2269 // higher id gets constrained in terms of
2270 // the one with lower id
2271 (primary_neigh->level() == primary_elem->level() &&
2272 primary_neigh->id() > primary_elem->id()) ||
2273
2274 // On a one-element-thick mesh, we compare
2275 // points to see what side gets constrained
2276 (primary_neigh == primary_elem &&
2277 (neigh_pt > primary_pt)))
2278 continue;
2279
2280 primary_elem = primary_neigh;
2281 primary_pt = neigh_pt;
2282 }
2283
2284 if (!primary_elem ||
2285 primary_elem != main_neigh ||
2286 primary_pt != main_pt)
2287 continue;
2288 }
2289 else if (elem->is_edge(n))
2290 {
2291 // Find which edge we're on
2292 unsigned int e=0, ne = elem->n_edges();
2293 for (; e != ne; ++e)
2294 {
2295 if (elem->is_node_on_edge(n,e))
2296 break;
2297 }
2298 libmesh_assert_less (e, elem->n_edges());
2299
2300 // Find the edge end nodes
2301 const Node
2302 * e1 = nullptr,
2303 * e2 = nullptr;
2304 for (auto nn : elem->node_index_range())
2305 {
2306 if (nn == n)
2307 continue;
2308
2309 if (elem->is_node_on_edge(nn, e))
2310 {
2311 if (e1 == nullptr)
2312 {
2313 e1 = elem->node_ptr(nn);
2314 }
2315 else
2316 {
2317 e2 = elem->node_ptr(nn);
2318 break;
2319 }
2320 }
2321 }
2322 libmesh_assert (e1 && e2);
2323
2324 // Find all boundary ids that include this
2325 // edge and have periodic boundary
2326 // conditions for this variable
2327 std::set<boundary_id_type> edge_bcids;
2328
2329 for (unsigned int new_s = 0;
2330 new_s != max_ns; ++new_s)
2331 {
2332 if (!elem->is_node_on_side(n,new_s))
2333 continue;
2334
2335 // We're reusing the new_bc_ids vector created outside the loop over nodes.
2336 mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
2337
2338 for (const auto & new_boundary_id : new_bc_ids)
2339 {
2340 const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2341 if (new_periodic && new_periodic->is_my_variable(variable_number))
2342 edge_bcids.insert(new_boundary_id);
2343 }
2344 }
2345
2346
2347 // See if this edge has neighbors to defer to
2348 if (primary_boundary_edge_neighbor
2349 (elem, *e1, *e2, mesh.get_boundary_info(), edge_bcids)
2350 != elem)
2351 continue;
2352
2353 // Find the complementary boundary id set
2354 std::set<boundary_id_type> edge_pairedids;
2355 for (const auto & new_boundary_id : edge_bcids)
2356 {
2357 const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2358 edge_pairedids.insert(new_periodic->pairedboundary);
2359 }
2360
2361 // What do we want to constrain against?
2362 const Elem * primary_elem = nullptr;
2363 const Elem * main_neigh = nullptr;
2364 Point main_pt1 = *e1,
2365 main_pt2 = *e2,
2366 primary_pt1 = *e1,
2367 primary_pt2 = *e2;
2368
2369 for (const auto & new_boundary_id : edge_bcids)
2370 {
2371 // Find the corresponding periodic edge and
2372 // its primary neighbor
2373 const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2374
2375 Point neigh_pt1 = new_periodic->get_corresponding_pos(*e1),
2376 neigh_pt2 = new_periodic->get_corresponding_pos(*e2);
2377
2378 // If the edge is getting constrained
2379 // to itself by this PBC then we don't
2380 // generate any constraints
2381 if (neigh_pt1.absolute_fuzzy_equals
2382 (*e1, primary_hmin*TOLERANCE) &&
2383 neigh_pt2.absolute_fuzzy_equals
2384 (*e2, primary_hmin*TOLERANCE))
2385 continue;
2386
2387 // Otherwise we'll have a constraint in
2388 // one direction or another
2389 if (!primary_elem)
2390 primary_elem = elem;
2391
2392 const Elem * primary_neigh = primary_boundary_edge_neighbor
2393 (neigh, neigh_pt1, neigh_pt2,
2394 mesh.get_boundary_info(), edge_pairedids);
2395
2396 libmesh_assert(primary_neigh);
2397
2398 if (new_boundary_id == boundary_id)
2399 {
2400 main_neigh = primary_neigh;
2401 main_pt1 = neigh_pt1;
2402 main_pt2 = neigh_pt2;
2403 }
2404
2405 // If we have a one-element thick mesh,
2406 // we'll need to sort our points to get a
2407 // consistent ordering rule
2408 //
2409 // Use >= in this test to make sure that,
2410 // for angular constraints, no node gets
2411 // constrained to itself.
2412 if (primary_neigh == primary_elem)
2413 {
2414 if (primary_pt1 > primary_pt2)
2415 std::swap(primary_pt1, primary_pt2);
2416 if (neigh_pt1 > neigh_pt2)
2417 std::swap(neigh_pt1, neigh_pt2);
2418
2419 if (neigh_pt2 >= primary_pt2)
2420 continue;
2421 }
2422
2423 // Otherwise:
2424 // Finer elements will get constrained in
2425 // terms of coarser ones, not the other way
2426 // around
2427 if ((primary_neigh->level() > primary_elem->level()) ||
2428
2429 // For equal-level elements, the one with
2430 // higher id gets constrained in terms of
2431 // the one with lower id
2432 (primary_neigh->level() == primary_elem->level() &&
2433 primary_neigh->id() > primary_elem->id()))
2434 continue;
2435
2436 primary_elem = primary_neigh;
2437 primary_pt1 = neigh_pt1;
2438 primary_pt2 = neigh_pt2;
2439 }
2440
2441 if (!primary_elem ||
2442 primary_elem != main_neigh ||
2443 primary_pt1 != main_pt1 ||
2444 primary_pt2 != main_pt2)
2445 continue;
2446 }
2447 else if (elem->is_face(n))
2448 {
2449 // If we have a one-element thick mesh,
2450 // use the ordering of the face node and its
2451 // periodic counterpart to determine what
2452 // gets constrained
2453 if (neigh == elem)
2454 {
2455 const Point neigh_pt =
2456 periodic->get_corresponding_pos(my_node);
2457 if (neigh_pt > my_node)
2458 continue;
2459 }
2460
2461 // Otherwise:
2462 // Finer elements will get constrained in
2463 // terms of coarser ones, not the other way
2464 // around
2465 if ((neigh->level() > elem->level()) ||
2466
2467 // For equal-level elements, the one with
2468 // higher id gets constrained in terms of
2469 // the one with lower id
2470 (neigh->level() == elem->level() &&
2471 neigh->id() > elem->id()))
2472 continue;
2473 }
2474
2475 // If we made it here without hitting a continue
2476 // statement, then we're at a node whose dofs
2477 // should be constrained by this element's
2478 // calculations.
2479 const unsigned int n_comp =
2480 my_node.n_comp(sys_number, variable_number);
2481
2482 for (unsigned int i=0; i != n_comp; ++i)
2483 my_constrained_dofs.insert
2484 (my_node.dof_number
2485 (sys_number, variable_number, i));
2486 }
2487
2488 // FIXME: old code for disambiguating periodic BCs:
2489 // this is not threadsafe nor safe to run on a
2490 // non-serialized mesh.
2491 /*
2492 std::vector<bool> recursive_constraint(n_side_dofs, false);
2493
2494 for (unsigned int is = 0; is != n_side_dofs; ++is)
2495 {
2496 const unsigned int i = neigh_side_dofs[is];
2497 const dof_id_type their_dof_g = neigh_dof_indices[i];
2498 libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2499
2500 {
2501 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2502
2503 if (!dof_map.is_constrained_dof(their_dof_g))
2504 continue;
2505 }
2506
2507 DofConstraintRow & their_constraint_row =
2508 constraints[their_dof_g].first;
2509
2510 for (unsigned int js = 0; js != n_side_dofs; ++js)
2511 {
2512 const unsigned int j = my_side_dofs[js];
2513 const dof_id_type my_dof_g = my_dof_indices[j];
2514 libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2515
2516 if (their_constraint_row.count(my_dof_g))
2517 recursive_constraint[js] = true;
2518 }
2519 }
2520 */
2521
2522 for (unsigned int js = 0; js != n_side_dofs; ++js)
2523 {
2524 // FIXME: old code path
2525 // if (recursive_constraint[js])
2526 // continue;
2527
2528 const unsigned int j = my_side_dofs[js];
2529 const dof_id_type my_dof_g = my_dof_indices[j];
2530 libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2531
2532 // FIXME: new code path
2533 if (!my_constrained_dofs.count(my_dof_g))
2534 continue;
2535
2536 DofConstraintRow * constraint_row;
2537
2538 // we may be running constraint methods concurrently
2539 // on multiple threads, so we need a lock to
2540 // ensure that this constraint is "ours"
2541 {
2542 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2543
2544 if (dof_map.is_constrained_dof(my_dof_g))
2545 continue;
2546
2547 constraint_row = &(constraints[my_dof_g]);
2548 libmesh_assert(constraint_row->empty());
2549 }
2550
2551 for (unsigned int is = 0; is != n_side_dofs; ++is)
2552 {
2553 const unsigned int i = neigh_side_dofs[is];
2554 const dof_id_type their_dof_g = neigh_dof_indices[i];
2555 libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2556
2557 // Periodic constraints should never be
2558 // self-constraints
2559 // libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
2560
2561 const Real their_dof_value = Ue[is](js);
2562
2563 if (their_dof_g == my_dof_g)
2564 {
2565 libmesh_assert_less (std::abs(their_dof_value-1.), 1.e-5);
2566 for (unsigned int k = 0; k != n_side_dofs; ++k)
2567 libmesh_assert(k == is || std::abs(Ue[k](js)) < 1.e-5);
2568 continue;
2569 }
2570
2571 if (std::abs(their_dof_value) < 10*TOLERANCE)
2572 continue;
2573
2574 if(!periodic->has_transformation_matrix())
2575 {
2576 constraint_row->emplace(their_dof_g, their_dof_value);
2577 }
2578 else
2579 {
2580 // In this case the current variable is constrained in terms of other variables.
2581 // We assume that all variables in this constraint have the same FE type (this
2582 // is asserted below), and hence we can create the constraint row contribution
2583 // by multiplying their_dof_value by the corresponding row of the transformation
2584 // matrix.
2585
2586 const std::set<unsigned int> & variables = periodic->get_variables();
2587 neigh_dof_indices_all_variables.resize(variables.size());
2588 unsigned int index = 0;
2589 for(unsigned int other_var : variables)
2590 {
2591 libmesh_assert_msg(base_fe_type == dof_map.variable_type(other_var), "FE types must match for all variables involved in constraint");
2592
2593 Real var_weighting = periodic->get_transformation_matrix()(variable_number, other_var);
2594 constraint_row->emplace(neigh_dof_indices_all_variables[index][i],
2595 var_weighting*their_dof_value);
2596 index++;
2597 }
2598 }
2599
2600 }
2601 }
2602 }
2603 // p refinement constraints:
2604 // constrain dofs shared between
2605 // active elements and neighbors with
2606 // lower polynomial degrees
2607#ifdef LIBMESH_ENABLE_AMR
2608 const unsigned int min_p_level =
2609 neigh->min_p_level_by_neighbor(elem, elem->p_level());
2610 if (min_p_level < elem->p_level())
2611 {
2612 // Adaptive p refinement of non-hierarchic bases will
2613 // require more coding
2614 libmesh_assert(my_fe->is_hierarchic());
2615 dof_map.constrain_p_dofs(variable_number, elem,
2616 s, min_p_level);
2617 }
2618#endif // #ifdef LIBMESH_ENABLE_AMR
2619 }
2620 }
2621}
2622
2623#endif // LIBMESH_ENABLE_PERIODIC
2624
2625// ------------------------------------------------------------
2626// Explicit instantiations
2627template class LIBMESH_EXPORT FEGenericBase<Real>;
2628template class LIBMESH_EXPORT FEGenericBase<RealGradient>;
2629
2630} // namespace libMesh
dim
unsigned int dim
Definition adaptivity_ex3.C:114
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces,...
Definition boundary_info.h:58
libMesh::BoundaryInfo::boundary_ids
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
Definition boundary_info.C:1389
libMesh::DenseMatrix
Defines a dense matrix for use in Finite Element-type computations.
Definition dense_matrix.h:73
libMesh::DenseMatrix::cholesky_solve
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
Definition dense_matrix_impl.h:844
libMesh::DenseMatrix::resize
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition dense_matrix.h:953
libMesh::DenseMatrix::zero
virtual void zero() override final
Sets all elements of the matrix to 0 and resets any decomposition flag which may have been previously...
Definition dense_matrix.h:969
libMesh::DenseVector
Defines a dense vector for use in Finite Element-type computations.
Definition dense_vector.h:61
libMesh::DenseVector::resize
void resize(const unsigned int n)
Resize the vector.
Definition dense_vector.h:440
libMesh::DenseVector::zero
virtual void zero() override final
Set every element in the vector to 0.
Definition dense_vector.h:464
libMesh::DenseVector::append
void append(const DenseVector< T2 > &other_vector)
Append additional entries to (resizing, but unchanging) the vector.
Definition dense_vector.h:452
libMesh::DofConstraints
The constraint matrix storage format.
Definition dof_map.h:112
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition dof_map.h:181
libMesh::DofMap::constrain_p_dofs
void constrain_p_dofs(unsigned int var, const Elem *elem, unsigned int s, unsigned int p)
Constrains degrees of freedom on side s of element elem which correspond to variable number var and t...
Definition dof_map_constraints.C:5379
libMesh::DofMap::n_variables
unsigned int n_variables() const override
Definition dof_map.h:743
libMesh::DofMap::dof_indices
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition dof_map.C:2203
libMesh::DofMap::variable_type
const FEType & variable_type(const unsigned int i) const
Definition dof_map.h:2401
libMesh::DofMap::sys_number
unsigned int sys_number() const
Definition dof_map.h:2353
libMesh::DofMap::is_constrained_dof
bool is_constrained_dof(const dof_id_type dof) const
Definition dof_map.h:2439
libMesh::DofMap::variable
const Variable & variable(const unsigned int c) const override
Definition dof_map.h:2371
libMesh::DofMap::old_dof_indices
void old_dof_indices(const Elem &elem, unsigned int n, std::vector< dof_id_type > &di, const unsigned int vn) const
Appends to the vector di the old global degree of freedom indices for elem.node_ref(n),...
Definition dof_map.C:2480
libMesh::DofObject::n_comp
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition dof_object.h:992
libMesh::DofObject::dof_number
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition dof_object.h:1022
libMesh::DofObject::invalid_id
static constexpr dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition dof_object.h:473
libMesh::DofObject::id
dof_id_type id() const
Definition dof_object.h:833
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::is_node_on_side
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
libMesh::Elem::is_child_on_side
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const =0
libMesh::Elem::active
bool active() const
Definition elem.h:3079
libMesh::Elem::node_ref
const Node & node_ref(const unsigned int i) const
Definition elem.h:2659
libMesh::Elem::is_child_on_edge
virtual bool is_child_on_edge(const unsigned int c, const unsigned int e) const
Definition elem.C:2435
libMesh::Elem::n_nodes
virtual unsigned int n_nodes() const =0
libMesh::Elem::child_ptr
const Elem * child_ptr(unsigned int i) const
Definition elem.h:3301
libMesh::Elem::which_neighbor_am_i
unsigned int which_neighbor_am_i(const Elem *e) const
This function tells you which neighbor e is.
Definition elem.h:3057
libMesh::Elem::n_neighbors
unsigned int n_neighbors() const
Definition elem.h:812
libMesh::Elem::is_face
virtual bool is_face(const unsigned int i) const =0
libMesh::Elem::find_point_neighbors
void find_point_neighbors(const Point &p, std::set< const Elem * > &neighbor_set) const
This function finds all active elements (including this one) which are in the same manifold as this e...
Definition elem.C:907
libMesh::Elem::hmin
virtual Real hmin() const
Definition elem.C:623
libMesh::Elem::dim
virtual unsigned short dim() const =0
libMesh::Elem::find_edge_neighbors
void find_edge_neighbors(const Point &p1, const Point &p2, std::set< const Elem * > &neighbor_set) const
This function finds all active elements in the same manifold as this element which touch the current ...
Definition elem.C:997
libMesh::Elem::subdomain_id
subdomain_id_type subdomain_id() const
Definition elem.h:2712
libMesh::Elem::child_ref_range
SimpleRange< ChildRefIter > child_ref_range()
Returns a range with all children of a parent element, usable in range-based for loops.
Definition elem.h:2476
libMesh::Elem::is_edge
virtual bool is_edge(const unsigned int i) const =0
libMesh::Elem::level
unsigned int level() const
Definition elem.h:3212
libMesh::Elem::min_p_level_by_neighbor
unsigned int min_p_level_by_neighbor(const Elem *neighbor, unsigned int current_min) const
Definition elem.C:2452
libMesh::Elem::node_ptr
const Node * node_ptr(const unsigned int i) const
Definition elem.h:2637
libMesh::Elem::edge_index_range
IntRange< unsigned short > edge_index_range() const
Definition elem.h:2830
libMesh::Elem::is_vertex
virtual bool is_vertex(const unsigned int i) const =0
libMesh::Elem::n_children
virtual unsigned int n_children() const =0
libMesh::Elem::max_descendant_p_level
unsigned int max_descendant_p_level() const
Definition elem.h:3384
libMesh::Elem::node_index_range
IntRange< unsigned short > node_index_range() const
Definition elem.h:2821
libMesh::Elem::is_node_on_edge
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const =0
libMesh::Elem::n_edges
virtual unsigned int n_edges() const =0
libMesh::Elem::n_sides
virtual unsigned int n_sides() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3246
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition elem.h:2736
libMesh::Elem::infinite
virtual bool infinite() const =0
libMesh::Elem::side_index_range
IntRange< unsigned short > side_index_range() const
Definition elem.h:2848
libMesh::FEGenericBase
This class forms the foundation from which generic finite elements may be derived.
Definition fe_base.h:86
libMesh::FEGenericBase::determine_calculations
virtual_for_inffe void determine_calculations()
Determine which values are to be calculated, for both the FE itself and for the FEMap.
Definition fe_base.C:951
libMesh::FEGenericBase::build
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libMesh::FEGenericBase::compute_dual_shape_coeffs
void compute_dual_shape_coeffs(const std::vector< Real > &JxW, const std::vector< std::vector< OutputShape > > &phi)
Compute the dual basis coefficients dual_coeff we rely on the JxW (or weights) and the phi values,...
Definition fe_base.h:838
libMesh::FEGenericBase::print_dphi
virtual void print_dphi(std::ostream &os) const override
Prints the value of each shape function's derivative at each quadrature point.
Definition fe_base.C:933
libMesh::FEGenericBase::coarsened_dof_values
static void coarsened_dof_values(const NumericVector< Number > &global_vector, const DofMap &dof_map, const Elem *coarse_elem, DenseVector< Number > &coarse_dofs, const unsigned int var, const bool use_old_dof_indices=false)
Creates a local projection on coarse_elem, based on the DoF values in global_vector for it's children...
Definition fe_base.C:1014
libMesh::FEGenericBase::compute_proj_constraints
static void compute_proj_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...
Definition fe_base.C:1572
libMesh::FEGenericBase::build_InfFE
static std::unique_ptr< FEGenericBase > build_InfFE(const unsigned int dim, const FEType &type)
Builds a specific infinite element type.
libMesh::FEGenericBase::print_dual_d2phi
virtual void print_dual_d2phi(std::ostream &os) const override
Definition fe_base.C:999
libMesh::FEGenericBase::print_d2phi
virtual void print_d2phi(std::ostream &os) const override
Prints the value of each shape function's second derivatives at each quadrature point.
Definition fe_base.C:991
libMesh::FEGenericBase::print_phi
virtual void print_phi(std::ostream &os) const override
Prints the value of each shape function at each quadrature point.
Definition fe_base.C:914
libMesh::FEGenericBase::compute_periodic_constraints
static void compute_periodic_constraints(DofConstraints &constraints, DofMap &dof_map, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const unsigned int variable_number, const Elem *elem)
Computes the constraint matrix contributions (for meshes with periodic boundary conditions) correspon...
Definition fe_base.C:1880
libMesh::FEGenericBase::print_dual_phi
virtual void print_dual_phi(std::ostream &os) const override
Definition fe_base.C:922
libMesh::FEGenericBase::OutputNumber
TensorTools::MakeNumber< OutputShape >::type OutputNumber
Definition fe_base.h:123
libMesh::FEGenericBase::compute_shape_functions
virtual void compute_shape_functions(const Elem *elem, const std::vector< Point > &qp) override
After having updated the jacobian and the transformation from local to global coordinates in FEMap::c...
Definition fe_base.C:774
libMesh::FEGenericBase::print_dual_dphi
virtual void print_dual_dphi(std::ostream &os) const override
Definition fe_base.C:941
libMesh::FEGenericBase::compute_dual_shape_functions
void compute_dual_shape_functions()
Compute dual_phi, dual_dphi, dual_d2phi It is only valid for this to be called after reinit has occur...
Definition fe_base.h:830
libMesh::FEInterface::dofs_on_side
static void dofs_on_side(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int s, std::vector< unsigned int > &di, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
Definition fe_interface.C:601
libMesh::FEInterface::n_dofs
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition fe_interface.C:356
libMesh::FEInterface::n_dofs_at_node
static unsigned int n_dofs_at_node(const unsigned int dim, const FEType &fe_t, const ElemType t, const unsigned int n)
Definition fe_interface.C:440
libMesh::FEInterface::dofs_on_edge
static void dofs_on_edge(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int e, std::vector< unsigned int > &di, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
Definition fe_interface.C:615
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, const bool extra_checks=true)
Definition fe_map.C:1512
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::inf_map
InfMapType inf_map
The coordinate mapping type of the infinite element.
Definition fe_type.h:284
libMesh::FEType::default_quadrature_rule
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition fe_type.C:35
libMesh::FEType::default_quadrature_order
Order default_quadrature_order() const
Definition fe_type.h:415
libMesh::FEType::order
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
libMesh::FEType::p_refinement
bool p_refinement
Whether or not the finite elements for this type increase their p refinement level on geometric eleme...
Definition fe_type.h:292
libMesh::FEType::radial_family
FEFamily radial_family
The type of approximation in radial direction.
Definition fe_type.h:276
libMesh::FEType::family
FEFamily family
The type of finite element.
Definition fe_type.h:228
libMesh::MeshBase
This is the MeshBase class.
Definition mesh_base.h:81
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition mesh_base.h:170
libMesh::Node
A Node is like a Point, but with more information.
Definition node.h:55
libMesh::NumericVector
Provides a uniform interface to vector storage schemes for different linear algebra libraries.
Definition numeric_vector.h:70
libMesh::OrderWrapper::get_order
int get_order() const
Explicitly request the order as an int.
Definition fe_type.h:80
libMesh::PeriodicBoundaries
We're using a class instead of a typedef to allow forward declarations and future flexibility.
Definition periodic_boundaries.h:50
libMesh::PeriodicBoundaries::boundary
PeriodicBoundaryBase * boundary(boundary_id_type id)
Definition periodic_boundaries.C:38
libMesh::PeriodicBoundaries::neighbor
const Elem * neighbor(boundary_id_type boundary_id, const PointLocatorBase &point_locator, const Elem *e, unsigned int side, unsigned int *neigh_side=nullptr) const
Definition periodic_boundaries.C:59
libMesh::PeriodicBoundaryBase
The base class for defining periodic boundaries.
Definition periodic_boundary_base.h:49
libMesh::PeriodicBoundaryBase::get_transformation_matrix
const DenseMatrix< Real > & get_transformation_matrix() const
Get the transformation matrix, if it is defined.
Definition periodic_boundary_base.C:86
libMesh::PeriodicBoundaryBase::is_my_variable
bool is_my_variable(unsigned int var_num) const
Definition periodic_boundary_base.C:71
libMesh::PeriodicBoundaryBase::get_variables
const std::set< unsigned int > & get_variables() const
Get the set of variables for this periodic boundary condition.
Definition periodic_boundary_base.C:109
libMesh::PeriodicBoundaryBase::get_corresponding_pos
virtual Point get_corresponding_pos(const Point &pt) const =0
This function should be overridden by derived classes to define how one finds corresponding nodes on ...
libMesh::PointLocatorBase
This is the base class for point locators.
Definition point_locator_base.h:55
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:41
libMesh::QBase::n_points
unsigned int n_points() const
Definition quadrature.h:131
libMesh::QGauss
This class implements specific orders of Gauss quadrature.
Definition quadrature_gauss.h:40
libMesh::TypeVector::absolute_fuzzy_equals
bool absolute_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:1043
libMesh::Variable
This class defines the notion of a variable in the system.
Definition variable.h:51
libMesh::Variable::active_on_subdomain
bool active_on_subdomain(subdomain_id_type sid) const
Definition variable.h:167
libMesh::Variable::type
const FEType & type() const
Definition variable.h:144
libMesh::VectorValue
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
Definition vector_value.h:49
mesh
MeshBase & mesh
Definition mesh_communication.C:1473
libMesh::TensorTools::inner_product
LIBMESH_DEVICE_INLINE std::enable_if< ScalarTraits< T >::value &&ScalarTraits< T2 >::value, typenameCompareTypes< T, T2 >::supertype >::type inner_product(const T &a, const T2 &b)
Definition tensor_tools.h:51
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition threads.C:30
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::SIDE_DISCONTINUOUS
@ SIDE_DISCONTINUOUS
Definition enum_fe_family.h:93
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition int_range.h:153
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::zero
const Number zero
.
Definition libmesh.h:297
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:100
libMesh::is
PetscErrorCode PetscInt const PetscInt IS * is
Definition petsc_dm_wrapper.C:62
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition libmesh_common.h:155
libMesh::dof_id_type
uint8_t dof_id_type
Definition id_types.h:67
libMesh::L2_RAVIART_THOMAS
@ L2_RAVIART_THOMAS
Definition enum_fe_family.h:64
libMesh::L2_HIERARCHIC_VEC
@ L2_HIERARCHIC_VEC
Definition enum_fe_family.h:67
libMesh::RATIONAL_BERNSTEIN
@ RATIONAL_BERNSTEIN
Definition enum_fe_family.h:69
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:148
libMesh::make_range
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition int_range.h:176
n_nodes
const dof_id_type n_nodes
Definition tecplot_io.C:67
Generated on Thu Nov 12 2026 20:15:21 for libMesh by doxygen 1.9.8