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