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