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