Line data Source code
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 57362 : 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 9289 : const Elem * primary = elem;
64 :
65 : // Container to catch boundary IDs passed back by BoundaryInfo.
66 18578 : std::vector<boundary_id_type> bc_ids;
67 :
68 : // Pull object out of the loop to reduce heap operations
69 57362 : std::unique_ptr<const Elem> periodic_side;
70 :
71 9289 : std::set<const Elem *> point_neighbors;
72 57362 : elem->find_point_neighbors(p, point_neighbors);
73 218932 : 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 320590 : if ((pt_neighbor->level() > primary->level()) ||
80 275873 : (pt_neighbor->level() == primary->level() &&
81 157886 : pt_neighbor->id() < primary->id()))
82 62428 : 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 14067 : bool vertex_on_periodic_side = false;
88 245325 : for (auto ns : pt_neighbor->side_index_range())
89 : {
90 230634 : boundary_info.boundary_ids (pt_neighbor, ns, bc_ids);
91 :
92 33585 : bool on_relevant_boundary = false;
93 461268 : for (const auto & id : boundary_ids)
94 230634 : if (std::find(bc_ids.begin(), bc_ids.end(), id) != bc_ids.end())
95 14022 : on_relevant_boundary = true;
96 :
97 230634 : if (!on_relevant_boundary)
98 112315 : continue;
99 :
100 98740 : pt_neighbor->build_side_ptr(periodic_side, ns);
101 98740 : if (!periodic_side->contains_point(p))
102 208 : continue;
103 :
104 14004 : vertex_on_periodic_side = true;
105 14004 : break;
106 : }
107 :
108 99142 : if (vertex_on_periodic_side)
109 98519 : primary = pt_neighbor;
110 : }
111 :
112 66651 : return primary;
113 38784 : }
114 :
115 : // Find the "primary" element around a boundary edge:
116 480 : 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 40 : const Elem * primary = elem;
125 :
126 80 : std::set<const Elem *> edge_neighbors;
127 480 : elem->find_edge_neighbors(p1, p2, edge_neighbors);
128 :
129 : // Container to catch boundary IDs handed back by BoundaryInfo
130 80 : std::vector<boundary_id_type> bc_ids;
131 :
132 : // Pull object out of the loop to reduce heap operations
133 440 : std::unique_ptr<const Elem> periodic_side;
134 :
135 960 : 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 960 : if ((e_neighbor->level() > primary->level()) ||
142 880 : (e_neighbor->level() == primary->level() &&
143 480 : e_neighbor->id() < primary->id()))
144 0 : 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 40 : bool vertex_on_periodic_side = false;
150 808 : for (auto ns : e_neighbor->side_index_range())
151 : {
152 768 : boundary_info.boundary_ids (e_neighbor, ns, bc_ids);
153 :
154 64 : bool on_relevant_boundary = false;
155 1536 : for (const auto & id : boundary_ids)
156 768 : if (std::find(bc_ids.begin(), bc_ids.end(), id) != bc_ids.end())
157 58 : on_relevant_boundary = true;
158 :
159 768 : if (!on_relevant_boundary)
160 66 : continue;
161 :
162 696 : e_neighbor->build_side_ptr(periodic_side, ns);
163 1284 : if (!(periodic_side->contains_point(p1) &&
164 588 : periodic_side->contains_point(p2)))
165 216 : continue;
166 :
167 40 : vertex_on_periodic_side = true;
168 40 : break;
169 : }
170 :
171 480 : if (vertex_on_periodic_side)
172 480 : primary = e_neighbor;
173 : }
174 :
175 520 : return primary;
176 400 : }
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 5266814 : FEGenericBase<Real>::build (const unsigned int dim,
192 : const FEType & fet)
193 : {
194 5266814 : switch (dim)
195 : {
196 : // 0D
197 12 : case 0:
198 : {
199 12 : switch (fet.family)
200 : {
201 0 : case CLOUGH:
202 0 : return std::make_unique<FE<0,CLOUGH>>(fet);
203 :
204 0 : case HERMITE:
205 0 : return std::make_unique<FE<0,HERMITE>>(fet);
206 :
207 12 : case LAGRANGE:
208 12 : return std::make_unique<FE<0,LAGRANGE>>(fet);
209 :
210 0 : case L2_LAGRANGE:
211 0 : return std::make_unique<FE<0,L2_LAGRANGE>>(fet);
212 :
213 0 : case HIERARCHIC:
214 0 : return std::make_unique<FE<0,HIERARCHIC>>(fet);
215 :
216 0 : case L2_HIERARCHIC:
217 0 : return std::make_unique<FE<0,L2_HIERARCHIC>>(fet);
218 :
219 0 : case SIDE_HIERARCHIC:
220 0 : return std::make_unique<FE<0,SIDE_HIERARCHIC>>(fet);
221 :
222 0 : case MONOMIAL:
223 0 : return std::make_unique<FE<0,MONOMIAL>>(fet);
224 :
225 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
226 0 : case SZABAB:
227 0 : return std::make_unique<FE<0,SZABAB>>(fet);
228 :
229 0 : case BERNSTEIN:
230 0 : return std::make_unique<FE<0,BERNSTEIN>>(fet);
231 :
232 0 : case RATIONAL_BERNSTEIN:
233 0 : return std::make_unique<FE<0,RATIONAL_BERNSTEIN>>(fet);
234 : #endif
235 :
236 0 : case XYZ:
237 0 : return std::make_unique<FEXYZ<0>>(fet);
238 :
239 0 : case SCALAR:
240 0 : return std::make_unique<FEScalar<0>>(fet);
241 :
242 0 : default:
243 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
244 : }
245 : }
246 : // 1D
247 436876 : case 1:
248 : {
249 436876 : switch (fet.family)
250 : {
251 0 : case CLOUGH:
252 0 : return std::make_unique<FE<1,CLOUGH>>(fet);
253 :
254 398834 : case HERMITE:
255 398834 : return std::make_unique<FE<1,HERMITE>>(fet);
256 :
257 3563 : case LAGRANGE:
258 3563 : return std::make_unique<FE<1,LAGRANGE>>(fet);
259 :
260 2130 : case L2_LAGRANGE:
261 2130 : return std::make_unique<FE<1,L2_LAGRANGE>>(fet);
262 :
263 11337 : case HIERARCHIC:
264 11337 : return std::make_unique<FE<1,HIERARCHIC>>(fet);
265 :
266 3550 : case L2_HIERARCHIC:
267 3550 : return std::make_unique<FE<1,L2_HIERARCHIC>>(fet);
268 :
269 3550 : case SIDE_HIERARCHIC:
270 3550 : return std::make_unique<FE<1,SIDE_HIERARCHIC>>(fet);
271 :
272 3592 : case MONOMIAL:
273 3592 : return std::make_unique<FE<1,MONOMIAL>>(fet);
274 :
275 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
276 2840 : case SZABAB:
277 2840 : return std::make_unique<FE<1,SZABAB>>(fet);
278 :
279 2840 : case BERNSTEIN:
280 2840 : return std::make_unique<FE<1,BERNSTEIN>>(fet);
281 :
282 1800 : case RATIONAL_BERNSTEIN:
283 1800 : return std::make_unique<FE<1,RATIONAL_BERNSTEIN>>(fet);
284 : #endif
285 :
286 2840 : case XYZ:
287 2840 : return std::make_unique<FEXYZ<1>>(fet);
288 :
289 0 : case SCALAR:
290 0 : return std::make_unique<FEScalar<1>>(fet);
291 :
292 0 : default:
293 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
294 : }
295 : }
296 :
297 :
298 : // 2D
299 1678524 : case 2:
300 : {
301 1678524 : switch (fet.family)
302 : {
303 69780 : case CLOUGH:
304 69780 : return std::make_unique<FE<2,CLOUGH>>(fet);
305 :
306 35446 : case HERMITE:
307 35446 : return std::make_unique<FE<2,HERMITE>>(fet);
308 :
309 539227 : case LAGRANGE:
310 539227 : return std::make_unique<FE<2,LAGRANGE>>(fet);
311 :
312 53964 : case L2_LAGRANGE:
313 53964 : return std::make_unique<FE<2,L2_LAGRANGE>>(fet);
314 :
315 538812 : case HIERARCHIC:
316 538812 : return std::make_unique<FE<2,HIERARCHIC>>(fet);
317 :
318 100174 : case L2_HIERARCHIC:
319 100174 : return std::make_unique<FE<2,L2_HIERARCHIC>>(fet);
320 :
321 147652 : case SIDE_HIERARCHIC:
322 147652 : return std::make_unique<FE<2,SIDE_HIERARCHIC>>(fet);
323 :
324 18337 : case MONOMIAL:
325 18337 : return std::make_unique<FE<2,MONOMIAL>>(fet);
326 :
327 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
328 12082 : case SZABAB:
329 12082 : return std::make_unique<FE<2,SZABAB>>(fet);
330 :
331 11698 : case BERNSTEIN:
332 11698 : return std::make_unique<FE<2,BERNSTEIN>>(fet);
333 :
334 117469 : case RATIONAL_BERNSTEIN:
335 117469 : return std::make_unique<FE<2,RATIONAL_BERNSTEIN>>(fet);
336 : #endif
337 :
338 33812 : case XYZ:
339 33812 : return std::make_unique<FEXYZ<2>>(fet);
340 :
341 0 : case SCALAR:
342 0 : return std::make_unique<FEScalar<2>>(fet);
343 :
344 71 : case SUBDIVISION:
345 71 : return std::make_unique<FESubdivision>(fet);
346 :
347 0 : default:
348 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
349 : }
350 : }
351 :
352 :
353 : // 3D
354 3151402 : case 3:
355 : {
356 3151402 : switch (fet.family)
357 : {
358 0 : case CLOUGH:
359 0 : libmesh_error_msg("ERROR: Clough-Tocher elements currently only support 1D and 2D");
360 :
361 926 : case HERMITE:
362 926 : return std::make_unique<FE<3,HERMITE>>(fet);
363 :
364 2462495 : case LAGRANGE:
365 2462495 : return std::make_unique<FE<3,LAGRANGE>>(fet);
366 :
367 85165 : case L2_LAGRANGE:
368 85165 : return std::make_unique<FE<3,L2_LAGRANGE>>(fet);
369 :
370 56554 : case HIERARCHIC:
371 56554 : return std::make_unique<FE<3,HIERARCHIC>>(fet);
372 :
373 89076 : case L2_HIERARCHIC:
374 89076 : return std::make_unique<FE<3,L2_HIERARCHIC>>(fet);
375 :
376 305448 : case SIDE_HIERARCHIC:
377 305448 : return std::make_unique<FE<3,SIDE_HIERARCHIC>>(fet);
378 :
379 25603 : case MONOMIAL:
380 25603 : return std::make_unique<FE<3,MONOMIAL>>(fet);
381 :
382 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
383 0 : case SZABAB:
384 0 : return std::make_unique<FE<3,SZABAB>>(fet);
385 :
386 21386 : case BERNSTEIN:
387 21386 : return std::make_unique<FE<3,BERNSTEIN>>(fet);
388 :
389 12419 : case RATIONAL_BERNSTEIN:
390 12419 : return std::make_unique<FE<3,RATIONAL_BERNSTEIN>>(fet);
391 : #endif
392 :
393 92330 : case XYZ:
394 92330 : return std::make_unique<FEXYZ<3>>(fet);
395 :
396 0 : case SCALAR:
397 0 : return std::make_unique<FEScalar<3>>(fet);
398 :
399 0 : default:
400 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
401 : }
402 : }
403 :
404 0 : default:
405 0 : libmesh_error_msg("Invalid dimension dim = " << dim);
406 : }
407 : }
408 :
409 :
410 :
411 : template <>
412 : std::unique_ptr<FEGenericBase<RealGradient>>
413 574854 : FEGenericBase<RealGradient>::build (const unsigned int dim,
414 : const FEType & fet)
415 : {
416 574854 : switch (dim)
417 : {
418 : // 0D
419 0 : case 0:
420 : {
421 0 : switch (fet.family)
422 : {
423 0 : case HIERARCHIC_VEC:
424 0 : return std::make_unique<FEHierarchicVec<0>>(fet);
425 :
426 0 : case L2_HIERARCHIC_VEC:
427 0 : return std::make_unique<FEL2HierarchicVec<0>>(fet);
428 :
429 0 : case LAGRANGE_VEC:
430 0 : return std::make_unique<FELagrangeVec<0>>(fet);
431 :
432 0 : case L2_LAGRANGE_VEC:
433 0 : return std::make_unique<FEL2LagrangeVec<0>>(fet);
434 :
435 0 : case MONOMIAL_VEC:
436 0 : return std::make_unique<FEMonomialVec<0>>(fet);
437 :
438 0 : default:
439 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
440 : }
441 : }
442 0 : case 1:
443 : {
444 0 : switch (fet.family)
445 : {
446 0 : case HIERARCHIC_VEC:
447 0 : return std::make_unique<FEHierarchicVec<1>>(fet);
448 :
449 0 : case L2_HIERARCHIC_VEC:
450 0 : return std::make_unique<FEL2HierarchicVec<1>>(fet);
451 :
452 0 : case LAGRANGE_VEC:
453 0 : return std::make_unique<FELagrangeVec<1>>(fet);
454 :
455 0 : case L2_LAGRANGE_VEC:
456 0 : return std::make_unique<FEL2LagrangeVec<1>>(fet);
457 :
458 0 : case MONOMIAL_VEC:
459 0 : return std::make_unique<FEMonomialVec<1>>(fet);
460 :
461 0 : default:
462 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
463 : }
464 : }
465 212152 : case 2:
466 : {
467 212152 : switch (fet.family)
468 : {
469 16572 : case HIERARCHIC_VEC:
470 16572 : return std::make_unique<FEHierarchicVec<2>>(fet);
471 :
472 2840 : case L2_HIERARCHIC_VEC:
473 2840 : return std::make_unique<FEL2HierarchicVec<2>>(fet);
474 :
475 16121 : case LAGRANGE_VEC:
476 16121 : return std::make_unique<FELagrangeVec<2>>(fet);
477 :
478 3250 : case L2_LAGRANGE_VEC:
479 3250 : return std::make_unique<FEL2LagrangeVec<2>>(fet);
480 :
481 1470 : case MONOMIAL_VEC:
482 1470 : return std::make_unique<FEMonomialVec<2>>(fet);
483 :
484 55937 : case NEDELEC_ONE:
485 55937 : return std::make_unique<FENedelecOne<2>>(fet);
486 :
487 114542 : case RAVIART_THOMAS:
488 114542 : return std::make_unique<FERaviartThomas<2>>(fet);
489 :
490 1420 : case L2_RAVIART_THOMAS:
491 1420 : return std::make_unique<FEL2RaviartThomas<2>>(fet);
492 :
493 0 : default:
494 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
495 : }
496 : }
497 362702 : case 3:
498 : {
499 362702 : switch (fet.family)
500 : {
501 0 : case HIERARCHIC_VEC:
502 0 : return std::make_unique<FEHierarchicVec<3>>(fet);
503 :
504 710 : case L2_HIERARCHIC_VEC:
505 710 : return std::make_unique<FEL2HierarchicVec<3>>(fet);
506 :
507 23093 : case LAGRANGE_VEC:
508 23093 : return std::make_unique<FELagrangeVec<3>>(fet);
509 :
510 1420 : case L2_LAGRANGE_VEC:
511 1420 : return std::make_unique<FEL2LagrangeVec<3>>(fet);
512 :
513 0 : case MONOMIAL_VEC:
514 0 : return std::make_unique<FEMonomialVec<3>>(fet);
515 :
516 33874 : case NEDELEC_ONE:
517 33874 : return std::make_unique<FENedelecOne<3>>(fet);
518 :
519 302895 : case RAVIART_THOMAS:
520 302895 : return std::make_unique<FERaviartThomas<3>>(fet);
521 :
522 710 : case L2_RAVIART_THOMAS:
523 710 : return std::make_unique<FEL2RaviartThomas<3>>(fet);
524 :
525 0 : default:
526 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
527 : }
528 : }
529 :
530 0 : default:
531 0 : 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 474 : FEGenericBase<Real>::build_InfFE (const unsigned int dim,
547 : const FEType & fet)
548 : {
549 474 : switch (dim)
550 : {
551 :
552 : // 1D
553 0 : case 1:
554 : {
555 0 : switch (fet.radial_family)
556 : {
557 0 : case INFINITE_MAP:
558 0 : libmesh_error_msg("ERROR: Can't build an infinite element with FEFamily = " << Utility::enum_to_string(fet.radial_family));
559 :
560 0 : case JACOBI_20_00:
561 : {
562 0 : switch (fet.inf_map)
563 : {
564 0 : case CARTESIAN:
565 0 : return std::make_unique<InfFE<1,JACOBI_20_00,CARTESIAN>>(fet);
566 :
567 0 : default:
568 0 : libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
569 : }
570 : }
571 :
572 0 : case JACOBI_30_00:
573 : {
574 0 : switch (fet.inf_map)
575 : {
576 0 : case CARTESIAN:
577 0 : return std::make_unique<InfFE<1,JACOBI_30_00,CARTESIAN>>(fet);
578 :
579 0 : default:
580 0 : libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
581 : }
582 : }
583 :
584 0 : case LEGENDRE:
585 : {
586 0 : switch (fet.inf_map)
587 : {
588 0 : case CARTESIAN:
589 0 : return std::make_unique<InfFE<1,LEGENDRE,CARTESIAN>>(fet);
590 :
591 0 : default:
592 0 : libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
593 : }
594 : }
595 :
596 0 : case LAGRANGE:
597 : {
598 0 : switch (fet.inf_map)
599 : {
600 0 : case CARTESIAN:
601 0 : return std::make_unique<InfFE<1,LAGRANGE,CARTESIAN>>(fet);
602 :
603 0 : default:
604 0 : libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
605 : }
606 : }
607 :
608 0 : default:
609 0 : libmesh_error_msg("ERROR: Bad FEType.radial_family= " << Utility::enum_to_string(fet.radial_family));
610 : }
611 : }
612 :
613 :
614 :
615 :
616 : // 2D
617 0 : case 2:
618 : {
619 0 : switch (fet.radial_family)
620 : {
621 0 : case INFINITE_MAP:
622 0 : libmesh_error_msg("ERROR: Can't build an infinite element with FEFamily = " << Utility::enum_to_string(fet.radial_family));
623 :
624 0 : case JACOBI_20_00:
625 : {
626 0 : switch (fet.inf_map)
627 : {
628 0 : case CARTESIAN:
629 0 : return std::make_unique<InfFE<2,JACOBI_20_00,CARTESIAN>>(fet);
630 :
631 0 : default:
632 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
633 : }
634 : }
635 :
636 0 : case JACOBI_30_00:
637 : {
638 0 : switch (fet.inf_map)
639 : {
640 0 : case CARTESIAN:
641 0 : return std::make_unique<InfFE<2,JACOBI_30_00,CARTESIAN>>(fet);
642 :
643 0 : default:
644 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
645 : }
646 : }
647 :
648 0 : case LEGENDRE:
649 : {
650 0 : switch (fet.inf_map)
651 : {
652 0 : case CARTESIAN:
653 0 : return std::make_unique<InfFE<2,LEGENDRE,CARTESIAN>>(fet);
654 :
655 0 : default:
656 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
657 : }
658 : }
659 :
660 0 : case LAGRANGE:
661 : {
662 0 : switch (fet.inf_map)
663 : {
664 0 : case CARTESIAN:
665 0 : return std::make_unique<InfFE<2,LAGRANGE,CARTESIAN>>(fet);
666 :
667 0 : default:
668 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
669 : }
670 : }
671 :
672 0 : default:
673 0 : libmesh_error_msg("ERROR: Bad FEType.radial_family= " << Utility::enum_to_string(fet.radial_family));
674 : }
675 : }
676 :
677 :
678 :
679 :
680 : // 3D
681 474 : case 3:
682 : {
683 474 : switch (fet.radial_family)
684 : {
685 0 : case INFINITE_MAP:
686 0 : libmesh_error_msg("ERROR: Don't build an infinite element with FEFamily = " << Utility::enum_to_string(fet.radial_family));
687 :
688 419 : case JACOBI_20_00:
689 : {
690 419 : switch (fet.inf_map)
691 : {
692 419 : case CARTESIAN:
693 419 : return std::make_unique<InfFE<3,JACOBI_20_00,CARTESIAN>>(fet);
694 :
695 0 : default:
696 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
697 : }
698 : }
699 :
700 20 : case JACOBI_30_00:
701 : {
702 20 : switch (fet.inf_map)
703 : {
704 20 : case CARTESIAN:
705 20 : return std::make_unique<InfFE<3,JACOBI_30_00,CARTESIAN>>(fet);
706 :
707 0 : default:
708 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
709 : }
710 : }
711 :
712 20 : case LEGENDRE:
713 : {
714 20 : switch (fet.inf_map)
715 : {
716 20 : case CARTESIAN:
717 20 : return std::make_unique<InfFE<3,LEGENDRE,CARTESIAN>>(fet);
718 :
719 0 : default:
720 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
721 : }
722 : }
723 :
724 15 : case LAGRANGE:
725 : {
726 15 : switch (fet.inf_map)
727 : {
728 15 : case CARTESIAN:
729 15 : return std::make_unique<InfFE<3,LAGRANGE,CARTESIAN>>(fet);
730 :
731 0 : default:
732 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
733 : }
734 : }
735 :
736 0 : default:
737 0 : libmesh_error_msg("ERROR: Bad FEType.radial_family= " << Utility::enum_to_string(fet.radial_family));
738 : }
739 : }
740 :
741 0 : default:
742 0 : libmesh_error_msg("Invalid dimension dim = " << dim);
743 : }
744 : }
745 :
746 :
747 :
748 : template <>
749 : std::unique_ptr<FEGenericBase<RealGradient>>
750 0 : FEGenericBase<RealGradient>::build_InfFE (const unsigned int,
751 : const FEType & )
752 : {
753 : // No vector types defined... YET.
754 0 : 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 122869529 : 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 22032914 : LOG_SCOPE("compute_shape_functions()", "FE");
772 :
773 122869529 : this->determine_calculations();
774 :
775 122869529 : if (calculate_phi)
776 111419034 : this->_fe_trans->map_phi(this->dim, elem, qp, (*this), this->phi, this->_add_p_level_in_reinit);
777 :
778 122869529 : if (calculate_dphi)
779 88615807 : this->_fe_trans->map_dphi(this->dim, elem, qp, (*this), this->dphi,
780 81153692 : this->dphidx, this->dphidy, this->dphidz);
781 :
782 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
783 122869529 : if (calculate_d2phi)
784 7132611 : this->_fe_trans->map_d2phi(this->dim, qp, (*this), this->d2phi,
785 6608208 : this->d2phidx2, this->d2phidxdy, this->d2phidxdz,
786 6608208 : this->d2phidy2, this->d2phidydz, this->d2phidz2);
787 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
788 :
789 : // Only compute curl for vector-valued elements
790 9947971 : if (calculate_curl_phi && TypesEqual<OutputType,RealGradient>::value)
791 859569 : this->_fe_trans->map_curl(this->dim, elem, qp, (*this), this->curl_phi);
792 :
793 : // Only compute div for vector-valued elements
794 9947971 : if (calculate_div_phi && TypesEqual<OutputType,RealGradient>::value)
795 1476125 : this->_fe_trans->map_div(this->dim, elem, qp, (*this), this->div_phi);
796 122869529 : }
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 6384 : 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 856 : LOG_SCOPE("compute_dual_shape_coeffs()", "FE");
808 :
809 6384 : const unsigned int sz=phi_vals.size();
810 6384 : libmesh_error_msg_if(!sz, "ERROR: cannot compute dual shape coefficients with empty phi values");
811 :
812 : //compute dual basis coefficient (dual_coeff)
813 5956 : dual_coeff.resize(sz, sz);
814 7240 : DenseMatrix<Real> A(sz, sz), D(sz, sz);
815 :
816 108314 : for (const auto i : index_range(phi_vals))
817 6044160 : for (const auto qp : index_range(phi_vals[i]))
818 : {
819 7198867 : D(i,i) += JxW[qp]*phi_vals[i][qp];
820 306293128 : for (const auto j : index_range(phi_vals))
821 368484603 : A(i,j) += JxW[qp]*phi_vals[i][qp]*phi_vals[j][qp];
822 : }
823 :
824 : // dual_coeff = A^-1*D
825 108314 : for (const auto j : index_range(phi_vals))
826 : {
827 101930 : DenseVector<Real> Dcol(sz), coeffcol(sz);
828 3505680 : for (const auto i : index_range(phi_vals))
829 3902270 : Dcol(i) = D(i, j);
830 101930 : A.cholesky_solve(Dcol, coeffcol);
831 :
832 3505680 : for (const auto row : index_range(phi_vals))
833 3902270 : dual_coeff(row, j)=coeffcol(row);
834 : }
835 6384 : }
836 :
837 : template <>
838 6384 : void FEGenericBase<Real>::compute_dual_shape_functions ()
839 : {
840 : // Start logging the shape function computation
841 856 : LOG_SCOPE("compute_dual_shape_functions()", "FE");
842 :
843 : // The dual coeffs matrix should have the same size as phi
844 428 : libmesh_assert(dual_coeff.m() == phi.size());
845 428 : libmesh_assert(dual_coeff.n() == phi.size());
846 :
847 : // initialize dual basis
848 108314 : for (const auto j : index_range(phi))
849 6044160 : for (const auto qp : index_range(phi[j]))
850 : {
851 6361109 : dual_phi[j][qp] = 0;
852 5942230 : if (calculate_dphi)
853 1256625 : dual_dphi[j][qp] = 0;
854 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
855 5942230 : if (calculate_d2phi)
856 1248093 : dual_d2phi[j][qp] = 0;
857 : #endif
858 : }
859 :
860 : // compute dual basis
861 108314 : for (const auto j : index_range(phi))
862 3505680 : for (const auto i : index_range(phi))
863 303754648 : for (const auto qp : index_range(phi[j]))
864 : {
865 391195838 : dual_phi[j][qp] += dual_coeff(i, j) * phi[i][qp];
866 300350898 : if (calculate_dphi)
867 136267362 : dual_dphi[j][qp] += dual_coeff(i, j) * dphi[i][qp];
868 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
869 300350898 : if (calculate_d2phi)
870 667700859 : dual_d2phi[j][qp] += dual_coeff(i, j) * d2phi[i][qp];
871 : #endif
872 : }
873 6384 : }
874 :
875 : template <typename OutputType>
876 0 : void FEGenericBase<OutputType>::print_phi(std::ostream & os) const
877 : {
878 0 : for (auto i : index_range(phi))
879 0 : for (auto j : index_range(phi[i]))
880 0 : os << " phi[" << i << "][" << j << "]=" << phi[i][j] << std::endl;
881 0 : }
882 :
883 : template <typename OutputType>
884 0 : void FEGenericBase<OutputType>::print_dual_phi(std::ostream & os) const
885 : {
886 0 : for (auto i : index_range(dual_phi))
887 0 : for (auto j : index_range(dual_phi[i]))
888 0 : os << " dual_phi[" << i << "][" << j << "]=" << dual_phi[i][j] << std::endl;
889 0 : }
890 :
891 :
892 :
893 :
894 : template <typename OutputType>
895 0 : void FEGenericBase<OutputType>::print_dphi(std::ostream & os) const
896 : {
897 0 : for (auto i : index_range(dphi))
898 0 : for (auto j : index_range(dphi[i]))
899 0 : os << " dphi[" << i << "][" << j << "]=" << dphi[i][j];
900 0 : }
901 :
902 : template <typename OutputType>
903 0 : void FEGenericBase<OutputType>::print_dual_dphi(std::ostream & os) const
904 : {
905 0 : for (auto i : index_range(dphi))
906 0 : for (auto j : index_range(dphi[i]))
907 0 : os << " dual_dphi[" << i << "][" << j << "]=" << dual_dphi[i][j];
908 0 : }
909 :
910 :
911 :
912 : template <typename OutputType>
913 344118291 : void FEGenericBase<OutputType>::determine_calculations()
914 : {
915 344118291 : 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 29891376 : bool requested_ok =
920 334783029 : this->calculate_nothing || this->calculate_phi || this->calculate_dphi ||
921 375755589 : this->calculate_dphiref || this->calculate_curl_phi || this->calculate_div_phi ||
922 366149 : this->calculate_map;
923 :
924 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
925 29891376 : requested_ok = requested_ok || this->calculate_d2phi;
926 : #endif
927 :
928 29891376 : 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 344118291 : if (this->calculate_phi)
936 320375498 : this->_fe_trans->init_map_phi(*this);
937 :
938 344118291 : if (this->calculate_dphiref)
939 225612355 : this->_fe_trans->init_map_dphi(*this);
940 :
941 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
942 344118291 : if (this->calculate_d2phi)
943 58343657 : this->_fe_trans->init_map_d2phi(*this);
944 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
945 344118291 : }
946 :
947 :
948 :
949 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
950 :
951 :
952 : template <typename OutputType>
953 0 : void FEGenericBase<OutputType>::print_d2phi(std::ostream & os) const
954 : {
955 0 : for (auto i : index_range(dphi))
956 0 : for (auto j : index_range(dphi[i]))
957 0 : os << " d2phi[" << i << "][" << j << "]=" << d2phi[i][j];
958 0 : }
959 :
960 : template <typename OutputType>
961 0 : void FEGenericBase<OutputType>::print_dual_d2phi(std::ostream & os) const
962 : {
963 0 : for (auto i : index_range(dual_d2phi))
964 0 : for (auto j : index_range(dual_d2phi[i]))
965 0 : os << " dual_d2phi[" << i << "][" << j << "]=" << dual_d2phi[i][j];
966 0 : }
967 :
968 : #endif
969 :
970 :
971 :
972 : #ifdef LIBMESH_ENABLE_AMR
973 :
974 : template <typename OutputType>
975 : void
976 0 : 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 0 : std::vector<unsigned int> new_side_dofs, old_side_dofs;
985 :
986 : // FIXME: what about 2D shells in 3D space?
987 0 : unsigned int dim = elem->dim();
988 :
989 : // Cache n_children(); it's a virtual call but it's const.
990 0 : 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 0 : const FEType & base_fe_type = dof_map.variable_type(var);
995 0 : std::unique_ptr<FEGenericBase<OutputShape>> fe
996 : (FEGenericBase<OutputShape>::build(dim, base_fe_type));
997 0 : std::unique_ptr<FEGenericBase<OutputShape>> fe_coarse
998 : (FEGenericBase<OutputShape>::build(dim, base_fe_type));
999 :
1000 0 : std::unique_ptr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
1001 0 : std::unique_ptr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
1002 0 : std::unique_ptr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
1003 0 : std::vector<Point> coarse_qpoints;
1004 :
1005 : // The values of the shape functions at the quadrature
1006 : // points
1007 0 : const std::vector<std::vector<OutputShape>> & phi_values =
1008 : fe->get_phi();
1009 0 : 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 0 : const std::vector<std::vector<OutputGradient>> * dphi_values =
1015 : nullptr;
1016 0 : const std::vector<std::vector<OutputGradient>> * dphi_coarse =
1017 : nullptr;
1018 :
1019 0 : const FEContinuity cont = fe->get_continuity();
1020 :
1021 0 : if (cont == C_ONE)
1022 : {
1023 : const std::vector<std::vector<OutputGradient>> &
1024 0 : ref_dphi_values = fe->get_dphi();
1025 0 : dphi_values = &ref_dphi_values;
1026 : const std::vector<std::vector<OutputGradient>> &
1027 0 : ref_dphi_coarse = fe_coarse->get_dphi();
1028 0 : dphi_coarse = &ref_dphi_coarse;
1029 : }
1030 :
1031 : // The Jacobian * quadrature weight at the quadrature points
1032 0 : const std::vector<Real> & JxW =
1033 0 : fe->get_JxW();
1034 :
1035 : // The XYZ locations of the quadrature points on the
1036 : // child element
1037 0 : const std::vector<Point> & xyz_values =
1038 0 : fe->get_xyz();
1039 :
1040 : // Number of nodes on parent element
1041 0 : const unsigned int n_nodes = elem->n_nodes();
1042 :
1043 : // Number of dofs on parent element
1044 : const unsigned int new_n_dofs =
1045 0 : FEInterface::n_dofs(base_fe_type, elem->max_descendant_p_level(), elem);
1046 :
1047 : // Fixed vs. free DoFs on edge/face projections
1048 0 : std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
1049 0 : std::vector<int> free_dof(new_n_dofs, 0);
1050 :
1051 0 : DenseMatrix<Real> Ke;
1052 0 : DenseVector<Number> Fe;
1053 0 : 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 0 : std::vector<dof_id_type> node_dof_indices;
1066 0 : if (use_old_dof_indices)
1067 0 : dof_map.old_dof_indices (elem, node_dof_indices, var);
1068 : else
1069 0 : dof_map.dof_indices (elem, node_dof_indices, var);
1070 :
1071 0 : unsigned int current_dof = 0;
1072 0 : 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 0 : FEInterface::n_dofs_at_node (base_fe_type, elem->max_descendant_p_level(), elem, n);
1078 0 : if (!elem->is_vertex(n))
1079 : {
1080 0 : current_dof += my_nc;
1081 0 : 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 0 : 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 0 : FEInterface::n_dofs_at_node (base_fe_type, extra_order, elem, n);
1094 0 : for (unsigned int i=0; i!= nc; ++i)
1095 : {
1096 0 : Ue(current_dof) =
1097 0 : old_vector(node_dof_indices[current_dof]);
1098 0 : dof_is_fixed[current_dof] = true;
1099 0 : current_dof++;
1100 : }
1101 : }
1102 : }
1103 :
1104 0 : FEType fe_type = base_fe_type, temp_fe_type;
1105 0 : fe_type.order = fe_type.order + elem->max_descendant_p_level();
1106 :
1107 : // In 3D, project any edge values next
1108 0 : if (dim > 2 && cont != DISCONTINUOUS)
1109 0 : for (auto e : elem->edge_index_range())
1110 : {
1111 0 : FEInterface::dofs_on_edge(elem, dim, fe_type,
1112 : e, new_side_dofs);
1113 :
1114 : const unsigned int n_new_side_dofs =
1115 0 : 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 0 : unsigned int free_dofs = 0;
1120 0 : for (unsigned int i=0; i != n_new_side_dofs; ++i)
1121 0 : if (!dof_is_fixed[new_side_dofs[i]])
1122 0 : free_dof[free_dofs++] = i;
1123 0 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1124 0 : Fe.resize (free_dofs); Fe.zero();
1125 : // The new edge coefficients
1126 0 : DenseVector<Number> Uedge(free_dofs);
1127 :
1128 : // Add projection terms from each child sharing
1129 : // this edge
1130 0 : for (unsigned int c=0; c != n_children; ++c)
1131 : {
1132 0 : if (!elem->is_child_on_edge(c,e))
1133 0 : continue;
1134 0 : const Elem * child = elem->child_ptr(c);
1135 :
1136 0 : std::vector<dof_id_type> child_dof_indices;
1137 0 : if (use_old_dof_indices)
1138 0 : dof_map.old_dof_indices (child,
1139 : child_dof_indices, var);
1140 : else
1141 0 : dof_map.dof_indices (child,
1142 : child_dof_indices, var);
1143 : const unsigned int child_n_dofs =
1144 : cast_int<unsigned int>
1145 0 : (child_dof_indices.size());
1146 :
1147 0 : temp_fe_type = base_fe_type;
1148 0 : temp_fe_type.order = temp_fe_type.order + child->p_level();
1149 :
1150 0 : 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 0 : fe->attach_quadrature_rule (qedgerule.get());
1156 0 : fe->edge_reinit (child, e);
1157 0 : const unsigned int n_qp = qedgerule->n_points();
1158 :
1159 0 : FEMap::inverse_map (dim, elem, xyz_values,
1160 : coarse_qpoints);
1161 :
1162 0 : fe_coarse->reinit(elem, &coarse_qpoints);
1163 :
1164 : // Loop over the quadrature points
1165 0 : for (unsigned int qp=0; qp<n_qp; qp++)
1166 : {
1167 : // solution value at the quadrature point
1168 0 : OutputNumber fineval = libMesh::zero;
1169 : // solution grad at the quadrature point
1170 0 : 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 0 : for (unsigned int i=0; i<child_n_dofs;
1176 : i++)
1177 : {
1178 0 : fineval +=
1179 0 : (old_vector(child_dof_indices[i])*
1180 0 : phi_values[i][qp]);
1181 0 : if (cont == C_ONE)
1182 0 : finegrad += (*dphi_values)[i][qp] *
1183 0 : old_vector(child_dof_indices[i]);
1184 : }
1185 :
1186 : // Form edge projection matrix
1187 0 : for (unsigned int sidei=0, freei=0; sidei != n_new_side_dofs; ++sidei)
1188 : {
1189 0 : unsigned int i = new_side_dofs[sidei];
1190 : // fixed DoFs aren't test functions
1191 0 : if (dof_is_fixed[i])
1192 0 : continue;
1193 0 : for (unsigned int sidej=0, freej=0; sidej != n_new_side_dofs; ++sidej)
1194 : {
1195 0 : unsigned int j =
1196 : new_side_dofs[sidej];
1197 0 : if (dof_is_fixed[j])
1198 0 : Fe(freei) -=
1199 0 : TensorTools::inner_product(phi_coarse[i][qp],
1200 0 : phi_coarse[j][qp]) *
1201 0 : JxW[qp] * Ue(j);
1202 : else
1203 0 : Ke(freei,freej) +=
1204 0 : TensorTools::inner_product(phi_coarse[i][qp],
1205 0 : phi_coarse[j][qp]) *
1206 : JxW[qp];
1207 0 : if (cont == C_ONE)
1208 : {
1209 0 : if (dof_is_fixed[j])
1210 0 : Fe(freei) -=
1211 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1212 0 : (*dphi_coarse)[j][qp]) *
1213 0 : JxW[qp] * Ue(j);
1214 : else
1215 0 : Ke(freei,freej) +=
1216 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1217 0 : (*dphi_coarse)[j][qp]) *
1218 : JxW[qp];
1219 : }
1220 0 : if (!dof_is_fixed[j])
1221 0 : freej++;
1222 : }
1223 0 : Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp],
1224 0 : fineval) * JxW[qp];
1225 0 : if (cont == C_ONE)
1226 0 : Fe(freei) +=
1227 0 : TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1228 0 : freei++;
1229 : }
1230 : }
1231 : }
1232 0 : Ke.cholesky_solve(Fe, Uedge);
1233 :
1234 : // Transfer new edge solutions to element
1235 0 : for (unsigned int i=0; i != free_dofs; ++i)
1236 : {
1237 0 : Number & ui = Ue(new_side_dofs[free_dof[i]]);
1238 0 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1239 : std::abs(ui - Uedge(i)) < TOLERANCE);
1240 0 : ui = Uedge(i);
1241 0 : 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 0 : if (dim > 1 && cont != DISCONTINUOUS)
1247 0 : for (auto s : elem->side_index_range())
1248 : {
1249 0 : FEInterface::dofs_on_side(elem, dim, fe_type,
1250 : s, new_side_dofs);
1251 :
1252 : const unsigned int n_new_side_dofs =
1253 0 : 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 0 : unsigned int free_dofs = 0;
1258 0 : for (unsigned int i=0; i != n_new_side_dofs; ++i)
1259 0 : if (!dof_is_fixed[new_side_dofs[i]])
1260 0 : free_dof[free_dofs++] = i;
1261 0 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1262 0 : Fe.resize (free_dofs); Fe.zero();
1263 : // The new side coefficients
1264 0 : DenseVector<Number> Uside(free_dofs);
1265 :
1266 : // Add projection terms from each child sharing
1267 : // this side
1268 0 : for (unsigned int c=0; c != n_children; ++c)
1269 : {
1270 0 : if (!elem->is_child_on_side(c,s))
1271 0 : continue;
1272 0 : const Elem * child = elem->child_ptr(c);
1273 :
1274 0 : std::vector<dof_id_type> child_dof_indices;
1275 0 : if (use_old_dof_indices)
1276 0 : dof_map.old_dof_indices (child,
1277 : child_dof_indices, var);
1278 : else
1279 0 : dof_map.dof_indices (child,
1280 : child_dof_indices, var);
1281 : const unsigned int child_n_dofs =
1282 : cast_int<unsigned int>
1283 0 : (child_dof_indices.size());
1284 :
1285 0 : temp_fe_type = base_fe_type;
1286 0 : temp_fe_type.order = temp_fe_type.order + child->p_level();
1287 :
1288 0 : 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 0 : fe->attach_quadrature_rule (qsiderule.get());
1294 0 : fe->reinit (child, s);
1295 0 : const unsigned int n_qp = qsiderule->n_points();
1296 :
1297 0 : FEMap::inverse_map (dim, elem, xyz_values,
1298 : coarse_qpoints);
1299 :
1300 0 : fe_coarse->reinit(elem, &coarse_qpoints);
1301 :
1302 : // Loop over the quadrature points
1303 0 : for (unsigned int qp=0; qp<n_qp; qp++)
1304 : {
1305 : // solution value at the quadrature point
1306 0 : OutputNumber fineval = libMesh::zero;
1307 : // solution grad at the quadrature point
1308 0 : 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 0 : for (unsigned int i=0; i<child_n_dofs;
1314 : i++)
1315 : {
1316 0 : fineval +=
1317 0 : old_vector(child_dof_indices[i]) *
1318 0 : phi_values[i][qp];
1319 0 : if (cont == C_ONE)
1320 0 : finegrad += (*dphi_values)[i][qp] *
1321 0 : old_vector(child_dof_indices[i]);
1322 : }
1323 :
1324 : // Form side projection matrix
1325 0 : for (unsigned int sidei=0, freei=0; sidei != n_new_side_dofs; ++sidei)
1326 : {
1327 0 : unsigned int i = new_side_dofs[sidei];
1328 : // fixed DoFs aren't test functions
1329 0 : if (dof_is_fixed[i])
1330 0 : continue;
1331 0 : for (unsigned int sidej=0, freej=0; sidej != n_new_side_dofs; ++sidej)
1332 : {
1333 0 : unsigned int j =
1334 : new_side_dofs[sidej];
1335 0 : if (dof_is_fixed[j])
1336 0 : Fe(freei) -=
1337 0 : TensorTools::inner_product(phi_coarse[i][qp],
1338 0 : phi_coarse[j][qp]) *
1339 0 : JxW[qp] * Ue(j);
1340 : else
1341 0 : Ke(freei,freej) +=
1342 0 : TensorTools::inner_product(phi_coarse[i][qp],
1343 0 : phi_coarse[j][qp]) *
1344 : JxW[qp];
1345 0 : if (cont == C_ONE)
1346 : {
1347 0 : if (dof_is_fixed[j])
1348 0 : Fe(freei) -=
1349 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1350 0 : (*dphi_coarse)[j][qp]) *
1351 0 : JxW[qp] * Ue(j);
1352 : else
1353 0 : Ke(freei,freej) +=
1354 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1355 0 : (*dphi_coarse)[j][qp]) *
1356 : JxW[qp];
1357 : }
1358 0 : if (!dof_is_fixed[j])
1359 0 : freej++;
1360 : }
1361 0 : Fe(freei) += TensorTools::inner_product(fineval, phi_coarse[i][qp]) * JxW[qp];
1362 0 : if (cont == C_ONE)
1363 0 : Fe(freei) +=
1364 0 : TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1365 0 : freei++;
1366 : }
1367 : }
1368 : }
1369 0 : Ke.cholesky_solve(Fe, Uside);
1370 :
1371 : // Transfer new side solutions to element
1372 0 : for (unsigned int i=0; i != free_dofs; ++i)
1373 : {
1374 0 : Number & ui = Ue(new_side_dofs[free_dof[i]]);
1375 0 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1376 : std::abs(ui - Uside(i)) < TOLERANCE);
1377 0 : ui = Uside(i);
1378 0 : 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 0 : unsigned int free_dofs = 0;
1387 0 : for (unsigned int i=0; i != new_n_dofs; ++i)
1388 0 : if (!dof_is_fixed[i])
1389 0 : free_dof[free_dofs++] = i;
1390 0 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1391 0 : Fe.resize (free_dofs); Fe.zero();
1392 : // The new interior coefficients
1393 0 : DenseVector<Number> Uint(free_dofs);
1394 :
1395 : // Add projection terms from each child
1396 0 : for (auto & child : elem->child_ref_range())
1397 : {
1398 0 : std::vector<dof_id_type> child_dof_indices;
1399 0 : if (use_old_dof_indices)
1400 0 : dof_map.old_dof_indices (&child,
1401 : child_dof_indices, var);
1402 : else
1403 0 : dof_map.dof_indices (&child,
1404 : child_dof_indices, var);
1405 : const unsigned int child_n_dofs =
1406 : cast_int<unsigned int>
1407 0 : (child_dof_indices.size());
1408 :
1409 : // Initialize both child and parent FE data
1410 : // on the child's quadrature points
1411 0 : fe->attach_quadrature_rule (qrule.get());
1412 0 : fe->reinit (&child);
1413 0 : const unsigned int n_qp = qrule->n_points();
1414 :
1415 0 : FEMap::inverse_map (dim, elem, xyz_values, coarse_qpoints);
1416 :
1417 0 : fe_coarse->reinit(elem, &coarse_qpoints);
1418 :
1419 : // Loop over the quadrature points
1420 0 : for (unsigned int qp=0; qp<n_qp; qp++)
1421 : {
1422 : // solution value at the quadrature point
1423 0 : OutputNumber fineval = libMesh::zero;
1424 : // solution grad at the quadrature point
1425 0 : 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 0 : for (unsigned int i=0; i<child_n_dofs; i++)
1431 : {
1432 0 : fineval +=
1433 0 : (old_vector(child_dof_indices[i]) *
1434 0 : phi_values[i][qp]);
1435 0 : if (cont == C_ONE)
1436 0 : finegrad += (*dphi_values)[i][qp] *
1437 0 : old_vector(child_dof_indices[i]);
1438 : }
1439 :
1440 : // Form interior projection matrix
1441 0 : for (unsigned int i=0, freei=0;
1442 0 : i != new_n_dofs; ++i)
1443 : {
1444 : // fixed DoFs aren't test functions
1445 0 : if (dof_is_fixed[i])
1446 0 : continue;
1447 0 : for (unsigned int j=0, freej=0; j !=
1448 : new_n_dofs; ++j)
1449 : {
1450 0 : if (dof_is_fixed[j])
1451 0 : Fe(freei) -=
1452 0 : TensorTools::inner_product(phi_coarse[i][qp],
1453 0 : phi_coarse[j][qp]) *
1454 0 : JxW[qp] * Ue(j);
1455 : else
1456 0 : Ke(freei,freej) +=
1457 0 : TensorTools::inner_product(phi_coarse[i][qp],
1458 0 : phi_coarse[j][qp]) *
1459 : JxW[qp];
1460 0 : if (cont == C_ONE)
1461 : {
1462 0 : if (dof_is_fixed[j])
1463 0 : Fe(freei) -=
1464 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1465 0 : (*dphi_coarse)[j][qp]) *
1466 0 : JxW[qp] * Ue(j);
1467 : else
1468 0 : Ke(freei,freej) +=
1469 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1470 0 : (*dphi_coarse)[j][qp]) *
1471 : JxW[qp];
1472 : }
1473 0 : if (!dof_is_fixed[j])
1474 0 : freej++;
1475 : }
1476 0 : Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp], fineval) *
1477 : JxW[qp];
1478 0 : if (cont == C_ONE)
1479 0 : Fe(freei) += TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1480 0 : freei++;
1481 : }
1482 : }
1483 : }
1484 0 : Ke.cholesky_solve(Fe, Uint);
1485 :
1486 : // Transfer new interior solutions to element
1487 0 : for (unsigned int i=0; i != free_dofs; ++i)
1488 : {
1489 0 : Number & ui = Ue(free_dof[i]);
1490 0 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1491 : std::abs(ui - Uint(i)) < TOLERANCE);
1492 0 : 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 0 : dof_is_fixed[free_dof[i]] = true;
1497 : #endif
1498 : }
1499 :
1500 : #ifndef NDEBUG
1501 : // Make sure every DoF got reached!
1502 0 : for (unsigned int i=0; i != new_n_dofs; ++i)
1503 0 : libmesh_assert(dof_is_fixed[i]);
1504 : #endif
1505 0 : }
1506 :
1507 :
1508 :
1509 : template <typename OutputType>
1510 : void
1511 0 : 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 0 : Ue.resize(0);
1518 :
1519 0 : for (auto v : make_range(dof_map.n_variables()))
1520 : {
1521 0 : DenseVector<Number> Usub;
1522 :
1523 0 : coarsened_dof_values(old_vector, dof_map, elem, Usub,
1524 : v, use_old_dof_indices);
1525 :
1526 0 : Ue.append (Usub);
1527 : }
1528 0 : }
1529 :
1530 :
1531 :
1532 : template <typename OutputType>
1533 : void
1534 671882 : FEGenericBase<OutputType>::compute_proj_constraints (DofConstraints & constraints,
1535 : DofMap & dof_map,
1536 : const unsigned int variable_number,
1537 : const Elem * elem)
1538 : {
1539 56462 : libmesh_assert(elem);
1540 :
1541 671882 : const unsigned int Dim = elem->dim();
1542 :
1543 : // Only constrain elements in 2,3D.
1544 671882 : if (Dim == 1)
1545 108039 : return;
1546 :
1547 : // Only constrain active elements with this method
1548 56462 : if (!elem->active())
1549 9226 : return;
1550 :
1551 563843 : const Variable & var = dof_map.variable(variable_number);
1552 47236 : const FEType & base_fe_type = var.type();
1553 563843 : const bool add_p_level = base_fe_type.p_refinement;
1554 :
1555 : // Construct FE objects for this element and its neighbors.
1556 563843 : std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1557 : (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1558 47236 : my_fe->add_p_level_in_reinit(add_p_level);
1559 563843 : const FEContinuity cont = my_fe->get_continuity();
1560 :
1561 : // We don't need to constrain discontinuous elements
1562 563843 : if (cont == DISCONTINUOUS)
1563 0 : return;
1564 47236 : 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 126075 : if (elem->infinite())
1572 0 : libmesh_not_implemented();
1573 :
1574 611079 : std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1575 : (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1576 47236 : neigh_fe->add_p_level_in_reinit(add_p_level);
1577 :
1578 611079 : QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1579 563843 : my_fe->attach_quadrature_rule (&my_qface);
1580 94472 : std::vector<Point> neigh_qface;
1581 :
1582 126075 : const std::vector<Real> & JxW = my_fe->get_JxW();
1583 126075 : const std::vector<Point> & q_point = my_fe->get_xyz();
1584 47236 : const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1585 47236 : const std::vector<std::vector<OutputShape>> & neigh_phi =
1586 : neigh_fe->get_phi();
1587 47236 : const std::vector<Point> * face_normals = nullptr;
1588 47236 : const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1589 47236 : const std::vector<std::vector<OutputGradient>> * neigh_dphi = nullptr;
1590 :
1591 94472 : std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1592 94472 : std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1593 :
1594 563843 : if (cont == C_ONE)
1595 : {
1596 7218 : const std::vector<Point> & ref_face_normals =
1597 7434 : my_fe->get_normals();
1598 3609 : face_normals = &ref_face_normals;
1599 3609 : const std::vector<std::vector<OutputGradient>> & ref_dphi =
1600 : my_fe->get_dphi();
1601 3609 : dphi = &ref_dphi;
1602 3609 : const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1603 : neigh_fe->get_dphi();
1604 3609 : neigh_dphi = &ref_neigh_dphi;
1605 : }
1606 :
1607 658315 : DenseMatrix<Real> Ke;
1608 563843 : DenseVector<Real> Fe;
1609 141708 : std::vector<DenseVector<Real>> Ue;
1610 :
1611 : // Look at the element faces. Check to see if we need to
1612 : // build constraints.
1613 2683416 : for (auto s : elem->side_index_range())
1614 : {
1615 : // Get pointers to the element's neighbor.
1616 2119573 : const Elem * neigh = elem->neighbor_ptr(s);
1617 :
1618 2119573 : if (!neigh)
1619 202457 : continue;
1620 :
1621 1898594 : if (!var.active_on_subdomain(neigh->subdomain_id()))
1622 1636 : continue;
1623 :
1624 : // h refinement constraints:
1625 : // constrain dofs shared between
1626 : // this element and ones coarser
1627 : // than this element.
1628 1896814 : if (neigh->level() < elem->level())
1629 : {
1630 119142 : unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
1631 10120 : 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 10120 : libmesh_assert(neigh->active());
1637 139382 : const unsigned int min_p_level = add_p_level *
1638 139382 : 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 119142 : const unsigned int old_elem_level = add_p_level * elem->p_level();
1642 119142 : if (old_elem_level != min_p_level)
1643 360 : my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + min_p_level - old_elem_level);
1644 119142 : const unsigned int old_neigh_level = add_p_level * neigh->p_level();
1645 119142 : if (old_neigh_level != min_p_level)
1646 0 : neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + min_p_level - old_neigh_level);
1647 :
1648 119142 : 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 119142 : my_dof_indices.reserve (elem->n_nodes());
1657 119142 : neigh_dof_indices.reserve (neigh->n_nodes());
1658 :
1659 119142 : dof_map.dof_indices (elem, my_dof_indices,
1660 : variable_number,
1661 : min_p_level);
1662 119142 : dof_map.dof_indices (neigh, neigh_dof_indices,
1663 : variable_number,
1664 : min_p_level);
1665 :
1666 10120 : const unsigned int n_qp = my_qface.n_points();
1667 :
1668 119142 : FEMap::inverse_map (Dim, neigh, q_point, neigh_qface);
1669 :
1670 119142 : 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 119142 : FEType elem_fe_type = base_fe_type;
1675 119142 : if (old_elem_level != min_p_level)
1676 360 : elem_fe_type.order = base_fe_type.order.get_order() + min_p_level - old_elem_level;
1677 119142 : FEType neigh_fe_type = base_fe_type;
1678 119142 : if (old_neigh_level != min_p_level)
1679 0 : neigh_fe_type.order = base_fe_type.order.get_order() + min_p_level - old_neigh_level;
1680 119142 : FEInterface::dofs_on_side(elem, Dim, elem_fe_type, s, my_side_dofs);
1681 119142 : FEInterface::dofs_on_side(neigh, Dim, neigh_fe_type, s_neigh, neigh_side_dofs);
1682 :
1683 10120 : const unsigned int n_side_dofs =
1684 20240 : cast_int<unsigned int>(my_side_dofs.size());
1685 10120 : libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
1686 :
1687 : #ifndef NDEBUG
1688 55500 : for (auto i : my_side_dofs)
1689 45380 : libmesh_assert_less(i, my_dof_indices.size());
1690 55500 : for (auto i : neigh_side_dofs)
1691 45380 : libmesh_assert_less(i, neigh_dof_indices.size());
1692 : #endif
1693 :
1694 109022 : Ke.resize (n_side_dofs, n_side_dofs);
1695 119142 : Ue.resize(n_side_dofs);
1696 :
1697 : // Form the projection matrix, (inner product of fine basis
1698 : // functions against fine test functions)
1699 655720 : for (unsigned int is = 0; is != n_side_dofs; ++is)
1700 : {
1701 581958 : const unsigned int i = my_side_dofs[is];
1702 3220000 : for (unsigned int js = 0; js != n_side_dofs; ++js)
1703 : {
1704 2908306 : const unsigned int j = my_side_dofs[js];
1705 14928726 : for (unsigned int qp = 0; qp != n_qp; ++qp)
1706 : {
1707 18281592 : Ke(is,js) += JxW[qp] * TensorTools::inner_product(phi[i][qp], phi[j][qp]);
1708 12245304 : if (cont == C_ONE)
1709 4220984 : Ke(is,js) += JxW[qp] *
1710 646400 : TensorTools::inner_product((*dphi)[i][qp] *
1711 : (*face_normals)[qp],
1712 1616000 : (*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 655720 : for (unsigned int is = 0; is != n_side_dofs; ++is)
1721 : {
1722 581958 : const unsigned int i = neigh_side_dofs[is];
1723 491198 : Fe.resize (n_side_dofs);
1724 3220000 : for (unsigned int js = 0; js != n_side_dofs; ++js)
1725 : {
1726 2908306 : const unsigned int j = my_side_dofs[js];
1727 14928726 : for (unsigned int qp = 0; qp != n_qp; ++qp)
1728 : {
1729 14257400 : Fe(js) += JxW[qp] *
1730 14257400 : TensorTools::inner_product(neigh_phi[i][qp],
1731 13251352 : phi[j][qp]);
1732 12245304 : if (cont == C_ONE)
1733 4220984 : Fe(js) += JxW[qp] *
1734 969600 : TensorTools::inner_product((*neigh_dphi)[i][qp] *
1735 : (*face_normals)[qp],
1736 1616000 : (*dphi)[j][qp] *
1737 : (*face_normals)[qp]);
1738 : }
1739 : }
1740 581958 : Ke.cholesky_solve(Fe, Ue[is]);
1741 : }
1742 :
1743 655720 : for (unsigned int js = 0; js != n_side_dofs; ++js)
1744 : {
1745 536578 : const unsigned int j = my_side_dofs[js];
1746 581958 : const dof_id_type my_dof_g = my_dof_indices[j];
1747 45380 : 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 45380 : bool self_constraint = false;
1752 2592091 : for (unsigned int is = 0; is != n_side_dofs; ++is)
1753 : {
1754 2217915 : const unsigned int i = neigh_side_dofs[is];
1755 2217915 : const dof_id_type their_dof_g = neigh_dof_indices[i];
1756 185750 : libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1757 :
1758 2217915 : if (their_dof_g == my_dof_g)
1759 : {
1760 : #ifndef NDEBUG
1761 13708 : const Real their_dof_value = Ue[is](js);
1762 13708 : libmesh_assert_less (std::abs(their_dof_value-1.),
1763 : 10*TOLERANCE);
1764 :
1765 86448 : for (unsigned int k = 0; k != n_side_dofs; ++k)
1766 72740 : libmesh_assert(k == is ||
1767 : std::abs(Ue[k](js)) <
1768 : 10*TOLERANCE);
1769 : #endif
1770 :
1771 13708 : self_constraint = true;
1772 13708 : break;
1773 : }
1774 : }
1775 :
1776 536578 : if (self_constraint)
1777 239232 : 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 31672 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1786 :
1787 374176 : if (dof_map.is_constrained_dof(my_dof_g))
1788 6863 : continue;
1789 :
1790 297346 : constraint_row = &(constraints[my_dof_g]);
1791 24809 : libmesh_assert(constraint_row->empty());
1792 : }
1793 :
1794 1685324 : for (unsigned int is = 0; is != n_side_dofs; ++is)
1795 : {
1796 1387978 : const unsigned int i = neigh_side_dofs[is];
1797 1387978 : const dof_id_type their_dof_g = neigh_dof_indices[i];
1798 114227 : libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1799 114227 : libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
1800 :
1801 1502205 : const Real their_dof_value = Ue[is](js);
1802 :
1803 1387978 : if (std::abs(their_dof_value) < 10*TOLERANCE)
1804 760373 : continue;
1805 :
1806 51566 : constraint_row->emplace(their_dof_g, their_dof_value);
1807 : }
1808 : }
1809 :
1810 119142 : my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + old_elem_level - min_p_level);
1811 119142 : neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + old_neigh_level - min_p_level);
1812 : }
1813 :
1814 1896814 : 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 2055899 : neigh->min_p_level_by_neighbor(elem, elem->p_level());
1822 2055899 : if (min_p_level < elem->p_level())
1823 : {
1824 : // Adaptive p refinement of non-hierarchic bases will
1825 : // require more coding
1826 48 : libmesh_assert(my_fe->is_hierarchic());
1827 576 : dof_map.constrain_p_dofs(variable_number, elem,
1828 : s, min_p_level);
1829 : }
1830 : }
1831 : }
1832 1408113 : }
1833 :
1834 : #endif // #ifdef LIBMESH_ENABLE_AMR
1835 :
1836 :
1837 :
1838 : #ifdef LIBMESH_ENABLE_PERIODIC
1839 : template <typename OutputType>
1840 : void
1841 264164 : 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 264164 : if (boundaries.empty())
1852 58473 : return;
1853 :
1854 67082 : libmesh_assert(elem);
1855 :
1856 : // Only constrain active elements with this method
1857 67082 : if (!elem->active())
1858 19491 : return;
1859 :
1860 138677 : if (elem->infinite())
1861 0 : libmesh_not_implemented();
1862 :
1863 205691 : const unsigned int Dim = elem->dim();
1864 :
1865 : // We need sys_number and variable_number for DofObject methods
1866 : // later
1867 95182 : const unsigned int sys_number = dof_map.sys_number();
1868 :
1869 47591 : const FEType & base_fe_type = dof_map.variable_type(variable_number);
1870 :
1871 : // Construct FE objects for this element and its pseudo-neighbors.
1872 205691 : std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1873 : (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1874 205691 : const FEContinuity cont = my_fe->get_continuity();
1875 :
1876 : // We don't need to constrain discontinuous elements
1877 205691 : if (cont == DISCONTINUOUS)
1878 0 : return;
1879 47591 : libmesh_assert (cont == C_ZERO || cont == C_ONE);
1880 :
1881 : // We'll use element size to generate relative tolerances later
1882 205691 : const Real primary_hmin = elem->hmin();
1883 :
1884 253282 : std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1885 : (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1886 :
1887 253282 : QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1888 205691 : my_fe->attach_quadrature_rule (&my_qface);
1889 95182 : std::vector<Point> neigh_qface;
1890 :
1891 138677 : const std::vector<Real> & JxW = my_fe->get_JxW();
1892 138677 : const std::vector<Point> & q_point = my_fe->get_xyz();
1893 47591 : const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1894 47591 : const std::vector<std::vector<OutputShape>> & neigh_phi =
1895 : neigh_fe->get_phi();
1896 47591 : const std::vector<Point> * face_normals = nullptr;
1897 47591 : const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1898 47591 : const std::vector<std::vector<OutputGradient>> * neigh_dphi = nullptr;
1899 95182 : std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1900 95182 : std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1901 :
1902 205691 : if (cont != C_ZERO)
1903 : {
1904 8192 : const std::vector<Point> & ref_face_normals =
1905 4096 : my_fe->get_normals();
1906 4096 : face_normals = &ref_face_normals;
1907 4096 : const std::vector<std::vector<OutputGradient>> & ref_dphi =
1908 : my_fe->get_dphi();
1909 4096 : dphi = &ref_dphi;
1910 4096 : const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1911 : neigh_fe->get_dphi();
1912 4096 : neigh_dphi = &ref_neigh_dphi;
1913 : }
1914 :
1915 300873 : DenseMatrix<Real> Ke;
1916 205691 : DenseVector<Real> Fe;
1917 142773 : std::vector<DenseVector<Real>> Ue;
1918 :
1919 : // Container to catch the boundary ids that BoundaryInfo hands us.
1920 95182 : 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 205691 : const unsigned short int max_ns = elem->n_sides();
1925 1018023 : for (unsigned short int s = 0; s != max_ns; ++s)
1926 : {
1927 1001144 : if (elem->neighbor_ptr(s))
1928 589194 : continue;
1929 :
1930 39424 : mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
1931 :
1932 59984 : for (const auto & boundary_id : bc_ids)
1933 : {
1934 20560 : const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
1935 20560 : if (!periodic || !periodic->is_my_variable(variable_number))
1936 5784 : continue;
1937 :
1938 3044 : libmesh_assert(point_locator);
1939 :
1940 : // Get pointers to the element's neighbor.
1941 : unsigned int s_neigh;
1942 14776 : const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s, &s_neigh);
1943 :
1944 14776 : 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 14776 : 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 2596 : libmesh_assert(neigh->active());
1958 13432 : const unsigned int min_p_level =
1959 18624 : 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 2596 : const unsigned int old_elem_level = elem->p_level();
1967 13432 : if (old_elem_level != min_p_level)
1968 0 : (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
1969 5192 : const unsigned int old_neigh_level = neigh->p_level();
1970 13432 : if (old_neigh_level != min_p_level)
1971 0 : (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 13432 : my_fe->reinit(elem, s);
1981 :
1982 13432 : dof_map.dof_indices (elem, my_dof_indices,
1983 : variable_number);
1984 13432 : 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 7788 : std::vector<std::vector<dof_id_type>> neigh_dof_indices_all_variables;
1991 13432 : if(periodic->has_transformation_matrix())
1992 : {
1993 7200 : const std::set<unsigned int> & variables = periodic->get_variables();
1994 7200 : neigh_dof_indices_all_variables.resize(variables.size());
1995 600 : unsigned int index = 0;
1996 28800 : for(unsigned int var : variables)
1997 : {
1998 23400 : dof_map.dof_indices (neigh, neigh_dof_indices_all_variables[index],
1999 : var);
2000 21600 : index++;
2001 : }
2002 : }
2003 :
2004 2596 : const unsigned int n_qp = my_qface.n_points();
2005 :
2006 : // Translate the quadrature points over to the
2007 : // neighbor's boundary
2008 18624 : std::vector<Point> neigh_point(q_point.size());
2009 54848 : for (auto i : index_range(neigh_point))
2010 47820 : neigh_point[i] = periodic->get_corresponding_pos(q_point[i]);
2011 :
2012 13432 : FEMap::inverse_map (Dim, neigh, neigh_point,
2013 : neigh_qface);
2014 :
2015 13432 : 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 13432 : FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
2020 13432 : 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 16028 : if (elem->p_level() != old_elem_level)
2026 0 : (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
2027 16028 : if (neigh->p_level() != old_neigh_level)
2028 0 : (const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
2029 : #endif // #ifdef LIBMESH_ENABLE_AMR
2030 :
2031 2596 : const unsigned int n_side_dofs =
2032 : cast_int<unsigned int>
2033 5192 : (my_side_dofs.size());
2034 2596 : libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
2035 :
2036 10836 : Ke.resize (n_side_dofs, n_side_dofs);
2037 13432 : Ue.resize(n_side_dofs);
2038 :
2039 : // Form the projection matrix, (inner product of fine basis
2040 : // functions against fine test functions)
2041 54936 : for (unsigned int is = 0; is != n_side_dofs; ++is)
2042 : {
2043 47916 : const unsigned int i = my_side_dofs[is];
2044 182832 : for (unsigned int js = 0; js != n_side_dofs; ++js)
2045 : {
2046 159012 : const unsigned int j = my_side_dofs[js];
2047 656192 : for (unsigned int qp = 0; qp != n_qp; ++qp)
2048 : {
2049 624296 : Ke(is,js) += JxW[qp] *
2050 624296 : TensorTools::inner_product(phi[i][qp],
2051 569580 : phi[j][qp]);
2052 514864 : if (cont != C_ZERO)
2053 384 : Ke(is,js) += JxW[qp] *
2054 64 : TensorTools::inner_product((*dphi)[i][qp] *
2055 : (*face_normals)[qp],
2056 160 : (*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 54936 : for (unsigned int is = 0; is != n_side_dofs; ++is)
2065 : {
2066 47916 : const unsigned int i = neigh_side_dofs[is];
2067 35092 : Fe.resize (n_side_dofs);
2068 182832 : for (unsigned int js = 0; js != n_side_dofs; ++js)
2069 : {
2070 159012 : const unsigned int j = my_side_dofs[js];
2071 656192 : for (unsigned int qp = 0; qp != n_qp; ++qp)
2072 : {
2073 624296 : Fe(js) += JxW[qp] *
2074 624296 : TensorTools::inner_product(neigh_phi[i][qp],
2075 569580 : phi[j][qp]);
2076 514864 : if (cont != C_ZERO)
2077 384 : Fe(js) += JxW[qp] *
2078 96 : TensorTools::inner_product((*neigh_dphi)[i][qp] *
2079 : (*face_normals)[qp],
2080 160 : (*dphi)[j][qp] *
2081 : (*face_normals)[qp]);
2082 : }
2083 : }
2084 47916 : 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 5192 : std::set<dof_id_type> my_constrained_dofs;
2136 :
2137 : // Container to catch boundary IDs handed back by BoundaryInfo.
2138 5192 : std::vector<boundary_id_type> new_bc_ids;
2139 :
2140 96264 : for (auto n : elem->node_index_range())
2141 : {
2142 82832 : if (!elem->is_node_on_side(n,s))
2143 35012 : continue;
2144 :
2145 6404 : const Node & my_node = elem->node_ref(n);
2146 :
2147 41416 : 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 6384 : std::set<boundary_id_type> point_bcids;
2153 :
2154 256440 : for (unsigned int new_s = 0;
2155 262824 : new_s != max_ns; ++new_s)
2156 : {
2157 221648 : if (!elem->is_node_on_side(n,new_s))
2158 95464 : continue;
2159 :
2160 111064 : mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
2161 :
2162 222128 : for (const auto & new_boundary_id : new_bc_ids)
2163 : {
2164 111064 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2165 111064 : if (new_periodic && new_periodic->is_my_variable(variable_number))
2166 95904 : point_bcids.insert(new_boundary_id);
2167 : }
2168 : }
2169 :
2170 : // See if this vertex has point neighbors to
2171 : // defer to
2172 44655 : if (primary_boundary_point_neighbor
2173 41176 : (elem, my_node, mesh.get_boundary_info(), point_bcids)
2174 6384 : != elem)
2175 21511 : continue;
2176 :
2177 : // Find the complementary boundary id set
2178 2905 : std::set<boundary_id_type> point_pairedids;
2179 32372 : for (const auto & new_boundary_id : point_bcids)
2180 : {
2181 16186 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2182 16186 : point_pairedids.insert(new_periodic->pairedboundary);
2183 : }
2184 :
2185 : // What do we want to constrain against?
2186 2905 : const Elem * primary_elem = nullptr;
2187 2905 : const Elem * main_neigh = nullptr;
2188 16186 : Point main_pt = my_node,
2189 16186 : primary_pt = my_node;
2190 :
2191 32372 : for (const auto & new_boundary_id : point_bcids)
2192 : {
2193 : // Find the corresponding periodic point and
2194 : // its primary neighbor
2195 16186 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2196 :
2197 : const Point neigh_pt =
2198 16186 : 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 16186 : if (neigh_pt.absolute_fuzzy_equals
2204 16186 : (my_node, primary_hmin*TOLERANCE))
2205 6828 : continue;
2206 :
2207 : // Otherwise we'll have a constraint in
2208 : // one direction or another
2209 16186 : if (!primary_elem)
2210 2905 : primary_elem = elem;
2211 :
2212 2905 : const Elem * primary_neigh =
2213 16186 : primary_boundary_point_neighbor(neigh, neigh_pt,
2214 : mesh.get_boundary_info(),
2215 : point_pairedids);
2216 :
2217 2905 : libmesh_assert(primary_neigh);
2218 :
2219 16186 : if (new_boundary_id == boundary_id)
2220 : {
2221 2905 : main_neigh = primary_neigh;
2222 16186 : 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 29467 : 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 25646 : (primary_neigh->level() == primary_elem->level() &&
2234 28558 : 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 1880 : (primary_neigh == primary_elem &&
2239 0 : (neigh_pt > primary_pt)))
2240 6828 : continue;
2241 :
2242 1880 : primary_elem = primary_neigh;
2243 9358 : primary_pt = neigh_pt;
2244 : }
2245 :
2246 13493 : if (!primary_elem ||
2247 19091 : primary_elem != main_neigh ||
2248 1880 : primary_pt != main_pt)
2249 1025 : continue;
2250 : }
2251 240 : else if (elem->is_edge(n))
2252 : {
2253 : // Find which edge we're on
2254 240 : unsigned int e=0, ne = elem->n_edges();
2255 384 : for (; e != ne; ++e)
2256 : {
2257 384 : if (elem->is_node_on_edge(n,e))
2258 20 : break;
2259 : }
2260 20 : libmesh_assert_less (e, elem->n_edges());
2261 :
2262 : // Find the edge end nodes
2263 : const Node
2264 20 : * e1 = nullptr,
2265 20 : * e2 = nullptr;
2266 612 : for (auto nn : elem->node_index_range())
2267 : {
2268 612 : if (nn == n)
2269 0 : continue;
2270 :
2271 612 : if (elem->is_node_on_edge(nn, e))
2272 : {
2273 480 : if (e1 == nullptr)
2274 : {
2275 40 : e1 = elem->node_ptr(nn);
2276 : }
2277 : else
2278 : {
2279 40 : e2 = elem->node_ptr(nn);
2280 240 : break;
2281 : }
2282 : }
2283 : }
2284 20 : 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 20 : std::set<boundary_id_type> edge_bcids;
2290 :
2291 940 : for (unsigned int new_s = 0;
2292 960 : new_s != max_ns; ++new_s)
2293 : {
2294 720 : if (!elem->is_node_on_side(n,new_s))
2295 440 : continue;
2296 :
2297 : // We're reusing the new_bc_ids vector created outside the loop over nodes.
2298 240 : mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
2299 :
2300 480 : for (const auto & new_boundary_id : new_bc_ids)
2301 : {
2302 240 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2303 240 : if (new_periodic && new_periodic->is_my_variable(variable_number))
2304 220 : edge_bcids.insert(new_boundary_id);
2305 : }
2306 : }
2307 :
2308 :
2309 : // See if this edge has neighbors to defer to
2310 240 : if (primary_boundary_edge_neighbor
2311 240 : (elem, *e1, *e2, mesh.get_boundary_info(), edge_bcids)
2312 20 : != elem)
2313 0 : continue;
2314 :
2315 : // Find the complementary boundary id set
2316 20 : std::set<boundary_id_type> edge_pairedids;
2317 480 : for (const auto & new_boundary_id : edge_bcids)
2318 : {
2319 240 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2320 240 : edge_pairedids.insert(new_periodic->pairedboundary);
2321 : }
2322 :
2323 : // What do we want to constrain against?
2324 20 : const Elem * primary_elem = nullptr;
2325 20 : const Elem * main_neigh = nullptr;
2326 240 : Point main_pt1 = *e1,
2327 240 : main_pt2 = *e2,
2328 240 : primary_pt1 = *e1,
2329 240 : primary_pt2 = *e2;
2330 :
2331 480 : for (const auto & new_boundary_id : edge_bcids)
2332 : {
2333 : // Find the corresponding periodic edge and
2334 : // its primary neighbor
2335 240 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2336 :
2337 240 : Point neigh_pt1 = new_periodic->get_corresponding_pos(*e1),
2338 240 : 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 20 : if (neigh_pt1.absolute_fuzzy_equals
2344 260 : (*e1, primary_hmin*TOLERANCE) &&
2345 : neigh_pt2.absolute_fuzzy_equals
2346 0 : (*e2, primary_hmin*TOLERANCE))
2347 120 : continue;
2348 :
2349 : // Otherwise we'll have a constraint in
2350 : // one direction or another
2351 240 : if (!primary_elem)
2352 20 : primary_elem = elem;
2353 :
2354 20 : const Elem * primary_neigh = primary_boundary_edge_neighbor
2355 240 : (neigh, neigh_pt1, neigh_pt2,
2356 : mesh.get_boundary_info(), edge_pairedids);
2357 :
2358 20 : libmesh_assert(primary_neigh);
2359 :
2360 240 : if (new_boundary_id == boundary_id)
2361 : {
2362 20 : main_neigh = primary_neigh;
2363 240 : main_pt1 = neigh_pt1;
2364 240 : 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 240 : if (primary_neigh == primary_elem)
2375 : {
2376 0 : if (primary_pt1 > primary_pt2)
2377 0 : std::swap(primary_pt1, primary_pt2);
2378 0 : if (neigh_pt1 > neigh_pt2)
2379 0 : std::swap(neigh_pt1, neigh_pt2);
2380 :
2381 0 : if (neigh_pt2 >= primary_pt2)
2382 0 : continue;
2383 : }
2384 :
2385 : // Otherwise:
2386 : // Finer elements will get constrained in
2387 : // terms of coarser ones, not the other way
2388 : // around
2389 480 : 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 440 : (primary_neigh->level() == primary_elem->level() &&
2395 40 : primary_neigh->id() > primary_elem->id()))
2396 120 : continue;
2397 :
2398 10 : primary_elem = primary_neigh;
2399 120 : primary_pt1 = neigh_pt1;
2400 120 : primary_pt2 = neigh_pt2;
2401 : }
2402 :
2403 160 : if (!primary_elem ||
2404 230 : primary_elem != main_neigh ||
2405 270 : primary_pt1 != main_pt1 ||
2406 10 : primary_pt2 != main_pt2)
2407 10 : continue;
2408 : }
2409 0 : 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 0 : if (neigh == elem)
2416 : {
2417 : const Point neigh_pt =
2418 0 : periodic->get_corresponding_pos(my_node);
2419 0 : if (neigh_pt > my_node)
2420 0 : continue;
2421 : }
2422 :
2423 : // Otherwise:
2424 : // Finer elements will get constrained in
2425 : // terms of coarser ones, not the other way
2426 : // around
2427 0 : 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 0 : (neigh->level() == elem->level() &&
2433 0 : neigh->id() > elem->id()))
2434 0 : 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 9478 : my_node.n_comp(sys_number, variable_number);
2443 :
2444 19000 : for (unsigned int i=0; i != n_comp; ++i)
2445 : my_constrained_dofs.insert
2446 7628 : (my_node.dof_number
2447 9522 : (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 54936 : 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 41504 : const unsigned int j = my_side_dofs[js];
2491 47916 : const dof_id_type my_dof_g = my_dof_indices[j];
2492 6412 : libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2493 :
2494 : // FIXME: new code path
2495 29358 : if (!my_constrained_dofs.count(my_dof_g))
2496 32254 : 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 1894 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2505 :
2506 9522 : if (dof_map.is_constrained_dof(my_dof_g))
2507 81 : continue;
2508 :
2509 9250 : constraint_row = &(constraints[my_dof_g]);
2510 1813 : libmesh_assert(constraint_row->empty());
2511 : }
2512 :
2513 37506 : for (unsigned int is = 0; is != n_side_dofs; ++is)
2514 : {
2515 28256 : const unsigned int i = neigh_side_dofs[is];
2516 28256 : const dof_id_type their_dof_g = neigh_dof_indices[i];
2517 4439 : 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 28256 : const Real their_dof_value = Ue[is](js);
2524 :
2525 28256 : if (their_dof_g == my_dof_g)
2526 : {
2527 0 : libmesh_assert_less (std::abs(their_dof_value-1.), 1.e-5);
2528 0 : for (unsigned int k = 0; k != n_side_dofs; ++k)
2529 0 : libmesh_assert(k == is || std::abs(Ue[k](js)) < 1.e-5);
2530 17662 : continue;
2531 0 : }
2532 :
2533 28256 : if (std::abs(their_dof_value) < 10*TOLERANCE)
2534 15484 : continue;
2535 :
2536 10594 : if(!periodic->has_transformation_matrix())
2537 : {
2538 1865 : 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 4752 : const std::set<unsigned int> & variables = periodic->get_variables();
2549 4752 : neigh_dof_indices_all_variables.resize(variables.size());
2550 396 : unsigned int index = 0;
2551 19008 : for(unsigned int other_var : variables)
2552 : {
2553 1188 : libmesh_assert_msg(base_fe_type == dof_map.variable_type(other_var), "FE types must match for all variables involved in constraint");
2554 :
2555 14256 : Real var_weighting = periodic->get_transformation_matrix()(variable_number, other_var);
2556 14256 : constraint_row->emplace(neigh_dof_indices_all_variables[index][i],
2557 14256 : var_weighting*their_dof_value);
2558 14256 : index++;
2559 : }
2560 : }
2561 :
2562 : }
2563 : }
2564 8240 : }
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 17820 : neigh->min_p_level_by_neighbor(elem, elem->p_level());
2572 17820 : if (min_p_level < elem->p_level())
2573 : {
2574 : // Adaptive p refinement of non-hierarchic bases will
2575 : // require more coding
2576 0 : libmesh_assert(my_fe->is_hierarchic());
2577 0 : dof_map.constrain_p_dofs(variable_number, elem,
2578 : s, min_p_level);
2579 : }
2580 : #endif // #ifdef LIBMESH_ENABLE_AMR
2581 : }
2582 : }
2583 331527 : }
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
|
