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 320615 : if ((pt_neighbor->level() > primary->level()) ||
80 275898 : (pt_neighbor->level() == primary->level() &&
81 157911 : pt_neighbor->id() < primary->id()))
82 62370 : 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 14148 : bool vertex_on_periodic_side = false;
88 245533 : for (auto ns : pt_neighbor->side_index_range())
89 : {
90 230808 : boundary_info.boundary_ids (pt_neighbor, ns, bc_ids);
91 :
92 33807 : bool on_relevant_boundary = false;
93 461616 : for (const auto & id : boundary_ids)
94 230808 : if (std::find(bc_ids.begin(), bc_ids.end(), id) != bc_ids.end())
95 14091 : on_relevant_boundary = true;
96 :
97 230808 : if (!on_relevant_boundary)
98 112410 : continue;
99 :
100 98771 : pt_neighbor->build_side_ptr(periodic_side, ns);
101 98771 : if (!periodic_side->contains_point(p))
102 193 : continue;
103 :
104 14080 : vertex_on_periodic_side = true;
105 14080 : break;
106 : }
107 :
108 99200 : if (vertex_on_periodic_side)
109 98564 : 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 5554315 : FEGenericBase<Real>::build (const unsigned int dim,
192 : const FEType & fet)
193 : {
194 5554315 : 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 L2_LAGRANGE_GLL:
214 0 : return std::make_unique<FE<0,L2_LAGRANGE_GLL>>(fet);
215 :
216 0 : case HIERARCHIC:
217 0 : return std::make_unique<FE<0,HIERARCHIC>>(fet);
218 :
219 0 : case L2_HIERARCHIC:
220 0 : return std::make_unique<FE<0,L2_HIERARCHIC>>(fet);
221 :
222 0 : case SIDE_HIERARCHIC:
223 0 : return std::make_unique<FE<0,SIDE_HIERARCHIC>>(fet);
224 :
225 0 : case MONOMIAL:
226 0 : return std::make_unique<FE<0,MONOMIAL>>(fet);
227 :
228 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
229 0 : case SZABAB:
230 0 : return std::make_unique<FE<0,SZABAB>>(fet);
231 :
232 0 : case BERNSTEIN:
233 0 : return std::make_unique<FE<0,BERNSTEIN>>(fet);
234 :
235 0 : case RATIONAL_BERNSTEIN:
236 0 : return std::make_unique<FE<0,RATIONAL_BERNSTEIN>>(fet);
237 : #endif
238 :
239 0 : case XYZ:
240 0 : return std::make_unique<FEXYZ<0>>(fet);
241 :
242 0 : case SCALAR:
243 0 : return std::make_unique<FEScalar<0>>(fet);
244 :
245 0 : default:
246 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
247 : }
248 : }
249 : // 1D
250 458974 : case 1:
251 : {
252 458974 : switch (fet.family)
253 : {
254 0 : case CLOUGH:
255 0 : return std::make_unique<FE<1,CLOUGH>>(fet);
256 :
257 399080 : case HERMITE:
258 399080 : return std::make_unique<FE<1,HERMITE>>(fet);
259 :
260 4380 : case LAGRANGE:
261 4380 : return std::make_unique<FE<1,LAGRANGE>>(fet);
262 :
263 2982 : case L2_LAGRANGE:
264 2982 : return std::make_unique<FE<1,L2_LAGRANGE>>(fet);
265 :
266 12354 : case L2_LAGRANGE_GLL:
267 12354 : return std::make_unique<FE<1,L2_LAGRANGE_GLL>>(fet);
268 :
269 12350 : case HIERARCHIC:
270 12350 : return std::make_unique<FE<1,HIERARCHIC>>(fet);
271 :
272 4970 : case L2_HIERARCHIC:
273 4970 : return std::make_unique<FE<1,L2_HIERARCHIC>>(fet);
274 :
275 3550 : case SIDE_HIERARCHIC:
276 3550 : return std::make_unique<FE<1,SIDE_HIERARCHIC>>(fet);
277 :
278 5012 : case MONOMIAL:
279 5012 : return std::make_unique<FE<1,MONOMIAL>>(fet);
280 :
281 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
282 3976 : case SZABAB:
283 3976 : return std::make_unique<FE<1,SZABAB>>(fet);
284 :
285 3976 : case BERNSTEIN:
286 3976 : return std::make_unique<FE<1,BERNSTEIN>>(fet);
287 :
288 2368 : case RATIONAL_BERNSTEIN:
289 2368 : return std::make_unique<FE<1,RATIONAL_BERNSTEIN>>(fet);
290 : #endif
291 :
292 3976 : case XYZ:
293 3976 : return std::make_unique<FEXYZ<1>>(fet);
294 :
295 0 : case SCALAR:
296 0 : return std::make_unique<FEScalar<1>>(fet);
297 :
298 0 : default:
299 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
300 : }
301 : }
302 :
303 :
304 : // 2D
305 1726013 : case 2:
306 : {
307 1726013 : switch (fet.family)
308 : {
309 70656 : case CLOUGH:
310 70656 : return std::make_unique<FE<2,CLOUGH>>(fet);
311 :
312 36038 : case HERMITE:
313 36038 : return std::make_unique<FE<2,HERMITE>>(fet);
314 :
315 540455 : case LAGRANGE:
316 540455 : return std::make_unique<FE<2,LAGRANGE>>(fet);
317 :
318 56860 : case L2_LAGRANGE:
319 56860 : return std::make_unique<FE<2,L2_LAGRANGE>>(fet);
320 :
321 12496 : case L2_LAGRANGE_GLL:
322 12496 : return std::make_unique<FE<2,L2_LAGRANGE_GLL>>(fet);
323 :
324 543558 : case HIERARCHIC:
325 543558 : return std::make_unique<FE<2,HIERARCHIC>>(fet);
326 :
327 103866 : case L2_HIERARCHIC:
328 103866 : return std::make_unique<FE<2,L2_HIERARCHIC>>(fet);
329 :
330 147652 : case SIDE_HIERARCHIC:
331 147652 : return std::make_unique<FE<2,SIDE_HIERARCHIC>>(fet);
332 :
333 23816 : case MONOMIAL:
334 23816 : return std::make_unique<FE<2,MONOMIAL>>(fet);
335 :
336 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
337 16934 : case SZABAB:
338 16934 : return std::make_unique<FE<2,SZABAB>>(fet);
339 :
340 16550 : case BERNSTEIN:
341 16550 : return std::make_unique<FE<2,BERNSTEIN>>(fet);
342 :
343 118705 : case RATIONAL_BERNSTEIN:
344 118705 : return std::make_unique<FE<2,RATIONAL_BERNSTEIN>>(fet);
345 : #endif
346 :
347 38356 : case XYZ:
348 38356 : return std::make_unique<FEXYZ<2>>(fet);
349 :
350 0 : case SCALAR:
351 0 : return std::make_unique<FEScalar<2>>(fet);
352 :
353 71 : case SUBDIVISION:
354 71 : return std::make_unique<FESubdivision>(fet);
355 :
356 0 : default:
357 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
358 : }
359 : }
360 :
361 :
362 : // 3D
363 3369316 : case 3:
364 : {
365 3369316 : switch (fet.family)
366 : {
367 0 : case CLOUGH:
368 0 : libmesh_error_msg("ERROR: Clough-Tocher elements currently only support 1D and 2D");
369 :
370 1306 : case HERMITE:
371 1306 : return std::make_unique<FE<3,HERMITE>>(fet);
372 :
373 2603133 : case LAGRANGE:
374 2603133 : return std::make_unique<FE<3,LAGRANGE>>(fet);
375 :
376 89135 : case L2_LAGRANGE:
377 89135 : return std::make_unique<FE<3,L2_LAGRANGE>>(fet);
378 :
379 12354 : case L2_LAGRANGE_GLL:
380 12354 : return std::make_unique<FE<3,L2_LAGRANGE_GLL>>(fet);
381 :
382 83402 : case HIERARCHIC:
383 83402 : return std::make_unique<FE<3,HIERARCHIC>>(fet);
384 :
385 94460 : case L2_HIERARCHIC:
386 94460 : return std::make_unique<FE<3,L2_HIERARCHIC>>(fet);
387 :
388 305832 : case SIDE_HIERARCHIC:
389 305832 : return std::make_unique<FE<3,SIDE_HIERARCHIC>>(fet);
390 :
391 33827 : case MONOMIAL:
392 33827 : return std::make_unique<FE<3,MONOMIAL>>(fet);
393 :
394 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
395 0 : case SZABAB:
396 0 : return std::make_unique<FE<3,SZABAB>>(fet);
397 :
398 30670 : case BERNSTEIN:
399 30670 : return std::make_unique<FE<3,BERNSTEIN>>(fet);
400 :
401 15767 : case RATIONAL_BERNSTEIN:
402 15767 : return std::make_unique<FE<3,RATIONAL_BERNSTEIN>>(fet);
403 : #endif
404 :
405 99430 : case XYZ:
406 99430 : return std::make_unique<FEXYZ<3>>(fet);
407 :
408 0 : case SCALAR:
409 0 : return std::make_unique<FEScalar<3>>(fet);
410 :
411 0 : default:
412 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
413 : }
414 : }
415 :
416 0 : default:
417 0 : libmesh_error_msg("Invalid dimension dim = " << dim);
418 : }
419 : }
420 :
421 :
422 :
423 : template <>
424 : std::unique_ptr<FEGenericBase<RealGradient>>
425 574854 : FEGenericBase<RealGradient>::build (const unsigned int dim,
426 : const FEType & fet)
427 : {
428 574854 : switch (dim)
429 : {
430 : // 0D
431 0 : case 0:
432 : {
433 0 : switch (fet.family)
434 : {
435 0 : case HIERARCHIC_VEC:
436 0 : return std::make_unique<FEHierarchicVec<0>>(fet);
437 :
438 0 : case L2_HIERARCHIC_VEC:
439 0 : return std::make_unique<FEL2HierarchicVec<0>>(fet);
440 :
441 0 : case LAGRANGE_VEC:
442 0 : return std::make_unique<FELagrangeVec<0>>(fet);
443 :
444 0 : case L2_LAGRANGE_VEC:
445 0 : return std::make_unique<FEL2LagrangeVec<0>>(fet);
446 :
447 0 : case MONOMIAL_VEC:
448 0 : return std::make_unique<FEMonomialVec<0>>(fet);
449 :
450 0 : default:
451 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
452 : }
453 : }
454 0 : case 1:
455 : {
456 0 : switch (fet.family)
457 : {
458 0 : case HIERARCHIC_VEC:
459 0 : return std::make_unique<FEHierarchicVec<1>>(fet);
460 :
461 0 : case L2_HIERARCHIC_VEC:
462 0 : return std::make_unique<FEL2HierarchicVec<1>>(fet);
463 :
464 0 : case LAGRANGE_VEC:
465 0 : return std::make_unique<FELagrangeVec<1>>(fet);
466 :
467 0 : case L2_LAGRANGE_VEC:
468 0 : return std::make_unique<FEL2LagrangeVec<1>>(fet);
469 :
470 0 : case MONOMIAL_VEC:
471 0 : return std::make_unique<FEMonomialVec<1>>(fet);
472 :
473 0 : default:
474 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
475 : }
476 : }
477 212172 : case 2:
478 : {
479 212172 : switch (fet.family)
480 : {
481 16572 : case HIERARCHIC_VEC:
482 16572 : return std::make_unique<FEHierarchicVec<2>>(fet);
483 :
484 2840 : case L2_HIERARCHIC_VEC:
485 2840 : return std::make_unique<FEL2HierarchicVec<2>>(fet);
486 :
487 15977 : case LAGRANGE_VEC:
488 15977 : return std::make_unique<FELagrangeVec<2>>(fet);
489 :
490 3250 : case L2_LAGRANGE_VEC:
491 3250 : return std::make_unique<FEL2LagrangeVec<2>>(fet);
492 :
493 1470 : case MONOMIAL_VEC:
494 1470 : return std::make_unique<FEMonomialVec<2>>(fet);
495 :
496 56033 : case NEDELEC_ONE:
497 56033 : return std::make_unique<FENedelecOne<2>>(fet);
498 :
499 114610 : case RAVIART_THOMAS:
500 114610 : return std::make_unique<FERaviartThomas<2>>(fet);
501 :
502 1420 : case L2_RAVIART_THOMAS:
503 1420 : return std::make_unique<FEL2RaviartThomas<2>>(fet);
504 :
505 0 : default:
506 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
507 : }
508 : }
509 362682 : case 3:
510 : {
511 362682 : switch (fet.family)
512 : {
513 0 : case HIERARCHIC_VEC:
514 0 : return std::make_unique<FEHierarchicVec<3>>(fet);
515 :
516 710 : case L2_HIERARCHIC_VEC:
517 710 : return std::make_unique<FEL2HierarchicVec<3>>(fet);
518 :
519 22913 : case LAGRANGE_VEC:
520 22913 : return std::make_unique<FELagrangeVec<3>>(fet);
521 :
522 1420 : case L2_LAGRANGE_VEC:
523 1420 : return std::make_unique<FEL2LagrangeVec<3>>(fet);
524 :
525 0 : case MONOMIAL_VEC:
526 0 : return std::make_unique<FEMonomialVec<3>>(fet);
527 :
528 33874 : case NEDELEC_ONE:
529 33874 : return std::make_unique<FENedelecOne<3>>(fet);
530 :
531 303055 : case RAVIART_THOMAS:
532 303055 : return std::make_unique<FERaviartThomas<3>>(fet);
533 :
534 710 : case L2_RAVIART_THOMAS:
535 710 : return std::make_unique<FEL2RaviartThomas<3>>(fet);
536 :
537 0 : default:
538 0 : libmesh_error_msg("ERROR: Bad FEType.family == " << Utility::enum_to_string(fet.family));
539 : }
540 : }
541 :
542 0 : default:
543 0 : libmesh_error_msg("Invalid dimension dim = " << dim);
544 : } // switch(dim)
545 : }
546 :
547 :
548 :
549 :
550 :
551 :
552 :
553 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
554 :
555 :
556 : template <>
557 : std::unique_ptr<FEGenericBase<Real>>
558 452 : FEGenericBase<Real>::build_InfFE (const unsigned int dim,
559 : const FEType & fet)
560 : {
561 452 : switch (dim)
562 : {
563 :
564 : // 1D
565 0 : case 1:
566 : {
567 0 : switch (fet.radial_family)
568 : {
569 0 : case INFINITE_MAP:
570 0 : libmesh_error_msg("ERROR: Can't build an infinite element with FEFamily = " << Utility::enum_to_string(fet.radial_family));
571 :
572 0 : case JACOBI_20_00:
573 : {
574 0 : switch (fet.inf_map)
575 : {
576 0 : case CARTESIAN:
577 0 : return std::make_unique<InfFE<1,JACOBI_20_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 JACOBI_30_00:
585 : {
586 0 : switch (fet.inf_map)
587 : {
588 0 : case CARTESIAN:
589 0 : return std::make_unique<InfFE<1,JACOBI_30_00,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 LEGENDRE:
597 : {
598 0 : switch (fet.inf_map)
599 : {
600 0 : case CARTESIAN:
601 0 : return std::make_unique<InfFE<1,LEGENDRE,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 : case LAGRANGE:
609 : {
610 0 : switch (fet.inf_map)
611 : {
612 0 : case CARTESIAN:
613 0 : return std::make_unique<InfFE<1,LAGRANGE,CARTESIAN>>(fet);
614 :
615 0 : default:
616 0 : libmesh_error_msg("ERROR: Can't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
617 : }
618 : }
619 :
620 0 : default:
621 0 : libmesh_error_msg("ERROR: Bad FEType.radial_family= " << Utility::enum_to_string(fet.radial_family));
622 : }
623 : }
624 :
625 :
626 :
627 :
628 : // 2D
629 0 : case 2:
630 : {
631 0 : switch (fet.radial_family)
632 : {
633 0 : case INFINITE_MAP:
634 0 : libmesh_error_msg("ERROR: Can't build an infinite element with FEFamily = " << Utility::enum_to_string(fet.radial_family));
635 :
636 0 : case JACOBI_20_00:
637 : {
638 0 : switch (fet.inf_map)
639 : {
640 0 : case CARTESIAN:
641 0 : return std::make_unique<InfFE<2,JACOBI_20_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 JACOBI_30_00:
649 : {
650 0 : switch (fet.inf_map)
651 : {
652 0 : case CARTESIAN:
653 0 : return std::make_unique<InfFE<2,JACOBI_30_00,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 LEGENDRE:
661 : {
662 0 : switch (fet.inf_map)
663 : {
664 0 : case CARTESIAN:
665 0 : return std::make_unique<InfFE<2,LEGENDRE,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 : case LAGRANGE:
673 : {
674 0 : switch (fet.inf_map)
675 : {
676 0 : case CARTESIAN:
677 0 : return std::make_unique<InfFE<2,LAGRANGE,CARTESIAN>>(fet);
678 :
679 0 : default:
680 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
681 : }
682 : }
683 :
684 0 : default:
685 0 : libmesh_error_msg("ERROR: Bad FEType.radial_family= " << Utility::enum_to_string(fet.radial_family));
686 : }
687 : }
688 :
689 :
690 :
691 :
692 : // 3D
693 452 : case 3:
694 : {
695 452 : switch (fet.radial_family)
696 : {
697 0 : case INFINITE_MAP:
698 0 : libmesh_error_msg("ERROR: Don't build an infinite element with FEFamily = " << Utility::enum_to_string(fet.radial_family));
699 :
700 397 : case JACOBI_20_00:
701 : {
702 397 : switch (fet.inf_map)
703 : {
704 397 : case CARTESIAN:
705 397 : return std::make_unique<InfFE<3,JACOBI_20_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 JACOBI_30_00:
713 : {
714 20 : switch (fet.inf_map)
715 : {
716 20 : case CARTESIAN:
717 20 : return std::make_unique<InfFE<3,JACOBI_30_00,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 20 : case LEGENDRE:
725 : {
726 20 : switch (fet.inf_map)
727 : {
728 20 : case CARTESIAN:
729 20 : return std::make_unique<InfFE<3,LEGENDRE,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 15 : case LAGRANGE:
737 : {
738 15 : switch (fet.inf_map)
739 : {
740 15 : case CARTESIAN:
741 15 : return std::make_unique<InfFE<3,LAGRANGE,CARTESIAN>>(fet);
742 :
743 0 : default:
744 0 : libmesh_error_msg("ERROR: Don't build an infinite element with InfMapType = " << Utility::enum_to_string(fet.inf_map));
745 : }
746 : }
747 :
748 0 : default:
749 0 : libmesh_error_msg("ERROR: Bad FEType.radial_family= " << Utility::enum_to_string(fet.radial_family));
750 : }
751 : }
752 :
753 0 : default:
754 0 : libmesh_error_msg("Invalid dimension dim = " << dim);
755 : }
756 : }
757 :
758 :
759 :
760 : template <>
761 : std::unique_ptr<FEGenericBase<RealGradient>>
762 0 : FEGenericBase<RealGradient>::build_InfFE (const unsigned int,
763 : const FEType & )
764 : {
765 : // No vector types defined... YET.
766 0 : libmesh_not_implemented();
767 : return std::unique_ptr<FEVectorBase>();
768 : }
769 :
770 : #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
771 :
772 :
773 : template <typename OutputType>
774 124412404 : void FEGenericBase<OutputType>::compute_shape_functions (const Elem * elem,
775 : const std::vector<Point> & qp)
776 : {
777 : //-------------------------------------------------------------------------
778 : // Compute the shape function values (and derivatives)
779 : // at the Quadrature points. Note that the actual values
780 : // have already been computed via init_shape_functions
781 :
782 : // Start logging the shape function computation
783 33340651 : LOG_ONE_ENTITY_SCOPE("compute_shape_functions()", "FE");
784 :
785 124412404 : this->determine_calculations();
786 :
787 124412404 : if (calculate_phi)
788 112825371 : this->_fe_trans->map_phi(this->dim, elem, qp, (*this), this->phi, this->_add_p_level_in_reinit);
789 :
790 124412404 : if (calculate_dphi)
791 90016541 : this->_fe_trans->map_dphi(this->dim, elem, qp, (*this), this->dphi,
792 82459487 : this->dphidx, this->dphidy, this->dphidz);
793 :
794 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
795 124412404 : if (calculate_d2phi)
796 8480328 : this->_fe_trans->map_d2phi(this->dim, qp, (*this), this->d2phi,
797 7861374 : this->d2phidx2, this->d2phidxdy, this->d2phidxdz,
798 7861374 : this->d2phidy2, this->d2phidydz, this->d2phidz2);
799 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
800 :
801 : // Only compute curl for vector-valued elements
802 9954615 : if (calculate_curl_phi && TypesEqual<OutputType,RealGradient>::value)
803 859569 : this->_fe_trans->map_curl(this->dim, elem, qp, (*this), this->curl_phi);
804 :
805 : // Only compute div for vector-valued elements
806 9954615 : if (calculate_div_phi && TypesEqual<OutputType,RealGradient>::value)
807 1476125 : this->_fe_trans->map_div(this->dim, elem, qp, (*this), this->div_phi);
808 124412404 : }
809 :
810 :
811 : // Here, we rely on the input \p phi_vals for accurate integration of the mass matrix.
812 : // This is because in contact, we often have customized qrule for mortar segments
813 : // due to deformation of the element, and the size of \p phi_vals for the secondary
814 : // element changes accordingly.
815 : template <>
816 6974 : void FEGenericBase<Real>::compute_dual_shape_coeffs (const std::vector<Real> & JxW, const std::vector<std::vector<OutputShape>> & phi_vals)
817 : {
818 : // Start logging the dual coeff computation
819 1428 : LOG_SCOPE("compute_dual_shape_coeffs()", "FE");
820 :
821 6974 : const unsigned int sz=phi_vals.size();
822 6974 : libmesh_error_msg_if(!sz, "ERROR: cannot compute dual shape coefficients with empty phi values");
823 :
824 : //compute dual basis coefficient (dual_coeff)
825 6498 : dual_coeff.resize(sz, sz);
826 7926 : DenseMatrix<Real> A(sz, sz), D(sz, sz);
827 :
828 120780 : for (const auto i : index_range(phi_vals))
829 6854504 : for (const auto qp : index_range(phi_vals[i]))
830 : {
831 8194957 : D(i,i) += JxW[qp]*phi_vals[i][qp];
832 388594300 : for (const auto j : index_range(phi_vals))
833 470334381 : A(i,j) += JxW[qp]*phi_vals[i][qp]*phi_vals[j][qp];
834 : }
835 :
836 : // The integral positivity condition on the dual basis (Popp, Wohlmuth, Gee and Wall, SIAM J. Sci.
837 : // 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
838 : // QUAD8 corner and exactly 0 at a TRI6 vertex (their Eqs. (4.3) and (4.4)), the latter leaving those
839 : // dual shape functions identically zero. Biorthogonalize instead against the locally quadratic
840 : // transformed basis Ntilde = T N of their Sec. 4.4.1, in which each vertex absorbs a fraction alpha
841 : // of its adjacent mid-edge shapes. That replaces D by T^-1 diag(T d), which for this T is the sparse
842 : // update below. alpha = 1/5 is recommended there, and makes the weights strictly positive (QUAD8
843 : // 1/5 and 4/5, TRI6 1/15 and 1/10) while preserving the partition of unity. T = I elsewhere.
844 : // SECOND order on TRI7 is the TRI6 basis, with no degree of freedom at the centroid node, and
845 : // L2_LAGRANGE shares the LAGRANGE shape functions.
846 8688 : if (_elem && (get_family() == LAGRANGE || get_family() == L2_LAGRANGE) &&
847 4762 : (((_elem->type() == TRI6 || _elem->type() == TRI7) && sz == 6) ||
848 1571 : (_elem->type() == QUAD8 && sz == 8)))
849 : {
850 19 : const Real alpha = Real(1)/5;
851 : // Mid-edge nodes are the trailing indices, and only vertex entries are written, so D(m,m) here
852 : // is never a value an earlier iteration modified.
853 1058 : for (const auto m : make_range(_elem->n_vertices(), sz))
854 2568 : for (const auto v : make_range(_elem->n_second_order_adjacent_vertices(m)))
855 : {
856 1344 : const auto vertex = _elem->second_order_adjacent_vertex(m, v);
857 1344 : D(vertex, vertex) += alpha*D(m,m);
858 1464 : D(vertex, m) -= alpha*D(m,m);
859 : }
860 : }
861 :
862 : // dual_coeff = A^-1*D
863 120780 : for (const auto j : index_range(phi_vals))
864 : {
865 113806 : DenseVector<Real> Dcol(sz), coeffcol(sz);
866 4306932 : for (const auto i : index_range(phi_vals))
867 4822538 : Dcol(i) = D(i, j);
868 113806 : A.cholesky_solve(Dcol, coeffcol);
869 :
870 4306932 : for (const auto row : index_range(phi_vals))
871 4822538 : dual_coeff(row, j)=coeffcol(row);
872 : }
873 6974 : }
874 :
875 : template <>
876 6974 : void FEGenericBase<Real>::compute_dual_shape_functions ()
877 : {
878 : // Start logging the shape function computation
879 1428 : LOG_SCOPE("compute_dual_shape_functions()", "FE");
880 :
881 : // The dual coeffs matrix should have the same size as phi
882 476 : libmesh_assert(dual_coeff.m() == phi.size());
883 476 : libmesh_assert(dual_coeff.n() == phi.size());
884 :
885 : // initialize dual basis
886 120780 : for (const auto j : index_range(phi))
887 6854504 : for (const auto qp : index_range(phi[j]))
888 : {
889 7225451 : dual_phi[j][qp] = 0;
890 6740698 : if (calculate_dphi)
891 1452807 : dual_dphi[j][qp] = 0;
892 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
893 6740698 : if (calculate_d2phi)
894 1444275 : dual_d2phi[j][qp] = 0;
895 : #endif
896 : }
897 :
898 : // compute dual basis
899 120780 : for (const auto j : index_range(phi))
900 4306932 : for (const auto i : index_range(phi))
901 386046728 : for (const auto qp : index_range(phi[j]))
902 : {
903 499827974 : dual_phi[j][qp] += dual_coeff(i, j) * phi[i][qp];
904 381853602 : if (calculate_dphi)
905 176942502 : dual_dphi[j][qp] += dual_coeff(i, j) * dphi[i][qp];
906 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
907 381853602 : if (calculate_d2phi)
908 850971837 : dual_d2phi[j][qp] += dual_coeff(i, j) * d2phi[i][qp];
909 : #endif
910 : }
911 6974 : }
912 :
913 : template <typename OutputType>
914 0 : void FEGenericBase<OutputType>::print_phi(std::ostream & os) const
915 : {
916 0 : for (auto i : index_range(phi))
917 0 : for (auto j : index_range(phi[i]))
918 0 : os << " phi[" << i << "][" << j << "]=" << phi[i][j] << std::endl;
919 0 : }
920 :
921 : template <typename OutputType>
922 0 : void FEGenericBase<OutputType>::print_dual_phi(std::ostream & os) const
923 : {
924 0 : for (auto i : index_range(dual_phi))
925 0 : for (auto j : index_range(dual_phi[i]))
926 0 : os << " dual_phi[" << i << "][" << j << "]=" << dual_phi[i][j] << std::endl;
927 0 : }
928 :
929 :
930 :
931 :
932 : template <typename OutputType>
933 0 : void FEGenericBase<OutputType>::print_dphi(std::ostream & os) const
934 : {
935 0 : for (auto i : index_range(dphi))
936 0 : for (auto j : index_range(dphi[i]))
937 0 : os << " dphi[" << i << "][" << j << "]=" << dphi[i][j];
938 0 : }
939 :
940 : template <typename OutputType>
941 0 : void FEGenericBase<OutputType>::print_dual_dphi(std::ostream & os) const
942 : {
943 0 : for (auto i : index_range(dphi))
944 0 : for (auto j : index_range(dphi[i]))
945 0 : os << " dual_dphi[" << i << "][" << j << "]=" << dual_dphi[i][j];
946 0 : }
947 :
948 :
949 :
950 : template <typename OutputType>
951 347557510 : void FEGenericBase<OutputType>::determine_calculations()
952 : {
953 347557510 : this->calculations_started = true;
954 :
955 : // If the user did not explicitly pre-request something (or nothing)
956 : // to be computed, then we throw an error here.
957 30154093 : bool requested_ok =
958 338095617 : this->calculate_nothing || this->calculate_phi || this->calculate_dphi ||
959 379721291 : this->calculate_dphiref || this->calculate_curl_phi || this->calculate_div_phi ||
960 640437 : this->calculate_map;
961 :
962 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
963 30154093 : requested_ok = requested_ok || this->calculate_d2phi;
964 : #endif
965 :
966 30154093 : libmesh_error_msg_if(
967 : !requested_ok,
968 : "You must call one or more of the FE accessors "
969 : "(e.g. get_phi(), get_dphi(), get_nothing()) "
970 : "_before_ calling reinit()!");
971 :
972 : // Request whichever terms are necessary from the FEMap
973 347557510 : if (this->calculate_phi)
974 323541641 : this->_fe_trans->init_map_phi(*this);
975 :
976 347557510 : if (this->calculate_dphiref)
977 228463243 : this->_fe_trans->init_map_dphi(*this);
978 :
979 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
980 347557510 : if (this->calculate_d2phi)
981 61087551 : this->_fe_trans->init_map_d2phi(*this);
982 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
983 347557510 : }
984 :
985 :
986 :
987 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
988 :
989 :
990 : template <typename OutputType>
991 0 : void FEGenericBase<OutputType>::print_d2phi(std::ostream & os) const
992 : {
993 0 : for (auto i : index_range(dphi))
994 0 : for (auto j : index_range(dphi[i]))
995 0 : os << " d2phi[" << i << "][" << j << "]=" << d2phi[i][j];
996 0 : }
997 :
998 : template <typename OutputType>
999 0 : void FEGenericBase<OutputType>::print_dual_d2phi(std::ostream & os) const
1000 : {
1001 0 : for (auto i : index_range(dual_d2phi))
1002 0 : for (auto j : index_range(dual_d2phi[i]))
1003 0 : os << " dual_d2phi[" << i << "][" << j << "]=" << dual_d2phi[i][j];
1004 0 : }
1005 :
1006 : #endif
1007 :
1008 :
1009 :
1010 : #ifdef LIBMESH_ENABLE_AMR
1011 :
1012 : template <typename OutputType>
1013 : void
1014 0 : FEGenericBase<OutputType>::coarsened_dof_values(const NumericVector<Number> & old_vector,
1015 : const DofMap & dof_map,
1016 : const Elem * elem,
1017 : DenseVector<Number> & Ue,
1018 : const unsigned int var,
1019 : const bool use_old_dof_indices)
1020 : {
1021 : // Side/edge local DOF indices
1022 0 : std::vector<unsigned int> new_side_dofs, old_side_dofs;
1023 :
1024 : // FIXME: what about 2D shells in 3D space?
1025 0 : unsigned int dim = elem->dim();
1026 :
1027 : // Cache n_children(); it's a virtual call but it's const.
1028 0 : const unsigned int n_children = elem->n_children();
1029 :
1030 : // We use local FE objects for now
1031 : // FIXME: we should use more, external objects instead for efficiency
1032 0 : const FEType & base_fe_type = dof_map.variable_type(var);
1033 0 : std::unique_ptr<FEGenericBase<OutputShape>> fe
1034 : (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1035 0 : std::unique_ptr<FEGenericBase<OutputShape>> fe_coarse
1036 : (FEGenericBase<OutputShape>::build(dim, base_fe_type));
1037 :
1038 0 : std::unique_ptr<QBase> qrule (base_fe_type.default_quadrature_rule(dim));
1039 0 : std::unique_ptr<QBase> qedgerule (base_fe_type.default_quadrature_rule(1));
1040 0 : std::unique_ptr<QBase> qsiderule (base_fe_type.default_quadrature_rule(dim-1));
1041 0 : std::vector<Point> coarse_qpoints;
1042 :
1043 : // The values of the shape functions at the quadrature
1044 : // points
1045 0 : const std::vector<std::vector<OutputShape>> & phi_values =
1046 : fe->get_phi();
1047 0 : const std::vector<std::vector<OutputShape>> & phi_coarse =
1048 : fe_coarse->get_phi();
1049 :
1050 : // The gradients of the shape functions at the quadrature
1051 : // points on the child element.
1052 0 : const std::vector<std::vector<OutputGradient>> * dphi_values =
1053 : nullptr;
1054 0 : const std::vector<std::vector<OutputGradient>> * dphi_coarse =
1055 : nullptr;
1056 :
1057 0 : const FEContinuity cont = fe->get_continuity();
1058 :
1059 0 : if (cont == C_ONE)
1060 : {
1061 : const std::vector<std::vector<OutputGradient>> &
1062 0 : ref_dphi_values = fe->get_dphi();
1063 0 : dphi_values = &ref_dphi_values;
1064 : const std::vector<std::vector<OutputGradient>> &
1065 0 : ref_dphi_coarse = fe_coarse->get_dphi();
1066 0 : dphi_coarse = &ref_dphi_coarse;
1067 : }
1068 :
1069 : // The Jacobian * quadrature weight at the quadrature points
1070 0 : const std::vector<Real> & JxW =
1071 0 : fe->get_JxW();
1072 :
1073 : // The XYZ locations of the quadrature points on the
1074 : // child element
1075 0 : const std::vector<Point> & xyz_values =
1076 0 : fe->get_xyz();
1077 :
1078 : // Number of nodes on parent element
1079 0 : const unsigned int n_nodes = elem->n_nodes();
1080 :
1081 : // Number of dofs on parent element
1082 : const unsigned int new_n_dofs =
1083 0 : FEInterface::n_dofs(base_fe_type, elem->max_descendant_p_level(), elem);
1084 :
1085 : // Fixed vs. free DoFs on edge/face projections
1086 0 : std::vector<char> dof_is_fixed(new_n_dofs, false); // bools
1087 0 : std::vector<int> free_dof(new_n_dofs, 0);
1088 :
1089 0 : DenseMatrix<Real> Ke;
1090 0 : DenseVector<Number> Fe;
1091 0 : Ue.resize(new_n_dofs); Ue.zero();
1092 :
1093 :
1094 : // When coarsening, in general, we need a series of
1095 : // projections to ensure a unique and continuous
1096 : // solution. We start by interpolating nodes, then
1097 : // hold those fixed and project edges, then
1098 : // hold those fixed and project faces, then
1099 : // hold those fixed and project interiors
1100 :
1101 : // Copy node values first
1102 : {
1103 0 : std::vector<dof_id_type> node_dof_indices;
1104 0 : if (use_old_dof_indices)
1105 0 : dof_map.old_dof_indices (elem, node_dof_indices, var);
1106 : else
1107 0 : dof_map.dof_indices (elem, node_dof_indices, var);
1108 :
1109 0 : unsigned int current_dof = 0;
1110 0 : for (unsigned int n=0; n!= n_nodes; ++n)
1111 : {
1112 : // FIXME: this should go through the DofMap,
1113 : // not duplicate dof_indices code badly!
1114 : const unsigned int my_nc =
1115 0 : FEInterface::n_dofs_at_node (base_fe_type, elem->max_descendant_p_level(), elem, n);
1116 0 : if (!elem->is_vertex(n))
1117 : {
1118 0 : current_dof += my_nc;
1119 0 : continue;
1120 : }
1121 :
1122 : // We're assuming here that child n shares vertex n,
1123 : // which is wrong on non-simplices right now
1124 : // ... but this code isn't necessary except on elements
1125 : // where p refinement creates more vertex dofs; we have
1126 : // no such elements yet.
1127 0 : int extra_order = 0;
1128 : // if (elem->child_ptr(n)->p_level() < elem->p_level())
1129 : // extra_order = elem->child_ptr(n)->p_level();
1130 : const unsigned int nc =
1131 0 : FEInterface::n_dofs_at_node (base_fe_type, extra_order, elem, n);
1132 0 : for (unsigned int i=0; i!= nc; ++i)
1133 : {
1134 0 : Ue(current_dof) =
1135 0 : old_vector(node_dof_indices[current_dof]);
1136 0 : dof_is_fixed[current_dof] = true;
1137 0 : current_dof++;
1138 : }
1139 : }
1140 : }
1141 :
1142 0 : FEType fe_type = base_fe_type, temp_fe_type;
1143 0 : fe_type.order = fe_type.order + elem->max_descendant_p_level();
1144 :
1145 : // In 3D, project any edge values next
1146 0 : if (dim > 2 && cont != DISCONTINUOUS)
1147 0 : for (auto e : elem->edge_index_range())
1148 : {
1149 0 : FEInterface::dofs_on_edge(elem, dim, fe_type,
1150 : e, new_side_dofs);
1151 :
1152 : const unsigned int n_new_side_dofs =
1153 0 : cast_int<unsigned int>(new_side_dofs.size());
1154 :
1155 : // Some edge dofs are on nodes and already
1156 : // fixed, others are free to calculate
1157 0 : unsigned int free_dofs = 0;
1158 0 : for (unsigned int i=0; i != n_new_side_dofs; ++i)
1159 0 : if (!dof_is_fixed[new_side_dofs[i]])
1160 0 : free_dof[free_dofs++] = i;
1161 0 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1162 0 : Fe.resize (free_dofs); Fe.zero();
1163 : // The new edge coefficients
1164 0 : DenseVector<Number> Uedge(free_dofs);
1165 :
1166 : // Add projection terms from each child sharing
1167 : // this edge
1168 0 : for (unsigned int c=0; c != n_children; ++c)
1169 : {
1170 0 : if (!elem->is_child_on_edge(c,e))
1171 0 : continue;
1172 0 : const Elem * child = elem->child_ptr(c);
1173 :
1174 0 : std::vector<dof_id_type> child_dof_indices;
1175 0 : if (use_old_dof_indices)
1176 0 : dof_map.old_dof_indices (child,
1177 : child_dof_indices, var);
1178 : else
1179 0 : dof_map.dof_indices (child,
1180 : child_dof_indices, var);
1181 : const unsigned int child_n_dofs =
1182 : cast_int<unsigned int>
1183 0 : (child_dof_indices.size());
1184 :
1185 0 : temp_fe_type = base_fe_type;
1186 0 : temp_fe_type.order = temp_fe_type.order + child->p_level();
1187 :
1188 0 : FEInterface::dofs_on_edge(child, dim,
1189 : temp_fe_type, e, old_side_dofs);
1190 :
1191 : // Initialize both child and parent FE data
1192 : // on the child's edge
1193 0 : fe->attach_quadrature_rule (qedgerule.get());
1194 0 : fe->edge_reinit (child, e);
1195 0 : const unsigned int n_qp = qedgerule->n_points();
1196 :
1197 0 : FEMap::inverse_map (dim, elem, xyz_values,
1198 : coarse_qpoints);
1199 :
1200 0 : fe_coarse->reinit(elem, &coarse_qpoints);
1201 :
1202 : // Loop over the quadrature points
1203 0 : for (unsigned int qp=0; qp<n_qp; qp++)
1204 : {
1205 : // solution value at the quadrature point
1206 0 : OutputNumber fineval = libMesh::zero;
1207 : // solution grad at the quadrature point
1208 0 : OutputNumberGradient finegrad;
1209 :
1210 : // Sum the solution values * the DOF
1211 : // values at the quadrature point to
1212 : // get the solution value and gradient.
1213 0 : for (unsigned int i=0; i<child_n_dofs;
1214 : i++)
1215 : {
1216 0 : fineval +=
1217 0 : (old_vector(child_dof_indices[i])*
1218 0 : phi_values[i][qp]);
1219 0 : if (cont == C_ONE)
1220 0 : finegrad += (*dphi_values)[i][qp] *
1221 0 : old_vector(child_dof_indices[i]);
1222 : }
1223 :
1224 : // Form edge projection matrix
1225 0 : for (unsigned int sidei=0, freei=0; sidei != n_new_side_dofs; ++sidei)
1226 : {
1227 0 : unsigned int i = new_side_dofs[sidei];
1228 : // fixed DoFs aren't test functions
1229 0 : if (dof_is_fixed[i])
1230 0 : continue;
1231 0 : for (unsigned int sidej=0, freej=0; sidej != n_new_side_dofs; ++sidej)
1232 : {
1233 0 : unsigned int j =
1234 : new_side_dofs[sidej];
1235 0 : if (dof_is_fixed[j])
1236 0 : Fe(freei) -=
1237 0 : TensorTools::inner_product(phi_coarse[i][qp],
1238 0 : phi_coarse[j][qp]) *
1239 0 : JxW[qp] * Ue(j);
1240 : else
1241 0 : Ke(freei,freej) +=
1242 0 : TensorTools::inner_product(phi_coarse[i][qp],
1243 0 : phi_coarse[j][qp]) *
1244 : JxW[qp];
1245 0 : if (cont == C_ONE)
1246 : {
1247 0 : if (dof_is_fixed[j])
1248 0 : Fe(freei) -=
1249 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1250 0 : (*dphi_coarse)[j][qp]) *
1251 0 : JxW[qp] * Ue(j);
1252 : else
1253 0 : Ke(freei,freej) +=
1254 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1255 0 : (*dphi_coarse)[j][qp]) *
1256 : JxW[qp];
1257 : }
1258 0 : if (!dof_is_fixed[j])
1259 0 : freej++;
1260 : }
1261 0 : Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp],
1262 0 : fineval) * JxW[qp];
1263 0 : if (cont == C_ONE)
1264 0 : Fe(freei) +=
1265 0 : TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1266 0 : freei++;
1267 : }
1268 : }
1269 : }
1270 0 : Ke.cholesky_solve(Fe, Uedge);
1271 :
1272 : // Transfer new edge solutions to element
1273 0 : for (unsigned int i=0; i != free_dofs; ++i)
1274 : {
1275 0 : Number & ui = Ue(new_side_dofs[free_dof[i]]);
1276 0 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1277 : std::abs(ui - Uedge(i)) < TOLERANCE);
1278 0 : ui = Uedge(i);
1279 0 : dof_is_fixed[new_side_dofs[free_dof[i]]] = true;
1280 : }
1281 : }
1282 :
1283 : // Project any side values (edges in 2D, faces in 3D)
1284 0 : if (dim > 1 && cont != DISCONTINUOUS)
1285 0 : for (auto s : elem->side_index_range())
1286 : {
1287 0 : FEInterface::dofs_on_side(elem, dim, fe_type,
1288 : s, new_side_dofs);
1289 :
1290 : const unsigned int n_new_side_dofs =
1291 0 : cast_int<unsigned int>(new_side_dofs.size());
1292 :
1293 : // Some side dofs are on nodes/edges and already
1294 : // fixed, others are free to calculate
1295 0 : unsigned int free_dofs = 0;
1296 0 : for (unsigned int i=0; i != n_new_side_dofs; ++i)
1297 0 : if (!dof_is_fixed[new_side_dofs[i]])
1298 0 : free_dof[free_dofs++] = i;
1299 0 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1300 0 : Fe.resize (free_dofs); Fe.zero();
1301 : // The new side coefficients
1302 0 : DenseVector<Number> Uside(free_dofs);
1303 :
1304 : // Add projection terms from each child sharing
1305 : // this side
1306 0 : for (unsigned int c=0; c != n_children; ++c)
1307 : {
1308 0 : if (!elem->is_child_on_side(c,s))
1309 0 : continue;
1310 0 : const Elem * child = elem->child_ptr(c);
1311 :
1312 0 : std::vector<dof_id_type> child_dof_indices;
1313 0 : if (use_old_dof_indices)
1314 0 : dof_map.old_dof_indices (child,
1315 : child_dof_indices, var);
1316 : else
1317 0 : dof_map.dof_indices (child,
1318 : child_dof_indices, var);
1319 : const unsigned int child_n_dofs =
1320 : cast_int<unsigned int>
1321 0 : (child_dof_indices.size());
1322 :
1323 0 : temp_fe_type = base_fe_type;
1324 0 : temp_fe_type.order = temp_fe_type.order + child->p_level();
1325 :
1326 0 : FEInterface::dofs_on_side(child, dim,
1327 : temp_fe_type, s, old_side_dofs);
1328 :
1329 : // Initialize both child and parent FE data
1330 : // on the child's side
1331 0 : fe->attach_quadrature_rule (qsiderule.get());
1332 0 : fe->reinit (child, s);
1333 0 : const unsigned int n_qp = qsiderule->n_points();
1334 :
1335 0 : FEMap::inverse_map (dim, elem, xyz_values,
1336 : coarse_qpoints);
1337 :
1338 0 : fe_coarse->reinit(elem, &coarse_qpoints);
1339 :
1340 : // Loop over the quadrature points
1341 0 : for (unsigned int qp=0; qp<n_qp; qp++)
1342 : {
1343 : // solution value at the quadrature point
1344 0 : OutputNumber fineval = libMesh::zero;
1345 : // solution grad at the quadrature point
1346 0 : OutputNumberGradient finegrad;
1347 :
1348 : // Sum the solution values * the DOF
1349 : // values at the quadrature point to
1350 : // get the solution value and gradient.
1351 0 : for (unsigned int i=0; i<child_n_dofs;
1352 : i++)
1353 : {
1354 0 : fineval +=
1355 0 : old_vector(child_dof_indices[i]) *
1356 0 : phi_values[i][qp];
1357 0 : if (cont == C_ONE)
1358 0 : finegrad += (*dphi_values)[i][qp] *
1359 0 : old_vector(child_dof_indices[i]);
1360 : }
1361 :
1362 : // Form side projection matrix
1363 0 : for (unsigned int sidei=0, freei=0; sidei != n_new_side_dofs; ++sidei)
1364 : {
1365 0 : unsigned int i = new_side_dofs[sidei];
1366 : // fixed DoFs aren't test functions
1367 0 : if (dof_is_fixed[i])
1368 0 : continue;
1369 0 : for (unsigned int sidej=0, freej=0; sidej != n_new_side_dofs; ++sidej)
1370 : {
1371 0 : unsigned int j =
1372 : new_side_dofs[sidej];
1373 0 : if (dof_is_fixed[j])
1374 0 : Fe(freei) -=
1375 0 : TensorTools::inner_product(phi_coarse[i][qp],
1376 0 : phi_coarse[j][qp]) *
1377 0 : JxW[qp] * Ue(j);
1378 : else
1379 0 : Ke(freei,freej) +=
1380 0 : TensorTools::inner_product(phi_coarse[i][qp],
1381 0 : phi_coarse[j][qp]) *
1382 : JxW[qp];
1383 0 : if (cont == C_ONE)
1384 : {
1385 0 : if (dof_is_fixed[j])
1386 0 : Fe(freei) -=
1387 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1388 0 : (*dphi_coarse)[j][qp]) *
1389 0 : JxW[qp] * Ue(j);
1390 : else
1391 0 : Ke(freei,freej) +=
1392 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1393 0 : (*dphi_coarse)[j][qp]) *
1394 : JxW[qp];
1395 : }
1396 0 : if (!dof_is_fixed[j])
1397 0 : freej++;
1398 : }
1399 0 : Fe(freei) += TensorTools::inner_product(fineval, phi_coarse[i][qp]) * JxW[qp];
1400 0 : if (cont == C_ONE)
1401 0 : Fe(freei) +=
1402 0 : TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1403 0 : freei++;
1404 : }
1405 : }
1406 : }
1407 0 : Ke.cholesky_solve(Fe, Uside);
1408 :
1409 : // Transfer new side solutions to element
1410 0 : for (unsigned int i=0; i != free_dofs; ++i)
1411 : {
1412 0 : Number & ui = Ue(new_side_dofs[free_dof[i]]);
1413 0 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1414 : std::abs(ui - Uside(i)) < TOLERANCE);
1415 0 : ui = Uside(i);
1416 0 : dof_is_fixed[new_side_dofs[free_dof[i]]] = true;
1417 : }
1418 : }
1419 :
1420 : // Project the interior values, finally
1421 :
1422 : // Some interior dofs are on nodes/edges/sides and
1423 : // already fixed, others are free to calculate
1424 0 : unsigned int free_dofs = 0;
1425 0 : for (unsigned int i=0; i != new_n_dofs; ++i)
1426 0 : if (!dof_is_fixed[i])
1427 0 : free_dof[free_dofs++] = i;
1428 0 : Ke.resize (free_dofs, free_dofs); Ke.zero();
1429 0 : Fe.resize (free_dofs); Fe.zero();
1430 : // The new interior coefficients
1431 0 : DenseVector<Number> Uint(free_dofs);
1432 :
1433 : // Add projection terms from each child
1434 0 : for (auto & child : elem->child_ref_range())
1435 : {
1436 0 : std::vector<dof_id_type> child_dof_indices;
1437 0 : if (use_old_dof_indices)
1438 0 : dof_map.old_dof_indices (&child,
1439 : child_dof_indices, var);
1440 : else
1441 0 : dof_map.dof_indices (&child,
1442 : child_dof_indices, var);
1443 : const unsigned int child_n_dofs =
1444 : cast_int<unsigned int>
1445 0 : (child_dof_indices.size());
1446 :
1447 : // Initialize both child and parent FE data
1448 : // on the child's quadrature points
1449 0 : fe->attach_quadrature_rule (qrule.get());
1450 0 : fe->reinit (&child);
1451 0 : const unsigned int n_qp = qrule->n_points();
1452 :
1453 0 : FEMap::inverse_map (dim, elem, xyz_values, coarse_qpoints);
1454 :
1455 0 : fe_coarse->reinit(elem, &coarse_qpoints);
1456 :
1457 : // Loop over the quadrature points
1458 0 : for (unsigned int qp=0; qp<n_qp; qp++)
1459 : {
1460 : // solution value at the quadrature point
1461 0 : OutputNumber fineval = libMesh::zero;
1462 : // solution grad at the quadrature point
1463 0 : OutputNumberGradient finegrad;
1464 :
1465 : // Sum the solution values * the DOF
1466 : // values at the quadrature point to
1467 : // get the solution value and gradient.
1468 0 : for (unsigned int i=0; i<child_n_dofs; i++)
1469 : {
1470 0 : fineval +=
1471 0 : (old_vector(child_dof_indices[i]) *
1472 0 : phi_values[i][qp]);
1473 0 : if (cont == C_ONE)
1474 0 : finegrad += (*dphi_values)[i][qp] *
1475 0 : old_vector(child_dof_indices[i]);
1476 : }
1477 :
1478 : // Form interior projection matrix
1479 0 : for (unsigned int i=0, freei=0;
1480 0 : i != new_n_dofs; ++i)
1481 : {
1482 : // fixed DoFs aren't test functions
1483 0 : if (dof_is_fixed[i])
1484 0 : continue;
1485 0 : for (unsigned int j=0, freej=0; j !=
1486 : new_n_dofs; ++j)
1487 : {
1488 0 : if (dof_is_fixed[j])
1489 0 : Fe(freei) -=
1490 0 : TensorTools::inner_product(phi_coarse[i][qp],
1491 0 : phi_coarse[j][qp]) *
1492 0 : JxW[qp] * Ue(j);
1493 : else
1494 0 : Ke(freei,freej) +=
1495 0 : TensorTools::inner_product(phi_coarse[i][qp],
1496 0 : phi_coarse[j][qp]) *
1497 : JxW[qp];
1498 0 : if (cont == C_ONE)
1499 : {
1500 0 : if (dof_is_fixed[j])
1501 0 : Fe(freei) -=
1502 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1503 0 : (*dphi_coarse)[j][qp]) *
1504 0 : JxW[qp] * Ue(j);
1505 : else
1506 0 : Ke(freei,freej) +=
1507 0 : TensorTools::inner_product((*dphi_coarse)[i][qp],
1508 0 : (*dphi_coarse)[j][qp]) *
1509 : JxW[qp];
1510 : }
1511 0 : if (!dof_is_fixed[j])
1512 0 : freej++;
1513 : }
1514 0 : Fe(freei) += TensorTools::inner_product(phi_coarse[i][qp], fineval) *
1515 : JxW[qp];
1516 0 : if (cont == C_ONE)
1517 0 : Fe(freei) += TensorTools::inner_product(finegrad, (*dphi_coarse)[i][qp]) * JxW[qp];
1518 0 : freei++;
1519 : }
1520 : }
1521 : }
1522 0 : Ke.cholesky_solve(Fe, Uint);
1523 :
1524 : // Transfer new interior solutions to element
1525 0 : for (unsigned int i=0; i != free_dofs; ++i)
1526 : {
1527 0 : Number & ui = Ue(free_dof[i]);
1528 0 : libmesh_assert(std::abs(ui) < TOLERANCE ||
1529 : std::abs(ui - Uint(i)) < TOLERANCE);
1530 0 : ui = Uint(i);
1531 : // We should be fixing all dofs by now; no need to keep track of
1532 : // that unless we're debugging
1533 : #ifndef NDEBUG
1534 0 : dof_is_fixed[free_dof[i]] = true;
1535 : #endif
1536 : }
1537 :
1538 : #ifndef NDEBUG
1539 : // Make sure every DoF got reached!
1540 0 : for (unsigned int i=0; i != new_n_dofs; ++i)
1541 0 : libmesh_assert(dof_is_fixed[i]);
1542 : #endif
1543 0 : }
1544 :
1545 :
1546 :
1547 : template <typename OutputType>
1548 : void
1549 0 : FEGenericBase<OutputType>::coarsened_dof_values(const NumericVector<Number> & old_vector,
1550 : const DofMap & dof_map,
1551 : const Elem * elem,
1552 : DenseVector<Number> & Ue,
1553 : const bool use_old_dof_indices)
1554 : {
1555 0 : Ue.resize(0);
1556 :
1557 0 : for (auto v : make_range(dof_map.n_variables()))
1558 : {
1559 0 : DenseVector<Number> Usub;
1560 :
1561 0 : coarsened_dof_values(old_vector, dof_map, elem, Usub,
1562 : v, use_old_dof_indices);
1563 :
1564 0 : Ue.append (Usub);
1565 : }
1566 0 : }
1567 :
1568 :
1569 :
1570 : template <typename OutputType>
1571 : void
1572 691258 : FEGenericBase<OutputType>::compute_proj_constraints (DofConstraints & constraints,
1573 : DofMap & dof_map,
1574 : const unsigned int variable_number,
1575 : const Elem * elem)
1576 : {
1577 58078 : libmesh_assert(elem);
1578 :
1579 691258 : const unsigned int Dim = elem->dim();
1580 :
1581 : // Only constrain elements in 2,3D.
1582 691258 : if (Dim == 1)
1583 108035 : return;
1584 :
1585 : // Only constrain active elements with this method
1586 58078 : if (!elem->active())
1587 9226 : return;
1588 :
1589 583223 : const Variable & var = dof_map.variable(variable_number);
1590 48852 : const FEType & base_fe_type = var.type();
1591 583223 : const bool add_p_level = base_fe_type.p_refinement;
1592 :
1593 : // Construct FE objects for this element and its neighbors.
1594 583223 : std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1595 : (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1596 48852 : my_fe->add_p_level_in_reinit(add_p_level);
1597 583223 : const FEContinuity cont = my_fe->get_continuity();
1598 :
1599 : // We don't need to constrain discontinuous elements
1600 583223 : if (cont == DISCONTINUOUS)
1601 0 : return;
1602 48852 : libmesh_assert (cont == C_ZERO || cont == C_ONE ||
1603 : cont == SIDE_DISCONTINUOUS);
1604 :
1605 : // this would require some generalisation:
1606 : // - e.g. the 'my_fe'-object needs generalisation
1607 : // - due to lack of one-to-one correspondence of DOFs and nodes,
1608 : // this doesn't work easily.
1609 130923 : if (elem->infinite())
1610 0 : libmesh_not_implemented();
1611 :
1612 632075 : std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1613 : (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1614 48852 : neigh_fe->add_p_level_in_reinit(add_p_level);
1615 :
1616 632075 : QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1617 583223 : my_fe->attach_quadrature_rule (&my_qface);
1618 97704 : std::vector<Point> neigh_qface;
1619 :
1620 130923 : const std::vector<Real> & JxW = my_fe->get_JxW();
1621 130923 : const std::vector<Point> & q_point = my_fe->get_xyz();
1622 48852 : const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1623 48852 : const std::vector<std::vector<OutputShape>> & neigh_phi =
1624 : neigh_fe->get_phi();
1625 48852 : const std::vector<Point> * face_normals = nullptr;
1626 48852 : const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1627 48852 : const std::vector<std::vector<OutputGradient>> * neigh_dphi = nullptr;
1628 :
1629 97704 : std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1630 97704 : std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1631 :
1632 583223 : if (cont == C_ONE)
1633 : {
1634 7274 : const std::vector<Point> & ref_face_normals =
1635 7490 : my_fe->get_normals();
1636 3637 : face_normals = &ref_face_normals;
1637 3637 : const std::vector<std::vector<OutputGradient>> & ref_dphi =
1638 : my_fe->get_dphi();
1639 3637 : dphi = &ref_dphi;
1640 3637 : const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1641 : neigh_fe->get_dphi();
1642 3637 : neigh_dphi = &ref_neigh_dphi;
1643 : }
1644 :
1645 680927 : DenseMatrix<Real> Ke;
1646 583223 : DenseVector<Real> Fe;
1647 146556 : std::vector<DenseVector<Real>> Ue;
1648 :
1649 : // Look at the element faces. Check to see if we need to
1650 : // build constraints.
1651 2781288 : for (auto s : elem->side_index_range())
1652 : {
1653 : // Get pointers to the element's neighbor.
1654 2198065 : const Elem * neigh = elem->neighbor_ptr(s);
1655 :
1656 2198065 : if (!neigh)
1657 228637 : continue;
1658 :
1659 1948526 : if (!var.active_on_subdomain(neigh->subdomain_id()))
1660 1636 : continue;
1661 :
1662 : // h refinement constraints:
1663 : // constrain dofs shared between
1664 : // this element and ones coarser
1665 : // than this element.
1666 1946746 : if (neigh->level() < elem->level())
1667 : {
1668 119130 : unsigned int s_neigh = neigh->which_neighbor_am_i(elem);
1669 10120 : libmesh_assert_less (s_neigh, neigh->n_neighbors());
1670 :
1671 : // Find the minimum p level; we build the h constraint
1672 : // matrix with this and then constrain away all higher p
1673 : // DoFs.
1674 10120 : libmesh_assert(neigh->active());
1675 139370 : const unsigned int min_p_level = add_p_level *
1676 139370 : std::min(elem->p_level(), neigh->p_level());
1677 : // we may need to make the FE objects reinit with the
1678 : // minimum shared p_level
1679 119130 : const unsigned int old_elem_level = add_p_level * elem->p_level();
1680 119130 : if (old_elem_level != min_p_level)
1681 360 : my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + min_p_level - old_elem_level);
1682 119130 : const unsigned int old_neigh_level = add_p_level * neigh->p_level();
1683 119130 : if (old_neigh_level != min_p_level)
1684 0 : neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + min_p_level - old_neigh_level);
1685 :
1686 119130 : my_fe->reinit(elem, s);
1687 :
1688 : // This function gets called element-by-element, so there
1689 : // will be a lot of memory allocation going on. We can
1690 : // at least minimize this for the case of the dof indices
1691 : // by efficiently preallocating the requisite storage.
1692 : // n_nodes is not necessarily n_dofs, but it is better
1693 : // than nothing!
1694 119130 : my_dof_indices.reserve (elem->n_nodes());
1695 119130 : neigh_dof_indices.reserve (neigh->n_nodes());
1696 :
1697 119130 : dof_map.dof_indices (elem, my_dof_indices,
1698 : variable_number,
1699 : min_p_level);
1700 119130 : dof_map.dof_indices (neigh, neigh_dof_indices,
1701 : variable_number,
1702 : min_p_level);
1703 :
1704 10120 : const unsigned int n_qp = my_qface.n_points();
1705 :
1706 119130 : FEMap::inverse_map (Dim, neigh, q_point, neigh_qface);
1707 :
1708 119130 : neigh_fe->reinit(neigh, &neigh_qface);
1709 :
1710 : // We're only concerned with DOFs whose values (and/or first
1711 : // derivatives for C1 elements) are supported on side nodes
1712 119130 : FEType elem_fe_type = base_fe_type;
1713 119130 : if (old_elem_level != min_p_level)
1714 360 : elem_fe_type.order = base_fe_type.order.get_order() + min_p_level - old_elem_level;
1715 119130 : FEType neigh_fe_type = base_fe_type;
1716 119130 : if (old_neigh_level != min_p_level)
1717 0 : neigh_fe_type.order = base_fe_type.order.get_order() + min_p_level - old_neigh_level;
1718 119130 : FEInterface::dofs_on_side(elem, Dim, elem_fe_type, s, my_side_dofs);
1719 119130 : FEInterface::dofs_on_side(neigh, Dim, neigh_fe_type, s_neigh, neigh_side_dofs);
1720 :
1721 10120 : const unsigned int n_side_dofs =
1722 20240 : cast_int<unsigned int>(my_side_dofs.size());
1723 10120 : libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
1724 :
1725 : #ifndef NDEBUG
1726 55500 : for (auto i : my_side_dofs)
1727 45380 : libmesh_assert_less(i, my_dof_indices.size());
1728 55500 : for (auto i : neigh_side_dofs)
1729 45380 : libmesh_assert_less(i, neigh_dof_indices.size());
1730 : #endif
1731 :
1732 109010 : Ke.resize (n_side_dofs, n_side_dofs);
1733 119130 : Ue.resize(n_side_dofs);
1734 :
1735 : // Form the projection matrix, (inner product of fine basis
1736 : // functions against fine test functions)
1737 655624 : for (unsigned int is = 0; is != n_side_dofs; ++is)
1738 : {
1739 581874 : const unsigned int i = my_side_dofs[is];
1740 3219328 : for (unsigned int js = 0; js != n_side_dofs; ++js)
1741 : {
1742 2907718 : const unsigned int j = my_side_dofs[js];
1743 14925786 : for (unsigned int qp = 0; qp != n_qp; ++qp)
1744 : {
1745 18279240 : Ke(is,js) += JxW[qp] * TensorTools::inner_product(phi[i][qp], phi[j][qp]);
1746 12242952 : if (cont == C_ONE)
1747 4218632 : Ke(is,js) += JxW[qp] *
1748 646400 : TensorTools::inner_product((*dphi)[i][qp] *
1749 : (*face_normals)[qp],
1750 1616000 : (*dphi)[j][qp] *
1751 : (*face_normals)[qp]);
1752 : }
1753 : }
1754 : }
1755 :
1756 : // Form the right hand sides, (inner product of coarse basis
1757 : // functions against fine test functions)
1758 655624 : for (unsigned int is = 0; is != n_side_dofs; ++is)
1759 : {
1760 581874 : const unsigned int i = neigh_side_dofs[is];
1761 491114 : Fe.resize (n_side_dofs);
1762 3219328 : for (unsigned int js = 0; js != n_side_dofs; ++js)
1763 : {
1764 2907718 : const unsigned int j = my_side_dofs[js];
1765 14925786 : for (unsigned int qp = 0; qp != n_qp; ++qp)
1766 : {
1767 14255048 : Fe(js) += JxW[qp] *
1768 14255048 : TensorTools::inner_product(neigh_phi[i][qp],
1769 13249000 : phi[j][qp]);
1770 12242952 : if (cont == C_ONE)
1771 4218632 : Fe(js) += JxW[qp] *
1772 969600 : TensorTools::inner_product((*neigh_dphi)[i][qp] *
1773 : (*face_normals)[qp],
1774 1616000 : (*dphi)[j][qp] *
1775 : (*face_normals)[qp]);
1776 : }
1777 : }
1778 581874 : Ke.cholesky_solve(Fe, Ue[is]);
1779 : }
1780 :
1781 655624 : for (unsigned int js = 0; js != n_side_dofs; ++js)
1782 : {
1783 536494 : const unsigned int j = my_side_dofs[js];
1784 581874 : const dof_id_type my_dof_g = my_dof_indices[j];
1785 45380 : libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
1786 :
1787 : // Hunt for "constraining against myself" cases before
1788 : // we bother creating a constraint row
1789 45380 : bool self_constraint = false;
1790 2591581 : for (unsigned int is = 0; is != n_side_dofs; ++is)
1791 : {
1792 2217453 : const unsigned int i = neigh_side_dofs[is];
1793 2217453 : const dof_id_type their_dof_g = neigh_dof_indices[i];
1794 185750 : libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1795 :
1796 2217453 : if (their_dof_g == my_dof_g)
1797 : {
1798 : #ifndef NDEBUG
1799 13708 : const Real their_dof_value = Ue[is](js);
1800 13708 : libmesh_assert_less (std::abs(their_dof_value-1.),
1801 : 10*TOLERANCE);
1802 :
1803 86448 : for (unsigned int k = 0; k != n_side_dofs; ++k)
1804 72740 : libmesh_assert(k == is ||
1805 : std::abs(Ue[k](js)) <
1806 : 10*TOLERANCE);
1807 : #endif
1808 :
1809 13708 : self_constraint = true;
1810 13708 : break;
1811 : }
1812 : }
1813 :
1814 536494 : if (self_constraint)
1815 239175 : continue;
1816 :
1817 : DofConstraintRow * constraint_row;
1818 :
1819 : // we may be running constraint methods concurrently
1820 : // on multiple threads, so we need a lock to
1821 : // ensure that this constraint is "ours"
1822 : {
1823 31672 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1824 :
1825 374128 : if (dof_map.is_constrained_dof(my_dof_g))
1826 6863 : continue;
1827 :
1828 297319 : constraint_row = &(constraints[my_dof_g]);
1829 24809 : libmesh_assert(constraint_row->empty());
1830 : }
1831 :
1832 1685108 : for (unsigned int is = 0; is != n_side_dofs; ++is)
1833 : {
1834 1387789 : const unsigned int i = neigh_side_dofs[is];
1835 1387789 : const dof_id_type their_dof_g = neigh_dof_indices[i];
1836 114227 : libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
1837 114227 : libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
1838 :
1839 1502016 : const Real their_dof_value = Ue[is](js);
1840 :
1841 1387789 : if (std::abs(their_dof_value) < 10*TOLERANCE)
1842 760325 : continue;
1843 :
1844 51566 : constraint_row->emplace(their_dof_g, their_dof_value);
1845 : }
1846 : }
1847 :
1848 119130 : my_fe->set_fe_order(my_fe->get_fe_type().order.get_order() + old_elem_level - min_p_level);
1849 119130 : neigh_fe->set_fe_order(neigh_fe->get_fe_type().order.get_order() + old_neigh_level - min_p_level);
1850 : }
1851 :
1852 1946746 : if (add_p_level)
1853 : {
1854 : // p refinement constraints:
1855 : // constrain dofs shared between
1856 : // active elements and neighbors with
1857 : // lower polynomial degrees
1858 : const unsigned int min_p_level =
1859 2109995 : neigh->min_p_level_by_neighbor(elem, elem->p_level());
1860 2109995 : if (min_p_level < elem->p_level())
1861 : {
1862 : // Adaptive p refinement of non-hierarchic bases will
1863 : // require more coding
1864 48 : libmesh_assert(my_fe->is_hierarchic());
1865 576 : dof_map.constrain_p_dofs(variable_number, elem,
1866 : s, min_p_level);
1867 : }
1868 : }
1869 : }
1870 1456557 : }
1871 :
1872 : #endif // #ifdef LIBMESH_ENABLE_AMR
1873 :
1874 :
1875 :
1876 : #ifdef LIBMESH_ENABLE_PERIODIC
1877 : template <typename OutputType>
1878 : void
1879 264164 : FEGenericBase<OutputType>::
1880 : compute_periodic_constraints (DofConstraints & constraints,
1881 : DofMap & dof_map,
1882 : const PeriodicBoundaries & boundaries,
1883 : const MeshBase & mesh,
1884 : const PointLocatorBase * point_locator,
1885 : const unsigned int variable_number,
1886 : const Elem * elem)
1887 : {
1888 : // Only bother if we truly have periodic boundaries
1889 264164 : if (boundaries.empty())
1890 58473 : return;
1891 :
1892 67082 : libmesh_assert(elem);
1893 :
1894 : // Only constrain active elements with this method
1895 67082 : if (!elem->active())
1896 19491 : return;
1897 :
1898 138677 : if (elem->infinite())
1899 0 : libmesh_not_implemented();
1900 :
1901 205691 : const unsigned int Dim = elem->dim();
1902 :
1903 : // We need sys_number and variable_number for DofObject methods
1904 : // later
1905 95182 : const unsigned int sys_number = dof_map.sys_number();
1906 :
1907 47591 : const FEType & base_fe_type = dof_map.variable_type(variable_number);
1908 :
1909 : // Construct FE objects for this element and its pseudo-neighbors.
1910 205691 : std::unique_ptr<FEGenericBase<OutputShape>> my_fe
1911 : (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1912 205691 : const FEContinuity cont = my_fe->get_continuity();
1913 :
1914 : // We don't need to constrain discontinuous elements
1915 205691 : if (cont == DISCONTINUOUS)
1916 0 : return;
1917 47591 : libmesh_assert (cont == C_ZERO || cont == C_ONE);
1918 :
1919 : // We'll use element size to generate relative tolerances later
1920 205691 : const Real primary_hmin = elem->hmin();
1921 :
1922 253282 : std::unique_ptr<FEGenericBase<OutputShape>> neigh_fe
1923 : (FEGenericBase<OutputShape>::build(Dim, base_fe_type));
1924 :
1925 253282 : QGauss my_qface(Dim-1, base_fe_type.default_quadrature_order());
1926 205691 : my_fe->attach_quadrature_rule (&my_qface);
1927 95182 : std::vector<Point> neigh_qface;
1928 :
1929 138677 : const std::vector<Real> & JxW = my_fe->get_JxW();
1930 138677 : const std::vector<Point> & q_point = my_fe->get_xyz();
1931 47591 : const std::vector<std::vector<OutputShape>> & phi = my_fe->get_phi();
1932 47591 : const std::vector<std::vector<OutputShape>> & neigh_phi =
1933 : neigh_fe->get_phi();
1934 47591 : const std::vector<Point> * face_normals = nullptr;
1935 47591 : const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1936 47591 : const std::vector<std::vector<OutputGradient>> * neigh_dphi = nullptr;
1937 95182 : std::vector<dof_id_type> my_dof_indices, neigh_dof_indices;
1938 95182 : std::vector<unsigned int> my_side_dofs, neigh_side_dofs;
1939 :
1940 205691 : if (cont != C_ZERO)
1941 : {
1942 8192 : const std::vector<Point> & ref_face_normals =
1943 4096 : my_fe->get_normals();
1944 4096 : face_normals = &ref_face_normals;
1945 4096 : const std::vector<std::vector<OutputGradient>> & ref_dphi =
1946 : my_fe->get_dphi();
1947 4096 : dphi = &ref_dphi;
1948 4096 : const std::vector<std::vector<OutputGradient>> & ref_neigh_dphi =
1949 : neigh_fe->get_dphi();
1950 4096 : neigh_dphi = &ref_neigh_dphi;
1951 : }
1952 :
1953 300873 : DenseMatrix<Real> Ke;
1954 205691 : DenseVector<Real> Fe;
1955 142773 : std::vector<DenseVector<Real>> Ue;
1956 :
1957 : // Container to catch the boundary ids that BoundaryInfo hands us.
1958 95182 : std::vector<boundary_id_type> bc_ids;
1959 :
1960 : // Look at the element faces. Check to see if we need to
1961 : // build constraints.
1962 205691 : const unsigned short int max_ns = elem->n_sides();
1963 1018023 : for (unsigned short int s = 0; s != max_ns; ++s)
1964 : {
1965 1001144 : if (elem->neighbor_ptr(s))
1966 589194 : continue;
1967 :
1968 39424 : mesh.get_boundary_info().boundary_ids (elem, s, bc_ids);
1969 :
1970 59984 : for (const auto & boundary_id : bc_ids)
1971 : {
1972 20560 : const PeriodicBoundaryBase * periodic = boundaries.boundary(boundary_id);
1973 20560 : if (!periodic || !periodic->is_my_variable(variable_number))
1974 5784 : continue;
1975 :
1976 3044 : libmesh_assert(point_locator);
1977 :
1978 : // Get pointers to the element's neighbor.
1979 : unsigned int s_neigh;
1980 14776 : const Elem * neigh = boundaries.neighbor(boundary_id, *point_locator, elem, s, &s_neigh);
1981 :
1982 14776 : libmesh_error_msg_if(neigh == nullptr,
1983 : "PeriodicBoundaries point locator object returned nullptr!");
1984 :
1985 : // periodic (and possibly h refinement) constraints:
1986 : // constrain dofs shared between
1987 : // this element and ones as coarse
1988 : // as or coarser than this element.
1989 14776 : if (neigh->level() <= elem->level())
1990 : {
1991 : #ifdef LIBMESH_ENABLE_AMR
1992 : // Find the minimum p level; we build the h constraint
1993 : // matrix with this and then constrain away all higher p
1994 : // DoFs.
1995 2596 : libmesh_assert(neigh->active());
1996 13432 : const unsigned int min_p_level =
1997 18624 : std::min(elem->p_level(), neigh->p_level());
1998 :
1999 : // we may need to make the FE objects reinit with the
2000 : // minimum shared p_level
2001 : // FIXME - I hate using const_cast<> and avoiding
2002 : // accessor functions; there's got to be a
2003 : // better way to do this!
2004 2596 : const unsigned int old_elem_level = elem->p_level();
2005 13432 : if (old_elem_level != min_p_level)
2006 0 : (const_cast<Elem *>(elem))->hack_p_level(min_p_level);
2007 5192 : const unsigned int old_neigh_level = neigh->p_level();
2008 13432 : if (old_neigh_level != min_p_level)
2009 0 : (const_cast<Elem *>(neigh))->hack_p_level(min_p_level);
2010 : #endif // #ifdef LIBMESH_ENABLE_AMR
2011 :
2012 : // We can do a projection with a single integration,
2013 : // due to the assumption of nested finite element
2014 : // subspaces.
2015 : // FIXME: it might be more efficient to do nodes,
2016 : // then edges, then side, to reduce the size of the
2017 : // Cholesky factorization(s)
2018 13432 : my_fe->reinit(elem, s);
2019 :
2020 13432 : dof_map.dof_indices (elem, my_dof_indices,
2021 : variable_number);
2022 13432 : dof_map.dof_indices (neigh, neigh_dof_indices,
2023 : variable_number);
2024 :
2025 : // We use neigh_dof_indices_all_variables in the case that the
2026 : // periodic boundary condition involves mappings between multiple
2027 : // variables.
2028 7788 : std::vector<std::vector<dof_id_type>> neigh_dof_indices_all_variables;
2029 13432 : if(periodic->has_transformation_matrix())
2030 : {
2031 7200 : const std::set<unsigned int> & variables = periodic->get_variables();
2032 7200 : neigh_dof_indices_all_variables.resize(variables.size());
2033 600 : unsigned int index = 0;
2034 28800 : for(unsigned int var : variables)
2035 : {
2036 23400 : dof_map.dof_indices (neigh, neigh_dof_indices_all_variables[index],
2037 : var);
2038 21600 : index++;
2039 : }
2040 : }
2041 :
2042 2596 : const unsigned int n_qp = my_qface.n_points();
2043 :
2044 : // Translate the quadrature points over to the
2045 : // neighbor's boundary
2046 18624 : std::vector<Point> neigh_point(q_point.size());
2047 54848 : for (auto i : index_range(neigh_point))
2048 47820 : neigh_point[i] = periodic->get_corresponding_pos(q_point[i]);
2049 :
2050 13432 : FEMap::inverse_map (Dim, neigh, neigh_point,
2051 : neigh_qface);
2052 :
2053 13432 : neigh_fe->reinit(neigh, &neigh_qface);
2054 :
2055 : // We're only concerned with DOFs whose values (and/or first
2056 : // derivatives for C1 elements) are supported on side nodes
2057 13432 : FEInterface::dofs_on_side(elem, Dim, base_fe_type, s, my_side_dofs);
2058 13432 : FEInterface::dofs_on_side(neigh, Dim, base_fe_type, s_neigh, neigh_side_dofs);
2059 :
2060 : // We're done with functions that examine Elem::p_level(),
2061 : // so let's unhack those levels
2062 : #ifdef LIBMESH_ENABLE_AMR
2063 16028 : if (elem->p_level() != old_elem_level)
2064 0 : (const_cast<Elem *>(elem))->hack_p_level(old_elem_level);
2065 16028 : if (neigh->p_level() != old_neigh_level)
2066 0 : (const_cast<Elem *>(neigh))->hack_p_level(old_neigh_level);
2067 : #endif // #ifdef LIBMESH_ENABLE_AMR
2068 :
2069 2596 : const unsigned int n_side_dofs =
2070 : cast_int<unsigned int>
2071 5192 : (my_side_dofs.size());
2072 2596 : libmesh_assert_equal_to (n_side_dofs, neigh_side_dofs.size());
2073 :
2074 10836 : Ke.resize (n_side_dofs, n_side_dofs);
2075 13432 : Ue.resize(n_side_dofs);
2076 :
2077 : // Form the projection matrix, (inner product of fine basis
2078 : // functions against fine test functions)
2079 54936 : for (unsigned int is = 0; is != n_side_dofs; ++is)
2080 : {
2081 47916 : const unsigned int i = my_side_dofs[is];
2082 182832 : for (unsigned int js = 0; js != n_side_dofs; ++js)
2083 : {
2084 159012 : const unsigned int j = my_side_dofs[js];
2085 656192 : for (unsigned int qp = 0; qp != n_qp; ++qp)
2086 : {
2087 624296 : Ke(is,js) += JxW[qp] *
2088 624296 : TensorTools::inner_product(phi[i][qp],
2089 569580 : phi[j][qp]);
2090 514864 : if (cont != C_ZERO)
2091 384 : Ke(is,js) += JxW[qp] *
2092 64 : TensorTools::inner_product((*dphi)[i][qp] *
2093 : (*face_normals)[qp],
2094 160 : (*dphi)[j][qp] *
2095 : (*face_normals)[qp]);
2096 : }
2097 : }
2098 : }
2099 :
2100 : // Form the right hand sides, (inner product of coarse basis
2101 : // functions against fine test functions)
2102 54936 : for (unsigned int is = 0; is != n_side_dofs; ++is)
2103 : {
2104 47916 : const unsigned int i = neigh_side_dofs[is];
2105 35092 : Fe.resize (n_side_dofs);
2106 182832 : for (unsigned int js = 0; js != n_side_dofs; ++js)
2107 : {
2108 159012 : const unsigned int j = my_side_dofs[js];
2109 656192 : for (unsigned int qp = 0; qp != n_qp; ++qp)
2110 : {
2111 624296 : Fe(js) += JxW[qp] *
2112 624296 : TensorTools::inner_product(neigh_phi[i][qp],
2113 569580 : phi[j][qp]);
2114 514864 : if (cont != C_ZERO)
2115 384 : Fe(js) += JxW[qp] *
2116 96 : TensorTools::inner_product((*neigh_dphi)[i][qp] *
2117 : (*face_normals)[qp],
2118 160 : (*dphi)[j][qp] *
2119 : (*face_normals)[qp]);
2120 : }
2121 : }
2122 47916 : Ke.cholesky_solve(Fe, Ue[is]);
2123 : }
2124 :
2125 : // Make sure we're not adding recursive constraints
2126 : // due to the redundancy in the way we add periodic
2127 : // boundary constraints
2128 : //
2129 : // In order for this to work while threaded or on
2130 : // distributed meshes, we need a rigorous way to
2131 : // avoid recursive constraints. Here it is:
2132 : //
2133 : // For vertex DoFs, if there is a "prior" element
2134 : // (i.e. a coarser element or an equally refined
2135 : // element with a lower id) on this boundary which
2136 : // contains the vertex point, then we will avoid
2137 : // generating constraints; the prior element (or
2138 : // something prior to it) may do so. If we are the
2139 : // most prior (or "primary") element on this
2140 : // boundary sharing this point, then we look at the
2141 : // boundary periodic to us, we find the primary
2142 : // element there, and if that primary is coarser or
2143 : // equal-but-lower-id, then our vertex dofs are
2144 : // constrained in terms of that element.
2145 : //
2146 : // For edge DoFs, if there is a coarser element
2147 : // on this boundary sharing this edge, then we will
2148 : // avoid generating constraints (we will be
2149 : // constrained indirectly via AMR constraints
2150 : // connecting us to the coarser element's DoFs). If
2151 : // we are the coarsest element sharing this edge,
2152 : // then we generate constraints if and only if we
2153 : // are finer than the coarsest element on the
2154 : // boundary periodic to us sharing the corresponding
2155 : // periodic edge, or if we are at equal level but
2156 : // our edge nodes have higher ids than the periodic
2157 : // edge nodes (sorted from highest to lowest, then
2158 : // compared lexicographically)
2159 : //
2160 : // For face DoFs, we generate constraints if we are
2161 : // finer than our periodic neighbor, or if we are at
2162 : // equal level but our element id is higher than its
2163 : // element id.
2164 : //
2165 : // If the primary neighbor is also the current elem
2166 : // (a 1-element-thick mesh) then we choose which
2167 : // vertex dofs to constrain via lexicographic
2168 : // ordering on point locations
2169 :
2170 : // FIXME: This code doesn't yet properly handle
2171 : // cases where multiple different periodic BCs
2172 : // intersect.
2173 5192 : std::set<dof_id_type> my_constrained_dofs;
2174 :
2175 : // Container to catch boundary IDs handed back by BoundaryInfo.
2176 5192 : std::vector<boundary_id_type> new_bc_ids;
2177 :
2178 96264 : for (auto n : elem->node_index_range())
2179 : {
2180 82832 : if (!elem->is_node_on_side(n,s))
2181 35012 : continue;
2182 :
2183 6404 : const Node & my_node = elem->node_ref(n);
2184 :
2185 41416 : if (elem->is_vertex(n))
2186 : {
2187 : // Find all boundary ids that include this
2188 : // point and have periodic boundary
2189 : // conditions for this variable
2190 6384 : std::set<boundary_id_type> point_bcids;
2191 :
2192 256440 : for (unsigned int new_s = 0;
2193 262824 : new_s != max_ns; ++new_s)
2194 : {
2195 221648 : if (!elem->is_node_on_side(n,new_s))
2196 95464 : continue;
2197 :
2198 111064 : mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
2199 :
2200 222128 : for (const auto & new_boundary_id : new_bc_ids)
2201 : {
2202 111064 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2203 111064 : if (new_periodic && new_periodic->is_my_variable(variable_number))
2204 95904 : point_bcids.insert(new_boundary_id);
2205 : }
2206 : }
2207 :
2208 : // See if this vertex has point neighbors to
2209 : // defer to
2210 44655 : if (primary_boundary_point_neighbor
2211 41176 : (elem, my_node, mesh.get_boundary_info(), point_bcids)
2212 6384 : != elem)
2213 21511 : continue;
2214 :
2215 : // Find the complementary boundary id set
2216 2905 : std::set<boundary_id_type> point_pairedids;
2217 32372 : for (const auto & new_boundary_id : point_bcids)
2218 : {
2219 16186 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2220 16186 : point_pairedids.insert(new_periodic->pairedboundary);
2221 : }
2222 :
2223 : // What do we want to constrain against?
2224 2905 : const Elem * primary_elem = nullptr;
2225 2905 : const Elem * main_neigh = nullptr;
2226 16186 : Point main_pt = my_node,
2227 16186 : primary_pt = my_node;
2228 :
2229 32372 : for (const auto & new_boundary_id : point_bcids)
2230 : {
2231 : // Find the corresponding periodic point and
2232 : // its primary neighbor
2233 16186 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2234 :
2235 : const Point neigh_pt =
2236 16186 : new_periodic->get_corresponding_pos(my_node);
2237 :
2238 : // If the point is getting constrained
2239 : // to itself by this PBC then we don't
2240 : // generate any constraints
2241 16186 : if (neigh_pt.absolute_fuzzy_equals
2242 16186 : (my_node, primary_hmin*TOLERANCE))
2243 6828 : continue;
2244 :
2245 : // Otherwise we'll have a constraint in
2246 : // one direction or another
2247 16186 : if (!primary_elem)
2248 2905 : primary_elem = elem;
2249 :
2250 2905 : const Elem * primary_neigh =
2251 16186 : primary_boundary_point_neighbor(neigh, neigh_pt,
2252 : mesh.get_boundary_info(),
2253 : point_pairedids);
2254 :
2255 2905 : libmesh_assert(primary_neigh);
2256 :
2257 16186 : if (new_boundary_id == boundary_id)
2258 : {
2259 2905 : main_neigh = primary_neigh;
2260 16186 : main_pt = neigh_pt;
2261 : }
2262 :
2263 : // Finer elements will get constrained in
2264 : // terms of coarser neighbors, not the
2265 : // other way around
2266 29467 : if ((primary_neigh->level() > primary_elem->level()) ||
2267 :
2268 : // For equal-level elements, the one with
2269 : // higher id gets constrained in terms of
2270 : // the one with lower id
2271 25646 : (primary_neigh->level() == primary_elem->level() &&
2272 28558 : primary_neigh->id() > primary_elem->id()) ||
2273 :
2274 : // On a one-element-thick mesh, we compare
2275 : // points to see what side gets constrained
2276 1880 : (primary_neigh == primary_elem &&
2277 0 : (neigh_pt > primary_pt)))
2278 6828 : continue;
2279 :
2280 1880 : primary_elem = primary_neigh;
2281 9358 : primary_pt = neigh_pt;
2282 : }
2283 :
2284 13493 : if (!primary_elem ||
2285 19091 : primary_elem != main_neigh ||
2286 1880 : primary_pt != main_pt)
2287 1025 : continue;
2288 : }
2289 240 : else if (elem->is_edge(n))
2290 : {
2291 : // Find which edge we're on
2292 240 : unsigned int e=0, ne = elem->n_edges();
2293 384 : for (; e != ne; ++e)
2294 : {
2295 384 : if (elem->is_node_on_edge(n,e))
2296 20 : break;
2297 : }
2298 20 : libmesh_assert_less (e, elem->n_edges());
2299 :
2300 : // Find the edge end nodes
2301 : const Node
2302 20 : * e1 = nullptr,
2303 20 : * e2 = nullptr;
2304 612 : for (auto nn : elem->node_index_range())
2305 : {
2306 612 : if (nn == n)
2307 0 : continue;
2308 :
2309 612 : if (elem->is_node_on_edge(nn, e))
2310 : {
2311 480 : if (e1 == nullptr)
2312 : {
2313 40 : e1 = elem->node_ptr(nn);
2314 : }
2315 : else
2316 : {
2317 40 : e2 = elem->node_ptr(nn);
2318 240 : break;
2319 : }
2320 : }
2321 : }
2322 20 : libmesh_assert (e1 && e2);
2323 :
2324 : // Find all boundary ids that include this
2325 : // edge and have periodic boundary
2326 : // conditions for this variable
2327 20 : std::set<boundary_id_type> edge_bcids;
2328 :
2329 940 : for (unsigned int new_s = 0;
2330 960 : new_s != max_ns; ++new_s)
2331 : {
2332 720 : if (!elem->is_node_on_side(n,new_s))
2333 440 : continue;
2334 :
2335 : // We're reusing the new_bc_ids vector created outside the loop over nodes.
2336 240 : mesh.get_boundary_info().boundary_ids (elem, s, new_bc_ids);
2337 :
2338 480 : for (const auto & new_boundary_id : new_bc_ids)
2339 : {
2340 240 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2341 240 : if (new_periodic && new_periodic->is_my_variable(variable_number))
2342 220 : edge_bcids.insert(new_boundary_id);
2343 : }
2344 : }
2345 :
2346 :
2347 : // See if this edge has neighbors to defer to
2348 240 : if (primary_boundary_edge_neighbor
2349 240 : (elem, *e1, *e2, mesh.get_boundary_info(), edge_bcids)
2350 20 : != elem)
2351 0 : continue;
2352 :
2353 : // Find the complementary boundary id set
2354 20 : std::set<boundary_id_type> edge_pairedids;
2355 480 : for (const auto & new_boundary_id : edge_bcids)
2356 : {
2357 240 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2358 240 : edge_pairedids.insert(new_periodic->pairedboundary);
2359 : }
2360 :
2361 : // What do we want to constrain against?
2362 20 : const Elem * primary_elem = nullptr;
2363 20 : const Elem * main_neigh = nullptr;
2364 240 : Point main_pt1 = *e1,
2365 240 : main_pt2 = *e2,
2366 240 : primary_pt1 = *e1,
2367 240 : primary_pt2 = *e2;
2368 :
2369 480 : for (const auto & new_boundary_id : edge_bcids)
2370 : {
2371 : // Find the corresponding periodic edge and
2372 : // its primary neighbor
2373 240 : const PeriodicBoundaryBase * new_periodic = boundaries.boundary(new_boundary_id);
2374 :
2375 240 : Point neigh_pt1 = new_periodic->get_corresponding_pos(*e1),
2376 240 : neigh_pt2 = new_periodic->get_corresponding_pos(*e2);
2377 :
2378 : // If the edge is getting constrained
2379 : // to itself by this PBC then we don't
2380 : // generate any constraints
2381 20 : if (neigh_pt1.absolute_fuzzy_equals
2382 260 : (*e1, primary_hmin*TOLERANCE) &&
2383 : neigh_pt2.absolute_fuzzy_equals
2384 0 : (*e2, primary_hmin*TOLERANCE))
2385 120 : continue;
2386 :
2387 : // Otherwise we'll have a constraint in
2388 : // one direction or another
2389 240 : if (!primary_elem)
2390 20 : primary_elem = elem;
2391 :
2392 20 : const Elem * primary_neigh = primary_boundary_edge_neighbor
2393 240 : (neigh, neigh_pt1, neigh_pt2,
2394 : mesh.get_boundary_info(), edge_pairedids);
2395 :
2396 20 : libmesh_assert(primary_neigh);
2397 :
2398 240 : if (new_boundary_id == boundary_id)
2399 : {
2400 20 : main_neigh = primary_neigh;
2401 240 : main_pt1 = neigh_pt1;
2402 240 : main_pt2 = neigh_pt2;
2403 : }
2404 :
2405 : // If we have a one-element thick mesh,
2406 : // we'll need to sort our points to get a
2407 : // consistent ordering rule
2408 : //
2409 : // Use >= in this test to make sure that,
2410 : // for angular constraints, no node gets
2411 : // constrained to itself.
2412 240 : if (primary_neigh == primary_elem)
2413 : {
2414 0 : if (primary_pt1 > primary_pt2)
2415 0 : std::swap(primary_pt1, primary_pt2);
2416 0 : if (neigh_pt1 > neigh_pt2)
2417 0 : std::swap(neigh_pt1, neigh_pt2);
2418 :
2419 0 : if (neigh_pt2 >= primary_pt2)
2420 0 : continue;
2421 : }
2422 :
2423 : // Otherwise:
2424 : // Finer elements will get constrained in
2425 : // terms of coarser ones, not the other way
2426 : // around
2427 480 : if ((primary_neigh->level() > primary_elem->level()) ||
2428 :
2429 : // For equal-level elements, the one with
2430 : // higher id gets constrained in terms of
2431 : // the one with lower id
2432 440 : (primary_neigh->level() == primary_elem->level() &&
2433 40 : primary_neigh->id() > primary_elem->id()))
2434 120 : continue;
2435 :
2436 10 : primary_elem = primary_neigh;
2437 120 : primary_pt1 = neigh_pt1;
2438 120 : primary_pt2 = neigh_pt2;
2439 : }
2440 :
2441 160 : if (!primary_elem ||
2442 230 : primary_elem != main_neigh ||
2443 270 : primary_pt1 != main_pt1 ||
2444 10 : primary_pt2 != main_pt2)
2445 10 : continue;
2446 : }
2447 0 : else if (elem->is_face(n))
2448 : {
2449 : // If we have a one-element thick mesh,
2450 : // use the ordering of the face node and its
2451 : // periodic counterpart to determine what
2452 : // gets constrained
2453 0 : if (neigh == elem)
2454 : {
2455 : const Point neigh_pt =
2456 0 : periodic->get_corresponding_pos(my_node);
2457 0 : if (neigh_pt > my_node)
2458 0 : continue;
2459 : }
2460 :
2461 : // Otherwise:
2462 : // Finer elements will get constrained in
2463 : // terms of coarser ones, not the other way
2464 : // around
2465 0 : if ((neigh->level() > elem->level()) ||
2466 :
2467 : // For equal-level elements, the one with
2468 : // higher id gets constrained in terms of
2469 : // the one with lower id
2470 0 : (neigh->level() == elem->level() &&
2471 0 : neigh->id() > elem->id()))
2472 0 : continue;
2473 : }
2474 :
2475 : // If we made it here without hitting a continue
2476 : // statement, then we're at a node whose dofs
2477 : // should be constrained by this element's
2478 : // calculations.
2479 : const unsigned int n_comp =
2480 9478 : my_node.n_comp(sys_number, variable_number);
2481 :
2482 19000 : for (unsigned int i=0; i != n_comp; ++i)
2483 : my_constrained_dofs.insert
2484 7628 : (my_node.dof_number
2485 9522 : (sys_number, variable_number, i));
2486 : }
2487 :
2488 : // FIXME: old code for disambiguating periodic BCs:
2489 : // this is not threadsafe nor safe to run on a
2490 : // non-serialized mesh.
2491 : /*
2492 : std::vector<bool> recursive_constraint(n_side_dofs, false);
2493 :
2494 : for (unsigned int is = 0; is != n_side_dofs; ++is)
2495 : {
2496 : const unsigned int i = neigh_side_dofs[is];
2497 : const dof_id_type their_dof_g = neigh_dof_indices[i];
2498 : libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2499 :
2500 : {
2501 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2502 :
2503 : if (!dof_map.is_constrained_dof(their_dof_g))
2504 : continue;
2505 : }
2506 :
2507 : DofConstraintRow & their_constraint_row =
2508 : constraints[their_dof_g].first;
2509 :
2510 : for (unsigned int js = 0; js != n_side_dofs; ++js)
2511 : {
2512 : const unsigned int j = my_side_dofs[js];
2513 : const dof_id_type my_dof_g = my_dof_indices[j];
2514 : libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2515 :
2516 : if (their_constraint_row.count(my_dof_g))
2517 : recursive_constraint[js] = true;
2518 : }
2519 : }
2520 : */
2521 :
2522 54936 : for (unsigned int js = 0; js != n_side_dofs; ++js)
2523 : {
2524 : // FIXME: old code path
2525 : // if (recursive_constraint[js])
2526 : // continue;
2527 :
2528 41504 : const unsigned int j = my_side_dofs[js];
2529 47916 : const dof_id_type my_dof_g = my_dof_indices[j];
2530 6412 : libmesh_assert_not_equal_to (my_dof_g, DofObject::invalid_id);
2531 :
2532 : // FIXME: new code path
2533 29358 : if (!my_constrained_dofs.count(my_dof_g))
2534 32254 : continue;
2535 :
2536 : DofConstraintRow * constraint_row;
2537 :
2538 : // we may be running constraint methods concurrently
2539 : // on multiple threads, so we need a lock to
2540 : // ensure that this constraint is "ours"
2541 : {
2542 1894 : Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
2543 :
2544 9522 : if (dof_map.is_constrained_dof(my_dof_g))
2545 81 : continue;
2546 :
2547 9250 : constraint_row = &(constraints[my_dof_g]);
2548 1813 : libmesh_assert(constraint_row->empty());
2549 : }
2550 :
2551 37506 : for (unsigned int is = 0; is != n_side_dofs; ++is)
2552 : {
2553 28256 : const unsigned int i = neigh_side_dofs[is];
2554 28256 : const dof_id_type their_dof_g = neigh_dof_indices[i];
2555 4439 : libmesh_assert_not_equal_to (their_dof_g, DofObject::invalid_id);
2556 :
2557 : // Periodic constraints should never be
2558 : // self-constraints
2559 : // libmesh_assert_not_equal_to (their_dof_g, my_dof_g);
2560 :
2561 28256 : const Real their_dof_value = Ue[is](js);
2562 :
2563 28256 : if (their_dof_g == my_dof_g)
2564 : {
2565 0 : libmesh_assert_less (std::abs(their_dof_value-1.), 1.e-5);
2566 0 : for (unsigned int k = 0; k != n_side_dofs; ++k)
2567 0 : libmesh_assert(k == is || std::abs(Ue[k](js)) < 1.e-5);
2568 17662 : continue;
2569 0 : }
2570 :
2571 28256 : if (std::abs(their_dof_value) < 10*TOLERANCE)
2572 15484 : continue;
2573 :
2574 10594 : if(!periodic->has_transformation_matrix())
2575 : {
2576 1865 : constraint_row->emplace(their_dof_g, their_dof_value);
2577 : }
2578 : else
2579 : {
2580 : // In this case the current variable is constrained in terms of other variables.
2581 : // We assume that all variables in this constraint have the same FE type (this
2582 : // is asserted below), and hence we can create the constraint row contribution
2583 : // by multiplying their_dof_value by the corresponding row of the transformation
2584 : // matrix.
2585 :
2586 4752 : const std::set<unsigned int> & variables = periodic->get_variables();
2587 4752 : neigh_dof_indices_all_variables.resize(variables.size());
2588 396 : unsigned int index = 0;
2589 19008 : for(unsigned int other_var : variables)
2590 : {
2591 1188 : libmesh_assert_msg(base_fe_type == dof_map.variable_type(other_var), "FE types must match for all variables involved in constraint");
2592 :
2593 14256 : Real var_weighting = periodic->get_transformation_matrix()(variable_number, other_var);
2594 14256 : constraint_row->emplace(neigh_dof_indices_all_variables[index][i],
2595 14256 : var_weighting*their_dof_value);
2596 14256 : index++;
2597 : }
2598 : }
2599 :
2600 : }
2601 : }
2602 8240 : }
2603 : // p refinement constraints:
2604 : // constrain dofs shared between
2605 : // active elements and neighbors with
2606 : // lower polynomial degrees
2607 : #ifdef LIBMESH_ENABLE_AMR
2608 : const unsigned int min_p_level =
2609 17820 : neigh->min_p_level_by_neighbor(elem, elem->p_level());
2610 17820 : if (min_p_level < elem->p_level())
2611 : {
2612 : // Adaptive p refinement of non-hierarchic bases will
2613 : // require more coding
2614 0 : libmesh_assert(my_fe->is_hierarchic());
2615 0 : dof_map.constrain_p_dofs(variable_number, elem,
2616 : s, min_p_level);
2617 : }
2618 : #endif // #ifdef LIBMESH_ENABLE_AMR
2619 : }
2620 : }
2621 331527 : }
2622 :
2623 : #endif // LIBMESH_ENABLE_PERIODIC
2624 :
2625 : // ------------------------------------------------------------
2626 : // Explicit instantiations
2627 : template class LIBMESH_EXPORT FEGenericBase<Real>;
2628 : template class LIBMESH_EXPORT FEGenericBase<RealGradient>;
2629 :
2630 : } // namespace libMesh
|
