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