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