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