libMesh
fe_base.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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
49 namespace
50 {
51 using namespace libMesh;
52 
53 #ifdef LIBMESH_ENABLE_PERIODIC
54 
55 // Find the "primary" element around a boundary point:
56 const 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:
116 const 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 
182 namespace libMesh
183 {
184 
185 
186 
187 // ------------------------------------------------------------
188 // FEBase class members
189 template <>
190 std::unique_ptr<FEGenericBase<Real>>
191 FEGenericBase<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 
411 template <>
412 std::unique_ptr<FEGenericBase<RealGradient>>
413 FEGenericBase<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 
544 template <>
545 std::unique_ptr<FEGenericBase<Real>>
546 FEGenericBase<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 
748 template <>
749 std::unique_ptr<FEGenericBase<RealGradient>>
750 FEGenericBase<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 
761 template <typename OutputType>
762 void 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.
803 template <>
804 void 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 
837 template <>
838 void 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 
875 template <typename OutputType>
876 void 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 
883 template <typename OutputType>
884 void 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 
894 template <typename OutputType>
895 void 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 
902 template <typename OutputType>
903 void 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 
912 template <typename OutputType>
913 void FEGenericBase<OutputType>::determine_calculations()
914 {
915  this->calculations_started = true;
916 
917  // If the user forgot to request anything, but we're enabling
918  // deprecated backwards compatibility, then we'll be safe and
919  // calculate everything. If we haven't enable deprecated backwards
920  // compatibility then we'll scream and die.
921 #ifdef LIBMESH_ENABLE_DEPRECATED
922 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
923  if (!this->calculate_nothing &&
924  !this->calculate_phi && !this->calculate_dphi &&
925  !this->calculate_dphiref &&
926  !this->calculate_d2phi && !this->calculate_curl_phi &&
927  !this->calculate_div_phi && !this->calculate_map)
928  {
929  libmesh_deprecated();
930  this->calculate_phi = this->calculate_dphi = this->calculate_d2phi = this->calculate_dphiref = true;
931  if (FEInterface::field_type(fe_type.family) == TYPE_VECTOR)
932  {
933  this->calculate_curl_phi = true;
934  this->calculate_div_phi = true;
935  }
936  }
937 #else
938  if (!this->calculate_nothing &&
939  !this->calculate_phi && !this->calculate_dphi &&
940  !this->calculate_dphiref &&
941  !this->calculate_curl_phi && !this->calculate_div_phi &&
942  !this->calculate_map)
943  {
944  libmesh_deprecated();
945  this->calculate_phi = this->calculate_dphi = this->calculate_dphiref = true;
946  if (FEInterface::field_type(fe_type.family) == TYPE_VECTOR)
947  {
948  this->calculate_curl_phi = true;
949  this->calculate_div_phi = true;
950  }
951  }
952 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
953 #else
954 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
955  libmesh_assert (this->calculate_nothing ||
956  this->calculate_phi || this->calculate_dphi ||
957  this->calculate_d2phi ||
958  this->calculate_dphiref ||
959  this->calculate_curl_phi || this->calculate_div_phi ||
960  this->calculate_map);
961 #else
962  libmesh_assert (this->calculate_nothing ||
963  this->calculate_phi || this->calculate_dphi ||
964  this->calculate_dphiref ||
965  this->calculate_curl_phi || this->calculate_div_phi ||
966  this->calculate_map);
967 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
968 #endif // LIBMESH_ENABLE_DEPRECATED
969 
970  // Request whichever terms are necessary from the FEMap
971  if (this->calculate_phi)
972  this->_fe_trans->init_map_phi(*this);
973 
974  if (this->calculate_dphiref)
975  this->_fe_trans->init_map_dphi(*this);
976 
977 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
978  if (this->calculate_d2phi)
979  this->_fe_trans->init_map_d2phi(*this);
980 #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
981 }
982 
983 
984 
985 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
986 
987 
988 template <typename OutputType>
989 void FEGenericBase<OutputType>::print_d2phi(std::ostream & os) const
990 {
991  for (auto i : index_range(dphi))
992  for (auto j : index_range(dphi[i]))
993  os << " d2phi[" << i << "][" << j << "]=" << d2phi[i][j];
994 }
995 
996 template <typename OutputType>
997 void FEGenericBase<OutputType>::print_dual_d2phi(std::ostream & os) const
998 {
999  for (auto i : index_range(dual_d2phi))
1000  for (auto j : index_range(dual_d2phi[i]))
1001  os << " dual_d2phi[" << i << "][" << j << "]=" << dual_d2phi[i][j];
1002 }
1003 
1004 #endif
1005 
1006 
1007 
1008 #ifdef LIBMESH_ENABLE_AMR
1009 
1010 template <typename OutputType>
1011 void
1012 FEGenericBase<OutputType>::coarsened_dof_values(const NumericVector<Number> & old_vector,
1013  const DofMap & dof_map,
1014  const Elem * elem,
1015  DenseVector<Number> & Ue,
1016  const unsigned int var,
1017  const bool use_old_dof_indices)
1018 {
1019  // Side/edge local DOF indices
1020  std::vector<unsigned int> new_side_dofs, old_side_dofs;
1021 
1022  // FIXME: what about 2D shells in 3D space?
1023  unsigned int dim = elem->dim();
1024 
1025  // Cache n_children(); it's a virtual call but it's const.
1026  const unsigned int n_children = elem->n_children();
1027 
1028  // We use local FE objects for now
1029  // FIXME: we should use more, external objects instead for efficiency
1030  const FEType & base_fe_type = dof_map.variable_type(var);
1031  std::unique_ptr<FEGenericBase<OutputShape>> fe
1032  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1033  std::unique_ptr<FEGenericBase<OutputShape>> fe_coarse
1034  (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1035 
1036  std::unique_ptr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
1037  std::unique_ptr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
1038  std::unique_ptr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
1039  std::vector<Point> coarse_qpoints;
1040 
1041  // The values of the shape functions at the quadrature
1042  // points
1043  const std::vector<std::vector<OutputShape>> & phi_values =
1044  fe->get_phi();
1045  const std::vector<std::vector<OutputShape>> & phi_coarse =
1046  fe_coarse->get_phi();
1047 
1048  // The gradients of the shape functions at the quadrature
1049  // points on the child element.
1050  const std::vector<std::vector<OutputGradient>> * dphi_values =
1051  nullptr;
1052  const std::vector<std::vector<OutputGradient>> * dphi_coarse =
1053  nullptr;
1054 
1055  const FEContinuity cont = fe->get_continuity();
1056 
1057  if (cont == C_ONE)
1058  {
1059  const std::vector<std::vector<OutputGradient>> &
1060  ref_dphi_values = fe->get_dphi();
1061  dphi_values = &ref_dphi_values;
1062  const std::vector<std::vector<OutputGradient>> &
1063  ref_dphi_coarse = fe_coarse->get_dphi();
1064  dphi_coarse = &ref_dphi_coarse;
1065  }
1066 
1067  // The Jacobian * quadrature weight at the quadrature points
1068  const std::vector<Real> & JxW =
1069  fe->get_JxW();
1070 
1071  // The XYZ locations of the quadrature points on the
1072  // child element
1073  const std::vector<Point> & xyz_values =
1074  fe->get_xyz();
1075 
1076  // Number of nodes on parent element
1077  const unsigned int n_nodes = elem->n_nodes();
1078 
1079  // Number of dofs on parent element
1080  const unsigned int new_n_dofs =
1081  FEInterface::n_dofs(base_fe_type, elem->max_descendant_p_level(), elem);
1082 
1083  // Fixed vs. free DoFs on edge/face projections
1084  std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
1085  std::vector<int> free_dof(new_n_dofs, 0);
1086 
1087  DenseMatrix<Real> Ke;
1088  DenseVector<Number> Fe;
1089  Ue.resize(new_n_dofs); Ue.zero();
1090 
1091 
1092  // When coarsening, in general, we need a series of
1093  // projections to ensure a unique and continuous
1094  // solution. We start by interpolating nodes, then
1095  // hold those fixed and project edges, then
1096  // hold those fixed and project faces, then
1097  // hold those fixed and project interiors
1098 
1099  // Copy node values first
1100  {
1101  std::vector<dof_id_type> node_dof_indices;
1102  if (use_old_dof_indices)
1103  dof_map.old_dof_indices (elem, node_dof_indices, var);
1104  else
1105  dof_map.dof_indices (elem, node_dof_indices, var);
1106 
1107  unsigned int current_dof = 0;
1108  for (unsigned int n=0; n!= n_nodes; ++n)
1109  {
1110  // FIXME: this should go through the DofMap,
1111  // not duplicate dof_indices code badly!
1112  const unsigned int my_nc =
1113  FEInterface::n_dofs_at_node (base_fe_type, elem->max_descendant_p_level(), elem, n);
1114  if (!elem->is_vertex(n))
1115  {
1116  current_dof += my_nc;
1117  continue;
1118  }
1119 
1120  // We're assuming here that child n shares vertex n,
1121  // which is wrong on non-simplices right now
1122  // ... but this code isn't necessary except on elements
1123  // where p refinement creates more vertex dofs; we have
1124  // no such elements yet.
1125  int extra_order = 0;
1126  // if (elem->child_ptr(n)->p_level() < elem->p_level())
1127  // extra_order = elem->child_ptr(n)->p_level();
1128  const unsigned int nc =
1129  FEInterface::n_dofs_at_node (base_fe_type, extra_order, elem, n);
1130  for (unsigned int i=0; i!= nc; ++i)
1131  {
1132  Ue(current_dof) =
1133  old_vector(node_dof_indices[current_dof]);
1134  dof_is_fixed[current_dof] = true;
1135  current_dof++;
1136  }
1137  }
1138  }
1139 
1140  FEType fe_type = base_fe_type, temp_fe_type;
1141  fe_type.order = fe_type.order + elem->max_descendant_p_level();
1142 
1143  // In 3D, project any edge values next
1144  if (dim > 2 && cont != DISCONTINUOUS)
1145  for (auto e : elem->edge_index_range())
1146  {
1147  FEInterface::dofs_on_edge(elem, dim, fe_type,
1148  e, new_side_dofs);
1149 
1150  const unsigned int n_new_side_dofs =
1151  cast_int<unsigned int>(new_side_dofs.size());
1152 
1153  // Some edge dofs are on nodes and already
1154  // fixed, others are free to calculate
1155  unsigned int free_dofs = 0;
1156  for (unsigned int i=0; i != n_new_side_dofs; ++i)
1157  if (!dof_is_fixed[new_side_dofs[i]])
1158  free_dof[free_dofs++] = i;
1159  Ke.resize (free_dofs, free_dofs); Ke.zero();
1160  Fe.resize (free_dofs); Fe.zero();
1161  // The new edge coefficients
1162  DenseVector<Number> Uedge(free_dofs);
1163 
1164  // Add projection terms from each child sharing
1165  // this edge
1166  for (unsigned int c=0; c != n_children; ++c)
1167  {
1168  if (!elem->is_child_on_edge(c,e))
1169  continue;
1170  const Elem * child = elem->child_ptr(c);
1171 
1172  std::vector<dof_id_type> child_dof_indices;
1173  if (use_old_dof_indices)
1174  dof_map.old_dof_indices (child,
1175  child_dof_indices, var);
1176  else
1177  dof_map.dof_indices (child,
1178  child_dof_indices, var);
1179  const unsigned int child_n_dofs =
1180  cast_int<unsigned int>
1181  (child_dof_indices.size());
1182 
1183  temp_fe_type = base_fe_type;
1184  temp_fe_type.order = temp_fe_type.order + child->p_level();
1185 
1186  FEInterface::dofs_on_edge(child, dim,
1187  temp_fe_type, e, old_side_dofs);
1188 
1189  // Initialize both child and parent FE data
1190  // on the child's edge
1191  fe->attach_quadrature_rule (qedgerule.get());
1192  fe->edge_reinit (child, e);
1193  const unsigned int n_qp = qedgerule->n_points();
1194 
1195  FEMap::inverse_map (dim, elem, xyz_values,
1196  coarse_qpoints);
1197 
1198  fe_coarse->reinit(elem, &coarse_qpoints);
1199 
1200  // Loop over the quadrature points
1201  for (unsigned int qp=0; qp<n_qp; qp++)
1202  {
1203  // solution value at the quadrature point
1204  OutputNumber fineval = libMesh::zero;
1205  // solution grad at the quadrature point
1206  OutputNumberGradient finegrad;
1207 
1208  // Sum the solution values * the DOF
1209  // values at the quadrature point to
1210  // get the solution value and gradient.
1211  for (unsigned int i=0; i<child_n_dofs;
1212  i++)
1213  {
1214  fineval +=
1215  (old_vector(child_dof_indices[i])*
1216  phi_values[i][qp]);
1217  if (cont == C_ONE)
1218  finegrad += (*dphi_values)[i][qp] *
1219  old_vector(child_dof_indices[i]);
1220  }
1221 
1222  // Form edge projection matrix
1223  for (unsigned int sidei=0, freei=0; sidei != n_new_side_dofs; ++sidei)
1224  {
1225  unsigned int i = new_side_dofs[sidei];
1226  // fixed DoFs aren't test functions
1227  if (dof_is_fixed[i])
1228  continue;
1229  for (unsigned int sidej=0, freej=0; sidej != n_new_side_dofs; ++sidej)
1230  {
1231  unsigned int j =
1232  new_side_dofs[sidej];
1233  if (dof_is_fixed[j])
1234  Fe(freei) -=
1235  TensorTools::inner_product(phi_coarse[i][qp],
1236  phi_coarse[j][qp]) *
1237  JxW[qp] * Ue(j);
1238  else
1239  Ke(freei,freej) +=
1240  TensorTools::inner_product(phi_coarse[i][qp],
1241  phi_coarse[j][qp]) *
1242  JxW[qp];
1243  if (cont == C_ONE)
1244  {
1245  if (dof_is_fixed[j])
1246  Fe(freei) -=
1247  TensorTools::inner_product((*dphi_coarse)[i][qp],
1248  (*dphi_coarse)[j][qp]) *
1249  JxW[qp] * Ue(j);
1250  else
1251  Ke(freei,freej) +=
1252  TensorTools::inner_product((*dphi_coarse)[i][qp],
1253  (*dphi_coarse)[j][qp]) *
1254  JxW[qp];
1255  }
1256  if (!dof_is_fixed[j])
1257  freej++;
1258  }
1259  Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp],
1260  fineval) * JxW[qp];
1261  if (cont == C_ONE)
1262  Fe(freei) +=
1263  TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1264  freei++;
1265  }
1266  }
1267  }
1268  Ke.cholesky_solve(Fe, Uedge);
1269 
1270  // Transfer new edge solutions to element
1271  for (unsigned int i=0; i != free_dofs; ++i)
1272  {
1273  Number & ui = Ue(new_side_dofs[free_dof[i]]);
1274  libmesh_assert(std::abs(ui) < TOLERANCE ||
1275  std::abs(ui - Uedge(i)) < TOLERANCE);
1276  ui = Uedge(i);
1277  dof_is_fixed[new_side_dofs[free_dof[i]]] = true;
1278  }
1279  }
1280 
1281  // Project any side values (edges in 2D, faces in 3D)
1282  if (dim > 1 && cont != DISCONTINUOUS)
1283  for (auto s : elem->side_index_range())
1284  {
1285  FEInterface::dofs_on_side(elem, dim, fe_type,
1286  s, new_side_dofs);
1287 
1288  const unsigned int n_new_side_dofs =
1289  cast_int<unsigned int>(new_side_dofs.size());
1290 
1291  // Some side dofs are on nodes/edges and already
1292  // fixed, others are free to calculate
1293  unsigned int free_dofs = 0;
1294  for (unsigned int i=0; i != n_new_side_dofs; ++i)
1295  if (!dof_is_fixed[new_side_dofs[i]])
1296  free_dof[free_dofs++] = i;
1297  Ke.resize (free_dofs, free_dofs); Ke.zero();
1298  Fe.resize (free_dofs); Fe.zero();
1299  // The new side coefficients
1300  DenseVector<Number> Uside(free_dofs);
1301 
1302  // Add projection terms from each child sharing
1303  // this side
1304  for (unsigned int c=0; c != n_children; ++c)
1305  {
1306  if (!elem->is_child_on_side(c,s))
1307  continue;
1308  const Elem * child = elem->child_ptr(c);
1309 
1310  std::vector<dof_id_type> child_dof_indices;
1311  if (use_old_dof_indices)
1312  dof_map.old_dof_indices (child,
1313  child_dof_indices, var);
1314  else
1315  dof_map.dof_indices (child,
1316  child_dof_indices, var);
1317  const unsigned int child_n_dofs =
1318  cast_int<unsigned int>
1319  (child_dof_indices.size());
1320 
1321  temp_fe_type = base_fe_type;
1322  temp_fe_type.order = temp_fe_type.order + child->p_level();
1323 
1324  FEInterface::dofs_on_side(child, dim,
1325  temp_fe_type, s, old_side_dofs);
1326 
1327  // Initialize both child and parent FE data
1328  // on the child's side
1329  fe->attach_quadrature_rule (qsiderule.get());
1330  fe->reinit (child, s);
1331  const unsigned int n_qp = qsiderule->n_points();
1332 
1333  FEMap::inverse_map (dim, elem, xyz_values,
1334  coarse_qpoints);
1335 
1336  fe_coarse->reinit(elem, &coarse_qpoints);
1337 
1338  // Loop over the quadrature points
1339  for (unsigned int qp=0; qp<n_qp; qp++)
1340  {
1341  // solution value at the quadrature point
1342  OutputNumber fineval = libMesh::zero;
1343  // solution grad at the quadrature point
1344  OutputNumberGradient finegrad;
1345 
1346  // Sum the solution values * the DOF
1347  // values at the quadrature point to
1348  // get the solution value and gradient.
1349  for (unsigned int i=0; i<child_n_dofs;
1350  i++)
1351  {
1352  fineval +=
1353  old_vector(child_dof_indices[i]) *
1354  phi_values[i][qp];
1355  if (cont == C_ONE)
1356  finegrad += (*dphi_values)[i][qp] *
1357  old_vector(child_dof_indices[i]);
1358  }
1359 
1360  // Form side projection matrix
1361  for (unsigned int sidei=0, freei=0; sidei != n_new_side_dofs; ++sidei)
1362  {
1363  unsigned int i = new_side_dofs[sidei];
1364  // fixed DoFs aren't test functions
1365  if (dof_is_fixed[i])
1366  continue;
1367  for (unsigned int sidej=0, freej=0; sidej != n_new_side_dofs; ++sidej)
1368  {
1369  unsigned int j =
1370  new_side_dofs[sidej];
1371  if (dof_is_fixed[j])
1372  Fe(freei) -=
1373  TensorTools::inner_product(phi_coarse[i][qp],
1374  phi_coarse[j][qp]) *
1375  JxW[qp] * Ue(j);
1376  else
1377  Ke(freei,freej) +=
1378  TensorTools::inner_product(phi_coarse[i][qp],
1379  phi_coarse[j][qp]) *
1380  JxW[qp];
1381  if (cont == C_ONE)
1382  {
1383  if (dof_is_fixed[j])
1384  Fe(freei) -=
1385  TensorTools::inner_product((*dphi_coarse)[i][qp],
1386  (*dphi_coarse)[j][qp]) *
1387  JxW[qp] * Ue(j);
1388  else
1389  Ke(freei,freej) +=
1390  TensorTools::inner_product((*dphi_coarse)[i][qp],
1391  (*dphi_coarse)[j][qp]) *
1392  JxW[qp];
1393  }
1394  if (!dof_is_fixed[j])
1395  freej++;
1396  }
1397  Fe(freei) += TensorTools::inner_product(fineval, phi_coarse[i][qp]) * JxW[qp];
1398  if (cont == C_ONE)
1399  Fe(freei) +=
1400  TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1401  freei++;
1402  }
1403  }
1404  }
1405  Ke.cholesky_solve(Fe, Uside);
1406 
1407  // Transfer new side solutions to element
1408  for (unsigned int i=0; i != free_dofs; ++i)
1409  {
1410  Number & ui = Ue(new_side_dofs[free_dof[i]]);
1411  libmesh_assert(std::abs(ui) < TOLERANCE ||
1412  std::abs(ui - Uside(i)) < TOLERANCE);
1413  ui = Uside(i);
1414  dof_is_fixed[new_side_dofs[free_dof[i]]] = true;
1415  }
1416  }
1417 
1418  // Project the interior values, finally
1419 
1420  // Some interior dofs are on nodes/edges/sides and
1421  // already fixed, others are free to calculate
1422  unsigned int free_dofs = 0;
1423  for (unsigned int i=0; i != new_n_dofs; ++i)
1424  if (!dof_is_fixed[i])
1425  free_dof[free_dofs++] = i;
1426  Ke.resize (free_dofs, free_dofs); Ke.zero();
1427  Fe.resize (free_dofs); Fe.zero();
1428  // The new interior coefficients
1429  DenseVector<Number> Uint(free_dofs);
1430 
1431  // Add projection terms from each child
1432  for (auto & child : elem->child_ref_range())
1433  {
1434  std::vector<dof_id_type> child_dof_indices;
1435  if (use_old_dof_indices)
1436  dof_map.old_dof_indices (&child,
1437  child_dof_indices, var);
1438  else
1439  dof_map.dof_indices (&child,
1440  child_dof_indices, var);
1441  const unsigned int child_n_dofs =
1442  cast_int<unsigned int>
1443  (child_dof_indices.size());
1444 
1445  // Initialize both child and parent FE data
1446  // on the child's quadrature points
1447  fe->attach_quadrature_rule (qrule.get());
1448  fe->reinit (&child);
1449  const unsigned int n_qp = qrule->n_points();
1450 
1451  FEMap::inverse_map (dim, elem, xyz_values, coarse_qpoints);
1452 
1453  fe_coarse->reinit(elem, &coarse_qpoints);
1454 
1455  // Loop over the quadrature points
1456  for (unsigned int qp=0; qp<n_qp; qp++)
1457  {
1458  // solution value at the quadrature point
1459  OutputNumber fineval = libMesh::zero;
1460  // solution grad at the quadrature point
1461  OutputNumberGradient finegrad;
1462 
1463  // Sum the solution values * the DOF
1464  // values at the quadrature point to
1465  // get the solution value and gradient.
1466  for (unsigned int i=0; i<child_n_dofs; i++)
1467  {
1468  fineval +=
1469  (old_vector(child_dof_indices[i]) *
1470  phi_values[i][qp]);
1471  if (cont == C_ONE)
1472  finegrad += (*dphi_values)[i][qp] *
1473  old_vector(child_dof_indices[i]);
1474  }
1475 
1476  // Form interior projection matrix
1477  for (unsigned int i=0, freei=0;
1478  i != new_n_dofs; ++i)
1479  {
1480  // fixed DoFs aren't test functions
1481  if (dof_is_fixed[i])
1482  continue;
1483  for (unsigned int j=0, freej=0; j !=
1484  new_n_dofs; ++j)
1485  {
1486  if (dof_is_fixed[j])
1487  Fe(freei) -=
1488  TensorTools::inner_product(phi_coarse[i][qp],
1489  phi_coarse[j][qp]) *
1490  JxW[qp] * Ue(j);
1491  else
1492  Ke(freei,freej) +=
1493  TensorTools::inner_product(phi_coarse[i][qp],
1494  phi_coarse[j][qp]) *
1495  JxW[qp];
1496  if (cont == C_ONE)
1497  {
1498  if (dof_is_fixed[j])
1499  Fe(freei) -=
1500  TensorTools::inner_product((*dphi_coarse)[i][qp],
1501  (*dphi_coarse)[j][qp]) *
1502  JxW[qp] * Ue(j);
1503  else
1504  Ke(freei,freej) +=
1505  TensorTools::inner_product((*dphi_coarse)[i][qp],
1506  (*dphi_coarse)[j][qp]) *
1507  JxW[qp];
1508  }
1509  if (!dof_is_fixed[j])
1510  freej++;
1511  }
1512  Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp], fineval) *
1513  JxW[qp];
1514  if (cont == C_ONE)
1515  Fe(freei) += TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1516  freei++;
1517  }
1518  }
1519  }
1520  Ke.cholesky_solve(Fe, Uint);
1521 
1522  // Transfer new interior solutions to element
1523  for (unsigned int i=0; i != free_dofs; ++i)
1524  {
1525  Number & ui = Ue(free_dof[i]);
1526  libmesh_assert(std::abs(ui) < TOLERANCE ||
1527  std::abs(ui - Uint(i)) < TOLERANCE);
1528  ui = Uint(i);
1529  // We should be fixing all dofs by now; no need to keep track of
1530  // that unless we're debugging
1531 #ifndef NDEBUG
1532  dof_is_fixed[free_dof[i]] = true;
1533 #endif
1534  }
1535 
1536 #ifndef NDEBUG
1537  // Make sure every DoF got reached!
1538  for (unsigned int i=0; i != new_n_dofs; ++i)
1539  libmesh_assert(dof_is_fixed[i]);
1540 #endif
1541 }
1542 
1543 
1544 
1545 template <typename OutputType>
1546 void
1547 FEGenericBase<OutputType>::coarsened_dof_values(const NumericVector<Number> & old_vector,
1548  const DofMap & dof_map,
1549  const Elem * elem,
1550  DenseVector<Number> & Ue,
1551  const bool use_old_dof_indices)
1552 {
1553  Ue.resize(0);
1554 
1555  for (auto v : make_range(dof_map.n_variables()))
1556  {
1557  DenseVector<Number> Usub;
1558 
1559  coarsened_dof_values(old_vector, dof_map, elem, Usub,
1560  v, use_old_dof_indices);
1561 
1562  Ue.append (Usub);
1563  }
1564 }
1565 
1566 
1567 
1568 template <typename OutputType>
1569 void
1570 FEGenericBase<OutputType>::compute_proj_constraints (DofConstraints & constraints,
1571  DofMap & dof_map,
1572  const unsigned int variable_number,
1573  const Elem * elem)
1574 {
1575  libmesh_assert(elem);
1576 
1577  const unsigned int Dim = elem->dim();
1578 
1579  // Only constrain elements in 2,3D.
1580  if (Dim == 1)
1581  return;
1582 
1583  // Only constrain active elements with this method
1584  if (!elem->active())
1585  return;
1586 
1587  const Variable & var = dof_map.variable(variable_number);
1588  const FEType & base_fe_type = var.type();
1589  const bool add_p_level = dof_map.should_p_refine_var(variable_number);
1590 
1591  // Construct FE objects for this element and its neighbors.
1592  std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1593  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1594  my_fe->add_p_level_in_reinit(add_p_level);
1595  const FEContinuity cont = my_fe->get_continuity();
1596 
1597  // We don't need to constrain discontinuous elements
1598  if (cont == DISCONTINUOUS)
1599  return;
1600  libmesh_assert (cont == C_ZERO || cont == C_ONE ||
1601  cont == SIDE_DISCONTINUOUS);
1602 
1603  // this would require some generalisation:
1604  // - e.g. the 'my_fe'-object needs generalisation
1605  // - due to lack of one-to-one correspondence of DOFs and nodes,
1606  // this doesn't work easily.
1607  if (elem->infinite())
1608  libmesh_not_implemented();
1609 
1610  std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1611  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1612  neigh_fe->add_p_level_in_reinit(add_p_level);
1613 
1614  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1615  my_fe->attach_quadrature_rule (&my_qface);
1616  std::vector<Point> neigh_qface;
1617 
1618  const std::vector<Real> & JxW = my_fe->get_JxW();
1619  const std::vector<Point> & q_point = my_fe->get_xyz();
1620  const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1621  const std::vector<std::vector<OutputShape>> & neigh_phi =
1622  neigh_fe->get_phi();
1623  const std::vector<Point> * face_normals = nullptr;
1624  const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1625  const std::vector<std::vector<OutputGradient>> * neigh_dphi = nullptr;
1626 
1627  std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1628  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1629 
1630  if (cont == C_ONE)
1631  {
1632  const std::vector<Point> & ref_face_normals =
1633  my_fe->get_normals();
1634  face_normals = &ref_face_normals;
1635  const std::vector<std::vector<OutputGradient>> & ref_dphi =
1636  my_fe->get_dphi();
1637  dphi = &ref_dphi;
1638  const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1639  neigh_fe->get_dphi();
1640  neigh_dphi = &ref_neigh_dphi;
1641  }
1642 
1643  DenseMatrix<Real> Ke;
1644  DenseVector<Real> Fe;
1645  std::vector<DenseVector<Real>> Ue;
1646 
1647  // Look at the element faces. Check to see if we need to
1648  // build constraints.
1649  for (auto s : elem->side_index_range())
1650  {
1651  // Get pointers to the element's neighbor.
1652  const Elem * neigh = elem->neighbor_ptr(s);
1653 
1654  if (!neigh)
1655  continue;
1656 
1657  if (!var.active_on_subdomain(neigh->subdomain_id()))
1658  continue;
1659 
1660  // h refinement constraints:
1661  // constrain dofs shared between
1662  // this element and ones coarser
1663  // than this element.
1664  if (neigh->level() < elem->level())
1665  {
1666  unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
1667  libmesh_assert_less (s_neigh, neigh->n_neighbors());
1668 
1669  // Find the minimum p level; we build the h constraint
1670  // matrix with this and then constrain away all higher p
1671  // DoFs.
1672  libmesh_assert(neigh->active());
1673  const unsigned int min_p_level = add_p_level *
1674  std::min(elem->p_level(), neigh->p_level());
1675  // we may need to make the FE objects reinit with the
1676  // minimum shared p_level
1677  const unsigned int old_elem_level = add_p_level * elem->p_level();
1678  if (old_elem_level != min_p_level)
1679  my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + min_p_level - old_elem_level);
1680  const unsigned int old_neigh_level = add_p_level * neigh->p_level();
1681  if (old_neigh_level != min_p_level)
1682  neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + min_p_level - old_neigh_level);
1683 
1684  my_fe->reinit(elem, s);
1685 
1686  // This function gets called element-by-element, so there
1687  // will be a lot of memory allocation going on. We can
1688  // at least minimize this for the case of the dof indices
1689  // by efficiently preallocating the requisite storage.
1690  // n_nodes is not necessarily n_dofs, but it is better
1691  // than nothing!
1692  my_dof_indices.reserve (elem->n_nodes());
1693  neigh_dof_indices.reserve (neigh->n_nodes());
1694 
1695  dof_map.dof_indices (elem, my_dof_indices,
1696  variable_number,
1697  min_p_level);
1698  dof_map.dof_indices (neigh, neigh_dof_indices,
1699  variable_number,
1700  min_p_level);
1701 
1702  const unsigned int n_qp = my_qface.n_points();
1703 
1704  FEMap::inverse_map (Dim, neigh, q_point, neigh_qface);
1705 
1706  neigh_fe->reinit(neigh, &neigh_qface);
1707 
1708  // We're only concerned with DOFs whose values (and/or first
1709  // derivatives for C1 elements) are supported on side nodes
1710  FEType elem_fe_type = base_fe_type;
1711  if (old_elem_level != min_p_level)
1712  elem_fe_type.order = base_fe_type.order.get_order() + min_p_level - old_elem_level;
1713  FEType neigh_fe_type = base_fe_type;
1714  if (old_neigh_level != min_p_level)
1715  neigh_fe_type.order = base_fe_type.order.get_order() + min_p_level - old_neigh_level;
1716  FEInterface::dofs_on_side(elem, Dim, elem_fe_type, s, my_side_dofs);
1717  FEInterface::dofs_on_side(neigh, Dim, neigh_fe_type, s_neigh, neigh_side_dofs);
1718 
1719  const unsigned int n_side_dofs =
1720  cast_int<unsigned int>(my_side_dofs.size());
1721  libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
1722 
1723 #ifndef NDEBUG
1724  for (auto i : my_side_dofs)
1725  libmesh_assert_less(i, my_dof_indices.size());
1726  for (auto i : neigh_side_dofs)
1727  libmesh_assert_less(i, neigh_dof_indices.size());
1728 #endif
1729 
1730  Ke.resize (n_side_dofs, n_side_dofs);
1731  Ue.resize(n_side_dofs);
1732 
1733  // Form the projection matrix, (inner product of fine basis
1734  // functions against fine test functions)
1735  for (unsigned int is = 0; is != n_side_dofs; ++is)
1736  {
1737  const unsigned int i = my_side_dofs[is];
1738  for (unsigned int js = 0; js != n_side_dofs; ++js)
1739  {
1740  const unsigned int j = my_side_dofs[js];
1741  for (unsigned int qp = 0; qp != n_qp; ++qp)
1742  {
1743  Ke(is,js) += JxW[qp] * TensorTools::inner_product(phi[i][qp], phi[j][qp]);
1744  if (cont == C_ONE)
1745  Ke(is,js) += JxW[qp] *
1746  TensorTools::inner_product((*dphi)[i][qp] *
1747  (*face_normals)[qp],
1748  (*dphi)[j][qp] *
1749  (*face_normals)[qp]);
1750  }
1751  }
1752  }
1753 
1754  // Form the right hand sides, (inner product of coarse basis
1755  // functions against fine test functions)
1756  for (unsigned int is = 0; is != n_side_dofs; ++is)
1757  {
1758  const unsigned int i = neigh_side_dofs[is];
1759  Fe.resize (n_side_dofs);
1760  for (unsigned int js = 0; js != n_side_dofs; ++js)
1761  {
1762  const unsigned int j = my_side_dofs[js];
1763  for (unsigned int qp = 0; qp != n_qp; ++qp)
1764  {
1765  Fe(js) += JxW[qp] *
1766  TensorTools::inner_product(neigh_phi[i][qp],
1767  phi[j][qp]);
1768  if (cont == C_ONE)
1769  Fe(js) += JxW[qp] *
1770  TensorTools::inner_product((*neigh_dphi)[i][qp] *
1771  (*face_normals)[qp],
1772  (*dphi)[j][qp] *
1773  (*face_normals)[qp]);
1774  }
1775  }
1776  Ke.cholesky_solve(Fe, Ue[is]);
1777  }
1778 
1779  for (unsigned int js = 0; js != n_side_dofs; ++js)
1780  {
1781  const unsigned int j = my_side_dofs[js];
1782  const dof_id_type my_dof_g = my_dof_indices[j];
1783  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
1784 
1785  // Hunt for "constraining against myself" cases before
1786  // we bother creating a constraint row
1787  bool self_constraint = false;
1788  for (unsigned int is = 0; is != n_side_dofs; ++is)
1789  {
1790  const unsigned int i = neigh_side_dofs[is];
1791  const dof_id_type their_dof_g = neigh_dof_indices[i];
1792  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1793 
1794  if (their_dof_g == my_dof_g)
1795  {
1796 #ifndef NDEBUG
1797  const Real their_dof_value = Ue[is](js);
1798  libmesh_assert_less (std::abs(their_dof_value-1.),
1799  10*TOLERANCE);
1800 
1801  for (unsigned int k = 0; k != n_side_dofs; ++k)
1802  libmesh_assert(k == is ||
1803  std::abs(Ue[k](js)) <
1804  10*TOLERANCE);
1805 #endif
1806 
1807  self_constraint = true;
1808  break;
1809  }
1810  }
1811 
1812  if (self_constraint)
1813  continue;
1814 
1815  DofConstraintRow * constraint_row;
1816 
1817  // we may be running constraint methods concurrently
1818  // on multiple threads, so we need a lock to
1819  // ensure that this constraint is "ours"
1820  {
1821  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1822 
1823  if (dof_map.is_constrained_dof(my_dof_g))
1824  continue;
1825 
1826  constraint_row = &(constraints[my_dof_g]);
1827  libmesh_assert(constraint_row->empty());
1828  }
1829 
1830  for (unsigned int is = 0; is != n_side_dofs; ++is)
1831  {
1832  const unsigned int i = neigh_side_dofs[is];
1833  const dof_id_type their_dof_g = neigh_dof_indices[i];
1834  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1835  libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
1836 
1837  const Real their_dof_value = Ue[is](js);
1838 
1839  if (std::abs(their_dof_value) < 10*TOLERANCE)
1840  continue;
1841 
1842  constraint_row->emplace(their_dof_g, their_dof_value);
1843  }
1844  }
1845 
1846  my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + old_elem_level - min_p_level);
1847  neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + old_neigh_level - min_p_level);
1848  }
1849 
1850  if (add_p_level)
1851  {
1852  // p refinement constraints:
1853  // constrain dofs shared between
1854  // active elements and neighbors with
1855  // lower polynomial degrees
1856  const unsigned int min_p_level =
1857  neigh->min_p_level_by_neighbor(elem, elem->p_level());
1858  if (min_p_level < elem->p_level())
1859  {
1860  // Adaptive p refinement of non-hierarchic bases will
1861  // require more coding
1862  libmesh_assert(my_fe->is_hierarchic());
1863  dof_map.constrain_p_dofs(variable_number, elem,
1864  s, min_p_level);
1865  }
1866  }
1867  }
1868 }
1869 
1870 #endif // #ifdef LIBMESH_ENABLE_AMR
1871 
1872 
1873 
1874 #ifdef LIBMESH_ENABLE_PERIODIC
1875 template <typename OutputType>
1876 void
1877 FEGenericBase<OutputType>::
1878 compute_periodic_constraints (DofConstraints & constraints,
1879  DofMap & dof_map,
1880  const PeriodicBoundaries & boundaries,
1881  const MeshBase & mesh,
1882  const PointLocatorBase * point_locator,
1883  const unsigned int variable_number,
1884  const Elem * elem)
1885 {
1886  // Only bother if we truly have periodic boundaries
1887  if (boundaries.empty())
1888  return;
1889 
1890  libmesh_assert(elem);
1891 
1892  // Only constrain active elements with this method
1893  if (!elem->active())
1894  return;
1895 
1896  if (elem->infinite())
1897  libmesh_not_implemented();
1898 
1899  const unsigned int Dim = elem->dim();
1900 
1901  // We need sys_number and variable_number for DofObject methods
1902  // later
1903  const unsigned int sys_number = dof_map.sys_number();
1904 
1905  const FEType & base_fe_type = dof_map.variable_type(variable_number);
1906 
1907  // Construct FE objects for this element and its pseudo-neighbors.
1908  std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1909  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1910  const FEContinuity cont = my_fe->get_continuity();
1911 
1912  // We don't need to constrain discontinuous elements
1913  if (cont == DISCONTINUOUS)
1914  return;
1915  libmesh_assert (cont == C_ZERO || cont == C_ONE);
1916 
1917  // We'll use element size to generate relative tolerances later
1918  const Real primary_hmin = elem->hmin();
1919 
1920  std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1921  (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1922 
1923  QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1924  my_fe->attach_quadrature_rule (&my_qface);
1925  std::vector<Point> neigh_qface;
1926 
1927  const std::vector<Real> & JxW = my_fe->get_JxW();
1928  const std::vector<Point> & q_point = my_fe->get_xyz();
1929  const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1930  const std::vector<std::vector<OutputShape>> & neigh_phi =
1931  neigh_fe->get_phi();
1932  const std::vector<Point> * face_normals = nullptr;
1933  const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1934  const std::vector<std::vector<OutputGradient>> * neigh_dphi = nullptr;
1935  std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1936  std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1937 
1938  if (cont != C_ZERO)
1939  {
1940  const std::vector<Point> & ref_face_normals =
1941  my_fe->get_normals();
1942  face_normals = &ref_face_normals;
1943  const std::vector<std::vector<OutputGradient>> & ref_dphi =
1944  my_fe->get_dphi();
1945  dphi = &ref_dphi;
1946  const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1947  neigh_fe->get_dphi();
1948  neigh_dphi = &ref_neigh_dphi;
1949  }
1950 
1951  DenseMatrix<Real> Ke;
1952  DenseVector<Real> Fe;
1953  std::vector<DenseVector<Real>> Ue;
1954 
1955  // Container to catch the boundary ids that BoundaryInfo hands us.
1956  std::vector<boundary_id_type> bc_ids;
1957 
1958  // Look at the element faces. Check to see if we need to
1959  // build constraints.
1960  const unsigned short int max_ns = elem->n_sides();
1961  for (unsigned short int s = 0; s != max_ns; ++s)
1962  {
1963  if (elem->neighbor_ptr(s))
1964  continue;
1965 
1966  mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
1967 
1968  for (const auto & boundary_id : bc_ids)
1969  {
1970  const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
1971  if (!periodic || !periodic->is_my_variable(variable_number))
1972  continue;
1973 
1974  libmesh_assert(point_locator);
1975 
1976  // Get pointers to the element's neighbor.
1977  unsigned int s_neigh;
1978  const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s, &s_neigh);
1979 
1980  libmesh_error_msg_if(neigh == nullptr,
1981  "PeriodicBoundaries point locator object returned nullptr!");
1982 
1983  // periodic (and possibly h refinement) constraints:
1984  // constrain dofs shared between
1985  // this element and ones as coarse
1986  // as or coarser than this element.
1987  if (neigh->level() <= elem->level())
1988  {
1989 #ifdef LIBMESH_ENABLE_AMR
1990  // Find the minimum p level; we build the h constraint
1991  // matrix with this and then constrain away all higher p
1992  // DoFs.
1993  libmesh_assert(neigh->active());
1994  const unsigned int min_p_level =
1995  std::min(elem->p_level(), neigh->p_level());
1996 
1997  // we may need to make the FE objects reinit with the
1998  // minimum shared p_level
1999  // FIXME - I hate using const_cast<> and avoiding
2000  // accessor functions; there's got to be a
2001  // better way to do this!
2002  const unsigned int old_elem_level = elem->p_level();
2003  if (old_elem_level != min_p_level)
2004  (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
2005  const unsigned int old_neigh_level = neigh->p_level();
2006  if (old_neigh_level != min_p_level)
2007  (const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
2008 #endif // #ifdef LIBMESH_ENABLE_AMR
2009 
2010  // We can do a projection with a single integration,
2011  // due to the assumption of nested finite element
2012  // subspaces.
2013  // FIXME: it might be more efficient to do nodes,
2014  // then edges, then side, to reduce the size of the
2015  // Cholesky factorization(s)
2016  my_fe->reinit(elem, s);
2017 
2018  dof_map.dof_indices (elem, my_dof_indices,
2019  variable_number);
2020  dof_map.dof_indices (neigh, neigh_dof_indices,
2021  variable_number);
2022 
2023  // We use neigh_dof_indices_all_variables in the case that the
2024  // periodic boundary condition involves mappings between multiple
2025  // variables.
2026  std::vector<std::vector<dof_id_type>> neigh_dof_indices_all_variables;
2027  if(periodic->has_transformation_matrix())
2028  {
2029  const std::set<unsigned int> & variables = periodic->get_variables();
2030  neigh_dof_indices_all_variables.resize(variables.size());
2031  unsigned int index = 0;
2032  for(unsigned int var : variables)
2033  {
2034  dof_map.dof_indices (neigh, neigh_dof_indices_all_variables[index],
2035  var);
2036  index++;
2037  }
2038  }
2039 
2040  const unsigned int n_qp = my_qface.n_points();
2041 
2042  // Translate the quadrature points over to the
2043  // neighbor's boundary
2044  std::vector<Point> neigh_point(q_point.size());
2045  for (auto i : index_range(neigh_point))
2046  neigh_point[i] = periodic->get_corresponding_pos(q_point[i]);
2047 
2048  FEMap::inverse_map (Dim, neigh, neigh_point,
2049  neigh_qface);
2050 
2051  neigh_fe->reinit(neigh, &neigh_qface);
2052 
2053  // We're only concerned with DOFs whose values (and/or first
2054  // derivatives for C1 elements) are supported on side nodes
2055  FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
2056  FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
2057 
2058  // We're done with functions that examine Elem::p_level(),
2059  // so let's unhack those levels
2060 #ifdef LIBMESH_ENABLE_AMR
2061  if (elem->p_level() != old_elem_level)
2062  (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
2063  if (neigh->p_level() != old_neigh_level)
2064  (const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
2065 #endif // #ifdef LIBMESH_ENABLE_AMR
2066 
2067  const unsigned int n_side_dofs =
2068  cast_int<unsigned int>
2069  (my_side_dofs.size());
2070  libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
2071 
2072  Ke.resize (n_side_dofs, n_side_dofs);
2073  Ue.resize(n_side_dofs);
2074 
2075  // Form the projection matrix, (inner product of fine basis
2076  // functions against fine test functions)
2077  for (unsigned int is = 0; is != n_side_dofs; ++is)
2078  {
2079  const unsigned int i = my_side_dofs[is];
2080  for (unsigned int js = 0; js != n_side_dofs; ++js)
2081  {
2082  const unsigned int j = my_side_dofs[js];
2083  for (unsigned int qp = 0; qp != n_qp; ++qp)
2084  {
2085  Ke(is,js) += JxW[qp] *
2086  TensorTools::inner_product(phi[i][qp],
2087  phi[j][qp]);
2088  if (cont != C_ZERO)
2089  Ke(is,js) += JxW[qp] *
2090  TensorTools::inner_product((*dphi)[i][qp] *
2091  (*face_normals)[qp],
2092  (*dphi)[j][qp] *
2093  (*face_normals)[qp]);
2094  }
2095  }
2096  }
2097 
2098  // Form the right hand sides, (inner product of coarse basis
2099  // functions against fine test functions)
2100  for (unsigned int is = 0; is != n_side_dofs; ++is)
2101  {
2102  const unsigned int i = neigh_side_dofs[is];
2103  Fe.resize (n_side_dofs);
2104  for (unsigned int js = 0; js != n_side_dofs; ++js)
2105  {
2106  const unsigned int j = my_side_dofs[js];
2107  for (unsigned int qp = 0; qp != n_qp; ++qp)
2108  {
2109  Fe(js) += JxW[qp] *
2110  TensorTools::inner_product(neigh_phi[i][qp],
2111  phi[j][qp]);
2112  if (cont != C_ZERO)
2113  Fe(js) += JxW[qp] *
2114  TensorTools::inner_product((*neigh_dphi)[i][qp] *
2115  (*face_normals)[qp],
2116  (*dphi)[j][qp] *
2117  (*face_normals)[qp]);
2118  }
2119  }
2120  Ke.cholesky_solve(Fe, Ue[is]);
2121  }
2122 
2123  // Make sure we're not adding recursive constraints
2124  // due to the redundancy in the way we add periodic
2125  // boundary constraints
2126  //
2127  // In order for this to work while threaded or on
2128  // distributed meshes, we need a rigorous way to
2129  // avoid recursive constraints. Here it is:
2130  //
2131  // For vertex DoFs, if there is a "prior" element
2132  // (i.e. a coarser element or an equally refined
2133  // element with a lower id) on this boundary which
2134  // contains the vertex point, then we will avoid
2135  // generating constraints; the prior element (or
2136  // something prior to it) may do so. If we are the
2137  // most prior (or "primary") element on this
2138  // boundary sharing this point, then we look at the
2139  // boundary periodic to us, we find the primary
2140  // element there, and if that primary is coarser or
2141  // equal-but-lower-id, then our vertex dofs are
2142  // constrained in terms of that element.
2143  //
2144  // For edge DoFs, if there is a coarser element
2145  // on this boundary sharing this edge, then we will
2146  // avoid generating constraints (we will be
2147  // constrained indirectly via AMR constraints
2148  // connecting us to the coarser element's DoFs). If
2149  // we are the coarsest element sharing this edge,
2150  // then we generate constraints if and only if we
2151  // are finer than the coarsest element on the
2152  // boundary periodic to us sharing the corresponding
2153  // periodic edge, or if we are at equal level but
2154  // our edge nodes have higher ids than the periodic
2155  // edge nodes (sorted from highest to lowest, then
2156  // compared lexicographically)
2157  //
2158  // For face DoFs, we generate constraints if we are
2159  // finer than our periodic neighbor, or if we are at
2160  // equal level but our element id is higher than its
2161  // element id.
2162  //
2163  // If the primary neighbor is also the current elem
2164  // (a 1-element-thick mesh) then we choose which
2165  // vertex dofs to constrain via lexicographic
2166  // ordering on point locations
2167 
2168  // FIXME: This code doesn't yet properly handle
2169  // cases where multiple different periodic BCs
2170  // intersect.
2171  std::set<dof_id_type> my_constrained_dofs;
2172 
2173  // Container to catch boundary IDs handed back by BoundaryInfo.
2174  std::vector<boundary_id_type> new_bc_ids;
2175 
2176  for (auto n : elem->node_index_range())
2177  {
2178  if (!elem->is_node_on_side(n,s))
2179  continue;
2180 
2181  const Node & my_node = elem->node_ref(n);
2182 
2183  if (elem->is_vertex(n))
2184  {
2185  // Find all boundary ids that include this
2186  // point and have periodic boundary
2187  // conditions for this variable
2188  std::set<boundary_id_type> point_bcids;
2189 
2190  for (unsigned int new_s = 0;
2191  new_s != max_ns; ++new_s)
2192  {
2193  if (!elem->is_node_on_side(n,new_s))
2194  continue;
2195 
2196  mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
2197 
2198  for (const auto & new_boundary_id : new_bc_ids)
2199  {
2200  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2201  if (new_periodic && new_periodic->is_my_variable(variable_number))
2202  point_bcids.insert(new_boundary_id);
2203  }
2204  }
2205 
2206  // See if this vertex has point neighbors to
2207  // defer to
2208  if (primary_boundary_point_neighbor
2209  (elem, my_node, mesh.get_boundary_info(), point_bcids)
2210  != elem)
2211  continue;
2212 
2213  // Find the complementary boundary id set
2214  std::set<boundary_id_type> point_pairedids;
2215  for (const auto & new_boundary_id : point_bcids)
2216  {
2217  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2218  point_pairedids.insert(new_periodic->pairedboundary);
2219  }
2220 
2221  // What do we want to constrain against?
2222  const Elem * primary_elem = nullptr;
2223  const Elem * main_neigh = nullptr;
2224  Point main_pt = my_node,
2225  primary_pt = my_node;
2226 
2227  for (const auto & new_boundary_id : point_bcids)
2228  {
2229  // Find the corresponding periodic point and
2230  // its primary neighbor
2231  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2232 
2233  const Point neigh_pt =
2234  new_periodic->get_corresponding_pos(my_node);
2235 
2236  // If the point is getting constrained
2237  // to itself by this PBC then we don't
2238  // generate any constraints
2239  if (neigh_pt.absolute_fuzzy_equals
2240  (my_node, primary_hmin*TOLERANCE))
2241  continue;
2242 
2243  // Otherwise we'll have a constraint in
2244  // one direction or another
2245  if (!primary_elem)
2246  primary_elem = elem;
2247 
2248  const Elem * primary_neigh =
2249  primary_boundary_point_neighbor(neigh, neigh_pt,
2250  mesh.get_boundary_info(),
2251  point_pairedids);
2252 
2253  libmesh_assert(primary_neigh);
2254 
2255  if (new_boundary_id == boundary_id)
2256  {
2257  main_neigh = primary_neigh;
2258  main_pt = neigh_pt;
2259  }
2260 
2261  // Finer elements will get constrained in
2262  // terms of coarser neighbors, not the
2263  // other way around
2264  if ((primary_neigh->level() > primary_elem->level()) ||
2265 
2266  // For equal-level elements, the one with
2267  // higher id gets constrained in terms of
2268  // the one with lower id
2269  (primary_neigh->level() == primary_elem->level() &&
2270  primary_neigh->id() > primary_elem->id()) ||
2271 
2272  // On a one-element-thick mesh, we compare
2273  // points to see what side gets constrained
2274  (primary_neigh == primary_elem &&
2275  (neigh_pt > primary_pt)))
2276  continue;
2277 
2278  primary_elem = primary_neigh;
2279  primary_pt = neigh_pt;
2280  }
2281 
2282  if (!primary_elem ||
2283  primary_elem != main_neigh ||
2284  primary_pt != main_pt)
2285  continue;
2286  }
2287  else if (elem->is_edge(n))
2288  {
2289  // Find which edge we're on
2290  unsigned int e=0, ne = elem->n_edges();
2291  for (; e != ne; ++e)
2292  {
2293  if (elem->is_node_on_edge(n,e))
2294  break;
2295  }
2296  libmesh_assert_less (e, elem->n_edges());
2297 
2298  // Find the edge end nodes
2299  const Node
2300  * e1 = nullptr,
2301  * e2 = nullptr;
2302  for (auto nn : elem->node_index_range())
2303  {
2304  if (nn == n)
2305  continue;
2306 
2307  if (elem->is_node_on_edge(nn, e))
2308  {
2309  if (e1 == nullptr)
2310  {
2311  e1 = elem->node_ptr(nn);
2312  }
2313  else
2314  {
2315  e2 = elem->node_ptr(nn);
2316  break;
2317  }
2318  }
2319  }
2320  libmesh_assert (e1 && e2);
2321 
2322  // Find all boundary ids that include this
2323  // edge and have periodic boundary
2324  // conditions for this variable
2325  std::set<boundary_id_type> edge_bcids;
2326 
2327  for (unsigned int new_s = 0;
2328  new_s != max_ns; ++new_s)
2329  {
2330  if (!elem->is_node_on_side(n,new_s))
2331  continue;
2332 
2333  // We're reusing the new_bc_ids vector created outside the loop over nodes.
2334  mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
2335 
2336  for (const auto & new_boundary_id : new_bc_ids)
2337  {
2338  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2339  if (new_periodic && new_periodic->is_my_variable(variable_number))
2340  edge_bcids.insert(new_boundary_id);
2341  }
2342  }
2343 
2344 
2345  // See if this edge has neighbors to defer to
2346  if (primary_boundary_edge_neighbor
2347  (elem, *e1, *e2, mesh.get_boundary_info(), edge_bcids)
2348  != elem)
2349  continue;
2350 
2351  // Find the complementary boundary id set
2352  std::set<boundary_id_type> edge_pairedids;
2353  for (const auto & new_boundary_id : edge_bcids)
2354  {
2355  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2356  edge_pairedids.insert(new_periodic->pairedboundary);
2357  }
2358 
2359  // What do we want to constrain against?
2360  const Elem * primary_elem = nullptr;
2361  const Elem * main_neigh = nullptr;
2362  Point main_pt1 = *e1,
2363  main_pt2 = *e2,
2364  primary_pt1 = *e1,
2365  primary_pt2 = *e2;
2366 
2367  for (const auto & new_boundary_id : edge_bcids)
2368  {
2369  // Find the corresponding periodic edge and
2370  // its primary neighbor
2371  const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2372 
2373  Point neigh_pt1 = new_periodic->get_corresponding_pos(*e1),
2374  neigh_pt2 = new_periodic->get_corresponding_pos(*e2);
2375 
2376  // If the edge is getting constrained
2377  // to itself by this PBC then we don't
2378  // generate any constraints
2379  if (neigh_pt1.absolute_fuzzy_equals
2380  (*e1, primary_hmin*TOLERANCE) &&
2381  neigh_pt2.absolute_fuzzy_equals
2382  (*e2, primary_hmin*TOLERANCE))
2383  continue;
2384 
2385  // Otherwise we'll have a constraint in
2386  // one direction or another
2387  if (!primary_elem)
2388  primary_elem = elem;
2389 
2390  const Elem * primary_neigh = primary_boundary_edge_neighbor
2391  (neigh, neigh_pt1, neigh_pt2,
2392  mesh.get_boundary_info(), edge_pairedids);
2393 
2394  libmesh_assert(primary_neigh);
2395 
2396  if (new_boundary_id == boundary_id)
2397  {
2398  main_neigh = primary_neigh;
2399  main_pt1 = neigh_pt1;
2400  main_pt2 = neigh_pt2;
2401  }
2402 
2403  // If we have a one-element thick mesh,
2404  // we'll need to sort our points to get a
2405  // consistent ordering rule
2406  //
2407  // Use >= in this test to make sure that,
2408  // for angular constraints, no node gets
2409  // constrained to itself.
2410  if (primary_neigh == primary_elem)
2411  {
2412  if (primary_pt1 > primary_pt2)
2413  std::swap(primary_pt1, primary_pt2);
2414  if (neigh_pt1 > neigh_pt2)
2415  std::swap(neigh_pt1, neigh_pt2);
2416 
2417  if (neigh_pt2 >= primary_pt2)
2418  continue;
2419  }
2420 
2421  // Otherwise:
2422  // Finer elements will get constrained in
2423  // terms of coarser ones, not the other way
2424  // around
2425  if ((primary_neigh->level() > primary_elem->level()) ||
2426 
2427  // For equal-level elements, the one with
2428  // higher id gets constrained in terms of
2429  // the one with lower id
2430  (primary_neigh->level() == primary_elem->level() &&
2431  primary_neigh->id() > primary_elem->id()))
2432  continue;
2433 
2434  primary_elem = primary_neigh;
2435  primary_pt1 = neigh_pt1;
2436  primary_pt2 = neigh_pt2;
2437  }
2438 
2439  if (!primary_elem ||
2440  primary_elem != main_neigh ||
2441  primary_pt1 != main_pt1 ||
2442  primary_pt2 != main_pt2)
2443  continue;
2444  }
2445  else if (elem->is_face(n))
2446  {
2447  // If we have a one-element thick mesh,
2448  // use the ordering of the face node and its
2449  // periodic counterpart to determine what
2450  // gets constrained
2451  if (neigh == elem)
2452  {
2453  const Point neigh_pt =
2454  periodic->get_corresponding_pos(my_node);
2455  if (neigh_pt > my_node)
2456  continue;
2457  }
2458 
2459  // Otherwise:
2460  // Finer elements will get constrained in
2461  // terms of coarser ones, not the other way
2462  // around
2463  if ((neigh->level() > elem->level()) ||
2464 
2465  // For equal-level elements, the one with
2466  // higher id gets constrained in terms of
2467  // the one with lower id
2468  (neigh->level() == elem->level() &&
2469  neigh->id() > elem->id()))
2470  continue;
2471  }
2472 
2473  // If we made it here without hitting a continue
2474  // statement, then we're at a node whose dofs
2475  // should be constrained by this element's
2476  // calculations.
2477  const unsigned int n_comp =
2478  my_node.n_comp(sys_number, variable_number);
2479 
2480  for (unsigned int i=0; i != n_comp; ++i)
2481  my_constrained_dofs.insert
2482  (my_node.dof_number
2483  (sys_number, variable_number, i));
2484  }
2485 
2486  // FIXME: old code for disambiguating periodic BCs:
2487  // this is not threadsafe nor safe to run on a
2488  // non-serialized mesh.
2489  /*
2490  std::vector<bool> recursive_constraint(n_side_dofs, false);
2491 
2492  for (unsigned int is = 0; is != n_side_dofs; ++is)
2493  {
2494  const unsigned int i = neigh_side_dofs[is];
2495  const dof_id_type their_dof_g = neigh_dof_indices[i];
2496  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2497 
2498  {
2499  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2500 
2501  if (!dof_map.is_constrained_dof(their_dof_g))
2502  continue;
2503  }
2504 
2505  DofConstraintRow & their_constraint_row =
2506  constraints[their_dof_g].first;
2507 
2508  for (unsigned int js = 0; js != n_side_dofs; ++js)
2509  {
2510  const unsigned int j = my_side_dofs[js];
2511  const dof_id_type my_dof_g = my_dof_indices[j];
2512  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2513 
2514  if (their_constraint_row.count(my_dof_g))
2515  recursive_constraint[js] = true;
2516  }
2517  }
2518  */
2519 
2520  for (unsigned int js = 0; js != n_side_dofs; ++js)
2521  {
2522  // FIXME: old code path
2523  // if (recursive_constraint[js])
2524  // continue;
2525 
2526  const unsigned int j = my_side_dofs[js];
2527  const dof_id_type my_dof_g = my_dof_indices[j];
2528  libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2529 
2530  // FIXME: new code path
2531  if (!my_constrained_dofs.count(my_dof_g))
2532  continue;
2533 
2534  DofConstraintRow * constraint_row;
2535 
2536  // we may be running constraint methods concurrently
2537  // on multiple threads, so we need a lock to
2538  // ensure that this constraint is "ours"
2539  {
2540  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2541 
2542  if (dof_map.is_constrained_dof(my_dof_g))
2543  continue;
2544 
2545  constraint_row = &(constraints[my_dof_g]);
2546  libmesh_assert(constraint_row->empty());
2547  }
2548 
2549  for (unsigned int is = 0; is != n_side_dofs; ++is)
2550  {
2551  const unsigned int i = neigh_side_dofs[is];
2552  const dof_id_type their_dof_g = neigh_dof_indices[i];
2553  libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2554 
2555  // Periodic constraints should never be
2556  // self-constraints
2557  // libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
2558 
2559  const Real their_dof_value = Ue[is](js);
2560 
2561  if (their_dof_g == my_dof_g)
2562  {
2563  libmesh_assert_less (std::abs(their_dof_value-1.), 1.e-5);
2564  for (unsigned int k = 0; k != n_side_dofs; ++k)
2565  libmesh_assert(k == is || std::abs(Ue[k](js)) < 1.e-5);
2566  continue;
2567  }
2568 
2569  if (std::abs(their_dof_value) < 10*TOLERANCE)
2570  continue;
2571 
2572  if(!periodic->has_transformation_matrix())
2573  {
2574  constraint_row->emplace(their_dof_g, their_dof_value);
2575  }
2576  else
2577  {
2578  // In this case the current variable is constrained in terms of other variables.
2579  // We assume that all variables in this constraint have the same FE type (this
2580  // is asserted below), and hence we can create the constraint row contribution
2581  // by multiplying their_dof_value by the corresponding row of the transformation
2582  // matrix.
2583 
2584  const std::set<unsigned int> & variables = periodic->get_variables();
2585  neigh_dof_indices_all_variables.resize(variables.size());
2586  unsigned int index = 0;
2587  for(unsigned int other_var : variables)
2588  {
2589  libmesh_assert_msg(base_fe_type == dof_map.variable_type(other_var), "FE types must match for all variables involved in constraint");
2590 
2591  Real var_weighting = periodic->get_transformation_matrix()(variable_number, other_var);
2592  constraint_row->emplace(neigh_dof_indices_all_variables[index][i],
2593  var_weighting*their_dof_value);
2594  index++;
2595  }
2596  }
2597 
2598  }
2599  }
2600  }
2601  // p refinement constraints:
2602  // constrain dofs shared between
2603  // active elements and neighbors with
2604  // lower polynomial degrees
2605 #ifdef LIBMESH_ENABLE_AMR
2606  const unsigned int min_p_level =
2607  neigh->min_p_level_by_neighbor(elem, elem->p_level());
2608  if (min_p_level < elem->p_level())
2609  {
2610  // Adaptive p refinement of non-hierarchic bases will
2611  // require more coding
2612  libmesh_assert(my_fe->is_hierarchic());
2613  dof_map.constrain_p_dofs(variable_number, elem,
2614  s, min_p_level);
2615  }
2616 #endif // #ifdef LIBMESH_ENABLE_AMR
2617  }
2618  }
2619 }
2620 
2621 #endif // LIBMESH_ENABLE_PERIODIC
2622 
2623 // ------------------------------------------------------------
2624 // Explicit instantiations
2625 template class LIBMESH_EXPORT FEGenericBase<Real>;
2626 template class LIBMESH_EXPORT FEGenericBase<RealGradient>;
2627 
2628 } // namespace libMesh
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
libMesh::FEType::family
FEFamily family
The type of finite element.
Definition: fe_type.h:221
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:1032
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
libMesh::FEInterface::n_dofs
static unsigned int n_dofs(const unsigned int dim, const FEType &fe_t, const ElemType t)
Definition: fe_interface.C:355
libMesh::DenseVector::zero
virtual void zero() override final
Set every element in the vector to 0.
Definition: dense_vector.h:420
libMesh::DofObject::n_comp
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1002
libMesh::FEGenericBase::print_dphi
virtual void print_dphi(std::ostream &os) const override
Prints the value of each shape function&#39;s derivative at each quadrature point.
Definition: fe_base.C:895
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::Elem::is_face
virtual bool is_face(const unsigned int i) const =0
libMesh::Elem::side_index_range
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2710
libMesh::DofMap::dof_indices
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
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::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
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:1628
libMesh::PeriodicBoundaries
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
Definition: periodic_boundaries.h:49
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::Elem::is_child_on_side
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const =0
libMesh::DenseVector::resize
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
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:611
libMesh::DofMap::variable_type
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2220
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::PeriodicBoundaries::boundary
PeriodicBoundaryBase * boundary(boundary_id_type id)
Definition: periodic_boundaries.C:38
libMesh::FEInterface::field_type
static FEFieldType field_type(const FEType &fe_type)
Definition: fe_interface.C:2642
libMesh::Elem::is_node_on_side
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
libMesh::FEType::default_quadrature_order
Order default_quadrature_order() const
Definition: fe_type.h:371
libMesh::Elem::n_children
virtual unsigned int n_children() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:3108
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:215
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:1083
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:1333
libMesh::VectorValue
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
Definition: exact_solution.h:47
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::zero
const Number zero
.
Definition: libmesh.h:304
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
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:2323
libMesh::DofMap::sys_number
unsigned int sys_number() const
Definition: dof_map.h:2172
libMesh::MeshBase
This is the MeshBase class.
Definition: mesh_base.h:75
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:2346
libMesh::DofMap::variable
const Variable & variable(const unsigned int c) const override
Definition: dof_map.h:2190
libMesh::DofMap::should_p_refine_var
bool should_p_refine_var(unsigned int var) const
Whether the given variable should be p-refined.
Definition: dof_map.h:2437
libMesh::Elem::edge_index_range
IntRange< unsigned short > edge_index_range() const
Definition: elem.h:2692
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&#39;s children...
Definition: fe_base.C:1012
libMesh::DofMap
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
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::Elem::is_node_on_edge
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const =0
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 ...
n_nodes
const dof_id_type n_nodes
Definition: tecplot_io.C:67
libMesh::Variable
This class defines the notion of a variable in the system.
Definition: variable.h:50
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::Elem::node_ref
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2529
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:828
libMesh::Elem::hmin
virtual Real hmin() const
Definition: elem.C:702
libMesh::Elem::n_nodes
virtual unsigned int n_nodes() const =0
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:2919
libMesh::DofMap::n_variables
unsigned int n_variables() const override
Definition: dof_map.h:628
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:993
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::build
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::DofMap::is_constrained_dof
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:2258
libMesh::DofObject::invalid_id
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
libMesh::FEGenericBase::print_dual_dphi
virtual void print_dual_dphi(std::ostream &os) const override
Definition: fe_base.C:903
libMesh::is
PetscErrorCode PetscInt const PetscInt IS * is
Definition: petsc_dm_wrapper.C:62
libMesh::FEType::inf_map
InfMapType inf_map
The coordinate mapping type of the infinite element.
Definition: fe_type.h:275
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:1570
libMesh::PointLocatorBase
This is the base class for point locators.
Definition: point_locator_base.h:54
libMesh::Variable::active_on_subdomain
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:167
libMesh::Elem::n_edges
virtual unsigned int n_edges() const =0
libMesh::FEGenericBase< FEOutputType< T >::type >::OutputNumber
TensorTools::MakeNumber< OutputShape >::type OutputNumber
Definition: fe_base.h:123
libMesh::TypeVector::absolute_fuzzy_equals
bool absolute_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:972
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::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::Elem::max_descendant_p_level
unsigned int max_descendant_p_level() const
Definition: elem.h:3246
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:1878
libMesh::FEType::radial_family
FEFamily radial_family
The type of approximation in radial direction.
Definition: fe_type.h:267
libMesh::OrderWrapper::get_order
int get_order() const
Explicitly request the order as an int.
Definition: fe_type.h:80
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
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:436
libMesh::Elem::n_sides
virtual unsigned int n_sides() const =0
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
libMesh::Elem::level
unsigned int level() const
Definition: elem.h:3074
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::Elem::subdomain_id
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
libMesh::Elem::dim
virtual unsigned short dim() const =0
libMesh::FEContinuity
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
Definition: enum_fe_family.h:84
libMesh::Elem::node_ptr
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
libMesh::Elem::is_vertex
virtual bool is_vertex(const unsigned int i) const =0
libMesh::Elem::n_neighbors
unsigned int n_neighbors() const
Definition: elem.h:714
libMesh::PeriodicBoundaryBase::is_my_variable
bool is_my_variable(unsigned int var_num) const
Definition: periodic_boundary_base.C:71
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::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::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::DenseVector::append
void append(const DenseVector< T2 > &other_vector)
Append additional entries to (resizing, but unchanging) the vector.
Definition: dense_vector.h:408
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:140
libMesh::FEGenericBase::print_dual_phi
virtual void print_dual_phi(std::ostream &os) const override
Definition: fe_base.C:884
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::QGauss
This class implements specific orders of Gauss quadrature.
Definition: quadrature_gauss.h:39
libMesh::Elem::node_index_range
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2683
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::PeriodicBoundaryBase
The base class for defining periodic boundaries.
Definition: periodic_boundary_base.h:48
libMesh::DenseVector
Defines a dense vector for use in Finite Element-type computations.
Definition: dof_map.h:74
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:841
libMesh::Elem::infinite
virtual bool infinite() const =0
libMesh::FEGenericBase::print_d2phi
virtual void print_d2phi(std::ostream &os) const override
Prints the value of each shape function&#39;s second derivatives at each quadrature point.
Definition: fe_base.C:989
libMesh::Elem::active
bool active() const
Definition: elem.h:2941
libMesh::FEGenericBase::print_dual_d2phi
virtual void print_dual_d2phi(std::ostream &os) const override
Definition: fe_base.C:997
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:597
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::DofConstraints
The constraint matrix storage format.
Definition: dof_map.h:108
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:117
libMesh::Elem::is_child_on_edge
virtual bool is_child_on_edge(const unsigned int c, const unsigned int e) const
Definition: elem.C:2306
libMesh::Elem::is_edge
virtual bool is_edge(const unsigned int i) const =0
libMesh::FEGenericBase
This class forms the foundation from which generic finite elements may be derived.
Definition: exact_error_estimator.h:39
libMesh::TensorTools::inner_product
boostcopy::enable_if_c< ScalarTraits< T >::value &&ScalarTraits< T2 >::value, typename CompareTypes< T, T2 >::supertype >::type inner_product(const T &a, const T2 &b)
Definition: tensor_tools.h:51
libMesh::Elem::child_ptr
const Elem * child_ptr(unsigned int i) const
Definition: elem.h:3163
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:5352
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), for one variable vn.
Definition: dof_map.C:2418
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::Variable::type
const FEType & type() const
Definition: variable.h:144
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:30
Generated on Thu Jul 17 2025 01:29:05 for libMesh by   doxygen 1.8.14