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