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