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
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"
38 #include "libmesh/tensor_value.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 <>
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  {
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:
517  }
518  }
519
520
521
522
523  // 2D
524  case 2:
525  {
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:
581  }
582  }
583
584
585
586
587  // 3D
588  case 3:
589  {
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:
645  }
646  }
647
648  default:
649  libmesh_error_msg("Invalid dimension dim = " << dim);
650  }
651 }
652
653
654
655 template <>
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
698  this->_fe_trans->map_curl(this->dim, elem, qp, (*this), this->curl_phi);
699
700  // Only compute div for vector-valued elements
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
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
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  {
839  ref_dphi_values = fe->get_dphi();
840  dphi_values = &ref_dphi_values;
842  ref_dphi_coarse = fe_coarse->get_dphi();
843  dphi_coarse = &ref_dphi_coarse;
844  }
845
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
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;
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)
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) +=
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
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;
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)
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) +=
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
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;
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)
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
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  {
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
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  {
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  {
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>;
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
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
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
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