libMesh
dof_map_constraints.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 // Local includes
19 #include "libmesh/dof_map.h"
20 
21 // libMesh includes
22 #include "libmesh/boundary_info.h" // needed for dirichlet constraints
23 #include "libmesh/dense_matrix.h"
24 #include "libmesh/dense_vector.h"
25 #include "libmesh/dirichlet_boundaries.h"
26 #include "libmesh/elem.h"
27 #include "libmesh/elem_range.h"
28 #include "libmesh/fe_base.h"
29 #include "libmesh/fe_interface.h"
30 #include "libmesh/fe_type.h"
31 #include "libmesh/function_base.h"
32 #include "libmesh/int_range.h"
33 #include "libmesh/libmesh_logging.h"
34 #include "libmesh/linear_solver.h" // for spline Dirichlet projection solves
35 #include "libmesh/mesh_base.h"
36 #include "libmesh/null_output_iterator.h"
37 #include "libmesh/mesh_tools.h" // for libmesh_assert_valid_boundary_ids()
38 #include "libmesh/nonlinear_implicit_system.h"
39 #include "libmesh/numeric_vector.h" // for enforce_constraints_exactly()
40 #include "libmesh/parallel_algebra.h"
41 #include "libmesh/parallel_elem.h"
42 #include "libmesh/parallel_node.h"
43 #include "libmesh/periodic_boundaries.h"
44 #include "libmesh/periodic_boundary.h"
45 #include "libmesh/periodic_boundary_base.h"
46 #include "libmesh/point_locator_base.h"
47 #include "libmesh/quadrature.h" // for dirichlet constraints
48 #include "libmesh/raw_accessor.h"
49 #include "libmesh/sparse_matrix.h" // needed to constrain adjoint rhs
50 #include "libmesh/system.h" // needed by enforce_constraints_exactly()
51 #include "libmesh/tensor_tools.h"
52 #include "libmesh/threads.h"
53 #include "libmesh/enum_to_string.h"
54 #include "libmesh/coupling_matrix.h"
55 
56 // TIMPI includes
58 #include "timpi/parallel_sync.h"
59 
60 // C++ Includes
61 #include <set>
62 #include <algorithm> // for std::count, std::fill
63 #include <sstream>
64 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
65 #include <cmath>
66 #include <memory>
67 #include <numeric>
68 #include <unordered_set>
69 
70 // Anonymous namespace to hold helper classes
71 namespace {
72 
73 using namespace libMesh;
74 
75 class ComputeConstraints
76 {
77 public:
78  ComputeConstraints (DofConstraints & constraints,
79  DofMap & dof_map,
80 #ifdef LIBMESH_ENABLE_PERIODIC
81  PeriodicBoundaries & periodic_boundaries,
82 #endif
83  const MeshBase & mesh,
84  const unsigned int variable_number) :
85  _constraints(constraints),
86  _dof_map(dof_map),
87 #ifdef LIBMESH_ENABLE_PERIODIC
88  _periodic_boundaries(periodic_boundaries),
89 #endif
90  _mesh(mesh),
91  _variable_number(variable_number)
92  {}
93 
94  void operator()(const ConstElemRange & range) const
95  {
96  const Variable & var_description = _dof_map.variable(_variable_number);
97 
98 #ifdef LIBMESH_ENABLE_PERIODIC
99  std::unique_ptr<PointLocatorBase> point_locator;
100  const bool have_periodic_boundaries =
101  !_periodic_boundaries.empty();
102  if (have_periodic_boundaries && !range.empty())
103  point_locator = _mesh.sub_point_locator();
104 #endif
105 
106  for (const auto & elem : range)
107  if (var_description.active_on_subdomain(elem->subdomain_id()))
108  {
109 #ifdef LIBMESH_ENABLE_AMR
110  FEInterface::compute_constraints (_constraints,
111  _dof_map,
112  _variable_number,
113  elem);
114 #endif
115 #ifdef LIBMESH_ENABLE_PERIODIC
116  // FIXME: periodic constraints won't work on a non-serial
117  // mesh unless it's kept ghost elements from opposing
118  // boundaries!
119  if (have_periodic_boundaries)
121  _dof_map,
122  _periodic_boundaries,
123  _mesh,
124  point_locator.get(),
125  _variable_number,
126  elem);
127 #endif
128  }
129  }
130 
131 private:
132  DofConstraints & _constraints;
133  DofMap & _dof_map;
134 #ifdef LIBMESH_ENABLE_PERIODIC
135  PeriodicBoundaries & _periodic_boundaries;
136 #endif
137  const MeshBase & _mesh;
138  const unsigned int _variable_number;
139 };
140 
141 
142 
143 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
144 class ComputeNodeConstraints
145 {
146 public:
147  ComputeNodeConstraints (NodeConstraints & node_constraints,
148 #ifdef LIBMESH_ENABLE_PERIODIC
149  PeriodicBoundaries & periodic_boundaries,
150 #endif
151  const MeshBase & mesh) :
152  _node_constraints(node_constraints),
153 #ifdef LIBMESH_ENABLE_PERIODIC
154  _periodic_boundaries(periodic_boundaries),
155 #endif
156  _mesh(mesh)
157  {}
158 
159  void operator()(const ConstElemRange & range) const
160  {
161 #ifdef LIBMESH_ENABLE_PERIODIC
162  std::unique_ptr<PointLocatorBase> point_locator;
163  bool have_periodic_boundaries = !_periodic_boundaries.empty();
164  if (have_periodic_boundaries && !range.empty())
165  point_locator = _mesh.sub_point_locator();
166 #endif
167 
168  for (const auto & elem : range)
169  {
170 #ifdef LIBMESH_ENABLE_AMR
171  FEBase::compute_node_constraints (_node_constraints, elem);
172 #endif
173 #ifdef LIBMESH_ENABLE_PERIODIC
174  // FIXME: periodic constraints won't work on a non-serial
175  // mesh unless it's kept ghost elements from opposing
176  // boundaries!
177  if (have_periodic_boundaries)
179  _periodic_boundaries,
180  _mesh,
181  point_locator.get(),
182  elem);
183 #endif
184  }
185  }
186 
187 private:
188  NodeConstraints & _node_constraints;
189 #ifdef LIBMESH_ENABLE_PERIODIC
190  PeriodicBoundaries & _periodic_boundaries;
191 #endif
192  const MeshBase & _mesh;
193 };
194 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
195 
196 
197 #ifdef LIBMESH_ENABLE_DIRICHLET
198 
202 class AddConstraint
203 {
204 protected:
205  DofMap & dof_map;
206 
207 public:
208  AddConstraint(DofMap & dof_map_in) : dof_map(dof_map_in) {}
209 
210  virtual void operator()(dof_id_type dof_number,
211  const DofConstraintRow & constraint_row,
212  const Number constraint_rhs) const = 0;
213 };
214 
215 class AddPrimalConstraint : public AddConstraint
216 {
217 public:
218  AddPrimalConstraint(DofMap & dof_map_in) : AddConstraint(dof_map_in) {}
219 
220  virtual void operator()(dof_id_type dof_number,
221  const DofConstraintRow & constraint_row,
222  const Number constraint_rhs) const
223  {
224  if (!dof_map.is_constrained_dof(dof_number))
225  dof_map.add_constraint_row (dof_number, constraint_row,
226  constraint_rhs, true);
227  }
228 };
229 
230 class AddAdjointConstraint : public AddConstraint
231 {
232 private:
233  const unsigned int qoi_index;
234 
235 public:
236  AddAdjointConstraint(DofMap & dof_map_in, unsigned int qoi_index_in)
237  : AddConstraint(dof_map_in), qoi_index(qoi_index_in) {}
238 
239  virtual void operator()(dof_id_type dof_number,
240  const DofConstraintRow & constraint_row,
241  const Number constraint_rhs) const
242  {
243  dof_map.add_adjoint_constraint_row
244  (qoi_index, dof_number, constraint_row, constraint_rhs,
245  true);
246  }
247 };
248 
249 
250 
256 class ConstrainDirichlet
257 {
258 private:
259  DofMap & dof_map;
260  const MeshBase & mesh;
261  const Real time;
262  const DirichletBoundaries & dirichlets;
263 
264  const AddConstraint & add_fn;
265 
266  static Number f_component (FunctionBase<Number> * f,
267  FEMFunctionBase<Number> * f_fem,
268  const FEMContext * c,
269  unsigned int i,
270  const Point & p,
271  Real time)
272  {
273  if (f_fem)
274  {
275  if (c)
276  return f_fem->component(*c, i, p, time);
277  else
278  return std::numeric_limits<Real>::quiet_NaN();
279  }
280  return f->component(i, p, time);
281  }
282 
283  static Gradient g_component (FunctionBase<Gradient> * g,
285  const FEMContext * c,
286  unsigned int i,
287  const Point & p,
288  Real time)
289  {
290  if (g_fem)
291  {
292  if (c)
293  return g_fem->component(*c, i, p, time);
294  else
295  return std::numeric_limits<Number>::quiet_NaN();
296  }
297  return g->component(i, p, time);
298  }
299 
300 
301 
308  struct SingleElemBoundaryInfo
309  {
310  SingleElemBoundaryInfo(const BoundaryInfo & bi,
311  const std::map<boundary_id_type, std::set<std::pair<unsigned int, DirichletBoundary *>>> & ordered_map_in) :
312  boundary_info(bi),
313  boundary_id_to_ordered_dirichlet_boundaries(ordered_map_in),
314  elem(nullptr),
315  n_sides(0),
316  n_edges(0),
317  n_nodes(0)
318  {}
319 
320  const BoundaryInfo & boundary_info;
321  const std::map<boundary_id_type, std::set<std::pair<unsigned int, DirichletBoundary *>>> & boundary_id_to_ordered_dirichlet_boundaries;
322  const Elem * elem;
323 
324  unsigned short n_sides;
325  unsigned short n_edges;
326  unsigned short n_nodes;
327 
328  // Mapping from DirichletBoundary objects which are active on this
329  // element to sides/nodes/edges/shellfaces of this element which
330  // they are active on.
331  std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_node_map;
332  std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_side_map;
333  std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_edge_map;
334  std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_shellface_map;
335 
336  std::map<const DirichletBoundary *, std::vector<bool>> is_boundary_nodeset_map;
337 
338  // The set of (dirichlet_id, DirichletBoundary) pairs which have at least one boundary
339  // id related to this Elem.
340  std::set<std::pair<unsigned int, DirichletBoundary *>> ordered_dbs;
341 
349  bool reinit(const Elem * elem_in)
350  {
351  elem = elem_in;
352 
353  n_sides = elem->n_sides();
354  n_edges = elem->n_edges();
355  n_nodes = elem->n_nodes();
356 
357  // objects and node/side/edge/shellface ids.
358  is_boundary_node_map.clear();
359  is_boundary_side_map.clear();
360  is_boundary_edge_map.clear();
361  is_boundary_shellface_map.clear();
362  is_boundary_nodeset_map.clear();
363 
364  // Clear any DirichletBoundaries from the previous Elem
365  ordered_dbs.clear();
366 
367  // Update has_dirichlet_constraint below, and if it remains false then
368  // we can skip this element since there are not constraints to impose.
369  bool has_dirichlet_constraint = false;
370 
371  // Container to catch boundary ids handed back for sides,
372  // nodes, and edges in the loops below.
373  std::vector<boundary_id_type> ids_vec;
374 
375  for (unsigned char s = 0; s != n_sides; ++s)
376  {
377  // First see if this side has been requested
378  boundary_info.boundary_ids (elem, s, ids_vec);
379 
380  bool do_this_side = false;
381  for (const auto & bc_id : ids_vec)
382  {
383  if (auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
384  it != boundary_id_to_ordered_dirichlet_boundaries.end())
385  {
386  do_this_side = true;
387 
388  // Associate every DirichletBoundary object that has this bc_id with the current Elem
389  ordered_dbs.insert(it->second.begin(), it->second.end());
390 
391  // Turn on the flag for the current side for each DirichletBoundary
392  for (const auto & db_pair : it->second)
393  {
394  // Attempt to emplace an empty vector. If there
395  // is already an entry, the insertion will fail
396  // and we'll get an iterator back to the
397  // existing entry. Either way, we'll then set
398  // index s of that vector to "true".
399  auto pr = is_boundary_side_map.emplace(db_pair.second, std::vector<bool>(n_sides, false));
400  pr.first->second[s] = true;
401  }
402  }
403  }
404 
405  if (!do_this_side)
406  continue;
407 
408  has_dirichlet_constraint = true;
409 
410  // Then determine what nodes are on this side
411  for (unsigned int n = 0; n != n_nodes; ++n)
412  if (elem->is_node_on_side(n,s))
413  {
414  // Attempt to emplace an empty vector. If there is
415  // already an entry, the insertion will fail and we'll
416  // get an iterator back to the existing entry. Either
417  // way, we'll then set index n of that vector to
418  // "true".
419  for (const auto & db_pair : ordered_dbs)
420  {
421  // Only add this as a boundary node for this db if
422  // it is also a boundary side for this db.
423  if (auto side_it = is_boundary_side_map.find(db_pair.second);
424  side_it != is_boundary_side_map.end() && side_it->second[s])
425  {
426  auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
427  pr.first->second[n] = true;
428  }
429  }
430  }
431 
432  // Finally determine what edges are on this side
433  for (unsigned int e = 0; e != n_edges; ++e)
434  if (elem->is_edge_on_side(e,s))
435  {
436  // Attempt to emplace an empty vector. If there is
437  // already an entry, the insertion will fail and we'll
438  // get an iterator back to the existing entry. Either
439  // way, we'll then set index e of that vector to
440  // "true".
441  for (const auto & db_pair : ordered_dbs)
442  {
443  // Only add this as a boundary edge for this db if
444  // it is also a boundary side for this db.
445  if (auto side_it = is_boundary_side_map.find(db_pair.second);
446  side_it != is_boundary_side_map.end() && side_it->second[s])
447  {
448  auto pr = is_boundary_edge_map.emplace(db_pair.second, std::vector<bool>(n_edges, false));
449  pr.first->second[e] = true;
450  }
451  }
452  }
453  } // for (s = 0..n_sides)
454 
455  // We can also impose Dirichlet boundary conditions on nodes, so we should
456  // also independently check whether the nodes have been requested
457  for (unsigned int n=0; n != n_nodes; ++n)
458  {
459  boundary_info.boundary_ids (elem->node_ptr(n), ids_vec);
460 
461  for (const auto & bc_id : ids_vec)
462  {
463  if (auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
464  it != boundary_id_to_ordered_dirichlet_boundaries.end())
465  {
466  // Associate every DirichletBoundary object that has this bc_id with the current Elem
467  ordered_dbs.insert(it->second.begin(), it->second.end());
468 
469  // Turn on the flag for the current node for each DirichletBoundary
470  for (const auto & db_pair : it->second)
471  {
472  auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
473  pr.first->second[n] = true;
474 
475  auto pr2 = is_boundary_nodeset_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
476  pr2.first->second[n] = true;
477  }
478 
479  has_dirichlet_constraint = true;
480  }
481  }
482  } // for (n = 0..n_nodes)
483 
484  // We can also impose Dirichlet boundary conditions on edges, so we should
485  // also independently check whether the edges have been requested
486  for (unsigned short e=0; e != n_edges; ++e)
487  {
488  boundary_info.edge_boundary_ids (elem, e, ids_vec);
489 
490  bool do_this_side = false;
491  for (const auto & bc_id : ids_vec)
492  {
493  if (auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
494  it != boundary_id_to_ordered_dirichlet_boundaries.end())
495  {
496  do_this_side = true;
497 
498  // We need to loop over all DirichletBoundary objects associated with bc_id
499  ordered_dbs.insert(it->second.begin(), it->second.end());
500 
501  // Turn on the flag for the current edge for each DirichletBoundary
502  for (const auto & db_pair : it->second)
503  {
504  auto pr = is_boundary_edge_map.emplace(db_pair.second, std::vector<bool>(n_edges, false));
505  pr.first->second[e] = true;
506  }
507  }
508  }
509 
510  if (!do_this_side)
511  continue;
512 
513  has_dirichlet_constraint = true;
514 
515  // Then determine what nodes are on this edge
516  for (unsigned int n = 0; n != n_nodes; ++n)
517  if (elem->is_node_on_edge(n,e))
518  {
519  // Attempt to emplace an empty vector. If there is
520  // already an entry, the insertion will fail and we'll
521  // get an iterator back to the existing entry. Either
522  // way, we'll then set index n of that vector to
523  // "true".
524  for (const auto & db_pair : ordered_dbs)
525  {
526  // Only add this as a boundary node for this db if
527  // it is also a boundary edge for this db.
528  if (auto edge_it = is_boundary_edge_map.find(db_pair.second);
529  edge_it != is_boundary_edge_map.end() && edge_it->second[e])
530  {
531  auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
532  pr.first->second[n] = true;
533  }
534  }
535  }
536  }
537 
538  // We can also impose Dirichlet boundary conditions on shellfaces, so we should
539  // also independently check whether the shellfaces have been requested
540  for (unsigned short shellface=0; shellface != 2; ++shellface)
541  {
542  boundary_info.shellface_boundary_ids (elem, shellface, ids_vec);
543  bool do_this_shellface = false;
544 
545  for (const auto & bc_id : ids_vec)
546  {
547  if (auto it = boundary_id_to_ordered_dirichlet_boundaries.find(bc_id);
548  it != boundary_id_to_ordered_dirichlet_boundaries.end())
549  {
550  has_dirichlet_constraint = true;
551  do_this_shellface = true;
552 
553  // We need to loop over all DirichletBoundary objects associated with bc_id
554  ordered_dbs.insert(it->second.begin(), it->second.end());
555 
556  // Turn on the flag for the current shellface for each DirichletBoundary
557  for (const auto & db_pair : it->second)
558  {
559  auto pr = is_boundary_shellface_map.emplace(db_pair.second, std::vector<bool>(/*n_shellfaces=*/2, false));
560  pr.first->second[shellface] = true;
561  }
562  }
563  }
564 
565  if (do_this_shellface)
566  {
567  // Shellface BCs induce BCs on all the nodes of a shell Elem
568  for (unsigned int n = 0; n != n_nodes; ++n)
569  for (const auto & db_pair : ordered_dbs)
570  {
571  // Only add this as a boundary node for this db if
572  // it is also a boundary shellface for this db.
573  if (auto side_it = is_boundary_shellface_map.find(db_pair.second);
574  side_it != is_boundary_shellface_map.end() && side_it->second[shellface])
575  {
576  auto pr = is_boundary_node_map.emplace(db_pair.second, std::vector<bool>(n_nodes, false));
577  pr.first->second[n] = true;
578  }
579  }
580  }
581  } // for (shellface = 0..2)
582 
583  return has_dirichlet_constraint;
584  } // SingleElemBoundaryInfo::reinit()
585 
586  }; // struct SingleElemBoundaryInfo
587 
588 
589 
590  template<typename OutputType>
591  void apply_lagrange_dirichlet_impl(const SingleElemBoundaryInfo & sebi,
592  const Variable & variable,
593  const DirichletBoundary & dirichlet,
594  FEMContext & fem_context) const
595  {
596  // Get pointer to the Elem we are currently working on
597  const Elem * elem = sebi.elem;
598 
599  // Per-subdomain variables don't need to be projected on
600  // elements where they're not active
601  if (!variable.active_on_subdomain(elem->subdomain_id()))
602  return;
603 
604  FunctionBase<Number> * f = dirichlet.f.get();
605  FEMFunctionBase<Number> * f_fem = dirichlet.f_fem.get();
606 
607  const System * f_system = dirichlet.f_system;
608 
609  // We need data to project
610  libmesh_assert(f || f_fem);
611  libmesh_assert(!(f && f_fem));
612 
613  // Iff our data depends on a system, we should have it.
614  libmesh_assert(!(f && f_system));
615  libmesh_assert(!(f_fem && !f_system));
616 
617  // The new element coefficients. For Lagrange FEs, these are the
618  // nodal values.
620 
621  // Get a reference to the fe_type associated with this variable
622  const FEType & fe_type = variable.type();
623 
624  // Dimension of the vector-valued FE (1 for scalar-valued FEs)
625  unsigned int n_vec_dim = FEInterface::n_vec_dim(mesh, fe_type);
626 
627  const unsigned int var_component =
628  variable.first_scalar_number();
629 
630  // Get this Variable's number, as determined by the System.
631  const unsigned int var = variable.number();
632 
633  // If our supplied functions require a FEMContext, and if we have
634  // an initialized solution to use with that FEMContext, then
635  // create one. We're not going to use elem_jacobian or
636  // subjacobians here so don't allocate them.
637  std::unique_ptr<FEMContext> context;
638  if (f_fem)
639  {
640  libmesh_assert(f_system);
641  if (f_system->current_local_solution->initialized())
642  {
643  context = std::make_unique<FEMContext>
644  (*f_system, nullptr,
645  /* allocate local_matrices = */ false);
646  f_fem->init_context(*context);
647  }
648  }
649 
650  if (f_system && context.get())
651  context->pre_fe_reinit(*f_system, elem);
652 
653  // Also pre-init the fem_context() we were passed on the current Elem.
654  fem_context.pre_fe_reinit(fem_context.get_system(), elem);
655 
656  // Get a reference to the DOF indices for the current element
657  const std::vector<dof_id_type> & dof_indices =
658  fem_context.get_dof_indices(var);
659 
660  // The number of DOFs on the element
661  const unsigned int n_dofs =
662  cast_int<unsigned int>(dof_indices.size());
663 
664  // Fixed vs. free DoFs on edge/face projections
665  std::vector<char> dof_is_fixed(n_dofs, false); // bools
666 
667  // Zero the interpolated values
668  Ue.resize (n_dofs); Ue.zero();
669 
670  // For Lagrange elements, side, edge, and shellface BCs all
671  // "induce" boundary conditions on the nodes of those entities.
672  // In SingleElemBoundaryInfo::reinit(), we therefore set entries
673  // in the "is_boundary_node_map" container based on side and
674  // shellface BCs, Then, when we actually apply constraints, we
675  // only have to check whether any Nodes are in this container, and
676  // compute values as necessary.
677  unsigned int current_dof = 0;
678  for (unsigned int n=0; n!= sebi.n_nodes; ++n)
679  {
680  // For Lagrange this can return 0 (in case of a lower-order FE
681  // on a higher-order Elem) or 1. This function accounts for the
682  // Elem::p_level() internally.
683  const unsigned int nc =
684  FEInterface::n_dofs_at_node (fe_type, elem, n);
685 
686  // If there are no DOFs at this node, then it doesn't matter
687  // if it's technically a boundary node or not, there's nothing
688  // to constrain.
689  if (!nc)
690  continue;
691 
692  // Check whether the current node is a boundary node
693  auto is_boundary_node_it = sebi.is_boundary_node_map.find(&dirichlet);
694  const bool is_boundary_node =
695  (is_boundary_node_it != sebi.is_boundary_node_map.end() &&
696  is_boundary_node_it->second[n]);
697 
698  // Check whether the current node is in a boundary nodeset
699  auto is_boundary_nodeset_it = sebi.is_boundary_nodeset_map.find(&dirichlet);
700  const bool is_boundary_nodeset =
701  (is_boundary_nodeset_it != sebi.is_boundary_nodeset_map.end() &&
702  is_boundary_nodeset_it->second[n]);
703 
704  // If node is neither a boundary node or from a boundary nodeset, go to the next one.
705  if ( !(is_boundary_node || is_boundary_nodeset) )
706  {
707  current_dof += nc;
708  continue;
709  }
710 
711  // Compute function values, storing them in Ue
712  libmesh_assert_equal_to (nc, n_vec_dim);
713  for (unsigned int c = 0; c < n_vec_dim; c++)
714  {
715  Ue(current_dof+c) =
716  f_component(f, f_fem, context.get(), var_component+c,
717  elem->point(n), time);
718  dof_is_fixed[current_dof+c] = true;
719  }
720  current_dof += n_vec_dim;
721  } // end for (n=0..n_nodes)
722 
723  // Lock the DofConstraints since it is shared among threads.
724  {
725  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
726 
727  for (unsigned int i = 0; i < n_dofs; i++)
728  {
729  DofConstraintRow empty_row;
730  if (dof_is_fixed[i] && !libmesh_isnan(Ue(i)))
731  add_fn (dof_indices[i], empty_row, Ue(i));
732  }
733  }
734 
735  } // apply_lagrange_dirichlet_impl
736 
737 
738 
739  template<typename OutputType>
740  void apply_dirichlet_impl(const SingleElemBoundaryInfo & sebi,
741  const Variable & variable,
742  const DirichletBoundary & dirichlet,
743  FEMContext & fem_context) const
744  {
745  // Get pointer to the Elem we are currently working on
746  const Elem * elem = sebi.elem;
747 
748  // Per-subdomain variables don't need to be projected on
749  // elements where they're not active
750  if (!variable.active_on_subdomain(elem->subdomain_id()))
751  return;
752 
753  typedef OutputType OutputShape;
754  typedef typename TensorTools::IncrementRank<OutputShape>::type OutputGradient;
755  //typedef typename TensorTools::IncrementRank<OutputGradient>::type OutputTensor;
756  typedef typename TensorTools::MakeNumber<OutputShape>::type OutputNumber;
757  typedef typename TensorTools::IncrementRank<OutputNumber>::type OutputNumberGradient;
758  //typedef typename TensorTools::IncrementRank<OutputNumberGradient>::type OutputNumberTensor;
759 
760  FunctionBase<Number> * f = dirichlet.f.get();
761  FunctionBase<Gradient> * g = dirichlet.g.get();
762 
763  FEMFunctionBase<Number> * f_fem = dirichlet.f_fem.get();
764  FEMFunctionBase<Gradient> * g_fem = dirichlet.g_fem.get();
765 
766  const System * f_system = dirichlet.f_system;
767 
768  // We need data to project
769  libmesh_assert(f || f_fem);
770  libmesh_assert(!(f && f_fem));
771 
772  // Iff our data depends on a system, we should have it.
773  libmesh_assert(!(f && f_system));
774  libmesh_assert(!(f_fem && !f_system));
775 
776  // The element matrix and RHS for projections.
777  // Note that Ke is always real-valued, whereas
778  // Fe may be complex valued if complex number
779  // support is enabled
782  // The new element coefficients
784 
785  // The dimensionality of the current mesh
786  const unsigned int dim = mesh.mesh_dimension();
787 
788  // Get a reference to the fe_type associated with this variable
789  const FEType & fe_type = variable.type();
790 
791  // Dimension of the vector-valued FE (1 for scalar-valued FEs)
792  unsigned int n_vec_dim = FEInterface::n_vec_dim(mesh, fe_type);
793 
794  const unsigned int var_component =
795  variable.first_scalar_number();
796 
797  // Get this Variable's number, as determined by the System.
798  const unsigned int var = variable.number();
799 
800  // The type of projections done depend on the FE's continuity.
802 
803  // Make sure we have the right data available for C1 projections
804  if ((cont == C_ONE) && (fe_type.family != SUBDIVISION))
805  {
806  // We'll need gradient data for a C1 projection
807  libmesh_assert(g || g_fem);
808 
809  // We currently demand that either neither nor both function
810  // object depend on current FEM data.
811  libmesh_assert(!(g && g_fem));
812  libmesh_assert(!(f && g_fem));
813  libmesh_assert(!(f_fem && g));
814  }
815 
816  // If our supplied functions require a FEMContext, and if we have
817  // an initialized solution to use with that FEMContext, then
818  // create one. We're not going to use elem_jacobian or
819  // subjacobians here so don't allocate them.
820  std::unique_ptr<FEMContext> context;
821  if (f_fem)
822  {
823  libmesh_assert(f_system);
824  if (f_system->current_local_solution->initialized())
825  {
826  context = std::make_unique<FEMContext>
827  (*f_system, nullptr,
828  /* allocate local_matrices = */ false);
829  f_fem->init_context(*context);
830  if (g_fem)
831  g_fem->init_context(*context);
832  }
833  }
834 
835  // There's a chicken-and-egg problem with FEMFunction-based
836  // Dirichlet constraints: we can't evaluate the FEMFunction
837  // until we have an initialized local solution vector, we
838  // can't initialize local solution vectors until we have a
839  // send list, and we can't generate a send list until we know
840  // all our constraints
841  //
842  // We don't generate constraints on uninitialized systems;
843  // currently user code will have to reinit() before any
844  // FEMFunction-based constraints will be correct. This should
845  // be fine, since user code would want to reinit() after
846  // setting initial conditions anyway.
847  if (f_system && context.get())
848  context->pre_fe_reinit(*f_system, elem);
849 
850  // Also pre-init the fem_context() we were passed on the current Elem.
851  fem_context.pre_fe_reinit(fem_context.get_system(), elem);
852 
853  // Get a reference to the DOF indices for the current element
854  const std::vector<dof_id_type> & dof_indices =
855  fem_context.get_dof_indices(var);
856 
857  // The number of DOFs on the element
858  const unsigned int n_dofs =
859  cast_int<unsigned int>(dof_indices.size());
860 
861  // Fixed vs. free DoFs on edge/face projections
862  std::vector<char> dof_is_fixed(n_dofs, false); // bools
863  std::vector<int> free_dof(n_dofs, 0);
864 
865  // Zero the interpolated values
866  Ue.resize (n_dofs); Ue.zero();
867 
868  // In general, we need a series of
869  // projections to ensure a unique and continuous
870  // solution. We start by interpolating boundary nodes, then
871  // hold those fixed and project boundary edges, then hold
872  // those fixed and project boundary faces,
873 
874  // Interpolate node values first.
875  //
876  // If we have non-vertex nodes that have a boundary nodeset, we
877  // need to interpolate them directly, but that had better be
878  // happening in the apply_lagrange_dirichlet_impl code path,
879  // because you *can't* interpolate to evaluate non-nodal FE
880  // coefficients.
881  unsigned int current_dof = 0;
882  for (unsigned int n=0; n!= sebi.n_nodes; ++n)
883  {
884  // FIXME: this should go through the DofMap,
885  // not duplicate dof_indices code badly!
886 
887  // Get the number of DOFs at this node, accounting for
888  // Elem::p_level() internally.
889  const unsigned int nc =
890  FEInterface::n_dofs_at_node (fe_type, elem, n);
891 
892  // Get a reference to the "is_boundary_node" flags for the
893  // current DirichletBoundary object. In case the map does not
894  // contain an entry for this DirichletBoundary object, it
895  // means there are no boundary nodes active.
896  auto is_boundary_node_it = sebi.is_boundary_node_map.find(&dirichlet);
897 
898  // The current n is not a boundary node if either there is no
899  // boundary_node_map for this DirichletBoundary object, or if
900  // there is but the entry in the corresponding vector is
901  // false.
902  const bool not_boundary_node =
903  (is_boundary_node_it == sebi.is_boundary_node_map.end() ||
904  !is_boundary_node_it->second[n]);
905 
906  if (!elem->is_vertex(n) || not_boundary_node)
907  {
908  current_dof += nc;
909  continue;
910  }
911  if (cont == DISCONTINUOUS)
912  {
913  libmesh_assert_equal_to (nc, 0);
914  }
915  // Assume that C_ZERO elements have a single nodal
916  // value shape function
917  else if ((cont == C_ZERO) || (fe_type.family == SUBDIVISION))
918  {
919  libmesh_assert_equal_to (nc, n_vec_dim);
920  for (unsigned int c = 0; c < n_vec_dim; c++)
921  {
922  Ue(current_dof+c) =
923  f_component(f, f_fem, context.get(), var_component+c,
924  elem->point(n), time);
925  dof_is_fixed[current_dof+c] = true;
926  }
927  current_dof += n_vec_dim;
928  }
929  // The hermite element vertex shape functions are weird
930  else if (fe_type.family == HERMITE)
931  {
932  Ue(current_dof) =
933  f_component(f, f_fem, context.get(), var_component,
934  elem->point(n), time);
935  dof_is_fixed[current_dof] = true;
936  current_dof++;
937  Gradient grad =
938  g_component(g, g_fem, context.get(), var_component,
939  elem->point(n), time);
940  // x derivative
941  Ue(current_dof) = grad(0);
942  dof_is_fixed[current_dof] = true;
943  current_dof++;
944  if (dim > 1)
945  {
946  // We'll finite difference mixed derivatives
947  Point nxminus = elem->point(n),
948  nxplus = elem->point(n);
949  nxminus(0) -= TOLERANCE;
950  nxplus(0) += TOLERANCE;
951  Gradient gxminus =
952  g_component(g, g_fem, context.get(), var_component,
953  nxminus, time);
954  Gradient gxplus =
955  g_component(g, g_fem, context.get(), var_component,
956  nxplus, time);
957  // y derivative
958  Ue(current_dof) = grad(1);
959  dof_is_fixed[current_dof] = true;
960  current_dof++;
961  // xy derivative
962  Ue(current_dof) = (gxplus(1) - gxminus(1))
963  / 2. / TOLERANCE;
964  dof_is_fixed[current_dof] = true;
965  current_dof++;
966 
967  if (dim > 2)
968  {
969  // z derivative
970  Ue(current_dof) = grad(2);
971  dof_is_fixed[current_dof] = true;
972  current_dof++;
973  // xz derivative
974  Ue(current_dof) = (gxplus(2) - gxminus(2))
975  / 2. / TOLERANCE;
976  dof_is_fixed[current_dof] = true;
977  current_dof++;
978  // We need new points for yz
979  Point nyminus = elem->point(n),
980  nyplus = elem->point(n);
981  nyminus(1) -= TOLERANCE;
982  nyplus(1) += TOLERANCE;
983  Gradient gyminus =
984  g_component(g, g_fem, context.get(), var_component,
985  nyminus, time);
986  Gradient gyplus =
987  g_component(g, g_fem, context.get(), var_component,
988  nyplus, time);
989  // xz derivative
990  Ue(current_dof) = (gyplus(2) - gyminus(2))
991  / 2. / TOLERANCE;
992  dof_is_fixed[current_dof] = true;
993  current_dof++;
994  // Getting a 2nd order xyz is more tedious
995  Point nxmym = elem->point(n),
996  nxmyp = elem->point(n),
997  nxpym = elem->point(n),
998  nxpyp = elem->point(n);
999  nxmym(0) -= TOLERANCE;
1000  nxmym(1) -= TOLERANCE;
1001  nxmyp(0) -= TOLERANCE;
1002  nxmyp(1) += TOLERANCE;
1003  nxpym(0) += TOLERANCE;
1004  nxpym(1) -= TOLERANCE;
1005  nxpyp(0) += TOLERANCE;
1006  nxpyp(1) += TOLERANCE;
1007  Gradient gxmym =
1008  g_component(g, g_fem, context.get(), var_component,
1009  nxmym, time);
1010  Gradient gxmyp =
1011  g_component(g, g_fem, context.get(), var_component,
1012  nxmyp, time);
1013  Gradient gxpym =
1014  g_component(g, g_fem, context.get(), var_component,
1015  nxpym, time);
1016  Gradient gxpyp =
1017  g_component(g, g_fem, context.get(), var_component,
1018  nxpyp, time);
1019  Number gxzplus = (gxpyp(2) - gxmyp(2))
1020  / 2. / TOLERANCE;
1021  Number gxzminus = (gxpym(2) - gxmym(2))
1022  / 2. / TOLERANCE;
1023  // xyz derivative
1024  Ue(current_dof) = (gxzplus - gxzminus)
1025  / 2. / TOLERANCE;
1026  dof_is_fixed[current_dof] = true;
1027  current_dof++;
1028  }
1029  }
1030  }
1031  // Assume that other C_ONE elements have a single nodal
1032  // value shape function and nodal gradient component
1033  // shape functions
1034  else if (cont == C_ONE)
1035  {
1036  libmesh_assert_equal_to (nc, 1 + dim);
1037  Ue(current_dof) =
1038  f_component(f, f_fem, context.get(), var_component,
1039  elem->point(n), time);
1040  dof_is_fixed[current_dof] = true;
1041  current_dof++;
1042  Gradient grad =
1043  g_component(g, g_fem, context.get(), var_component,
1044  elem->point(n), time);
1045  for (unsigned int i=0; i!= dim; ++i)
1046  {
1047  Ue(current_dof) = grad(i);
1048  dof_is_fixed[current_dof] = true;
1049  current_dof++;
1050  }
1051  }
1052  else
1053  libmesh_error_msg("Unknown continuity cont = " << cont);
1054  } // end for (n=0..n_nodes)
1055 
1056  // In 3D, project any edge values next
1057  if (dim > 2 && cont != DISCONTINUOUS)
1058  {
1059  // Get a pointer to the 1 dimensional (edge) FE for the current
1060  // var which is stored in the fem_context. This will only be
1061  // different from side_fe in 3D.
1062  FEGenericBase<OutputType> * edge_fe = nullptr;
1063  fem_context.get_edge_fe(var, edge_fe);
1064 
1065  // Set tolerance on underlying FEMap object. This will allow us to
1066  // avoid spurious negative Jacobian errors while imposing BCs by
1067  // simply ignoring them. This should only be required in certain
1068  // special cases, see the DirichletBoundaries comments on this
1069  // parameter for more information.
1071 
1072  // Pre-request FE data
1073  const std::vector<std::vector<OutputShape>> & phi = edge_fe->get_phi();
1074  const std::vector<Point> & xyz_values = edge_fe->get_xyz();
1075  const std::vector<Real> & JxW = edge_fe->get_JxW();
1076 
1077  // Only pre-request gradients for C1 projections
1078  const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1079  if ((cont == C_ONE) && (fe_type.family != SUBDIVISION))
1080  {
1081  const std::vector<std::vector<OutputGradient>> & ref_dphi = edge_fe->get_dphi();
1082  dphi = &ref_dphi;
1083  }
1084 
1085  // Vector to hold edge local DOF indices
1086  std::vector<unsigned int> edge_dofs;
1087 
1088  // Get a reference to the "is_boundary_edge" flags for the
1089  // current DirichletBoundary object. In case the map does not
1090  // contain an entry for this DirichletBoundary object, it
1091  // means there are no boundary edges active.
1092  if (const auto is_boundary_edge_it = sebi.is_boundary_edge_map.find(&dirichlet);
1093  is_boundary_edge_it != sebi.is_boundary_edge_map.end())
1094  {
1095  for (unsigned int e=0; e != sebi.n_edges; ++e)
1096  {
1097  if (!is_boundary_edge_it->second[e])
1098  continue;
1099 
1100  FEInterface::dofs_on_edge(elem, dim, fe_type, e,
1101  edge_dofs);
1102 
1103  const unsigned int n_edge_dofs =
1104  cast_int<unsigned int>(edge_dofs.size());
1105 
1106  // Some edge dofs are on nodes and already
1107  // fixed, others are free to calculate
1108  unsigned int free_dofs = 0;
1109  for (unsigned int i=0; i != n_edge_dofs; ++i)
1110  if (!dof_is_fixed[edge_dofs[i]])
1111  free_dof[free_dofs++] = i;
1112 
1113  // There may be nothing to project
1114  if (!free_dofs)
1115  continue;
1116 
1117  Ke.resize (free_dofs, free_dofs); Ke.zero();
1118  Fe.resize (free_dofs); Fe.zero();
1119  // The new edge coefficients
1120  DenseVector<Number> Uedge(free_dofs);
1121 
1122  // Initialize FE data on the edge
1123  edge_fe->edge_reinit(elem, e);
1124  const unsigned int n_qp = fem_context.get_edge_qrule().n_points();
1125 
1126  // Loop over the quadrature points
1127  for (unsigned int qp=0; qp<n_qp; qp++)
1128  {
1129  // solution at the quadrature point
1130  OutputNumber fineval(0);
1131  libMesh::RawAccessor<OutputNumber> f_accessor( fineval, dim );
1132 
1133  for (unsigned int c = 0; c < n_vec_dim; c++)
1134  f_accessor(c) =
1135  f_component(f, f_fem, context.get(), var_component+c,
1136  xyz_values[qp], time);
1137 
1138  // solution grad at the quadrature point
1139  OutputNumberGradient finegrad;
1140  libMesh::RawAccessor<OutputNumberGradient> g_accessor( finegrad, dim );
1141 
1142  unsigned int g_rank;
1143  switch( FEInterface::field_type( fe_type ) )
1144  {
1145  case TYPE_SCALAR:
1146  {
1147  g_rank = 1;
1148  break;
1149  }
1150  case TYPE_VECTOR:
1151  {
1152  g_rank = 2;
1153  break;
1154  }
1155  default:
1156  libmesh_error_msg("Unknown field type!");
1157  }
1158 
1159  if (cont == C_ONE)
1160  for (unsigned int c = 0; c < n_vec_dim; c++)
1161  for (unsigned int d = 0; d < g_rank; d++)
1162  g_accessor(c + d*dim ) =
1163  g_component(g, g_fem, context.get(), var_component,
1164  xyz_values[qp], time)(c);
1165 
1166  // Form edge projection matrix
1167  for (unsigned int sidei=0, freei=0; sidei != n_edge_dofs; ++sidei)
1168  {
1169  unsigned int i = edge_dofs[sidei];
1170  // fixed DoFs aren't test functions
1171  if (dof_is_fixed[i])
1172  continue;
1173  for (unsigned int sidej=0, freej=0; sidej != n_edge_dofs; ++sidej)
1174  {
1175  unsigned int j = edge_dofs[sidej];
1176  if (dof_is_fixed[j])
1177  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1178  JxW[qp] * Ue(j);
1179  else
1180  Ke(freei,freej) += phi[i][qp] *
1181  phi[j][qp] * JxW[qp];
1182  if (cont == C_ONE)
1183  {
1184  if (dof_is_fixed[j])
1185  Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp]) ) *
1186  JxW[qp] * Ue(j);
1187  else
1188  Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1189  * JxW[qp];
1190  }
1191  if (!dof_is_fixed[j])
1192  freej++;
1193  }
1194  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
1195  if (cont == C_ONE)
1196  Fe(freei) += (finegrad.contract( (*dphi)[i][qp]) ) *
1197  JxW[qp];
1198  freei++;
1199  }
1200  }
1201 
1202  Ke.cholesky_solve(Fe, Uedge);
1203 
1204  // Transfer new edge solutions to element
1205  for (unsigned int i=0; i != free_dofs; ++i)
1206  {
1207  Number & ui = Ue(edge_dofs[free_dof[i]]);
1208  libmesh_assert(std::abs(ui) < TOLERANCE ||
1209  std::abs(ui - Uedge(i)) < TOLERANCE);
1210  ui = Uedge(i);
1211  dof_is_fixed[edge_dofs[free_dof[i]]] = true;
1212  }
1213  } // end for (e = 0..n_edges)
1214  } // end if (is_boundary_edge_it != sebi.is_boundary_edge_map.end())
1215  } // end if (dim > 2 && cont != DISCONTINUOUS)
1216 
1217  // Project any side values (edges in 2D, faces in 3D)
1218  if (dim > 1 && cont != DISCONTINUOUS)
1219  {
1220  FEGenericBase<OutputType> * side_fe = nullptr;
1221  fem_context.get_side_fe(var, side_fe);
1222 
1223  // Set tolerance on underlying FEMap object. This will allow us to
1224  // avoid spurious negative Jacobian errors while imposing BCs by
1225  // simply ignoring them. This should only be required in certain
1226  // special cases, see the DirichletBoundaries comments on this
1227  // parameter for more information.
1229 
1230  // Pre-request FE data
1231  const std::vector<std::vector<OutputShape>> & phi = side_fe->get_phi();
1232  const std::vector<Point> & xyz_values = side_fe->get_xyz();
1233  const std::vector<Real> & JxW = side_fe->get_JxW();
1234 
1235  // Only pre-request gradients for C1 projections
1236  const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1237  if ((cont == C_ONE) && (fe_type.family != SUBDIVISION))
1238  {
1239  const std::vector<std::vector<OutputGradient>> & ref_dphi = side_fe->get_dphi();
1240  dphi = &ref_dphi;
1241  }
1242 
1243  // Vector to hold side local DOF indices
1244  std::vector<unsigned int> side_dofs;
1245 
1246  // Get a reference to the "is_boundary_side" flags for the
1247  // current DirichletBoundary object. In case the map does not
1248  // contain an entry for this DirichletBoundary object, it
1249  // means there are no boundary sides active.
1250  if (const auto is_boundary_side_it = sebi.is_boundary_side_map.find(&dirichlet);
1251  is_boundary_side_it != sebi.is_boundary_side_map.end())
1252  {
1253  for (unsigned int s=0; s != sebi.n_sides; ++s)
1254  {
1255  if (!is_boundary_side_it->second[s])
1256  continue;
1257 
1258  FEInterface::dofs_on_side(elem, dim, fe_type, s,
1259  side_dofs);
1260 
1261  const unsigned int n_side_dofs =
1262  cast_int<unsigned int>(side_dofs.size());
1263 
1264  // Some side dofs are on nodes/edges and already
1265  // fixed, others are free to calculate
1266  unsigned int free_dofs = 0;
1267  for (unsigned int i=0; i != n_side_dofs; ++i)
1268  if (!dof_is_fixed[side_dofs[i]])
1269  free_dof[free_dofs++] = i;
1270 
1271  // There may be nothing to project
1272  if (!free_dofs)
1273  continue;
1274 
1275  Ke.resize (free_dofs, free_dofs); Ke.zero();
1276  Fe.resize (free_dofs); Fe.zero();
1277  // The new side coefficients
1278  DenseVector<Number> Uside(free_dofs);
1279 
1280  // Initialize FE data on the side
1281  side_fe->reinit(elem, s);
1282  const unsigned int n_qp = fem_context.get_side_qrule().n_points();
1283 
1284  // Loop over the quadrature points
1285  for (unsigned int qp=0; qp<n_qp; qp++)
1286  {
1287  // solution at the quadrature point
1288  OutputNumber fineval(0);
1289  libMesh::RawAccessor<OutputNumber> f_accessor( fineval, dim );
1290 
1291  for (unsigned int c = 0; c < n_vec_dim; c++)
1292  f_accessor(c) =
1293  f_component(f, f_fem, context.get(), var_component+c,
1294  xyz_values[qp], time);
1295 
1296  // solution grad at the quadrature point
1297  OutputNumberGradient finegrad;
1298  libMesh::RawAccessor<OutputNumberGradient> g_accessor( finegrad, dim );
1299 
1300  unsigned int g_rank;
1301  switch( FEInterface::field_type( fe_type ) )
1302  {
1303  case TYPE_SCALAR:
1304  {
1305  g_rank = 1;
1306  break;
1307  }
1308  case TYPE_VECTOR:
1309  {
1310  g_rank = 2;
1311  break;
1312  }
1313  default:
1314  libmesh_error_msg("Unknown field type!");
1315  }
1316 
1317  if (cont == C_ONE)
1318  for (unsigned int c = 0; c < n_vec_dim; c++)
1319  for (unsigned int d = 0; d < g_rank; d++)
1320  g_accessor(c + d*dim ) =
1321  g_component(g, g_fem, context.get(), var_component,
1322  xyz_values[qp], time)(c);
1323 
1324  // Form side projection matrix
1325  for (unsigned int sidei=0, freei=0; sidei != n_side_dofs; ++sidei)
1326  {
1327  unsigned int i = side_dofs[sidei];
1328  // fixed DoFs aren't test functions
1329  if (dof_is_fixed[i])
1330  continue;
1331  for (unsigned int sidej=0, freej=0; sidej != n_side_dofs; ++sidej)
1332  {
1333  unsigned int j = side_dofs[sidej];
1334  if (dof_is_fixed[j])
1335  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1336  JxW[qp] * Ue(j);
1337  else
1338  Ke(freei,freej) += phi[i][qp] *
1339  phi[j][qp] * JxW[qp];
1340  if (cont == C_ONE)
1341  {
1342  if (dof_is_fixed[j])
1343  Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp])) *
1344  JxW[qp] * Ue(j);
1345  else
1346  Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1347  * JxW[qp];
1348  }
1349  if (!dof_is_fixed[j])
1350  freej++;
1351  }
1352  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1353  if (cont == C_ONE)
1354  Fe(freei) += (finegrad.contract((*dphi)[i][qp])) *
1355  JxW[qp];
1356  freei++;
1357  }
1358  }
1359 
1360  Ke.cholesky_solve(Fe, Uside);
1361 
1362  // Transfer new side solutions to element
1363  for (unsigned int i=0; i != free_dofs; ++i)
1364  {
1365  Number & ui = Ue(side_dofs[free_dof[i]]);
1366 
1367  libmesh_assert(std::abs(ui) < TOLERANCE ||
1368  std::abs(ui - Uside(i)) < TOLERANCE);
1369  ui = Uside(i);
1370 
1371  dof_is_fixed[side_dofs[free_dof[i]]] = true;
1372  }
1373  } // end for (s = 0..n_sides)
1374  } // end if (is_boundary_side_it != sebi.is_boundary_side_map.end())
1375  } // end if (dim > 1 && cont != DISCONTINUOUS)
1376 
1377  // Project any shellface values
1378  if (dim == 2 && cont != DISCONTINUOUS)
1379  {
1380  FEGenericBase<OutputType> * fe = nullptr;
1381  fem_context.get_element_fe(var, fe, dim);
1382 
1383  // Set tolerance on underlying FEMap object. This will allow us to
1384  // avoid spurious negative Jacobian errors while imposing BCs by
1385  // simply ignoring them. This should only be required in certain
1386  // special cases, see the DirichletBoundaries comments on this
1387  // parameter for more information.
1389 
1390  // Pre-request FE data
1391  const std::vector<std::vector<OutputShape>> & phi = fe->get_phi();
1392  const std::vector<Point> & xyz_values = fe->get_xyz();
1393  const std::vector<Real> & JxW = fe->get_JxW();
1394 
1395  // Only pre-request gradients for C1 projections
1396  const std::vector<std::vector<OutputGradient>> * dphi = nullptr;
1397  if ((cont == C_ONE) && (fe_type.family != SUBDIVISION))
1398  {
1399  const std::vector<std::vector<OutputGradient>> & ref_dphi = fe->get_dphi();
1400  dphi = &ref_dphi;
1401  }
1402 
1403  // Get a reference to the "is_boundary_shellface" flags for the
1404  // current DirichletBoundary object. In case the map does not
1405  // contain an entry for this DirichletBoundary object, it
1406  // means there are no boundary shellfaces active.
1407  if (const auto is_boundary_shellface_it = sebi.is_boundary_shellface_map.find(&dirichlet);
1408  is_boundary_shellface_it != sebi.is_boundary_shellface_map.end())
1409  {
1410  for (unsigned int shellface=0; shellface != 2; ++shellface)
1411  {
1412  if (!is_boundary_shellface_it->second[shellface])
1413  continue;
1414 
1415  // A shellface has the same dof indices as the element itself
1416  std::vector<unsigned int> shellface_dofs(n_dofs);
1417  std::iota(shellface_dofs.begin(), shellface_dofs.end(), 0);
1418 
1419  // Some shellface dofs are on nodes/edges and already
1420  // fixed, others are free to calculate
1421  unsigned int free_dofs = 0;
1422  for (unsigned int i=0; i != n_dofs; ++i)
1423  if (!dof_is_fixed[shellface_dofs[i]])
1424  free_dof[free_dofs++] = i;
1425 
1426  // There may be nothing to project
1427  if (!free_dofs)
1428  continue;
1429 
1430  Ke.resize (free_dofs, free_dofs); Ke.zero();
1431  Fe.resize (free_dofs); Fe.zero();
1432  // The new shellface coefficients
1433  DenseVector<Number> Ushellface(free_dofs);
1434 
1435  // Initialize FE data on the element
1436  fe->reinit (elem);
1437  const unsigned int n_qp = fem_context.get_element_qrule().n_points();
1438 
1439  // Loop over the quadrature points
1440  for (unsigned int qp=0; qp<n_qp; qp++)
1441  {
1442  // solution at the quadrature point
1443  OutputNumber fineval(0);
1444  libMesh::RawAccessor<OutputNumber> f_accessor( fineval, dim );
1445 
1446  for (unsigned int c = 0; c < n_vec_dim; c++)
1447  f_accessor(c) =
1448  f_component(f, f_fem, context.get(), var_component+c,
1449  xyz_values[qp], time);
1450 
1451  // solution grad at the quadrature point
1452  OutputNumberGradient finegrad;
1453  libMesh::RawAccessor<OutputNumberGradient> g_accessor( finegrad, dim );
1454 
1455  unsigned int g_rank;
1456  switch( FEInterface::field_type( fe_type ) )
1457  {
1458  case TYPE_SCALAR:
1459  {
1460  g_rank = 1;
1461  break;
1462  }
1463  case TYPE_VECTOR:
1464  {
1465  g_rank = 2;
1466  break;
1467  }
1468  default:
1469  libmesh_error_msg("Unknown field type!");
1470  }
1471 
1472  if (cont == C_ONE)
1473  for (unsigned int c = 0; c < n_vec_dim; c++)
1474  for (unsigned int d = 0; d < g_rank; d++)
1475  g_accessor(c + d*dim ) =
1476  g_component(g, g_fem, context.get(), var_component,
1477  xyz_values[qp], time)(c);
1478 
1479  // Form shellface projection matrix
1480  for (unsigned int shellfacei=0, freei=0;
1481  shellfacei != n_dofs; ++shellfacei)
1482  {
1483  unsigned int i = shellface_dofs[shellfacei];
1484  // fixed DoFs aren't test functions
1485  if (dof_is_fixed[i])
1486  continue;
1487  for (unsigned int shellfacej=0, freej=0;
1488  shellfacej != n_dofs; ++shellfacej)
1489  {
1490  unsigned int j = shellface_dofs[shellfacej];
1491  if (dof_is_fixed[j])
1492  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1493  JxW[qp] * Ue(j);
1494  else
1495  Ke(freei,freej) += phi[i][qp] *
1496  phi[j][qp] * JxW[qp];
1497  if (cont == C_ONE)
1498  {
1499  if (dof_is_fixed[j])
1500  Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp])) *
1501  JxW[qp] * Ue(j);
1502  else
1503  Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1504  * JxW[qp];
1505  }
1506  if (!dof_is_fixed[j])
1507  freej++;
1508  }
1509  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1510  if (cont == C_ONE)
1511  Fe(freei) += (finegrad.contract((*dphi)[i][qp])) *
1512  JxW[qp];
1513  freei++;
1514  }
1515  }
1516 
1517  Ke.cholesky_solve(Fe, Ushellface);
1518 
1519  // Transfer new shellface solutions to element
1520  for (unsigned int i=0; i != free_dofs; ++i)
1521  {
1522  Number & ui = Ue(shellface_dofs[free_dof[i]]);
1523  libmesh_assert(std::abs(ui) < TOLERANCE ||
1524  std::abs(ui - Ushellface(i)) < TOLERANCE);
1525  ui = Ushellface(i);
1526  dof_is_fixed[shellface_dofs[free_dof[i]]] = true;
1527  }
1528  } // end for (shellface = 0..2)
1529  } // end if (is_boundary_shellface_it != sebi.is_boundary_shellface_map.end())
1530  } // end if (dim == 2 && cont != DISCONTINUOUS)
1531 
1532  // Lock the DofConstraints since it is shared among threads.
1533  {
1534  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1535 
1536  for (unsigned int i = 0; i < n_dofs; i++)
1537  {
1538  DofConstraintRow empty_row;
1539  if (dof_is_fixed[i] && !libmesh_isnan(Ue(i)))
1540  add_fn (dof_indices[i], empty_row, Ue(i));
1541  }
1542  }
1543  } // apply_dirichlet_impl
1544 
1545 public:
1546  ConstrainDirichlet (DofMap & dof_map_in,
1547  const MeshBase & mesh_in,
1548  const Real time_in,
1549  const DirichletBoundaries & dirichlets_in,
1550  const AddConstraint & add_in) :
1551  dof_map(dof_map_in),
1552  mesh(mesh_in),
1553  time(time_in),
1554  dirichlets(dirichlets_in),
1555  add_fn(add_in) { }
1556 
1557  // This class can be default copy/move constructed.
1558  ConstrainDirichlet (ConstrainDirichlet &&) = default;
1559  ConstrainDirichlet (const ConstrainDirichlet &) = default;
1560 
1561  // This class cannot be default copy/move assigned because it
1562  // contains reference members.
1563  ConstrainDirichlet & operator= (const ConstrainDirichlet &) = delete;
1564  ConstrainDirichlet & operator= (ConstrainDirichlet &&) = delete;
1565 
1566  void operator()(const ConstElemRange & range) const
1567  {
1574  // Figure out which System the DirichletBoundary objects are
1575  // defined for. We break out of the loop as soon as we encounter a
1576  // valid System pointer, the assumption is thus that all Variables
1577  // are defined on the same System.
1578  System * system = nullptr;
1579 
1580  // Map from boundary_id -> set<pair<id,DirichletBoundary*>> objects which
1581  // are active on that boundary_id. Later we will use this to determine
1582  // which DirichletBoundary objects to loop over for each Elem.
1583  std::map<boundary_id_type, std::set<std::pair<unsigned int, DirichletBoundary *>>>
1584  boundary_id_to_ordered_dirichlet_boundaries;
1585 
1586  for (auto dirichlet_id : index_range(dirichlets))
1587  {
1588  // Pointer to the current DirichletBoundary object
1589  const auto & dirichlet = dirichlets[dirichlet_id];
1590 
1591  // Construct mapping from boundary_id -> (dirichlet_id, DirichletBoundary)
1592  for (const auto & b_id : dirichlet->b)
1593  boundary_id_to_ordered_dirichlet_boundaries[b_id].emplace(dirichlet_id, dirichlet.get());
1594 
1595  for (const auto & var : dirichlet->variables)
1596  {
1597  const Variable & variable = dof_map.variable(var);
1598  auto current_system = variable.system();
1599 
1600  if (!system)
1601  system = current_system;
1602  else
1603  libmesh_error_msg_if(current_system != system,
1604  "All variables should be defined on the same System");
1605  }
1606  }
1607 
1608  // If we found no System, it could be because none of the
1609  // Variables have one defined, or because there are
1610  // DirichletBoundary objects with no Variables defined on
1611  // them. These situations both indicate a likely error in the
1612  // setup of a problem, so let's throw an error in this case.
1613  libmesh_error_msg_if(!system, "Valid System not found for any Variables.");
1614 
1615  // Construct a FEMContext object for the System on which the
1616  // Variables in our DirichletBoundary objects are defined. This
1617  // will be used in the apply_dirichlet_impl() function.
1618  // We're not going to use elem_jacobian or subjacobians here so
1619  // don't allocate them.
1620  auto fem_context = std::make_unique<FEMContext>
1621  (*system, nullptr, /* allocate local_matrices = */ false);
1622 
1623  // At the time we are using this FEMContext, the current_local_solution
1624  // vector is not initialized, but also we don't need it, so set
1625  // the algebraic_type flag to DOFS_ONLY.
1627 
1628  // Boundary info for the current mesh
1629  const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1630 
1631  // This object keeps track of the BoundaryInfo for a single Elem
1632  SingleElemBoundaryInfo sebi(boundary_info, boundary_id_to_ordered_dirichlet_boundaries);
1633 
1634  // Iterate over all the elements in the range
1635  for (const auto & elem : range)
1636  {
1637  // We only calculate Dirichlet constraints on active
1638  // elements
1639  if (!elem->active())
1640  continue;
1641 
1642  // Reinitialize BoundaryInfo data structures for the current elem
1643  bool has_dirichlet_constraint = sebi.reinit(elem);
1644 
1645  // If this Elem has no boundary ids, go to the next one.
1646  if (!has_dirichlet_constraint)
1647  continue;
1648 
1649  for (const auto & db_pair : sebi.ordered_dbs)
1650  {
1651  // Get pointer to the DirichletBoundary object
1652  const auto & dirichlet = db_pair.second;
1653 
1654  // Loop over all the variables which this DirichletBoundary object is responsible for
1655  for (const auto & var : dirichlet->variables)
1656  {
1657  const Variable & variable = dof_map.variable(var);
1658 
1659  // Make sure that the Variable and the DofMap agree on
1660  // what number this variable is.
1661  libmesh_assert_equal_to(variable.number(), var);
1662 
1663  const FEType & fe_type = variable.type();
1664 
1665  if (fe_type.family == SCALAR)
1666  continue;
1667 
1668  switch( FEInterface::field_type( fe_type ) )
1669  {
1670  case TYPE_SCALAR:
1671  {
1672  // For Lagrange FEs we don't need to do a full
1673  // blown projection, we can just interpolate
1674  // values directly.
1675  if (fe_type.family == LAGRANGE)
1676  this->apply_lagrange_dirichlet_impl<Real>( sebi, variable, *dirichlet, *fem_context );
1677  else
1678  this->apply_dirichlet_impl<Real>( sebi, variable, *dirichlet, *fem_context );
1679  break;
1680  }
1681  case TYPE_VECTOR:
1682  {
1683  this->apply_dirichlet_impl<RealGradient>( sebi, variable, *dirichlet, *fem_context );
1684  break;
1685  }
1686  default:
1687  libmesh_error_msg("Unknown field type!");
1688  }
1689  } // for (var : variables)
1690  } // for (db_pair : ordered_dbs)
1691  } // for (elem : range)
1692  } // operator()
1693 
1694 }; // class ConstrainDirichlet
1695 
1696 
1697 #endif // LIBMESH_ENABLE_DIRICHLET
1698 
1699 
1700 } // anonymous namespace
1701 
1702 
1703 
1704 namespace libMesh
1705 {
1706 
1707 // ------------------------------------------------------------
1708 // DofMap member functions
1709 
1710 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1711 
1712 
1714 {
1715  parallel_object_only();
1716 
1717  dof_id_type nc_dofs = this->n_local_constrained_dofs();
1718  this->comm().sum(nc_dofs);
1719  return nc_dofs;
1720 }
1721 
1722 
1724 {
1725  const DofConstraints::const_iterator lower =
1726  _dof_constraints.lower_bound(this->first_dof()),
1727  upper =
1728  _dof_constraints.lower_bound(this->end_dof());
1729 
1730  return cast_int<dof_id_type>(std::distance(lower, upper));
1731 }
1732 
1733 
1734 
1736 {
1737  parallel_object_only();
1738 
1739  LOG_SCOPE("create_dof_constraints()", "DofMap");
1740 
1742 
1743  // The user might have set boundary conditions after the mesh was
1744  // prepared; we should double-check that those boundary conditions
1745  // are still consistent.
1746 #ifdef DEBUG
1748 #endif
1749 
1750  // In a distributed mesh we might have constraint rows on some
1751  // processors but not all; if we have constraint rows on *any*
1752  // processor then we need to process them.
1753  bool constraint_rows_empty = mesh.get_constraint_rows().empty();
1754  this->comm().min(constraint_rows_empty);
1755 
1756  // We might get constraint equations from AMR hanging nodes in
1757  // 2D/3D, or from spline constraint rows or boundary conditions in
1758  // any dimension
1759  const bool possible_local_constraints = false
1760  || !mesh.n_elem()
1761  || !constraint_rows_empty
1762 #ifdef LIBMESH_ENABLE_AMR
1763  || mesh.mesh_dimension() > 1
1764 #endif
1765 #ifdef LIBMESH_ENABLE_PERIODIC
1766  || !_periodic_boundaries->empty()
1767 #endif
1768 #ifdef LIBMESH_ENABLE_DIRICHLET
1769  || !_dirichlet_boundaries->empty()
1770 #endif
1771  ;
1772 
1773  // Even if we don't have constraints, another processor might.
1774  bool possible_global_constraints = possible_local_constraints;
1775 #if defined(LIBMESH_ENABLE_PERIODIC) || defined(LIBMESH_ENABLE_DIRICHLET) || defined(LIBMESH_ENABLE_AMR)
1776  libmesh_assert(this->comm().verify(mesh.is_serial()));
1777 
1778  this->comm().max(possible_global_constraints);
1779 #endif
1780 
1781  // Recalculate dof constraints from scratch. (Or just clear them,
1782  // if the user has just deleted their last dirichlet/periodic/user
1783  // constraint)
1784  // Note: any _stashed_dof_constraints are not cleared as it
1785  // may be the user's intention to restore them later.
1786 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1787  _dof_constraints.clear();
1788  _primal_constraint_values.clear();
1790 #endif
1791 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1792  _node_constraints.clear();
1793 #endif
1794 
1795  if (!possible_global_constraints)
1796  return;
1797 
1798  // Here we build the hanging node constraints. This is done
1799  // by enforcing the condition u_a = u_b along hanging sides.
1800  // u_a = u_b is collocated at the nodes of side a, which gives
1801  // one row of the constraint matrix.
1802 
1803  // Processors only compute their local constraints
1804  ConstElemRange range (mesh.local_elements_begin(),
1805  mesh.local_elements_end());
1806 
1807  // Global computation fails if we're using a FEMFunctionBase BC on a
1808  // ReplicatedMesh in parallel
1809  // ConstElemRange range (mesh.elements_begin(),
1810  // mesh.elements_end());
1811 
1812  // compute_periodic_constraints requires a point_locator() from our
1813  // Mesh, but point_locator() construction is parallel and threaded.
1814  // Rather than nest threads within threads we'll make sure it's
1815  // preconstructed.
1816 #ifdef LIBMESH_ENABLE_PERIODIC
1817  bool need_point_locator = !_periodic_boundaries->empty() && !range.empty();
1818 
1819  this->comm().max(need_point_locator);
1820 
1821  if (need_point_locator)
1823 #endif
1824 
1825 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1826  Threads::parallel_for (range,
1827  ComputeNodeConstraints (_node_constraints,
1828 #ifdef LIBMESH_ENABLE_PERIODIC
1830 #endif
1831  mesh));
1832 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1833 
1834 
1835  // Look at all the variables in the system. Reset the element
1836  // range at each iteration -- there is no need to reconstruct it.
1837  const auto n_vars = this->n_variables();
1838  for (unsigned int variable_number=0; variable_number<n_vars;
1839  ++variable_number, range.reset())
1840  Threads::parallel_for (range,
1841  ComputeConstraints (_dof_constraints,
1842  *this,
1843 #ifdef LIBMESH_ENABLE_PERIODIC
1845 #endif
1846  mesh,
1847  variable_number));
1848 
1849 #ifdef LIBMESH_ENABLE_DIRICHLET
1850 
1851  if (!_dirichlet_boundaries->empty())
1852  {
1853  // Sanity check that the boundary ids associated with the
1854  // DirichletBoundary objects are actually present in the
1855  // mesh. We do this check by default, but in cases where you
1856  // intentionally add "inconsistent but valid" DirichletBoundary
1857  // objects in parallel, this check can deadlock since it does a
1858  // collective communication internally. In that case it is
1859  // possible to disable this check by setting the flag to false.
1861  for (const auto & dirichlet : *_dirichlet_boundaries)
1862  this->check_dirichlet_bcid_consistency(mesh, *dirichlet);
1863 
1864  // Threaded loop over local over elems applying all Dirichlet BCs
1866  (range,
1867  ConstrainDirichlet(*this, mesh, time, *_dirichlet_boundaries,
1868  AddPrimalConstraint(*this)));
1869 
1870  // Threaded loop over local over elems per QOI applying all adjoint
1871  // Dirichlet BCs. Note that the ConstElemRange is reset before each
1872  // execution of Threads::parallel_for().
1873 
1874  for (auto qoi_index : index_range(_adjoint_dirichlet_boundaries))
1875  {
1876  const DirichletBoundaries & adb_q =
1877  *(_adjoint_dirichlet_boundaries[qoi_index]);
1878 
1879  if (!adb_q.empty())
1881  (range.reset(),
1882  ConstrainDirichlet(*this, mesh, time, adb_q,
1883  AddAdjointConstraint(*this, qoi_index)));
1884  }
1885  }
1886 
1887 #endif // LIBMESH_ENABLE_DIRICHLET
1888 
1889  // Handle spline node constraints last, so we can try to move
1890  // existing constraints onto the spline basis if necessary.
1891  if (!constraint_rows_empty)
1892  this->process_mesh_constraint_rows(mesh);
1893 }
1894 
1895 
1896 
1898 {
1899  // If we already have simple Dirichlet constraints (with right hand
1900  // sides but with no coupling between DoFs) on spline-constrained FE
1901  // nodes, then we'll need a solve to compute the corresponding
1902  // constraints on the relevant spline nodes. (If we already have
1903  // constraints with coupling between DoFs on spline-constrained FE
1904  // nodes, then we'll need to go sit down and cry until we figure out
1905  // how to handle that.)
1906 
1907  const auto & constraint_rows = mesh.get_constraint_rows();
1908 
1909  // This routine is too expensive to use unless we really might
1910  // need it
1911 #ifdef DEBUG
1912  bool constraint_rows_empty = constraint_rows.empty();
1913  this->comm().min(constraint_rows_empty);
1914  libmesh_assert(!constraint_rows_empty);
1915 #endif
1916 
1917  // We can't handle periodic boundary conditions on spline meshes
1918  // yet.
1919 #ifdef LIBMESH_ENABLE_PERIODIC
1920  libmesh_error_msg_if (!_periodic_boundaries->empty(),
1921  "Periodic boundary conditions are not yet implemented for spline meshes");
1922 #endif
1923 
1924  // We can handle existing Dirichlet constraints, but we'll need
1925  // to do solves to project them down onto the spline basis.
1926  std::unique_ptr<SparsityPattern::Build> sp;
1927  std::unique_ptr<SparseMatrix<Number>> mat;
1928 
1929  const unsigned int n_adjoint_rhs =
1931 
1932  // [0] for primal rhs, [q+1] for adjoint qoi q
1933  std::vector<std::unique_ptr<NumericVector<Number>>>
1934  solve_rhs(n_adjoint_rhs+1);
1935 
1936  // Keep track of which spline DoFs will be Dirichlet.
1937  // We use a set here to make it easier to find what processors
1938  // to send the DoFs to later.
1939  std::set<dof_id_type> my_dirichlet_spline_dofs;
1940 
1941  // And keep track of which non-spline Dofs were Dirichlet
1942  std::unordered_set<dof_id_type> was_previously_constrained;
1943 
1944  const unsigned int sys_num = this->sys_number();
1945  for (auto & node_row : constraint_rows)
1946  {
1947  const Node * node = node_row.first;
1948  libmesh_assert(node == mesh.node_ptr(node->id()));
1949 
1950  // Each processor only computes its own (and in distributed
1951  // cases, is only guaranteed to have the dependency data to
1952  // compute its own) constraints here.
1953  if (node->processor_id() != mesh.processor_id())
1954  continue;
1955 
1956  for (auto var_num : IntRange<unsigned int>(0, this->n_variables()))
1957  {
1958  const FEFamily & fe_family = this->variable_type(var_num).family;
1959 
1960  // constraint_rows only applies to nodal variables
1961  if (fe_family != LAGRANGE &&
1962  fe_family != RATIONAL_BERNSTEIN)
1963  continue;
1964 
1965  DofConstraintRow dc_row;
1966 
1967  const dof_id_type constrained_id =
1968  node->dof_number(sys_num, var_num, 0);
1969  for (const auto & [pr, val] : node_row.second)
1970  {
1971  const Elem * spline_elem = pr.first;
1972  libmesh_assert(spline_elem == mesh.elem_ptr(spline_elem->id()));
1973 
1974  const Node & spline_node =
1975  spline_elem->node_ref(pr.second);
1976 
1977  const dof_id_type spline_dof_id =
1978  spline_node.dof_number(sys_num, var_num, 0);
1979  dc_row[spline_dof_id] = val;
1980  }
1981 
1982  // See if we already have a constraint here.
1983  if (this->is_constrained_dof(constrained_id))
1984  {
1985  was_previously_constrained.insert(constrained_id);
1986 
1987  // Keep track of which spline DoFs will be
1988  // inheriting this non-spline DoF's constraints
1989  for (auto & row_entry : dc_row)
1990  my_dirichlet_spline_dofs.insert(row_entry.first);
1991 
1992  // If it wasn't a simple Dirichlet-type constraint
1993  // then I don't know what to do with it. We'll make
1994  // this an assertion only because this should only
1995  // crop up with periodic boundary conditions, which
1996  // we've already made sure we don't have.
1997  libmesh_assert(_dof_constraints[constrained_id].empty());
1998  }
1999 
2000  // Add the constraint, replacing any previous, so we can
2001  // use the new constraint in setting up the solve below
2002  this->add_constraint_row(constrained_id, dc_row, false);
2003  }
2004  }
2005 
2006  // my_dirichlet_spline_dofs may now include DoFs whose owners
2007  // don't know they need to become spline DoFs! We need to push
2008  // this data to them.
2009  if (this->comm().size() > 1)
2010  {
2011  std::unordered_map
2012  <processor_id_type, std::vector<dof_id_type>>
2013  their_dirichlet_spline_dofs;
2014 
2015  // If we ever change the underlying container here then we'd
2016  // better do some kind of sort before using it; we'll rely
2017  // on sorting to make the processor id lookup efficient.
2018  libmesh_assert(std::is_sorted(my_dirichlet_spline_dofs.begin(),
2019  my_dirichlet_spline_dofs.end()));
2020  processor_id_type destination_pid = 0;
2021  for (auto d : my_dirichlet_spline_dofs)
2022  {
2023  libmesh_assert_less(d, this->end_dof(this->comm().size()-1));
2024  while (d >= this->end_dof(destination_pid))
2025  destination_pid++;
2026 
2027  if (destination_pid != this->processor_id())
2028  their_dirichlet_spline_dofs[destination_pid].push_back(d);
2029  }
2030 
2031  auto receive_dof_functor =
2032  [& my_dirichlet_spline_dofs]
2034  const std::vector<dof_id_type> & dofs)
2035  {
2036  my_dirichlet_spline_dofs.insert(dofs.begin(), dofs.end());
2037  };
2038 
2039  Parallel::push_parallel_vector_data
2040  (this->comm(), their_dirichlet_spline_dofs, receive_dof_functor);
2041  }
2042 
2043 
2044  // If anyone had any prior constraints in effect, then we need
2045  // to convert them to constraints on the spline nodes.
2046  //
2047  // NOT simply testing prior_constraints here; maybe it turned
2048  // out that all our constraints were on non-spline-constrained
2049  // parts of a hybrid mesh?
2050  bool important_prior_constraints =
2051  !was_previously_constrained.empty();
2052  this->comm().max(important_prior_constraints);
2053 
2054  if (important_prior_constraints)
2055  {
2056  // Now that we have the spline constraints added, we can
2057  // finally construct a sparsity pattern that correctly
2058  // accounts for those constraints!
2059  mat = SparseMatrix<Number>::build(this->comm());
2060  for (auto q : IntRange<unsigned int>(0, n_adjoint_rhs+1))
2061  {
2062  solve_rhs[q] = NumericVector<Number>::build(this->comm());
2063  solve_rhs[q]->init(this->n_dofs(), this->n_local_dofs(),
2064  false, PARALLEL);
2065  }
2066 
2067  // We need to compute our own sparsity pattern, to take into
2068  // account the particularly non-sparse rows that can be
2069  // created by the spline constraints we just added.
2070  mat->attach_dof_map(*this);
2071  sp = this->build_sparsity(mesh);
2072  mat->attach_sparsity_pattern(*sp);
2073  mat->init();
2074 
2075  for (auto & node_row : constraint_rows)
2076  {
2077  const Node * node = node_row.first;
2078  libmesh_assert(node == mesh.node_ptr(node->id()));
2079 
2080  for (auto var_num : IntRange<unsigned int>(0, this->n_variables()))
2081  {
2082  const FEFamily & fe_family = this->variable_type(var_num).family;
2083 
2084  // constraint_rows only applies to nodal variables
2085  if (fe_family != LAGRANGE &&
2086  fe_family != RATIONAL_BERNSTEIN)
2087  continue;
2088 
2089  const dof_id_type constrained_id =
2090  node->dof_number(sys_num, var_num, 0);
2091 
2092  if (was_previously_constrained.count(constrained_id))
2093  {
2094  for (auto q : IntRange<int>(0, n_adjoint_rhs+1))
2095  {
2096  DenseMatrix<Number> K(1,1);
2097  DenseVector<Number> F(1);
2098  std::vector<dof_id_type> dof_indices(1, constrained_id);
2099 
2100  K(0,0) = 1;
2101 
2102  DofConstraintValueMap & vals = q ?
2105 
2106  DofConstraintValueMap::const_iterator rhsit =
2107  vals.find(constrained_id);
2108  F(0) = (rhsit == vals.end()) ? 0 : rhsit->second;
2109 
2110  // We no longer need any rhs values here directly.
2111  if (rhsit != vals.end())
2112  vals.erase(rhsit);
2113 
2115  (K, F, dof_indices, false, q ? (q-1) : -1);
2116  if (!q)
2117  mat->add_matrix(K, dof_indices);
2118  solve_rhs[q]->add_vector(F, dof_indices);
2119  }
2120  }
2121  }
2122  }
2123 
2124  // Any DoFs that aren't part of any constraint, directly or
2125  // indirectly, need a diagonal term to make the matrix
2126  // here invertible.
2127  for (dof_id_type d : IntRange<dof_id_type>(this->first_dof(),
2128  this->end_dof()))
2129  if (!was_previously_constrained.count(d) &&
2130  !my_dirichlet_spline_dofs.count(d))
2131  mat->add(d,d,1);
2132 
2133  // At this point, we're finally ready to solve for Dirichlet
2134  // constraint values on spline nodes.
2135  std::unique_ptr<LinearSolver<Number>> linear_solver =
2137 
2138  std::unique_ptr<NumericVector<Number>> projected_vals =
2140 
2141  projected_vals->init(this->n_dofs(), this->n_local_dofs(),
2142  false, PARALLEL);
2143 
2144  DofConstraintRow empty_row;
2145  for (auto sd : my_dirichlet_spline_dofs)
2146  if (this->local_index(sd))
2147  this->add_constraint_row(sd, empty_row);
2148 
2149  for (auto q : IntRange<unsigned int>(0, n_adjoint_rhs+1))
2150  {
2151  // FIXME: we don't have an EquationSystems here, but I'd
2152  // rather not hardcode these...
2153  const double tol = double(TOLERANCE * TOLERANCE);
2154  const unsigned int max_its = 5000;
2155 
2156  linear_solver->solve(*mat, *projected_vals,
2157  *(solve_rhs[q]), tol, max_its);
2158 
2159  DofConstraintValueMap & vals = q ?
2162 
2163  for (auto sd : my_dirichlet_spline_dofs)
2164  if (this->local_index(sd))
2165  {
2166  Number constraint_rhs = (*projected_vals)(sd);
2167 
2168  std::pair<DofConstraintValueMap::iterator, bool> rhs_it =
2169  vals.emplace(sd, constraint_rhs);
2170  if (!rhs_it.second)
2171  rhs_it.first->second = constraint_rhs;
2172  }
2173  }
2174  }
2175 }
2176 
2177 
2178 
2180  const DofConstraintRow & constraint_row,
2181  const Number constraint_rhs,
2182  const bool forbid_constraint_overwrite)
2183 {
2184  // Optionally allow the user to overwrite constraints. Defaults to false.
2185  libmesh_error_msg_if(forbid_constraint_overwrite && this->is_constrained_dof(dof_number),
2186  "ERROR: DOF " << dof_number << " was already constrained!");
2187 
2188  libmesh_assert_less(dof_number, this->n_dofs());
2189 
2190  // There is an implied "1" on the diagonal of the constraint row, and the user
2191  // should not try to manually set _any_ value on the diagonal.
2192  libmesh_assert_msg(!constraint_row.count(dof_number),
2193  "Error: constraint_row for dof " << dof_number <<
2194  " should not contain an entry for dof " << dof_number);
2195 
2196 #ifndef NDEBUG
2197  for (const auto & pr : constraint_row)
2198  libmesh_assert_less(pr.first, this->n_dofs());
2199 #endif
2200 
2201  // Store the constraint_row in the map
2202  _dof_constraints.insert_or_assign(dof_number, constraint_row);
2203 
2204  std::pair<DofConstraintValueMap::iterator, bool> rhs_it =
2205  _primal_constraint_values.emplace(dof_number, constraint_rhs);
2206  if (!rhs_it.second)
2207  rhs_it.first->second = constraint_rhs;
2208 }
2209 
2210 
2211 void DofMap::add_adjoint_constraint_row (const unsigned int qoi_index,
2212  const dof_id_type dof_number,
2213  const DofConstraintRow & /*constraint_row*/,
2214  const Number constraint_rhs,
2215  const bool forbid_constraint_overwrite)
2216 {
2217  // Optionally allow the user to overwrite constraints. Defaults to false.
2218  if (forbid_constraint_overwrite)
2219  {
2220  libmesh_error_msg_if(!this->is_constrained_dof(dof_number),
2221  "ERROR: DOF " << dof_number << " has no corresponding primal constraint!");
2222 #ifndef NDEBUG
2223  // No way to do this without a non-normalized tolerance?
2224 
2225  // // If the user passed in more than just the rhs, let's check the
2226  // // coefficients for consistency
2227  // if (!constraint_row.empty())
2228  // {
2229  // DofConstraintRow row = _dof_constraints[dof_number];
2230  // for (const auto & [dof, val] : row)
2231  // libmesh_assert(constraint_row.find(dof)->second == val);
2232  // }
2233  //
2234  // if (_adjoint_constraint_values[qoi_index].find(dof_number) !=
2235  // _adjoint_constraint_values[qoi_index].end())
2236  // libmesh_assert_equal_to(_adjoint_constraint_values[qoi_index][dof_number],
2237  // constraint_rhs);
2238 
2239 #endif
2240  }
2241 
2242  // Creates the map of rhs values if it doesn't already exist; then
2243  // adds the current value to that map
2244 
2245  // Store the rhs value in the map
2246  _adjoint_constraint_values[qoi_index].insert_or_assign(dof_number, constraint_rhs);
2247 }
2248 
2249 
2250 
2251 
2252 void DofMap::print_dof_constraints(std::ostream & os,
2253  bool print_nonlocal) const
2254 {
2255  parallel_object_only();
2256 
2257  std::string local_constraints =
2258  this->get_local_constraints(print_nonlocal);
2259 
2260  if (this->processor_id())
2261  {
2262  this->comm().send(0, local_constraints);
2263  }
2264  else
2265  {
2266  os << "Processor 0:\n";
2267  os << local_constraints;
2268 
2269  for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
2270  {
2271  this->comm().receive(p, local_constraints);
2272  os << "Processor " << p << ":\n";
2273  os << local_constraints;
2274  }
2275  }
2276 }
2277 
2278 std::string DofMap::get_local_constraints(bool print_nonlocal) const
2279 {
2280  std::ostringstream os;
2281 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2282  if (print_nonlocal)
2283  os << "All ";
2284  else
2285  os << "Local ";
2286 
2287  os << "Node Constraints:"
2288  << std::endl;
2289 
2290  for (const auto & [node, pr] : _node_constraints)
2291  {
2292  // Skip non-local nodes if requested
2293  if (!print_nonlocal &&
2294  node->processor_id() != this->processor_id())
2295  continue;
2296 
2297  const NodeConstraintRow & row = pr.first;
2298  const Point & offset = pr.second;
2299 
2300  os << "Constraints for Node id " << node->id()
2301  << ": \t";
2302 
2303  for (const auto & [cnode, val] : row)
2304  os << " (" << cnode->id() << "," << val << ")\t";
2305 
2306  os << "rhs: " << offset;
2307 
2308  os << std::endl;
2309  }
2310 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2311 
2312  if (print_nonlocal)
2313  os << "All ";
2314  else
2315  os << "Local ";
2316 
2317  os << "DoF Constraints:"
2318  << std::endl;
2319 
2320  for (const auto & [i, row] : _dof_constraints)
2321  {
2322  // Skip non-local dofs if requested
2323  if (!print_nonlocal && !this->local_index(i))
2324  continue;
2325 
2326  DofConstraintValueMap::const_iterator rhsit =
2327  _primal_constraint_values.find(i);
2328  const Number rhs = (rhsit == _primal_constraint_values.end()) ?
2329  0 : rhsit->second;
2330 
2331  os << "Constraints for DoF " << i
2332  << ": \t";
2333 
2334  for (const auto & item : row)
2335  os << " (" << item.first << "," << item.second << ")\t";
2336 
2337  os << "rhs: " << rhs;
2338  os << std::endl;
2339  }
2340 
2341  for (unsigned int qoi_index = 0,
2342  n_qois = cast_int<unsigned int>(_adjoint_dirichlet_boundaries.size());
2343  qoi_index != n_qois; ++qoi_index)
2344  {
2345  os << "Adjoint " << qoi_index << " DoF rhs values:"
2346  << std::endl;
2347 
2348  if (auto adjoint_map_it = _adjoint_constraint_values.find(qoi_index);
2349  adjoint_map_it != _adjoint_constraint_values.end())
2350  for (const auto & [i, rhs] : adjoint_map_it->second)
2351  {
2352  // Skip non-local dofs if requested
2353  if (!print_nonlocal && !this->local_index(i))
2354  continue;
2355 
2356  os << "RHS for DoF " << i
2357  << ": " << rhs;
2358 
2359  os << std::endl;
2360  }
2361  }
2362 
2363  return os.str();
2364 }
2365 
2366 
2367 
2369  std::vector<dof_id_type> & elem_dofs,
2370  bool asymmetric_constraint_rows) const
2371 {
2372  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2373  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2374 
2375  // check for easy return
2376  if (this->_dof_constraints.empty())
2377  return;
2378 
2379  // The constrained matrix is built up as C^T K C.
2381 
2382 
2383  this->build_constraint_matrix (C, elem_dofs);
2384 
2385  LOG_SCOPE("constrain_elem_matrix()", "DofMap");
2386 
2387  // It is possible that the matrix is not constrained at all.
2388  if ((C.m() == matrix.m()) &&
2389  (C.n() == elem_dofs.size())) // It the matrix is constrained
2390  {
2391  // Compute the matrix-matrix-matrix product C^T K C
2392  matrix.left_multiply_transpose (C);
2393  matrix.right_multiply (C);
2394 
2395 
2396  libmesh_assert_equal_to (matrix.m(), matrix.n());
2397  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
2398  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
2399 
2400 
2401  for (unsigned int i=0,
2402  n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2403  i != n_elem_dofs; i++)
2404  // If the DOF is constrained
2405  if (this->is_constrained_dof(elem_dofs[i]))
2406  {
2407  for (auto j : make_range(matrix.n()))
2408  matrix(i,j) = 0.;
2409 
2410  matrix(i,i) = 1.;
2411 
2412  if (asymmetric_constraint_rows)
2413  {
2414  DofConstraints::const_iterator
2415  pos = _dof_constraints.find(elem_dofs[i]);
2416 
2417  libmesh_assert (pos != _dof_constraints.end());
2418 
2419  const DofConstraintRow & constraint_row = pos->second;
2420 
2421  // This is an overzealous assertion in the presence of
2422  // heterogeneous constraints: we now can constrain "u_i = c"
2423  // with no other u_j terms involved.
2424  //
2425  // libmesh_assert (!constraint_row.empty());
2426 
2427  for (const auto & item : constraint_row)
2428  for (unsigned int j=0; j != n_elem_dofs; j++)
2429  if (elem_dofs[j] == item.first)
2430  matrix(i,j) = -item.second;
2431  }
2432  }
2433  } // end if is constrained...
2434 }
2435 
2436 
2437 
2439  DenseVector<Number> & rhs,
2440  std::vector<dof_id_type> & elem_dofs,
2441  bool asymmetric_constraint_rows) const
2442 {
2443  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2444  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2445  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2446 
2447  // check for easy return
2448  if (this->_dof_constraints.empty())
2449  return;
2450 
2451  // The constrained matrix is built up as C^T K C.
2452  // The constrained RHS is built up as C^T F
2454 
2455  this->build_constraint_matrix (C, elem_dofs);
2456 
2457  LOG_SCOPE("cnstrn_elem_mat_vec()", "DofMap");
2458 
2459  // It is possible that the matrix is not constrained at all.
2460  if ((C.m() == matrix.m()) &&
2461  (C.n() == elem_dofs.size())) // It the matrix is constrained
2462  {
2463  // Compute the matrix-matrix-matrix product C^T K C
2464  matrix.left_multiply_transpose (C);
2465  matrix.right_multiply (C);
2466 
2467 
2468  libmesh_assert_equal_to (matrix.m(), matrix.n());
2469  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
2470  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
2471 
2472 
2473  for (unsigned int i=0,
2474  n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2475  i != n_elem_dofs; i++)
2476  if (this->is_constrained_dof(elem_dofs[i]))
2477  {
2478  for (auto j : make_range(matrix.n()))
2479  matrix(i,j) = 0.;
2480 
2481  // If the DOF is constrained
2482  matrix(i,i) = 1.;
2483 
2484  // This will put a nonsymmetric entry in the constraint
2485  // row to ensure that the linear system produces the
2486  // correct value for the constrained DOF.
2487  if (asymmetric_constraint_rows)
2488  {
2489  DofConstraints::const_iterator
2490  pos = _dof_constraints.find(elem_dofs[i]);
2491 
2492  libmesh_assert (pos != _dof_constraints.end());
2493 
2494  const DofConstraintRow & constraint_row = pos->second;
2495 
2496  // p refinement creates empty constraint rows
2497  // libmesh_assert (!constraint_row.empty());
2498 
2499  for (const auto & item : constraint_row)
2500  for (unsigned int j=0; j != n_elem_dofs; j++)
2501  if (elem_dofs[j] == item.first)
2502  matrix(i,j) = -item.second;
2503  }
2504  }
2505 
2506 
2507  // Compute the matrix-vector product C^T F
2508  DenseVector<Number> old_rhs(rhs);
2509 
2510  // compute matrix/vector product
2511  C.vector_mult_transpose(rhs, old_rhs);
2512  } // end if is constrained...
2513 }
2514 
2515 
2516 
2518  DenseVector<Number> & rhs,
2519  std::vector<dof_id_type> & elem_dofs,
2520  bool asymmetric_constraint_rows,
2521  int qoi_index) const
2522 {
2523  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2524  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2525  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2526 
2527  // check for easy return
2528  if (this->_dof_constraints.empty())
2529  return;
2530 
2531  // The constrained matrix is built up as C^T K C.
2532  // The constrained RHS is built up as C^T (F - K H)
2535 
2536  this->build_constraint_matrix_and_vector (C, H, elem_dofs, qoi_index);
2537 
2538  LOG_SCOPE("hetero_cnstrn_elem_mat_vec()", "DofMap");
2539 
2540  // It is possible that the matrix is not constrained at all.
2541  if ((C.m() == matrix.m()) &&
2542  (C.n() == elem_dofs.size())) // It the matrix is constrained
2543  {
2544  // We may have rhs values to use later
2545  const DofConstraintValueMap * rhs_values = nullptr;
2546  if (qoi_index < 0)
2547  rhs_values = &_primal_constraint_values;
2548  else if (auto it = _adjoint_constraint_values.find(qoi_index);
2549  it != _adjoint_constraint_values.end())
2550  rhs_values = &it->second;
2551 
2552  // Compute matrix/vector product K H
2554  matrix.vector_mult(KH, H);
2555 
2556  // Compute the matrix-vector product C^T (F - KH)
2557  DenseVector<Number> F_minus_KH(rhs);
2558  F_minus_KH -= KH;
2559  C.vector_mult_transpose(rhs, F_minus_KH);
2560 
2561  // Compute the matrix-matrix-matrix product C^T K C
2562  matrix.left_multiply_transpose (C);
2563  matrix.right_multiply (C);
2564 
2565  libmesh_assert_equal_to (matrix.m(), matrix.n());
2566  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
2567  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
2568 
2569  for (unsigned int i=0,
2570  n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2571  i != n_elem_dofs; i++)
2572  {
2573  const dof_id_type dof_id = elem_dofs[i];
2574 
2575  if (this->is_constrained_dof(dof_id))
2576  {
2577  for (auto j : make_range(matrix.n()))
2578  matrix(i,j) = 0.;
2579 
2580  // If the DOF is constrained
2581  matrix(i,i) = 1.;
2582 
2583  // This will put a nonsymmetric entry in the constraint
2584  // row to ensure that the linear system produces the
2585  // correct value for the constrained DOF.
2586  if (asymmetric_constraint_rows)
2587  {
2588  DofConstraints::const_iterator
2589  pos = _dof_constraints.find(dof_id);
2590 
2591  libmesh_assert (pos != _dof_constraints.end());
2592 
2593  const DofConstraintRow & constraint_row = pos->second;
2594 
2595  for (const auto & item : constraint_row)
2596  for (unsigned int j=0; j != n_elem_dofs; j++)
2597  if (elem_dofs[j] == item.first)
2598  matrix(i,j) = -item.second;
2599 
2600  if (rhs_values)
2601  {
2602  const DofConstraintValueMap::const_iterator valpos =
2603  rhs_values->find(dof_id);
2604 
2605  rhs(i) = (valpos == rhs_values->end()) ?
2606  0 : valpos->second;
2607  }
2608  }
2609  else
2610  rhs(i) = 0.;
2611  }
2612  }
2613 
2614  } // end if is constrained...
2615 }
2616 
2617 
2620  DenseVector<Number> & rhs,
2621  std::vector<dof_id_type> & elem_dofs,
2622  NumericVector<Number> & solution_local) const
2623 {
2624  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2625  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2626  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2627 
2628  libmesh_assert (solution_local.type() == SERIAL ||
2629  solution_local.type() == GHOSTED);
2630 
2631  // check for easy return
2632  if (this->_dof_constraints.empty())
2633  return;
2634 
2635  // The constrained matrix is built up as C^T K C.
2636  // The constrained RHS is built up as C^T F
2637  // Asymmetric residual terms are added if we do not have x = Cx+h
2640 
2641  this->build_constraint_matrix_and_vector (C, H, elem_dofs);
2642 
2643  LOG_SCOPE("hetero_cnstrn_elem_jac_res()", "DofMap");
2644 
2645  // It is possible that the matrix is not constrained at all.
2646  if ((C.m() != matrix.m()) ||
2647  (C.n() != elem_dofs.size()))
2648  return;
2649 
2650  // Compute the matrix-vector product C^T F
2651  DenseVector<Number> old_rhs(rhs);
2652  C.vector_mult_transpose(rhs, old_rhs);
2653 
2654  // Compute the matrix-matrix-matrix product C^T K C
2655  matrix.left_multiply_transpose (C);
2656  matrix.right_multiply (C);
2657 
2658  libmesh_assert_equal_to (matrix.m(), matrix.n());
2659  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
2660  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
2661 
2662  for (unsigned int i=0,
2663  n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2664  i != n_elem_dofs; i++)
2665  {
2666  const dof_id_type dof_id = elem_dofs[i];
2667 
2668  if (auto pos = _dof_constraints.find(dof_id);
2669  pos != _dof_constraints.end())
2670  {
2671  for (auto j : make_range(matrix.n()))
2672  matrix(i,j) = 0.;
2673 
2674  // If the DOF is constrained
2675  matrix(i,i) = 1.;
2676 
2677  // This will put a nonsymmetric entry in the constraint
2678  // row to ensure that the linear system produces the
2679  // correct value for the constrained DOF.
2680  const DofConstraintRow & constraint_row = pos->second;
2681 
2682  for (const auto & item : constraint_row)
2683  for (unsigned int j=0; j != n_elem_dofs; j++)
2684  if (elem_dofs[j] == item.first)
2685  matrix(i,j) = -item.second;
2686 
2687  const DofConstraintValueMap::const_iterator valpos =
2688  _primal_constraint_values.find(dof_id);
2689 
2690  Number & rhs_val = rhs(i);
2691  rhs_val = (valpos == _primal_constraint_values.end()) ?
2692  0 : -valpos->second;
2693  for (const auto & [constraining_dof, coef] : constraint_row)
2694  rhs_val -= coef * solution_local(constraining_dof);
2695  rhs_val += solution_local(dof_id);
2696  }
2697  }
2698 }
2699 
2700 
2703  std::vector<dof_id_type> & elem_dofs,
2704  NumericVector<Number> & solution_local) const
2705 {
2706  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2707 
2708  libmesh_assert (solution_local.type() == SERIAL ||
2709  solution_local.type() == GHOSTED);
2710 
2711  // check for easy return
2712  if (this->_dof_constraints.empty())
2713  return;
2714 
2715  // The constrained RHS is built up as C^T F
2716  // Asymmetric residual terms are added if we do not have x = Cx+h
2719 
2720  this->build_constraint_matrix_and_vector (C, H, elem_dofs);
2721 
2722  LOG_SCOPE("hetero_cnstrn_elem_res()", "DofMap");
2723 
2724  // It is possible that the element is not constrained at all.
2725  if ((C.m() != rhs.size()) ||
2726  (C.n() != elem_dofs.size()))
2727  return;
2728 
2729  // Compute the matrix-vector product C^T F
2730  DenseVector<Number> old_rhs(rhs);
2731  C.vector_mult_transpose(rhs, old_rhs);
2732 
2733  for (unsigned int i=0,
2734  n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2735  i != n_elem_dofs; i++)
2736  {
2737  const dof_id_type dof_id = elem_dofs[i];
2738 
2739  if (auto pos = _dof_constraints.find(dof_id);
2740  pos != _dof_constraints.end())
2741  {
2742  // This will put a nonsymmetric entry in the constraint
2743  // row to ensure that the linear system produces the
2744  // correct value for the constrained DOF.
2745  const DofConstraintRow & constraint_row = pos->second;
2746 
2747  const DofConstraintValueMap::const_iterator valpos =
2748  _primal_constraint_values.find(dof_id);
2749 
2750  Number & rhs_val = rhs(i);
2751  rhs_val = (valpos == _primal_constraint_values.end()) ?
2752  0 : -valpos->second;
2753  for (const auto & [constraining_dof, coef] : constraint_row)
2754  rhs_val -= coef * solution_local(constraining_dof);
2755  rhs_val += solution_local(dof_id);
2756  }
2757  }
2758 }
2759 
2760 
2763  std::vector<dof_id_type> & elem_dofs,
2764  NumericVector<Number> & solution_local) const
2765 {
2766  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2767 
2768  libmesh_assert (solution_local.type() == SERIAL ||
2769  solution_local.type() == GHOSTED);
2770 
2771  // check for easy return
2772  if (this->_dof_constraints.empty())
2773  return;
2774 
2775  // The constrained RHS is built up as C^T F
2777 
2778  this->build_constraint_matrix (C, elem_dofs);
2779 
2780  LOG_SCOPE("cnstrn_elem_residual()", "DofMap");
2781 
2782  // It is possible that the matrix is not constrained at all.
2783  if (C.n() != elem_dofs.size())
2784  return;
2785 
2786  // Compute the matrix-vector product C^T F
2787  DenseVector<Number> old_rhs(rhs);
2788  C.vector_mult_transpose(rhs, old_rhs);
2789 
2790  for (unsigned int i=0,
2791  n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2792  i != n_elem_dofs; i++)
2793  {
2794  const dof_id_type dof_id = elem_dofs[i];
2795 
2796  if (auto pos = _dof_constraints.find(dof_id);
2797  pos != _dof_constraints.end())
2798  {
2799  // This will put a nonsymmetric entry in the constraint
2800  // row to ensure that the linear system produces the
2801  // correct value for the constrained DOF.
2802  const DofConstraintRow & constraint_row = pos->second;
2803 
2804  Number & rhs_val = rhs(i);
2805  rhs_val = 0;
2806  for (const auto & [constraining_dof, coef] : constraint_row)
2807  rhs_val -= coef * solution_local(constraining_dof);
2808  rhs_val += solution_local(dof_id);
2809  }
2810  }
2811 }
2812 
2813 
2815  DenseVector<Number> & rhs,
2816  std::vector<dof_id_type> & elem_dofs,
2817  bool asymmetric_constraint_rows,
2818  int qoi_index) const
2819 {
2820  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
2821  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
2822  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
2823 
2824  // check for easy return
2825  if (this->_dof_constraints.empty())
2826  return;
2827 
2828  // The constrained matrix is built up as C^T K C.
2829  // The constrained RHS is built up as C^T (F - K H)
2832 
2833  this->build_constraint_matrix_and_vector (C, H, elem_dofs, qoi_index);
2834 
2835  LOG_SCOPE("hetero_cnstrn_elem_vec()", "DofMap");
2836 
2837  // It is possible that the matrix is not constrained at all.
2838  if ((C.m() == matrix.m()) &&
2839  (C.n() == elem_dofs.size())) // It the matrix is constrained
2840  {
2841  // We may have rhs values to use later
2842  const DofConstraintValueMap * rhs_values = nullptr;
2843  if (qoi_index < 0)
2844  rhs_values = &_primal_constraint_values;
2845  else if (auto it = _adjoint_constraint_values.find(qoi_index);
2846  it != _adjoint_constraint_values.end())
2847  rhs_values = &it->second;
2848 
2849  // Compute matrix/vector product K H
2851  matrix.vector_mult(KH, H);
2852 
2853  // Compute the matrix-vector product C^T (F - KH)
2854  DenseVector<Number> F_minus_KH(rhs);
2855  F_minus_KH -= KH;
2856  C.vector_mult_transpose(rhs, F_minus_KH);
2857 
2858  for (unsigned int i=0,
2859  n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
2860  i != n_elem_dofs; i++)
2861  {
2862  const dof_id_type dof_id = elem_dofs[i];
2863 
2864  if (this->is_constrained_dof(dof_id))
2865  {
2866  // This will put a nonsymmetric entry in the constraint
2867  // row to ensure that the linear system produces the
2868  // correct value for the constrained DOF.
2869  if (asymmetric_constraint_rows && rhs_values)
2870  {
2871  const DofConstraintValueMap::const_iterator valpos =
2872  rhs_values->find(dof_id);
2873 
2874  rhs(i) = (valpos == rhs_values->end()) ?
2875  0 : valpos->second;
2876  }
2877  else
2878  rhs(i) = 0.;
2879  }
2880  }
2881 
2882  } // end if is constrained...
2883 }
2884 
2885 
2886 
2887 
2889  std::vector<dof_id_type> & row_dofs,
2890  std::vector<dof_id_type> & col_dofs,
2891  bool asymmetric_constraint_rows) const
2892 {
2893  libmesh_assert_equal_to (row_dofs.size(), matrix.m());
2894  libmesh_assert_equal_to (col_dofs.size(), matrix.n());
2895 
2896  // check for easy return
2897  if (this->_dof_constraints.empty())
2898  return;
2899 
2900  // The constrained matrix is built up as R^T K C.
2903 
2904  // Safeguard against the user passing us the same
2905  // object for row_dofs and col_dofs. If that is done
2906  // the calls to build_matrix would fail
2907  std::vector<dof_id_type> orig_row_dofs(row_dofs);
2908  std::vector<dof_id_type> orig_col_dofs(col_dofs);
2909 
2910  this->build_constraint_matrix (R, orig_row_dofs);
2911  this->build_constraint_matrix (C, orig_col_dofs);
2912 
2913  LOG_SCOPE("constrain_elem_matrix()", "DofMap");
2914 
2915  row_dofs = orig_row_dofs;
2916  col_dofs = orig_col_dofs;
2917 
2918  bool constraint_found = false;
2919 
2920  // K_constrained = R^T K C
2921 
2922  if ((R.m() == matrix.m()) &&
2923  (R.n() == row_dofs.size()))
2924  {
2925  matrix.left_multiply_transpose (R);
2926  constraint_found = true;
2927  }
2928 
2929  if ((C.m() == matrix.n()) &&
2930  (C.n() == col_dofs.size()))
2931  {
2932  matrix.right_multiply (C);
2933  constraint_found = true;
2934  }
2935 
2936  // It is possible that the matrix is not constrained at all.
2937  if (constraint_found)
2938  {
2939  libmesh_assert_equal_to (matrix.m(), row_dofs.size());
2940  libmesh_assert_equal_to (matrix.n(), col_dofs.size());
2941 
2942 
2943  for (unsigned int i=0,
2944  n_row_dofs = cast_int<unsigned int>(row_dofs.size());
2945  i != n_row_dofs; i++)
2946  if (this->is_constrained_dof(row_dofs[i]))
2947  {
2948  for (auto j : make_range(matrix.n()))
2949  {
2950  if (row_dofs[i] != col_dofs[j])
2951  matrix(i,j) = 0.;
2952  else // If the DOF is constrained
2953  matrix(i,j) = 1.;
2954  }
2955 
2956  if (asymmetric_constraint_rows)
2957  {
2958  DofConstraints::const_iterator
2959  pos = _dof_constraints.find(row_dofs[i]);
2960 
2961  libmesh_assert (pos != _dof_constraints.end());
2962 
2963  const DofConstraintRow & constraint_row = pos->second;
2964 
2965  libmesh_assert (!constraint_row.empty());
2966 
2967  for (const auto & item : constraint_row)
2968  for (unsigned int j=0,
2969  n_col_dofs = cast_int<unsigned int>(col_dofs.size());
2970  j != n_col_dofs; j++)
2971  if (col_dofs[j] == item.first)
2972  matrix(i,j) = -item.second;
2973  }
2974  }
2975  } // end if is constrained...
2976 }
2977 
2978 
2979 
2981  std::vector<dof_id_type> & row_dofs,
2982  bool) const
2983 {
2984  libmesh_assert_equal_to (rhs.size(), row_dofs.size());
2985 
2986  // check for easy return
2987  if (this->_dof_constraints.empty())
2988  return;
2989 
2990  // The constrained RHS is built up as R^T F.
2992 
2993  this->build_constraint_matrix (R, row_dofs);
2994 
2995  LOG_SCOPE("constrain_elem_vector()", "DofMap");
2996 
2997  // It is possible that the vector is not constrained at all.
2998  if ((R.m() == rhs.size()) &&
2999  (R.n() == row_dofs.size())) // if the RHS is constrained
3000  {
3001  // Compute the matrix-vector product
3002  DenseVector<Number> old_rhs(rhs);
3003  R.vector_mult_transpose(rhs, old_rhs);
3004 
3005  libmesh_assert_equal_to (row_dofs.size(), rhs.size());
3006 
3007  for (unsigned int i=0,
3008  n_row_dofs = cast_int<unsigned int>(row_dofs.size());
3009  i != n_row_dofs; i++)
3010  if (this->is_constrained_dof(row_dofs[i]))
3011  {
3012  // If the DOF is constrained
3013  libmesh_assert (_dof_constraints.find(row_dofs[i]) != _dof_constraints.end());
3014 
3015  rhs(i) = 0;
3016  }
3017  } // end if the RHS is constrained.
3018 }
3019 
3020 
3021 
3023  DenseVector<Number> & w,
3024  std::vector<dof_id_type> & row_dofs,
3025  bool) const
3026 {
3027  libmesh_assert_equal_to (v.size(), row_dofs.size());
3028  libmesh_assert_equal_to (w.size(), row_dofs.size());
3029 
3030  // check for easy return
3031  if (this->_dof_constraints.empty())
3032  return;
3033 
3034  // The constrained RHS is built up as R^T F.
3036 
3037  this->build_constraint_matrix (R, row_dofs);
3038 
3039  LOG_SCOPE("cnstrn_elem_dyad_mat()", "DofMap");
3040 
3041  // It is possible that the vector is not constrained at all.
3042  if ((R.m() == v.size()) &&
3043  (R.n() == row_dofs.size())) // if the RHS is constrained
3044  {
3045  // Compute the matrix-vector products
3046  DenseVector<Number> old_v(v);
3047  DenseVector<Number> old_w(w);
3048 
3049  // compute matrix/vector product
3050  R.vector_mult_transpose(v, old_v);
3051  R.vector_mult_transpose(w, old_w);
3052 
3053  libmesh_assert_equal_to (row_dofs.size(), v.size());
3054  libmesh_assert_equal_to (row_dofs.size(), w.size());
3055 
3056  /* Constrain only v, not w. */
3057 
3058  for (unsigned int i=0,
3059  n_row_dofs = cast_int<unsigned int>(row_dofs.size());
3060  i != n_row_dofs; i++)
3061  if (this->is_constrained_dof(row_dofs[i]))
3062  {
3063  // If the DOF is constrained
3064  libmesh_assert (_dof_constraints.find(row_dofs[i]) != _dof_constraints.end());
3065 
3066  v(i) = 0;
3067  }
3068  } // end if the RHS is constrained.
3069 }
3070 
3071 
3072 
3073 void DofMap::constrain_nothing (std::vector<dof_id_type> & dofs) const
3074 {
3075  // check for easy return
3076  if (this->_dof_constraints.empty())
3077  return;
3078 
3079  // All the work is done by \p build_constraint_matrix. We just need
3080  // a dummy matrix.
3082  this->build_constraint_matrix (R, dofs);
3083 }
3084 
3085 
3086 
3087 void DofMap::enforce_constraints_exactly (const System & system,
3089  bool homogeneous) const
3090 {
3091  parallel_object_only();
3092 
3093  if (!this->n_constrained_dofs())
3094  return;
3095 
3096  LOG_SCOPE("enforce_constraints_exactly()","DofMap");
3097 
3098  if (!v)
3099  v = system.solution.get();
3100 
3101  if (!v->closed())
3102  v->close();
3103 
3104  NumericVector<Number> * v_local = nullptr; // will be initialized below
3105  NumericVector<Number> * v_global = nullptr; // will be initialized below
3106  std::unique_ptr<NumericVector<Number>> v_built;
3107  if (v->type() == SERIAL)
3108  {
3109  v_built = NumericVector<Number>::build(this->comm());
3110  v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
3111  v_built->close();
3112 
3113  for (dof_id_type i=v_built->first_local_index();
3114  i<v_built->last_local_index(); i++)
3115  v_built->set(i, (*v)(i));
3116  v_built->close();
3117  v_global = v_built.get();
3118 
3119  v_local = v;
3120  libmesh_assert (v_local->closed());
3121  }
3122  else if (v->type() == PARALLEL)
3123  {
3124  v_built = NumericVector<Number>::build(this->comm());
3125  v_built->init (v->size(), v->local_size(),
3126  this->get_send_list(), true,
3127  GHOSTED);
3128  v->localize(*v_built, this->get_send_list());
3129  v_built->close();
3130  v_local = v_built.get();
3131 
3132  v_global = v;
3133  }
3134  else if (v->type() == GHOSTED)
3135  {
3136  v_local = v;
3137  v_global = v;
3138  }
3139  else // unknown v->type()
3140  libmesh_error_msg("ERROR: Unsupported NumericVector type == " << Utility::enum_to_string(v->type()));
3141 
3142  // We should never hit these asserts because we should error-out in
3143  // else clause above. Just to be sure we don't try to use v_local
3144  // and v_global uninitialized...
3145  libmesh_assert(v_local);
3146  libmesh_assert(v_global);
3147  libmesh_assert_equal_to (this, &(system.get_dof_map()));
3148 
3149  for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
3150  {
3151  if (!this->local_index(constrained_dof))
3152  continue;
3153 
3154  Number exact_value = 0;
3155  if (!homogeneous)
3156  {
3157  if (auto rhsit = _primal_constraint_values.find(constrained_dof);
3158  rhsit != _primal_constraint_values.end())
3159  exact_value = rhsit->second;
3160  }
3161  for (const auto & [dof, val] : constraint_row)
3162  exact_value += val * (*v_local)(dof);
3163 
3164  v_global->set(constrained_dof, exact_value);
3165  }
3166 
3167  // If the old vector was serial, we probably need to send our values
3168  // to other processors
3169  if (v->type() == SERIAL)
3170  {
3171 #ifndef NDEBUG
3172  v_global->close();
3173 #endif
3174  v_global->localize (*v);
3175  }
3176  v->close();
3177 }
3178 
3180  NumericVector<Number> * rhs,
3181  NumericVector<Number> const * solution,
3182  bool homogeneous) const
3183 {
3184  parallel_object_only();
3185 
3186  if (!this->n_constrained_dofs())
3187  return;
3188 
3189  if (!rhs)
3190  rhs = system.rhs;
3191  if (!solution)
3192  solution = system.solution.get();
3193 
3194  NumericVector<Number> const * solution_local = nullptr; // will be initialized below
3195  std::unique_ptr<NumericVector<Number>> solution_built;
3196  if (solution->type() == SERIAL || solution->type() == GHOSTED)
3197  solution_local = solution;
3198  else if (solution->type() == PARALLEL)
3199  {
3200  solution_built = NumericVector<Number>::build(this->comm());
3201  solution_built->init (solution->size(), solution->local_size(),
3202  this->get_send_list(), true, GHOSTED);
3203  solution->localize(*solution_built, this->get_send_list());
3204  solution_built->close();
3205  solution_local = solution_built.get();
3206  }
3207  else // unknown solution->type()
3208  libmesh_error_msg("ERROR: Unsupported NumericVector type == " << Utility::enum_to_string(solution->type()));
3209 
3210  // We should never hit these asserts because we should error-out in
3211  // else clause above. Just to be sure we don't try to use solution_local
3212  libmesh_assert(solution_local);
3213  libmesh_assert_equal_to (this, &(system.get_dof_map()));
3214 
3215  for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
3216  {
3217  if (!this->local_index(constrained_dof))
3218  continue;
3219 
3220  Number exact_value = 0;
3221  for (const auto & [dof, val] : constraint_row)
3222  exact_value -= val * (*solution_local)(dof);
3223  exact_value += (*solution_local)(constrained_dof);
3224  if (!homogeneous)
3225  {
3226  if (auto rhsit = _primal_constraint_values.find(constrained_dof);
3227  rhsit != _primal_constraint_values.end())
3228  exact_value += rhsit->second;
3229  }
3230 
3231  rhs->set(constrained_dof, exact_value);
3232  }
3233 }
3234 
3236  SparseMatrix<Number> * jac) const
3237 {
3238  parallel_object_only();
3239 
3240  if (!this->n_constrained_dofs())
3241  return;
3242 
3243  if (!jac)
3244  jac = system.matrix;
3245 
3246  libmesh_assert_equal_to (this, &(system.get_dof_map()));
3247 
3248  for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
3249  {
3250  if (!this->local_index(constrained_dof))
3251  continue;
3252 
3253  for (const auto & j : constraint_row)
3254  jac->set(constrained_dof, j.first, -j.second);
3255  jac->set(constrained_dof, constrained_dof, 1);
3256  }
3257 }
3258 
3259 
3261  unsigned int q) const
3262 {
3263  parallel_object_only();
3264 
3265  if (!this->n_constrained_dofs())
3266  return;
3267 
3268  LOG_SCOPE("enforce_adjoint_constraints_exactly()", "DofMap");
3269 
3270  NumericVector<Number> * v_local = nullptr; // will be initialized below
3271  NumericVector<Number> * v_global = nullptr; // will be initialized below
3272  std::unique_ptr<NumericVector<Number>> v_built;
3273  if (v.type() == SERIAL)
3274  {
3275  v_built = NumericVector<Number>::build(this->comm());
3276  v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
3277  v_built->close();
3278 
3279  for (dof_id_type i=v_built->first_local_index();
3280  i<v_built->last_local_index(); i++)
3281  v_built->set(i, v(i));
3282  v_built->close();
3283  v_global = v_built.get();
3284 
3285  v_local = &v;
3286  libmesh_assert (v_local->closed());
3287  }
3288  else if (v.type() == PARALLEL)
3289  {
3290  v_built = NumericVector<Number>::build(this->comm());
3291  v_built->init (v.size(), v.local_size(),
3292  this->get_send_list(), true, GHOSTED);
3293  v.localize(*v_built, this->get_send_list());
3294  v_built->close();
3295  v_local = v_built.get();
3296 
3297  v_global = &v;
3298  }
3299  else if (v.type() == GHOSTED)
3300  {
3301  v_local = &v;
3302  v_global = &v;
3303  }
3304  else // unknown v.type()
3305  libmesh_error_msg("ERROR: Unknown v.type() == " << v.type());
3306 
3307  // We should never hit these asserts because we should error-out in
3308  // else clause above. Just to be sure we don't try to use v_local
3309  // and v_global uninitialized...
3310  libmesh_assert(v_local);
3311  libmesh_assert(v_global);
3312 
3313  // Do we have any non_homogeneous constraints?
3314  const AdjointDofConstraintValues::const_iterator
3315  adjoint_constraint_map_it = _adjoint_constraint_values.find(q);
3316  const DofConstraintValueMap * constraint_map =
3317  (adjoint_constraint_map_it == _adjoint_constraint_values.end()) ?
3318  nullptr : &adjoint_constraint_map_it->second;
3319 
3320  for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
3321  {
3322  if (!this->local_index(constrained_dof))
3323  continue;
3324 
3325  Number exact_value = 0;
3326  if (constraint_map)
3327  {
3328  if (const auto adjoint_constraint_it = constraint_map->find(constrained_dof);
3329  adjoint_constraint_it != constraint_map->end())
3330  exact_value = adjoint_constraint_it->second;
3331  }
3332 
3333  for (const auto & j : constraint_row)
3334  exact_value += j.second * (*v_local)(j.first);
3335 
3336  v_global->set(constrained_dof, exact_value);
3337  }
3338 
3339  // If the old vector was serial, we probably need to send our values
3340  // to other processors
3341  if (v.type() == SERIAL)
3342  {
3343 #ifndef NDEBUG
3344  v_global->close();
3345 #endif
3346  v_global->localize (v);
3347  }
3348  v.close();
3349 }
3350 
3351 
3352 
3353 std::pair<Real, Real>
3355  NumericVector<Number> * v) const
3356 {
3357  if (!v)
3358  v = system.solution.get();
3359  NumericVector<Number> & vec = *v;
3360 
3361  // We'll assume the vector is closed
3362  libmesh_assert (vec.closed());
3363 
3364  Real max_absolute_error = 0., max_relative_error = 0.;
3365 
3366  const MeshBase & mesh = system.get_mesh();
3367 
3368  libmesh_assert_equal_to (this, &(system.get_dof_map()));
3369 
3370  // indices on each element
3371  std::vector<dof_id_type> local_dof_indices;
3372 
3373  for (const auto & elem : mesh.active_local_element_ptr_range())
3374  {
3375  this->dof_indices(elem, local_dof_indices);
3376  std::vector<dof_id_type> raw_dof_indices = local_dof_indices;
3377 
3378  // Constraint matrix for each element
3380 
3381  this->build_constraint_matrix (C, local_dof_indices);
3382 
3383  // Continue if the element is unconstrained
3384  if (!C.m())
3385  continue;
3386 
3387  libmesh_assert_equal_to (C.m(), raw_dof_indices.size());
3388  libmesh_assert_equal_to (C.n(), local_dof_indices.size());
3389 
3390  for (auto i : make_range(C.m()))
3391  {
3392  // Recalculate any constrained dof owned by this processor
3393  dof_id_type global_dof = raw_dof_indices[i];
3394  if (this->is_constrained_dof(global_dof) &&
3395  global_dof >= vec.first_local_index() &&
3396  global_dof < vec.last_local_index())
3397  {
3398 #ifndef NDEBUG
3399  DofConstraints::const_iterator
3400  pos = _dof_constraints.find(global_dof);
3401 
3402  libmesh_assert (pos != _dof_constraints.end());
3403 #endif
3404 
3405  Number exact_value = 0;
3406  DofConstraintValueMap::const_iterator rhsit =
3407  _primal_constraint_values.find(global_dof);
3408  if (rhsit != _primal_constraint_values.end())
3409  exact_value = rhsit->second;
3410 
3411  for (auto j : make_range(C.n()))
3412  {
3413  if (local_dof_indices[j] != global_dof)
3414  exact_value += C(i,j) *
3415  vec(local_dof_indices[j]);
3416  }
3417 
3418  max_absolute_error = std::max(max_absolute_error,
3419  std::abs(vec(global_dof) - exact_value));
3420  max_relative_error = std::max(max_relative_error,
3421  std::abs(vec(global_dof) - exact_value)
3422  / std::abs(exact_value));
3423  }
3424  }
3425  }
3426 
3427  return std::pair<Real, Real>(max_absolute_error, max_relative_error);
3428 }
3429 
3430 
3431 
3433  std::vector<dof_id_type> & elem_dofs,
3434  const bool called_recursively) const
3435 {
3436  LOG_SCOPE_IF("build_constraint_matrix()", "DofMap", !called_recursively);
3437 
3438  // Create a set containing the DOFs we already depend on
3439  typedef std::set<dof_id_type> RCSet;
3440  RCSet dof_set;
3441 
3442  bool we_have_constraints = false;
3443 
3444  // Next insert any other dofs the current dofs might be constrained
3445  // in terms of. Note that in this case we may not be done: Those
3446  // may in turn depend on others. So, we need to repeat this process
3447  // in that case until the system depends only on unconstrained
3448  // degrees of freedom.
3449  for (const auto & dof : elem_dofs)
3450  if (this->is_constrained_dof(dof))
3451  {
3452  we_have_constraints = true;
3453 
3454  // If the DOF is constrained
3455  DofConstraints::const_iterator
3456  pos = _dof_constraints.find(dof);
3457 
3458  libmesh_assert (pos != _dof_constraints.end());
3459 
3460  const DofConstraintRow & constraint_row = pos->second;
3461 
3462  // Constraint rows in p refinement may be empty
3463  //libmesh_assert (!constraint_row.empty());
3464 
3465  for (const auto & item : constraint_row)
3466  dof_set.insert (item.first);
3467  }
3468 
3469  // May be safe to return at this point
3470  // (but remember to stop the perflog)
3471  if (!we_have_constraints)
3472  return;
3473 
3474  for (const auto & dof : elem_dofs)
3475  dof_set.erase (dof);
3476 
3477  // If we added any DOFS then we need to do this recursively.
3478  // It is possible that we just added a DOF that is also
3479  // constrained!
3480  //
3481  // Also, we need to handle the special case of an element having DOFs
3482  // constrained in terms of other, local DOFs
3483  if (!dof_set.empty() || // case 1: constrained in terms of other DOFs
3484  !called_recursively) // case 2: constrained in terms of our own DOFs
3485  {
3486  const unsigned int old_size =
3487  cast_int<unsigned int>(elem_dofs.size());
3488 
3489  // Add new dependency dofs to the end of the current dof set
3490  elem_dofs.insert(elem_dofs.end(),
3491  dof_set.begin(), dof_set.end());
3492 
3493  // Now we can build the constraint matrix.
3494  // Note that resize also zeros for a DenseMatrix<Number>.
3495  C.resize (old_size,
3496  cast_int<unsigned int>(elem_dofs.size()));
3497 
3498  // Create the C constraint matrix.
3499  for (unsigned int i=0; i != old_size; i++)
3500  if (this->is_constrained_dof(elem_dofs[i]))
3501  {
3502  // If the DOF is constrained
3503  DofConstraints::const_iterator
3504  pos = _dof_constraints.find(elem_dofs[i]);
3505 
3506  libmesh_assert (pos != _dof_constraints.end());
3507 
3508  const DofConstraintRow & constraint_row = pos->second;
3509 
3510  // p refinement creates empty constraint rows
3511  // libmesh_assert (!constraint_row.empty());
3512 
3513  for (const auto & item : constraint_row)
3514  for (unsigned int j=0,
3515  n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
3516  j != n_elem_dofs; j++)
3517  if (elem_dofs[j] == item.first)
3518  C(i,j) = item.second;
3519  }
3520  else
3521  {
3522  C(i,i) = 1.;
3523  }
3524 
3525  // May need to do this recursively. It is possible
3526  // that we just replaced a constrained DOF with another
3527  // constrained DOF.
3528  DenseMatrix<Number> Cnew;
3529 
3530  this->build_constraint_matrix (Cnew, elem_dofs, true);
3531 
3532  if ((C.n() == Cnew.m()) &&
3533  (Cnew.n() == elem_dofs.size())) // If the constraint matrix
3534  C.right_multiply(Cnew); // is constrained...
3535 
3536  libmesh_assert_equal_to (C.n(), elem_dofs.size());
3537  }
3538 }
3539 
3540 
3541 
3543  DenseVector<Number> & H,
3544  std::vector<dof_id_type> & elem_dofs,
3545  int qoi_index,
3546  const bool called_recursively) const
3547 {
3548  LOG_SCOPE_IF("build_constraint_matrix_and_vector()", "DofMap", !called_recursively);
3549 
3550  // Create a set containing the DOFs we already depend on
3551  typedef std::set<dof_id_type> RCSet;
3552  RCSet dof_set;
3553 
3554  bool we_have_constraints = false;
3555 
3556  // Next insert any other dofs the current dofs might be constrained
3557  // in terms of. Note that in this case we may not be done: Those
3558  // may in turn depend on others. So, we need to repeat this process
3559  // in that case until the system depends only on unconstrained
3560  // degrees of freedom.
3561  for (const auto & dof : elem_dofs)
3562  if (this->is_constrained_dof(dof))
3563  {
3564  we_have_constraints = true;
3565 
3566  // If the DOF is constrained
3567  DofConstraints::const_iterator
3568  pos = _dof_constraints.find(dof);
3569 
3570  libmesh_assert (pos != _dof_constraints.end());
3571 
3572  const DofConstraintRow & constraint_row = pos->second;
3573 
3574  // Constraint rows in p refinement may be empty
3575  //libmesh_assert (!constraint_row.empty());
3576 
3577  for (const auto & item : constraint_row)
3578  dof_set.insert (item.first);
3579  }
3580 
3581  // May be safe to return at this point
3582  // (but remember to stop the perflog)
3583  if (!we_have_constraints)
3584  return;
3585 
3586  for (const auto & dof : elem_dofs)
3587  dof_set.erase (dof);
3588 
3589  // If we added any DOFS then we need to do this recursively.
3590  // It is possible that we just added a DOF that is also
3591  // constrained!
3592  //
3593  // Also, we need to handle the special case of an element having DOFs
3594  // constrained in terms of other, local DOFs
3595  if (!dof_set.empty() || // case 1: constrained in terms of other DOFs
3596  !called_recursively) // case 2: constrained in terms of our own DOFs
3597  {
3598  const DofConstraintValueMap * rhs_values = nullptr;
3599  if (qoi_index < 0)
3600  rhs_values = &_primal_constraint_values;
3601  else if (auto it = _adjoint_constraint_values.find(qoi_index);
3602  it != _adjoint_constraint_values.end())
3603  rhs_values = &it->second;
3604 
3605  const unsigned int old_size =
3606  cast_int<unsigned int>(elem_dofs.size());
3607 
3608  // Add new dependency dofs to the end of the current dof set
3609  elem_dofs.insert(elem_dofs.end(),
3610  dof_set.begin(), dof_set.end());
3611 
3612  // Now we can build the constraint matrix and vector.
3613  // Note that resize also zeros for a DenseMatrix and DenseVector
3614  C.resize (old_size,
3615  cast_int<unsigned int>(elem_dofs.size()));
3616  H.resize (old_size);
3617 
3618  // Create the C constraint matrix.
3619  for (unsigned int i=0; i != old_size; i++)
3620  if (this->is_constrained_dof(elem_dofs[i]))
3621  {
3622  // If the DOF is constrained
3623  DofConstraints::const_iterator
3624  pos = _dof_constraints.find(elem_dofs[i]);
3625 
3626  libmesh_assert (pos != _dof_constraints.end());
3627 
3628  const DofConstraintRow & constraint_row = pos->second;
3629 
3630  // p refinement creates empty constraint rows
3631  // libmesh_assert (!constraint_row.empty());
3632 
3633  for (const auto & item : constraint_row)
3634  for (unsigned int j=0,
3635  n_elem_dofs = cast_int<unsigned int>(elem_dofs.size());
3636  j != n_elem_dofs; j++)
3637  if (elem_dofs[j] == item.first)
3638  C(i,j) = item.second;
3639 
3640  if (rhs_values)
3641  {
3642  if (const auto rhsit = rhs_values->find(elem_dofs[i]);
3643  rhsit != rhs_values->end())
3644  H(i) = rhsit->second;
3645  }
3646  }
3647  else
3648  {
3649  C(i,i) = 1.;
3650  }
3651 
3652  // May need to do this recursively. It is possible
3653  // that we just replaced a constrained DOF with another
3654  // constrained DOF.
3655  DenseMatrix<Number> Cnew;
3656  DenseVector<Number> Hnew;
3657 
3658  this->build_constraint_matrix_and_vector (Cnew, Hnew, elem_dofs,
3659  qoi_index, true);
3660 
3661  if ((C.n() == Cnew.m()) && // If the constraint matrix
3662  (Cnew.n() == elem_dofs.size())) // is constrained...
3663  {
3664  // If x = Cy + h and y = Dz + g
3665  // Then x = (CD)z + (Cg + h)
3666  C.vector_mult_add(H, 1, Hnew);
3667 
3668  C.right_multiply(Cnew);
3669  }
3670 
3671  libmesh_assert_equal_to (C.n(), elem_dofs.size());
3672  }
3673 }
3674 
3675 
3677 {
3678  // This function must be run on all processors at once
3679  parallel_object_only();
3680 
3681  // Return immediately if there's nothing to gather
3682  if (this->n_processors() == 1)
3683  return;
3684 
3685  // We might get to return immediately if none of the processors
3686  // found any constraints
3687  unsigned int has_constraints = !_dof_constraints.empty()
3688 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3689  || !_node_constraints.empty()
3690 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3691  ;
3692  this->comm().max(has_constraints);
3693  if (!has_constraints)
3694  return;
3695 
3696  // If we have heterogeneous adjoint constraints we need to
3697  // communicate those too.
3698  const unsigned int max_qoi_num =
3699  _adjoint_constraint_values.empty() ?
3700  0 : _adjoint_constraint_values.rbegin()->first+1;
3701 
3702 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3703  // We may need to send nodes ahead of data about them
3704  std::vector<Parallel::Request> packed_range_sends;
3705 
3706  // We may be receiving packed_range sends out of order with
3707  // parallel_sync tags, so make sure they're received correctly.
3708  Parallel::MessageTag range_tag = this->comm().get_unique_tag();
3709 
3710  // We only need to do these sends on a distributed mesh
3711  const bool dist_mesh = !mesh.is_serial();
3712 #endif
3713 
3714  // We might have calculated constraints for constrained dofs
3715  // which have support on other processors.
3716  // Push these out first.
3717  {
3718  std::map<processor_id_type, std::set<dof_id_type>> pushed_ids;
3719 
3720 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3721  std::map<processor_id_type, std::set<dof_id_type>> pushed_node_ids;
3722 #endif
3723 
3724  const unsigned int sys_num = this->sys_number();
3725 
3726  // Collect the constraints to push to each processor
3727  for (auto & elem : as_range(mesh.active_not_local_elements_begin(),
3728  mesh.active_not_local_elements_end()))
3729  {
3730  const unsigned short n_nodes = elem->n_nodes();
3731 
3732  // Just checking dof_indices on the foreign element isn't
3733  // enough. Consider a central hanging node between a coarse
3734  // Q2/Q1 element and its finer neighbors on a higher-ranked
3735  // processor. The coarse element's processor will own the node,
3736  // and will thereby own the pressure dof on that node, despite
3737  // the fact that that pressure dof doesn't directly exist on the
3738  // coarse element!
3739  //
3740  // So, we loop through dofs manually.
3741 
3742  {
3743  const unsigned int n_vars = elem->n_vars(sys_num);
3744  for (unsigned int v=0; v != n_vars; ++v)
3745  {
3746  const unsigned int n_comp = elem->n_comp(sys_num,v);
3747  for (unsigned int c=0; c != n_comp; ++c)
3748  {
3749  const unsigned int id =
3750  elem->dof_number(sys_num,v,c);
3751  if (this->is_constrained_dof(id))
3752  pushed_ids[elem->processor_id()].insert(id);
3753  }
3754  }
3755  }
3756 
3757  for (unsigned short n = 0; n != n_nodes; ++n)
3758  {
3759  const Node & node = elem->node_ref(n);
3760  const unsigned int n_vars = node.n_vars(sys_num);
3761  for (unsigned int v=0; v != n_vars; ++v)
3762  {
3763  const unsigned int n_comp = node.n_comp(sys_num,v);
3764  for (unsigned int c=0; c != n_comp; ++c)
3765  {
3766  const unsigned int id =
3767  node.dof_number(sys_num,v,c);
3768  if (this->is_constrained_dof(id))
3769  pushed_ids[elem->processor_id()].insert(id);
3770  }
3771  }
3772  }
3773 
3774 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3775  for (unsigned short n = 0; n != n_nodes; ++n)
3776  if (this->is_constrained_node(elem->node_ptr(n)))
3777  pushed_node_ids[elem->processor_id()].insert(elem->node_id(n));
3778 #endif
3779  }
3780 
3781  // Rewrite those id sets as vectors for sending and receiving,
3782  // then find the corresponding data for each id, then push it all.
3783  std::map<processor_id_type, std::vector<dof_id_type>>
3784  pushed_id_vecs, received_id_vecs;
3785  for (auto & p : pushed_ids)
3786  pushed_id_vecs[p.first].assign(p.second.begin(), p.second.end());
3787 
3788  std::map<processor_id_type, std::vector<std::vector<std::pair<dof_id_type,Real>>>>
3789  pushed_keys_vals, received_keys_vals;
3790  std::map<processor_id_type, std::vector<std::vector<Number>>> pushed_rhss, received_rhss;
3791  for (auto & p : pushed_id_vecs)
3792  {
3793  auto & keys_vals = pushed_keys_vals[p.first];
3794  keys_vals.reserve(p.second.size());
3795 
3796  auto & rhss = pushed_rhss[p.first];
3797  rhss.reserve(p.second.size());
3798  for (auto & pushed_id : p.second)
3799  {
3800  const DofConstraintRow & row = _dof_constraints[pushed_id];
3801  keys_vals.emplace_back(row.begin(), row.end());
3802 
3803  rhss.push_back(std::vector<Number>(max_qoi_num+1));
3804  std::vector<Number> & rhs = rhss.back();
3805  DofConstraintValueMap::const_iterator rhsit =
3806  _primal_constraint_values.find(pushed_id);
3807  rhs[max_qoi_num] =
3808  (rhsit == _primal_constraint_values.end()) ?
3809  0 : rhsit->second;
3810  for (unsigned int q = 0; q != max_qoi_num; ++q)
3811  {
3812  AdjointDofConstraintValues::const_iterator adjoint_map_it =
3814 
3815  if (adjoint_map_it == _adjoint_constraint_values.end())
3816  continue;
3817 
3818  const DofConstraintValueMap & constraint_map =
3819  adjoint_map_it->second;
3820 
3821  DofConstraintValueMap::const_iterator adj_rhsit =
3822  constraint_map.find(pushed_id);
3823 
3824  rhs[q] =
3825  (adj_rhsit == constraint_map.end()) ?
3826  0 : adj_rhsit->second;
3827  }
3828  }
3829  }
3830 
3831  auto ids_action_functor =
3832  [& received_id_vecs]
3833  (processor_id_type pid,
3834  const std::vector<dof_id_type> & data)
3835  {
3836  received_id_vecs[pid] = data;
3837  };
3838 
3839  Parallel::push_parallel_vector_data
3840  (this->comm(), pushed_id_vecs, ids_action_functor);
3841 
3842  auto keys_vals_action_functor =
3843  [& received_keys_vals]
3844  (processor_id_type pid,
3845  const std::vector<std::vector<std::pair<dof_id_type,Real>>> & data)
3846  {
3847  received_keys_vals[pid] = data;
3848  };
3849 
3850  Parallel::push_parallel_vector_data
3851  (this->comm(), pushed_keys_vals, keys_vals_action_functor);
3852 
3853  auto rhss_action_functor =
3854  [& received_rhss]
3855  (processor_id_type pid,
3856  const std::vector<std::vector<Number>> & data)
3857  {
3858  received_rhss[pid] = data;
3859  };
3860 
3861  Parallel::push_parallel_vector_data
3862  (this->comm(), pushed_rhss, rhss_action_functor);
3863 
3864  // Now we have all the DofConstraint rows and rhs values received
3865  // from others, so add the DoF constraints that we've been sent
3866 
3867 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3868  std::map<processor_id_type, std::vector<dof_id_type>>
3869  pushed_node_id_vecs, received_node_id_vecs;
3870  for (auto & p : pushed_node_ids)
3871  pushed_node_id_vecs[p.first].assign(p.second.begin(), p.second.end());
3872 
3873  std::map<processor_id_type, std::vector<std::vector<std::pair<dof_id_type,Real>>>>
3874  pushed_node_keys_vals, received_node_keys_vals;
3875  std::map<processor_id_type, std::vector<Point>> pushed_offsets, received_offsets;
3876 
3877  for (auto & p : pushed_node_id_vecs)
3878  {
3879  const processor_id_type pid = p.first;
3880 
3881  // FIXME - this could be an unordered set, given a
3882  // hash<pointers> specialization
3883  std::set<const Node *> nodes_requested;
3884 
3885  auto & node_keys_vals = pushed_node_keys_vals[pid];
3886  node_keys_vals.reserve(p.second.size());
3887 
3888  auto & offsets = pushed_offsets[pid];
3889  offsets.reserve(p.second.size());
3890 
3891  for (auto & pushed_node_id : p.second)
3892  {
3893  const Node * node = mesh.node_ptr(pushed_node_id);
3894  NodeConstraintRow & row = _node_constraints[node].first;
3895  const std::size_t row_size = row.size();
3896  node_keys_vals.push_back
3897  (std::vector<std::pair<dof_id_type,Real>>());
3898  std::vector<std::pair<dof_id_type,Real>> & this_node_kv =
3899  node_keys_vals.back();
3900  this_node_kv.reserve(row_size);
3901  for (const auto & j : row)
3902  {
3903  this_node_kv.emplace_back(j.first->id(), j.second);
3904 
3905  // If we're not sure whether our send
3906  // destination already has this node, let's give
3907  // it a copy.
3908  if (j.first->processor_id() != pid && dist_mesh)
3909  nodes_requested.insert(j.first);
3910  }
3911 
3912  offsets.push_back(_node_constraints[node].second);
3913 
3914  }
3915 
3916  // Constraining nodes might not even exist on our
3917  // correspondant's subset of a distributed mesh, so let's
3918  // make them exist.
3919  if (dist_mesh)
3920  {
3921  packed_range_sends.push_back(Parallel::Request());
3922  this->comm().send_packed_range
3923  (pid, &mesh, nodes_requested.begin(), nodes_requested.end(),
3924  packed_range_sends.back(), range_tag);
3925  }
3926  }
3927 
3928  auto node_ids_action_functor =
3929  [& received_node_id_vecs]
3930  (processor_id_type pid,
3931  const std::vector<dof_id_type> & data)
3932  {
3933  received_node_id_vecs[pid] = data;
3934  };
3935 
3936  Parallel::push_parallel_vector_data
3937  (this->comm(), pushed_node_id_vecs, node_ids_action_functor);
3938 
3939  auto node_keys_vals_action_functor =
3940  [& received_node_keys_vals]
3941  (processor_id_type pid,
3942  const std::vector<std::vector<std::pair<dof_id_type,Real>>> & data)
3943  {
3944  received_node_keys_vals[pid] = data;
3945  };
3946 
3947  Parallel::push_parallel_vector_data
3948  (this->comm(), pushed_node_keys_vals,
3949  node_keys_vals_action_functor);
3950 
3951  auto node_offsets_action_functor =
3952  [& received_offsets]
3953  (processor_id_type pid,
3954  const std::vector<Point> & data)
3955  {
3956  received_offsets[pid] = data;
3957  };
3958 
3959  Parallel::push_parallel_vector_data
3960  (this->comm(), pushed_offsets, node_offsets_action_functor);
3961 
3962 #endif
3963 
3964  // Add all the dof constraints that I've been sent
3965  for (auto & [pid, pushed_ids_to_me] : received_id_vecs)
3966  {
3967  libmesh_assert(received_keys_vals.count(pid));
3968  libmesh_assert(received_rhss.count(pid));
3969  const auto & pushed_keys_vals_to_me = received_keys_vals.at(pid);
3970  const auto & pushed_rhss_to_me = received_rhss.at(pid);
3971 
3972  libmesh_assert_equal_to (pushed_ids_to_me.size(),
3973  pushed_keys_vals_to_me.size());
3974  libmesh_assert_equal_to (pushed_ids_to_me.size(),
3975  pushed_rhss_to_me.size());
3976 
3977  for (auto i : index_range(pushed_ids_to_me))
3978  {
3979  dof_id_type constrained = pushed_ids_to_me[i];
3980 
3981  // If we don't already have a constraint for this dof,
3982  // add the one we were sent
3983  if (!this->is_constrained_dof(constrained))
3984  {
3985  DofConstraintRow & row = _dof_constraints[constrained];
3986  for (auto & kv : pushed_keys_vals_to_me[i])
3987  {
3988  libmesh_assert_less(kv.first, this->n_dofs());
3989  row[kv.first] = kv.second;
3990  }
3991 
3992  const Number primal_rhs = pushed_rhss_to_me[i][max_qoi_num];
3993 
3994  if (libmesh_isnan(primal_rhs))
3995  libmesh_assert(pushed_keys_vals_to_me[i].empty());
3996 
3997  if (primal_rhs != Number(0))
3998  _primal_constraint_values[constrained] = primal_rhs;
3999  else
4000  _primal_constraint_values.erase(constrained);
4001 
4002  for (unsigned int q = 0; q != max_qoi_num; ++q)
4003  {
4004  AdjointDofConstraintValues::iterator adjoint_map_it =
4006 
4007  const Number adj_rhs = pushed_rhss_to_me[i][q];
4008 
4009  if ((adjoint_map_it == _adjoint_constraint_values.end()) &&
4010  adj_rhs == Number(0))
4011  continue;
4012 
4013  if (adjoint_map_it == _adjoint_constraint_values.end())
4014  adjoint_map_it = _adjoint_constraint_values.emplace
4015  (q, DofConstraintValueMap()).first;
4016 
4017  DofConstraintValueMap & constraint_map =
4018  adjoint_map_it->second;
4019 
4020  if (adj_rhs != Number(0))
4021  constraint_map[constrained] = adj_rhs;
4022  else
4023  constraint_map.erase(constrained);
4024  }
4025  }
4026  }
4027  }
4028 
4029 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4030  // Add all the node constraints that I've been sent
4031  for (auto & [pid, pushed_node_ids_to_me] : received_node_id_vecs)
4032  {
4033  // Before we act on any new constraint rows, we may need to
4034  // make sure we have all the nodes involved!
4035  if (dist_mesh)
4036  this->comm().receive_packed_range
4037  (pid, &mesh, null_output_iterator<Node>(),
4038  (Node**)nullptr, range_tag);
4039 
4040  libmesh_assert(received_node_keys_vals.count(pid));
4041  libmesh_assert(received_offsets.count(pid));
4042  const auto & pushed_node_keys_vals_to_me = received_node_keys_vals.at(pid);
4043  const auto & pushed_offsets_to_me = received_offsets.at(pid);
4044 
4045  libmesh_assert_equal_to (pushed_node_ids_to_me.size(),
4046  pushed_node_keys_vals_to_me.size());
4047  libmesh_assert_equal_to (pushed_node_ids_to_me.size(),
4048  pushed_offsets_to_me.size());
4049 
4050  for (auto i : index_range(pushed_node_ids_to_me))
4051  {
4052  dof_id_type constrained_id = pushed_node_ids_to_me[i];
4053 
4054  // If we don't already have a constraint for this node,
4055  // add the one we were sent
4056  const Node * constrained = mesh.node_ptr(constrained_id);
4057  if (!this->is_constrained_node(constrained))
4058  {
4059  NodeConstraintRow & row = _node_constraints[constrained].first;
4060  for (auto & kv : pushed_node_keys_vals_to_me[i])
4061  {
4062  const Node * key_node = mesh.node_ptr(kv.first);
4063  libmesh_assert(key_node);
4064  row[key_node] = kv.second;
4065  }
4066  _node_constraints[constrained].second = pushed_offsets_to_me[i];
4067  }
4068  }
4069  }
4070 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4071  }
4072 
4073  // Now start checking for any other constraints we need
4074  // to know about, requesting them recursively.
4075 
4076  // Create sets containing the DOFs and nodes we already depend on
4077  typedef std::set<dof_id_type> DoF_RCSet;
4078  DoF_RCSet unexpanded_dofs;
4079 
4080  for (const auto & i : _dof_constraints)
4081  unexpanded_dofs.insert(i.first);
4082 
4083  // Gather all the dof constraints we need
4084  this->gather_constraints(mesh, unexpanded_dofs, false);
4085 
4086  // Gather all the node constraints we need
4087 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4088  typedef std::set<const Node *> Node_RCSet;
4089  Node_RCSet unexpanded_nodes;
4090 
4091  for (const auto & i : _node_constraints)
4092  unexpanded_nodes.insert(i.first);
4093 
4094  // We have to keep recursing while the unexpanded set is
4095  // nonempty on *any* processor
4096  bool unexpanded_set_nonempty = !unexpanded_nodes.empty();
4097  this->comm().max(unexpanded_set_nonempty);
4098 
4099  while (unexpanded_set_nonempty)
4100  {
4101  // Let's make sure we don't lose sync in this loop.
4102  parallel_object_only();
4103 
4104  // Request sets
4105  Node_RCSet node_request_set;
4106 
4107  // Request sets to send to each processor
4108  std::map<processor_id_type, std::vector<dof_id_type>>
4109  requested_node_ids;
4110 
4111  // And the sizes of each
4112  std::map<processor_id_type, dof_id_type> node_ids_on_proc;
4113 
4114  // Fill (and thereby sort and uniq!) the main request sets
4115  for (const auto & i : unexpanded_nodes)
4116  {
4117  NodeConstraintRow & row = _node_constraints[i].first;
4118  for (const auto & j : row)
4119  {
4120  const Node * const node = j.first;
4121  libmesh_assert(node);
4122 
4123  // If it's non-local and we haven't already got a
4124  // constraint for it, we might need to ask for one
4125  if ((node->processor_id() != this->processor_id()) &&
4126  !_node_constraints.count(node))
4127  node_request_set.insert(node);
4128  }
4129  }
4130 
4131  // Clear the unexpanded constraint sets; we're about to expand
4132  // them
4133  unexpanded_nodes.clear();
4134 
4135  // Count requests by processor
4136  for (const auto & node : node_request_set)
4137  {
4138  libmesh_assert(node);
4139  libmesh_assert_less (node->processor_id(), this->n_processors());
4140  node_ids_on_proc[node->processor_id()]++;
4141  }
4142 
4143  for (auto pair : node_ids_on_proc)
4144  requested_node_ids[pair.first].reserve(pair.second);
4145 
4146  // Prepare each processor's request set
4147  for (const auto & node : node_request_set)
4148  requested_node_ids[node->processor_id()].push_back(node->id());
4149 
4150  typedef std::vector<std::pair<dof_id_type, Real>> row_datum;
4151 
4152  auto node_row_gather_functor =
4153  [this,
4154  & mesh,
4155  dist_mesh,
4156  & packed_range_sends,
4157  & range_tag]
4158  (processor_id_type pid,
4159  const std::vector<dof_id_type> & ids,
4160  std::vector<row_datum> & data)
4161  {
4162  // FIXME - this could be an unordered set, given a
4163  // hash<pointers> specialization
4164  std::set<const Node *> nodes_requested;
4165 
4166  // Fill those requests
4167  const std::size_t query_size = ids.size();
4168 
4169  data.resize(query_size);
4170  for (std::size_t i=0; i != query_size; ++i)
4171  {
4172  dof_id_type constrained_id = ids[i];
4173  const Node * constrained_node = mesh.node_ptr(constrained_id);
4174  if (_node_constraints.count(constrained_node))
4175  {
4176  const NodeConstraintRow & row = _node_constraints[constrained_node].first;
4177  std::size_t row_size = row.size();
4178  data[i].reserve(row_size);
4179  for (const auto & j : row)
4180  {
4181  const Node * node = j.first;
4182  data[i].emplace_back(node->id(), j.second);
4183 
4184  // If we're not sure whether our send
4185  // destination already has this node, let's give
4186  // it a copy.
4187  if (node->processor_id() != pid && dist_mesh)
4188  nodes_requested.insert(node);
4189 
4190  // We can have 0 nodal constraint
4191  // coefficients, where no Lagrange constraint
4192  // exists but non-Lagrange basis constraints
4193  // might.
4194  // libmesh_assert(j.second);
4195  }
4196  }
4197  else
4198  {
4199  // We have to distinguish "constraint with no
4200  // constraining nodes" (e.g. due to user node
4201  // constraint equations) from "no constraint".
4202  // We'll use invalid_id for the latter.
4203  data[i].emplace_back(DofObject::invalid_id, Real(0));
4204  }
4205  }
4206 
4207  // Constraining nodes might not even exist on our
4208  // correspondant's subset of a distributed mesh, so let's
4209  // make them exist.
4210  if (dist_mesh)
4211  {
4212  packed_range_sends.push_back(Parallel::Request());
4213  this->comm().send_packed_range
4214  (pid, &mesh, nodes_requested.begin(), nodes_requested.end(),
4215  packed_range_sends.back(), range_tag);
4216  }
4217  };
4218 
4219  typedef Point node_rhs_datum;
4220 
4221  auto node_rhs_gather_functor =
4222  [this,
4223  & mesh]
4225  const std::vector<dof_id_type> & ids,
4226  std::vector<node_rhs_datum> & data)
4227  {
4228  // Fill those requests
4229  const std::size_t query_size = ids.size();
4230 
4231  data.resize(query_size);
4232  for (std::size_t i=0; i != query_size; ++i)
4233  {
4234  dof_id_type constrained_id = ids[i];
4235  const Node * constrained_node = mesh.node_ptr(constrained_id);
4236  if (_node_constraints.count(constrained_node))
4237  data[i] = _node_constraints[constrained_node].second;
4238  else
4239  data[i](0) = std::numeric_limits<Real>::quiet_NaN();
4240  }
4241  };
4242 
4243  auto node_row_action_functor =
4244  [this,
4245  & mesh,
4246  dist_mesh,
4247  & range_tag,
4248  & unexpanded_nodes]
4249  (processor_id_type pid,
4250  const std::vector<dof_id_type> & ids,
4251  const std::vector<row_datum> & data)
4252  {
4253  // Before we act on any new constraint rows, we may need to
4254  // make sure we have all the nodes involved!
4255  if (dist_mesh)
4256  this->comm().receive_packed_range
4257  (pid, &mesh, null_output_iterator<Node>(),
4258  (Node**)nullptr, range_tag);
4259 
4260  // Add any new constraint rows we've found
4261  const std::size_t query_size = ids.size();
4262 
4263  for (std::size_t i=0; i != query_size; ++i)
4264  {
4265  const dof_id_type constrained_id = ids[i];
4266 
4267  // An empty row is an constraint with an empty row; for
4268  // no constraint we use a "no row" placeholder
4269  if (data[i].empty())
4270  {
4271  const Node * constrained_node = mesh.node_ptr(constrained_id);
4272  NodeConstraintRow & row = _node_constraints[constrained_node].first;
4273  row.clear();
4274  }
4275  else if (data[i][0].first != DofObject::invalid_id)
4276  {
4277  const Node * constrained_node = mesh.node_ptr(constrained_id);
4278  NodeConstraintRow & row = _node_constraints[constrained_node].first;
4279  row.clear();
4280  for (auto & pair : data[i])
4281  {
4282  const Node * key_node =
4283  mesh.node_ptr(pair.first);
4284  libmesh_assert(key_node);
4285  row[key_node] = pair.second;
4286  }
4287 
4288  // And prepare to check for more recursive constraints
4289  unexpanded_nodes.insert(constrained_node);
4290  }
4291  }
4292  };
4293 
4294  auto node_rhs_action_functor =
4295  [this,
4296  & mesh]
4298  const std::vector<dof_id_type> & ids,
4299  const std::vector<node_rhs_datum> & data)
4300  {
4301  // Add rhs data for any new node constraint rows we've found
4302  const std::size_t query_size = ids.size();
4303 
4304  for (std::size_t i=0; i != query_size; ++i)
4305  {
4306  dof_id_type constrained_id = ids[i];
4307  const Node * constrained_node = mesh.node_ptr(constrained_id);
4308 
4309  if (!libmesh_isnan(data[i](0)))
4310  _node_constraints[constrained_node].second = data[i];
4311  else
4312  _node_constraints.erase(constrained_node);
4313  }
4314  };
4315 
4316  // Now request node constraint rows from other processors
4317  row_datum * node_row_ex = nullptr;
4318  Parallel::pull_parallel_vector_data
4319  (this->comm(), requested_node_ids, node_row_gather_functor,
4320  node_row_action_functor, node_row_ex);
4321 
4322  // And request node constraint right hand sides from other procesors
4323  node_rhs_datum * node_rhs_ex = nullptr;
4324  Parallel::pull_parallel_vector_data
4325  (this->comm(), requested_node_ids, node_rhs_gather_functor,
4326  node_rhs_action_functor, node_rhs_ex);
4327 
4328 
4329  // We have to keep recursing while the unexpanded set is
4330  // nonempty on *any* processor
4331  unexpanded_set_nonempty = !unexpanded_nodes.empty();
4332  this->comm().max(unexpanded_set_nonempty);
4333  }
4334  Parallel::wait(packed_range_sends);
4335 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4336 }
4337 
4338 
4339 
4341 {
4342  // We've computed our local constraints, but they may depend on
4343  // non-local constraints that we'll need to take into account.
4344  this->allgather_recursive_constraints(mesh);
4345 
4347  {
4348  // Optionally check for constraint loops and throw an error
4349  // if they're detected. We always do this check below in dbg/devel
4350  // mode but here we optionally do it in opt mode as well.
4352  }
4353 
4354  // Adjoints will be constrained where the primal is
4355  // Therefore, we will expand the adjoint_constraint_values
4356  // map whenever the primal_constraint_values map is expanded
4357 
4358  // First, figure out the total number of QoIs
4359  const unsigned int max_qoi_num =
4360  _adjoint_constraint_values.empty() ?
4361  0 : _adjoint_constraint_values.rbegin()->first+1;
4362 
4363  // Create a set containing the DOFs we already depend on
4364  typedef std::set<dof_id_type> RCSet;
4365  RCSet unexpanded_set;
4366 
4367  for (const auto & i : _dof_constraints)
4368  unexpanded_set.insert(i.first);
4369 
4370  while (!unexpanded_set.empty())
4371  for (RCSet::iterator i = unexpanded_set.begin();
4372  i != unexpanded_set.end(); /* nothing */)
4373  {
4374  // If the DOF is constrained
4375  DofConstraints::iterator
4376  pos = _dof_constraints.find(*i);
4377 
4378  libmesh_assert (pos != _dof_constraints.end());
4379 
4380  DofConstraintRow & constraint_row = pos->second;
4381 
4382  DofConstraintValueMap::iterator rhsit =
4383  _primal_constraint_values.find(*i);
4384  Number constraint_rhs = (rhsit == _primal_constraint_values.end()) ?
4385  0 : rhsit->second;
4386 
4387  // A vector of DofConstraintValueMaps for each adjoint variable
4388  std::vector<DofConstraintValueMap::iterator> adjoint_rhs_iterators;
4389  adjoint_rhs_iterators.resize(max_qoi_num);
4390 
4391  // Another to hold the adjoint constraint rhs
4392  std::vector<Number> adjoint_constraint_rhs(max_qoi_num, 0.0);
4393 
4394  // Find and gather recursive constraints for each adjoint variable
4395  for (auto & adjoint_map : _adjoint_constraint_values)
4396  {
4397  const std::size_t q = adjoint_map.first;
4398  adjoint_rhs_iterators[q] = adjoint_map.second.find(*i);
4399 
4400  adjoint_constraint_rhs[q] =
4401  (adjoint_rhs_iterators[q] == adjoint_map.second.end()) ?
4402  0 : adjoint_rhs_iterators[q]->second;
4403  }
4404 
4405  std::vector<dof_id_type> constraints_to_expand;
4406 
4407  for (const auto & item : constraint_row)
4408  if (item.first != *i && this->is_constrained_dof(item.first))
4409  {
4410  unexpanded_set.insert(item.first);
4411  constraints_to_expand.push_back(item.first);
4412  }
4413 
4414  for (const auto & expandable : constraints_to_expand)
4415  {
4416  const Real this_coef = constraint_row[expandable];
4417 
4418  DofConstraints::const_iterator
4419  subpos = _dof_constraints.find(expandable);
4420 
4421  libmesh_assert (subpos != _dof_constraints.end());
4422 
4423  const DofConstraintRow & subconstraint_row = subpos->second;
4424 
4425  for (const auto & item : subconstraint_row)
4426  {
4427  // Assert that the constraint does not form a cycle.
4428  libmesh_assert(item.first != expandable);
4429  constraint_row[item.first] += item.second * this_coef;
4430  }
4431 
4432  if (auto subrhsit = _primal_constraint_values.find(expandable);
4433  subrhsit != _primal_constraint_values.end())
4434  constraint_rhs += subrhsit->second * this_coef;
4435 
4436  // Find and gather recursive constraints for each adjoint variable
4437  for (const auto & adjoint_map : _adjoint_constraint_values)
4438  {
4439  if (auto adjoint_subrhsit = adjoint_map.second.find(expandable);
4440  adjoint_subrhsit != adjoint_map.second.end())
4441  adjoint_constraint_rhs[adjoint_map.first] += adjoint_subrhsit->second * this_coef;
4442  }
4443 
4444  constraint_row.erase(expandable);
4445  }
4446 
4447  if (rhsit == _primal_constraint_values.end())
4448  {
4449  if (constraint_rhs != Number(0))
4450  _primal_constraint_values[*i] = constraint_rhs;
4451  else
4452  _primal_constraint_values.erase(*i);
4453  }
4454  else
4455  {
4456  if (constraint_rhs != Number(0))
4457  rhsit->second = constraint_rhs;
4458  else
4459  _primal_constraint_values.erase(rhsit);
4460  }
4461 
4462  // Finally fill in the adjoint constraints for each adjoint variable if possible
4463  for (auto & adjoint_map : _adjoint_constraint_values)
4464  {
4465  const std::size_t q = adjoint_map.first;
4466 
4467  if(adjoint_rhs_iterators[q] == adjoint_map.second.end())
4468  {
4469  if (adjoint_constraint_rhs[q] != Number(0))
4470  (adjoint_map.second)[*i] = adjoint_constraint_rhs[q];
4471  else
4472  adjoint_map.second.erase(*i);
4473  }
4474  else
4475  {
4476  if (adjoint_constraint_rhs[q] != Number(0))
4477  adjoint_rhs_iterators[q]->second = adjoint_constraint_rhs[q];
4478  else
4479  adjoint_map.second.erase(adjoint_rhs_iterators[q]);
4480  }
4481  }
4482 
4483  if (constraints_to_expand.empty())
4484  i = unexpanded_set.erase(i);
4485  else
4486  ++i;
4487  }
4488 
4489  // In parallel we can't guarantee that nodes/dofs which constrain
4490  // others are on processors which are aware of that constraint, yet
4491  // we need such awareness for sparsity pattern generation. So send
4492  // other processors any constraints they might need to know about.
4493  this->scatter_constraints(mesh);
4494 
4495  // Now that we have our root constraint dependencies sorted out, add
4496  // them to the send_list
4498 }
4499 
4500 
4501 #ifdef LIBMESH_ENABLE_CONSTRAINTS
4503 {
4504  // Eventually make this officially libmesh_deprecated();
4506 }
4507 
4509 {
4510  // Create a set containing the DOFs we already depend on
4511  typedef std::set<dof_id_type> RCSet;
4512  RCSet unexpanded_set;
4513 
4514  // Use dof_constraints_copy in this method so that we don't
4515  // mess with _dof_constraints.
4516  DofConstraints dof_constraints_copy = _dof_constraints;
4517 
4518  for (const auto & i : dof_constraints_copy)
4519  unexpanded_set.insert(i.first);
4520 
4521  while (!unexpanded_set.empty())
4522  for (RCSet::iterator i = unexpanded_set.begin();
4523  i != unexpanded_set.end(); /* nothing */)
4524  {
4525  // If the DOF is constrained
4526  DofConstraints::iterator
4527  pos = dof_constraints_copy.find(*i);
4528 
4529  libmesh_assert (pos != dof_constraints_copy.end());
4530 
4531  DofConstraintRow & constraint_row = pos->second;
4532 
4533  // Comment out "rhs" parts of this method copied from process_constraints
4534  // DofConstraintValueMap::iterator rhsit =
4535  // _primal_constraint_values.find(*i);
4536  // Number constraint_rhs = (rhsit == _primal_constraint_values.end()) ?
4537  // 0 : rhsit->second;
4538 
4539  std::vector<dof_id_type> constraints_to_expand;
4540 
4541  for (const auto & item : constraint_row)
4542  if (item.first != *i && this->is_constrained_dof(item.first))
4543  {
4544  unexpanded_set.insert(item.first);
4545  constraints_to_expand.push_back(item.first);
4546  }
4547 
4548  for (const auto & expandable : constraints_to_expand)
4549  {
4550  const Real this_coef = constraint_row[expandable];
4551 
4552  DofConstraints::const_iterator
4553  subpos = dof_constraints_copy.find(expandable);
4554 
4555  libmesh_assert (subpos != dof_constraints_copy.end());
4556 
4557  const DofConstraintRow & subconstraint_row = subpos->second;
4558 
4559  for (const auto & item : subconstraint_row)
4560  {
4561  libmesh_error_msg_if(item.first == expandable, "Constraint loop detected");
4562 
4563  constraint_row[item.first] += item.second * this_coef;
4564  }
4565 
4566  // Comment out "rhs" parts of this method copied from process_constraints
4567  // DofConstraintValueMap::const_iterator subrhsit =
4568  // _primal_constraint_values.find(expandable);
4569  // if (subrhsit != _primal_constraint_values.end())
4570  // constraint_rhs += subrhsit->second * this_coef;
4571 
4572  constraint_row.erase(expandable);
4573  }
4574 
4575  // Comment out "rhs" parts of this method copied from process_constraints
4576  // if (rhsit == _primal_constraint_values.end())
4577  // {
4578  // if (constraint_rhs != Number(0))
4579  // _primal_constraint_values[*i] = constraint_rhs;
4580  // else
4581  // _primal_constraint_values.erase(*i);
4582  // }
4583  // else
4584  // {
4585  // if (constraint_rhs != Number(0))
4586  // rhsit->second = constraint_rhs;
4587  // else
4588  // _primal_constraint_values.erase(rhsit);
4589  // }
4590 
4591  if (constraints_to_expand.empty())
4592  i = unexpanded_set.erase(i);
4593  else
4594  ++i;
4595  }
4596 }
4597 #else
4600 {
4601  // Do nothing
4602 }
4603 #endif
4604 
4605 
4607 {
4608  // At this point each processor with a constrained node knows
4609  // the corresponding constraint row, but we also need each processor
4610  // with a constrainer node to know the corresponding row(s).
4611 
4612  // This function must be run on all processors at once
4613  parallel_object_only();
4614 
4615  // Return immediately if there's nothing to gather
4616  if (this->n_processors() == 1)
4617  return;
4618 
4619  // We might get to return immediately if none of the processors
4620  // found any constraints
4621  unsigned int has_constraints = !_dof_constraints.empty()
4622 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4623  || !_node_constraints.empty()
4624 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4625  ;
4626  this->comm().max(has_constraints);
4627  if (!has_constraints)
4628  return;
4629 
4630  // We may be receiving packed_range sends out of order with
4631  // parallel_sync tags, so make sure they're received correctly.
4632  Parallel::MessageTag range_tag = this->comm().get_unique_tag();
4633 
4634 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4635  std::map<processor_id_type, std::set<dof_id_type>> pushed_node_ids;
4636 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4637 
4638  std::map<processor_id_type, std::set<dof_id_type>> pushed_ids;
4639 
4640  // Collect the dof constraints I need to push to each processor
4641  dof_id_type constrained_proc_id = 0;
4642  for (const auto & [constrained, row] : _dof_constraints)
4643  {
4644  while (constrained >= _end_df[constrained_proc_id])
4645  constrained_proc_id++;
4646 
4647  if (constrained_proc_id != this->processor_id())
4648  continue;
4649 
4650  for (auto & j : row)
4651  {
4652  const dof_id_type constraining = j.first;
4653 
4654  processor_id_type constraining_proc_id = 0;
4655  while (constraining >= _end_df[constraining_proc_id])
4656  constraining_proc_id++;
4657 
4658  if (constraining_proc_id != this->processor_id() &&
4659  constraining_proc_id != constrained_proc_id)
4660  pushed_ids[constraining_proc_id].insert(constrained);
4661  }
4662  }
4663 
4664  // Pack the dof constraint rows and rhs's to push
4665 
4666  std::map<processor_id_type,
4667  std::vector<std::vector<std::pair<dof_id_type, Real>>>>
4668  pushed_keys_vals, pushed_keys_vals_to_me;
4669 
4670  std::map<processor_id_type, std::vector<std::pair<dof_id_type, Number>>>
4671  pushed_ids_rhss, pushed_ids_rhss_to_me;
4672 
4673  auto gather_ids =
4674  [this,
4675  & pushed_ids,
4676  & pushed_keys_vals,
4677  & pushed_ids_rhss]
4678  ()
4679  {
4680  for (const auto & [pid, pid_ids] : pushed_ids)
4681  {
4682  const std::size_t ids_size = pid_ids.size();
4683  std::vector<std::vector<std::pair<dof_id_type, Real>>> &
4684  keys_vals = pushed_keys_vals[pid];
4685  std::vector<std::pair<dof_id_type,Number>> &
4686  ids_rhss = pushed_ids_rhss[pid];
4687  keys_vals.resize(ids_size);
4688  ids_rhss.resize(ids_size);
4689 
4690  std::size_t push_i;
4691  std::set<dof_id_type>::const_iterator it;
4692  for (push_i = 0, it = pid_ids.begin();
4693  it != pid_ids.end(); ++push_i, ++it)
4694  {
4695  const dof_id_type constrained = *it;
4696  DofConstraintRow & row = _dof_constraints[constrained];
4697  keys_vals[push_i].assign(row.begin(), row.end());
4698 
4699  DofConstraintValueMap::const_iterator rhsit =
4700  _primal_constraint_values.find(constrained);
4701  ids_rhss[push_i].first = constrained;
4702  ids_rhss[push_i].second =
4703  (rhsit == _primal_constraint_values.end()) ?
4704  0 : rhsit->second;
4705  }
4706  }
4707  };
4708 
4709  gather_ids();
4710 
4711  auto ids_rhss_action_functor =
4712  [& pushed_ids_rhss_to_me]
4713  (processor_id_type pid,
4714  const std::vector<std::pair<dof_id_type, Number>> & data)
4715  {
4716  pushed_ids_rhss_to_me[pid] = data;
4717  };
4718 
4719  auto keys_vals_action_functor =
4720  [& pushed_keys_vals_to_me]
4721  (processor_id_type pid,
4722  const std::vector<std::vector<std::pair<dof_id_type, Real>>> & data)
4723  {
4724  pushed_keys_vals_to_me[pid] = data;
4725  };
4726 
4727  Parallel::push_parallel_vector_data
4728  (this->comm(), pushed_ids_rhss, ids_rhss_action_functor);
4729  Parallel::push_parallel_vector_data
4730  (this->comm(), pushed_keys_vals, keys_vals_action_functor);
4731 
4732  // Now work on traded dof constraint rows
4733  auto receive_dof_constraints =
4734  [this,
4735  & pushed_ids_rhss_to_me,
4736  & pushed_keys_vals_to_me]
4737  ()
4738  {
4739  for (const auto & [pid, ids_rhss] : pushed_ids_rhss_to_me)
4740  {
4741  const auto & keys_vals = pushed_keys_vals_to_me[pid];
4742 
4743  libmesh_assert_equal_to
4744  (ids_rhss.size(), keys_vals.size());
4745 
4746  // Add the dof constraints that I've been sent
4747  for (auto i : index_range(ids_rhss))
4748  {
4749  dof_id_type constrained = ids_rhss[i].first;
4750 
4751  // If we don't already have a constraint for this dof,
4752  // add the one we were sent
4753  if (!this->is_constrained_dof(constrained))
4754  {
4755  DofConstraintRow & row = _dof_constraints[constrained];
4756  for (auto & key_val : keys_vals[i])
4757  {
4758  libmesh_assert_less(key_val.first, this->n_dofs());
4759  row[key_val.first] = key_val.second;
4760  }
4761  if (ids_rhss[i].second != Number(0))
4762  _primal_constraint_values[constrained] =
4763  ids_rhss[i].second;
4764  else
4765  _primal_constraint_values.erase(constrained);
4766  }
4767  }
4768  }
4769  };
4770 
4771  receive_dof_constraints();
4772 
4773 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
4774  // Collect the node constraints to push to each processor
4775  for (auto & i : _node_constraints)
4776  {
4777  const Node * constrained = i.first;
4778 
4779  if (constrained->processor_id() != this->processor_id())
4780  continue;
4781 
4782  NodeConstraintRow & row = i.second.first;
4783  for (auto & j : row)
4784  {
4785  const Node * constraining = j.first;
4786 
4787  if (constraining->processor_id() != this->processor_id() &&
4788  constraining->processor_id() != constrained->processor_id())
4789  pushed_node_ids[constraining->processor_id()].insert(constrained->id());
4790  }
4791  }
4792 
4793  // Pack the node constraint rows and rhss to push
4794  std::map<processor_id_type,
4795  std::vector<std::vector<std::pair<dof_id_type,Real>>>>
4796  pushed_node_keys_vals, pushed_node_keys_vals_to_me;
4797  std::map<processor_id_type, std::vector<std::pair<dof_id_type, Point>>>
4798  pushed_node_ids_offsets, pushed_node_ids_offsets_to_me;
4799  std::map<processor_id_type, std::vector<const Node *>> pushed_node_vecs;
4800 
4801  for (const auto & [pid, pid_ids]: pushed_node_ids)
4802  {
4803  const std::size_t ids_size = pid_ids.size();
4804  std::vector<std::vector<std::pair<dof_id_type,Real>>> &
4805  keys_vals = pushed_node_keys_vals[pid];
4806  std::vector<std::pair<dof_id_type, Point>> &
4807  ids_offsets = pushed_node_ids_offsets[pid];
4808  keys_vals.resize(ids_size);
4809  ids_offsets.resize(ids_size);
4810  std::set<Node *> nodes;
4811 
4812  std::size_t push_i;
4813  std::set<dof_id_type>::const_iterator it;
4814  for (push_i = 0, it = pid_ids.begin();
4815  it != pid_ids.end(); ++push_i, ++it)
4816  {
4817  Node * constrained = mesh.node_ptr(*it);
4818 
4819  if (constrained->processor_id() != pid)
4820  nodes.insert(constrained);
4821 
4822  NodeConstraintRow & row = _node_constraints[constrained].first;
4823  std::size_t row_size = row.size();
4824  keys_vals[push_i].reserve(row_size);
4825  for (const auto & j : row)
4826  {
4827  Node * constraining = const_cast<Node *>(j.first);
4828 
4829  keys_vals[push_i].emplace_back(constraining->id(), j.second);
4830 
4831  if (constraining->processor_id() != pid)
4832  nodes.insert(constraining);
4833  }
4834 
4835  ids_offsets[push_i].first = *it;
4836  ids_offsets[push_i].second = _node_constraints[constrained].second;
4837  }
4838 
4839  if (!mesh.is_serial())
4840  {
4841  auto & pid_nodes = pushed_node_vecs[pid];
4842  pid_nodes.assign(nodes.begin(), nodes.end());
4843  }
4844  }
4845 
4846  auto node_ids_offsets_action_functor =
4847  [& pushed_node_ids_offsets_to_me]
4848  (processor_id_type pid,
4849  const std::vector<std::pair<dof_id_type, Point>> & data)
4850  {
4851  pushed_node_ids_offsets_to_me[pid] = data;
4852  };
4853 
4854  auto node_keys_vals_action_functor =
4855  [& pushed_node_keys_vals_to_me]
4856  (processor_id_type pid,
4857  const std::vector<std::vector<std::pair<dof_id_type, Real>>> & data)
4858  {
4859  pushed_node_keys_vals_to_me[pid] = data;
4860  };
4861 
4862  // Trade pushed node constraint rows
4863  Parallel::push_parallel_vector_data
4864  (this->comm(), pushed_node_ids_offsets, node_ids_offsets_action_functor);
4865  Parallel::push_parallel_vector_data
4866  (this->comm(), pushed_node_keys_vals, node_keys_vals_action_functor);
4867 
4868  // Constraining nodes might not even exist on our subset of a
4869  // distributed mesh, so let's make them exist.
4870 
4871  // Node unpack() now automatically adds them to the context mesh
4872  auto null_node_functor = [](processor_id_type, const std::vector<const Node *> &){};
4873 
4874  if (!mesh.is_serial())
4875  Parallel::push_parallel_packed_range
4876  (this->comm(), pushed_node_vecs, &mesh, null_node_functor);
4877 
4878  for (const auto & [pid, ids_offsets] : pushed_node_ids_offsets_to_me)
4879  {
4880  const auto & keys_vals = pushed_node_keys_vals_to_me[pid];
4881 
4882  libmesh_assert_equal_to
4883  (ids_offsets.size(), keys_vals.size());
4884 
4885  // Add the node constraints that I've been sent
4886  for (auto i : index_range(ids_offsets))
4887  {
4888  dof_id_type constrained_id = ids_offsets[i].first;
4889 
4890  // If we don't already have a constraint for this node,
4891  // add the one we were sent
4892  const Node * constrained = mesh.node_ptr(constrained_id);
4893  if (!this->is_constrained_node(constrained))
4894  {
4895  NodeConstraintRow & row = _node_constraints[constrained].first;
4896  for (auto & key_val : keys_vals[i])
4897  {
4898  const Node * key_node = mesh.node_ptr(key_val.first);
4899  row[key_node] = key_val.second;
4900  }
4901  _node_constraints[constrained].second =
4902  ids_offsets[i].second;
4903  }
4904  }
4905  }
4906 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
4907 
4908  // Next we need to push constraints to processors which don't own
4909  // the constrained dof, don't own the constraining dof, but own an
4910  // element supporting the constraining dof.
4911  //
4912  // We need to be able to quickly look up constrained dof ids by what
4913  // constrains them, so that we can handle the case where we see a
4914  // foreign element containing one of our constraining DoF ids and we
4915  // need to push that constraint.
4916  //
4917  // Getting distributed adaptive sparsity patterns right is hard.
4918 
4919  typedef std::map<dof_id_type, std::set<dof_id_type>> DofConstrainsMap;
4920  DofConstrainsMap dof_id_constrains;
4921 
4922  for (const auto & [constrained, row] : _dof_constraints)
4923  {
4924  for (const auto & j : row)
4925  {
4926  const dof_id_type constraining = j.first;
4927 
4928  dof_id_type constraining_proc_id = 0;
4929  while (constraining >= _end_df[constraining_proc_id])
4930  constraining_proc_id++;
4931 
4932  if (constraining_proc_id == this->processor_id())
4933  dof_id_constrains[constraining].insert(constrained);
4934  }
4935  }
4936 
4937  // Loop over all foreign elements, find any supporting our
4938  // constrained dof indices.
4939  pushed_ids.clear();
4940 
4941  for (const auto & elem : as_range(mesh.active_not_local_elements_begin(),
4942  mesh.active_not_local_elements_end()))
4943  {
4944  std::vector<dof_id_type> my_dof_indices;
4945  this->dof_indices (elem, my_dof_indices);
4946 
4947  for (const auto & dof : my_dof_indices)
4948  {
4949  if (auto dcmi = dof_id_constrains.find(dof);
4950  dcmi != dof_id_constrains.end())
4951  {
4952  for (const auto & constrained : dcmi->second)
4953  {
4954  dof_id_type the_constrained_proc_id = 0;
4955  while (constrained >= _end_df[the_constrained_proc_id])
4956  the_constrained_proc_id++;
4957 
4958  const processor_id_type elemproc = elem->processor_id();
4959  if (elemproc != the_constrained_proc_id)
4960  pushed_ids[elemproc].insert(constrained);
4961  }
4962  }
4963  }
4964  }
4965 
4966  pushed_ids_rhss.clear();
4967  pushed_ids_rhss_to_me.clear();
4968  pushed_keys_vals.clear();
4969  pushed_keys_vals_to_me.clear();
4970 
4971  gather_ids();
4972 
4973  // Trade pushed dof constraint rows
4974  Parallel::push_parallel_vector_data
4975  (this->comm(), pushed_ids_rhss, ids_rhss_action_functor);
4976  Parallel::push_parallel_vector_data
4977  (this->comm(), pushed_keys_vals, keys_vals_action_functor);
4978 
4979  receive_dof_constraints();
4980 
4981  // Finally, we need to handle the case of remote dof coupling. If a
4982  // processor's element is coupled to a ghost element, then the
4983  // processor needs to know about all constraints which affect the
4984  // dofs on that ghost element, so we'll have to query the ghost
4985  // element's owner.
4986 
4987  GhostingFunctor::map_type elements_to_couple;
4988  DofMap::CouplingMatricesSet temporary_coupling_matrices;
4989 
4991  (elements_to_couple,
4992  temporary_coupling_matrices,
4993  this->coupling_functors_begin(),
4994  this->coupling_functors_end(),
4995  mesh.active_local_elements_begin(),
4996  mesh.active_local_elements_end(),
4997  this->processor_id());
4998 
4999  // Each ghost-coupled element's owner should get a request for its dofs
5000  std::set<dof_id_type> requested_dofs;
5001 
5002  for (const auto & pr : elements_to_couple)
5003  {
5004  const Elem * elem = pr.first;
5005 
5006  // FIXME - optimize for the non-fully-coupled case?
5007  std::vector<dof_id_type> element_dofs;
5008  this->dof_indices(elem, element_dofs);
5009 
5010  for (auto dof : element_dofs)
5011  requested_dofs.insert(dof);
5012  }
5013 
5014  this->gather_constraints(mesh, requested_dofs, false);
5015 }
5016 
5017 
5019  std::set<dof_id_type> & unexpanded_dofs,
5020  bool /*look_for_constrainees*/)
5021 {
5022  typedef std::set<dof_id_type> DoF_RCSet;
5023 
5024  // If we have heterogeneous adjoint constraints we need to
5025  // communicate those too.
5026  const unsigned int max_qoi_num =
5027  _adjoint_constraint_values.empty() ?
5028  0 : _adjoint_constraint_values.rbegin()->first+1;
5029 
5030  // We have to keep recursing while the unexpanded set is
5031  // nonempty on *any* processor
5032  bool unexpanded_set_nonempty = !unexpanded_dofs.empty();
5033  this->comm().max(unexpanded_set_nonempty);
5034 
5035  while (unexpanded_set_nonempty)
5036  {
5037  // Let's make sure we don't lose sync in this loop.
5038  parallel_object_only();
5039 
5040  // Request sets
5041  DoF_RCSet dof_request_set;
5042 
5043  // Request sets to send to each processor
5044  std::map<processor_id_type, std::vector<dof_id_type>>
5045  requested_dof_ids;
5046 
5047  // And the sizes of each
5048  std::map<processor_id_type, dof_id_type>
5049  dof_ids_on_proc;
5050 
5051  // Fill (and thereby sort and uniq!) the main request sets
5052  for (const auto & unexpanded_dof : unexpanded_dofs)
5053  {
5054  // If we were asked for a DoF and we don't already have a
5055  // constraint for it, then we need to check for one.
5056  if (auto pos = _dof_constraints.find(unexpanded_dof);
5057  pos == _dof_constraints.end())
5058  {
5059  if (!this->local_index(unexpanded_dof) &&
5060  !_dof_constraints.count(unexpanded_dof) )
5061  dof_request_set.insert(unexpanded_dof);
5062  }
5063  // If we were asked for a DoF and we already have a
5064  // constraint for it, then we need to check if the
5065  // constraint is recursive.
5066  else
5067  {
5068  const DofConstraintRow & row = pos->second;
5069  for (const auto & j : row)
5070  {
5071  const dof_id_type constraining_dof = j.first;
5072 
5073  // If it's non-local and we haven't already got a
5074  // constraint for it, we might need to ask for one
5075  if (!this->local_index(constraining_dof) &&
5076  !_dof_constraints.count(constraining_dof))
5077  dof_request_set.insert(constraining_dof);
5078  }
5079  }
5080  }
5081 
5082  // Clear the unexpanded constraint set; we're about to expand it
5083  unexpanded_dofs.clear();
5084 
5085  // Count requests by processor
5086  processor_id_type proc_id = 0;
5087  for (const auto & i : dof_request_set)
5088  {
5089  while (i >= _end_df[proc_id])
5090  proc_id++;
5091  dof_ids_on_proc[proc_id]++;
5092  }
5093 
5094  for (auto & pair : dof_ids_on_proc)
5095  {
5096  requested_dof_ids[pair.first].reserve(pair.second);
5097  }
5098 
5099  // Prepare each processor's request set
5100  proc_id = 0;
5101  for (const auto & i : dof_request_set)
5102  {
5103  while (i >= _end_df[proc_id])
5104  proc_id++;
5105  requested_dof_ids[proc_id].push_back(i);
5106  }
5107 
5108  typedef std::vector<std::pair<dof_id_type, Real>> row_datum;
5109 
5110  typedef std::vector<Number> rhss_datum;
5111 
5112  auto row_gather_functor =
5113  [this]
5115  const std::vector<dof_id_type> & ids,
5116  std::vector<row_datum> & data)
5117  {
5118  // Fill those requests
5119  const std::size_t query_size = ids.size();
5120 
5121  data.resize(query_size);
5122  for (std::size_t i=0; i != query_size; ++i)
5123  {
5124  dof_id_type constrained = ids[i];
5125  if (_dof_constraints.count(constrained))
5126  {
5127  DofConstraintRow & row = _dof_constraints[constrained];
5128  std::size_t row_size = row.size();
5129  data[i].reserve(row_size);
5130  for (const auto & j : row)
5131  {
5132  data[i].push_back(j);
5133 
5134  // We should never have an invalid constraining
5135  // dof id
5137 
5138  // We should never have a 0 constraint
5139  // coefficient; that's implicit via sparse
5140  // constraint storage
5141  //
5142  // But we can't easily control how users add
5143  // constraints, so we can't safely assert that
5144  // we're being efficient here.
5145  //
5146  // libmesh_assert(j.second);
5147  }
5148  }
5149  else
5150  {
5151  // We have to distinguish "constraint with no
5152  // constraining dofs" (e.g. due to Dirichlet
5153  // constraint equations) from "no constraint".
5154  // We'll use invalid_id for the latter.
5155  data[i].emplace_back(DofObject::invalid_id, Real(0));
5156  }
5157  }
5158  };
5159 
5160  auto rhss_gather_functor =
5161  [this,
5162  max_qoi_num]
5164  const std::vector<dof_id_type> & ids,
5165  std::vector<rhss_datum> & data)
5166  {
5167  // Fill those requests
5168  const std::size_t query_size = ids.size();
5169 
5170  data.resize(query_size);
5171  for (std::size_t i=0; i != query_size; ++i)
5172  {
5173  dof_id_type constrained = ids[i];
5174  data[i].clear();
5175  if (_dof_constraints.count(constrained))
5176  {
5177  DofConstraintValueMap::const_iterator rhsit =
5178  _primal_constraint_values.find(constrained);
5179  data[i].push_back
5180  ((rhsit == _primal_constraint_values.end()) ?
5181  0 : rhsit->second);
5182 
5183  for (unsigned int q = 0; q != max_qoi_num; ++q)
5184  {
5185  AdjointDofConstraintValues::const_iterator adjoint_map_it =
5187 
5188  if (adjoint_map_it == _adjoint_constraint_values.end())
5189  {
5190  data[i].push_back(0);
5191  continue;
5192  }
5193 
5194  const DofConstraintValueMap & constraint_map =
5195  adjoint_map_it->second;
5196 
5197  DofConstraintValueMap::const_iterator adj_rhsit =
5198  constraint_map.find(constrained);
5199  data[i].push_back
5200  ((adj_rhsit == constraint_map.end()) ?
5201  0 : adj_rhsit->second);
5202  }
5203  }
5204  }
5205  };
5206 
5207  auto row_action_functor =
5208  [this,
5209  & unexpanded_dofs]
5211  const std::vector<dof_id_type> & ids,
5212  const std::vector<row_datum> & data)
5213  {
5214  // Add any new constraint rows we've found
5215  const std::size_t query_size = ids.size();
5216 
5217  for (std::size_t i=0; i != query_size; ++i)
5218  {
5219  const dof_id_type constrained = ids[i];
5220 
5221  // An empty row is an constraint with an empty row; for
5222  // no constraint we use a "no row" placeholder
5223  if (data[i].empty())
5224  {
5225  DofConstraintRow & row = _dof_constraints[constrained];
5226  row.clear();
5227  }
5228  else if (data[i][0].first != DofObject::invalid_id)
5229  {
5230  DofConstraintRow & row = _dof_constraints[constrained];
5231  row.clear();
5232  for (auto & pair : data[i])
5233  {
5234  libmesh_assert_less(pair.first, this->n_dofs());
5235  row[pair.first] = pair.second;
5236  }
5237 
5238  // And prepare to check for more recursive constraints
5239  unexpanded_dofs.insert(constrained);
5240  }
5241  }
5242  };
5243 
5244  auto rhss_action_functor =
5245  [this,
5246  max_qoi_num]
5248  const std::vector<dof_id_type> & ids,
5249  const std::vector<rhss_datum> & data)
5250  {
5251  // Add rhs data for any new constraint rows we've found
5252  const std::size_t query_size = ids.size();
5253 
5254  for (std::size_t i=0; i != query_size; ++i)
5255  {
5256  if (!data[i].empty())
5257  {
5258  dof_id_type constrained = ids[i];
5259  if (data[i][0] != Number(0))
5260  _primal_constraint_values[constrained] = data[i][0];
5261  else
5262  _primal_constraint_values.erase(constrained);
5263 
5264  for (unsigned int q = 0; q != max_qoi_num; ++q)
5265  {
5266  AdjointDofConstraintValues::iterator adjoint_map_it =
5268 
5269  if ((adjoint_map_it == _adjoint_constraint_values.end()) &&
5270  data[i][q+1] == Number(0))
5271  continue;
5272 
5273  if (adjoint_map_it == _adjoint_constraint_values.end())
5274  adjoint_map_it = _adjoint_constraint_values.emplace
5275  (q, DofConstraintValueMap()).first;
5276 
5277  DofConstraintValueMap & constraint_map =
5278  adjoint_map_it->second;
5279 
5280  if (data[i][q+1] != Number(0))
5281  constraint_map[constrained] =
5282  data[i][q+1];
5283  else
5284  constraint_map.erase(constrained);
5285  }
5286  }
5287  }
5288 
5289  };
5290 
5291  // Now request constraint rows from other processors
5292  row_datum * row_ex = nullptr;
5293  Parallel::pull_parallel_vector_data
5294  (this->comm(), requested_dof_ids, row_gather_functor,
5295  row_action_functor, row_ex);
5296 
5297  // And request constraint right hand sides from other procesors
5298  rhss_datum * rhs_ex = nullptr;
5299  Parallel::pull_parallel_vector_data
5300  (this->comm(), requested_dof_ids, rhss_gather_functor,
5301  rhss_action_functor, rhs_ex);
5302 
5303  // We have to keep recursing while the unexpanded set is
5304  // nonempty on *any* processor
5305  unexpanded_set_nonempty = !unexpanded_dofs.empty();
5306  this->comm().max(unexpanded_set_nonempty);
5307  }
5308 }
5309 
5311 {
5312  // This function must be run on all processors at once
5313  parallel_object_only();
5314 
5315  // Return immediately if there's nothing to gather
5316  if (this->n_processors() == 1)
5317  return;
5318 
5319  // We might get to return immediately if none of the processors
5320  // found any constraints
5321  unsigned int has_constraints = !_dof_constraints.empty();
5322  this->comm().max(has_constraints);
5323  if (!has_constraints)
5324  return;
5325 
5326  for (const auto & [constrained_dof, constraint_row] : _dof_constraints)
5327  {
5328  // We only need the dependencies of our own constrained dofs
5329  if (!this->local_index(constrained_dof))
5330  continue;
5331 
5332  for (const auto & j : constraint_row)
5333  {
5334  dof_id_type constraint_dependency = j.first;
5335 
5336  // No point in adding one of our own dofs to the send_list
5337  if (this->local_index(constraint_dependency))
5338  continue;
5339 
5340  _send_list.push_back(constraint_dependency);
5341  }
5342  }
5343 }
5344 
5345 
5346 
5347 #endif // LIBMESH_ENABLE_CONSTRAINTS
5348 
5349 
5350 #ifdef LIBMESH_ENABLE_AMR
5351 
5352 void DofMap::constrain_p_dofs (unsigned int var,
5353  const Elem * elem,
5354  unsigned int s,
5355  unsigned int p)
5356 {
5357  // We're constraining dofs on elem which correspond to p refinement
5358  // levels above p - this only makes sense if elem's p refinement
5359  // level is above p.
5360  libmesh_assert_greater (elem->p_level(), p);
5361  libmesh_assert_less (s, elem->n_sides());
5362 
5363  const unsigned int sys_num = this->sys_number();
5364  FEType fe_type = this->variable_type(var);
5365 
5366  const unsigned int n_nodes = elem->n_nodes();
5367  for (unsigned int n = 0; n != n_nodes; ++n)
5368  if (elem->is_node_on_side(n, s))
5369  {
5370  const Node & node = elem->node_ref(n);
5371  const unsigned int low_nc =
5372  FEInterface::n_dofs_at_node (fe_type, p, elem, n);
5373  const unsigned int high_nc =
5374  FEInterface::n_dofs_at_node (fe_type, elem, n);
5375 
5376  // since we may be running this method concurrently
5377  // on multiple threads we need to acquire a lock
5378  // before modifying the _dof_constraints object.
5379  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
5380 
5381  if (elem->is_vertex(n))
5382  {
5383  // Add "this is zero" constraint rows for high p vertex
5384  // dofs
5385  for (unsigned int i = low_nc; i != high_nc; ++i)
5386  {
5387  _dof_constraints[node.dof_number(sys_num,var,i)].clear();
5388  _primal_constraint_values.erase(node.dof_number(sys_num,var,i));
5389  }
5390  }
5391  else
5392  {
5393  const unsigned int total_dofs = node.n_comp(sys_num, var);
5394  libmesh_assert_greater_equal (total_dofs, high_nc);
5395  // Add "this is zero" constraint rows for high p
5396  // non-vertex dofs, which are numbered in reverse
5397  for (unsigned int j = low_nc; j != high_nc; ++j)
5398  {
5399  const unsigned int i = total_dofs - j - 1;
5400  _dof_constraints[node.dof_number(sys_num,var,i)].clear();
5401  _primal_constraint_values.erase(node.dof_number(sys_num,var,i));
5402  }
5403  }
5404  }
5405 }
5406 
5407 #endif // LIBMESH_ENABLE_AMR
5408 
5409 
5410 #ifdef LIBMESH_ENABLE_DIRICHLET
5411 void DofMap::add_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary)
5412 {
5413  _dirichlet_boundaries->push_back(std::make_unique<DirichletBoundary>(dirichlet_boundary));
5414 }
5415 
5416 
5418  unsigned int qoi_index)
5419 {
5420  unsigned int old_size = cast_int<unsigned int>
5422  for (unsigned int i = old_size; i <= qoi_index; ++i)
5423  _adjoint_dirichlet_boundaries.push_back(std::make_unique<DirichletBoundaries>());
5424 
5425  // Make copy of DirichletBoundary, owned by _adjoint_dirichlet_boundaries
5426  _adjoint_dirichlet_boundaries[qoi_index]->push_back
5427  (std::make_unique<DirichletBoundary>(dirichlet_boundary));
5428 }
5429 
5430 
5432 {
5433  if (_adjoint_dirichlet_boundaries.size() > q)
5434  return true;
5435 
5436  return false;
5437 }
5438 
5439 
5440 const DirichletBoundaries *
5442 {
5443  libmesh_assert_greater(_adjoint_dirichlet_boundaries.size(),q);
5444  return _adjoint_dirichlet_boundaries[q].get();
5445 }
5446 
5447 
5450 {
5451  unsigned int old_size = cast_int<unsigned int>
5453  for (unsigned int i = old_size; i <= q; ++i)
5454  _adjoint_dirichlet_boundaries.push_back(std::make_unique<DirichletBoundaries>());
5455 
5456  return _adjoint_dirichlet_boundaries[q].get();
5457 }
5458 
5459 
5460 void DofMap::remove_dirichlet_boundary (const DirichletBoundary & boundary_to_remove)
5461 {
5462  // Find a boundary condition matching the one to be removed
5463  auto lam = [&boundary_to_remove](const auto & bdy)
5464  {return bdy->b == boundary_to_remove.b && bdy->variables == boundary_to_remove.variables;};
5465 
5466  auto it = std::find_if(_dirichlet_boundaries->begin(), _dirichlet_boundaries->end(), lam);
5467 
5468  // Assert it was actually found and remove it from the vector
5469  libmesh_assert (it != _dirichlet_boundaries->end());
5470  _dirichlet_boundaries->erase(it);
5471 }
5472 
5473 
5475  unsigned int qoi_index)
5476 {
5477  libmesh_assert_greater(_adjoint_dirichlet_boundaries.size(),
5478  qoi_index);
5479 
5480  auto lam = [&boundary_to_remove](const auto & bdy)
5481  {return bdy->b == boundary_to_remove.b && bdy->variables == boundary_to_remove.variables;};
5482 
5483  auto it = std::find_if(_adjoint_dirichlet_boundaries[qoi_index]->begin(),
5484  _adjoint_dirichlet_boundaries[qoi_index]->end(),
5485  lam);
5486 
5487  // Assert it was actually found and remove it from the vector
5488  libmesh_assert (it != _adjoint_dirichlet_boundaries[qoi_index]->end());
5489  _adjoint_dirichlet_boundaries[qoi_index]->erase(it);
5490 }
5491 
5492 
5494  const DirichletBoundary & boundary) const
5495 {
5496  const std::set<boundary_id_type>& mesh_side_bcids =
5498  const std::set<boundary_id_type>& mesh_edge_bcids =
5500  const std::set<boundary_id_type>& mesh_node_bcids =
5502  const std::set<boundary_id_type>& dbc_bcids = boundary.b;
5503 
5504  // DirichletBoundary id sets should be consistent across all ranks
5505  libmesh_assert(mesh.comm().verify(dbc_bcids.size()));
5506 
5507  for (const auto & bc_id : dbc_bcids)
5508  {
5509  // DirichletBoundary id sets should be consistent across all ranks
5510  libmesh_assert(mesh.comm().verify(bc_id));
5511 
5512  bool found_bcid = (mesh_side_bcids.find(bc_id) != mesh_side_bcids.end() ||
5513  mesh_edge_bcids.find(bc_id) != mesh_edge_bcids.end() ||
5514  mesh_node_bcids.find(bc_id) != mesh_node_bcids.end());
5515 
5516  // On a distributed mesh, boundary id sets may *not* be
5517  // consistent across all ranks, since not all ranks see all
5518  // boundaries
5519  mesh.comm().max(found_bcid);
5520 
5521  libmesh_error_msg_if(!found_bcid,
5522  "Could not find Dirichlet boundary id " << bc_id << " in mesh!");
5523  }
5524 }
5525 
5526 #endif // LIBMESH_ENABLE_DIRICHLET
5527 
5528 
5529 #ifdef LIBMESH_ENABLE_PERIODIC
5530 
5531 void DofMap::add_periodic_boundary (const PeriodicBoundaryBase & periodic_boundary)
5532 {
5533  // See if we already have a periodic boundary associated myboundary...
5534  PeriodicBoundaryBase * existing_boundary = _periodic_boundaries->boundary(periodic_boundary.myboundary);
5535 
5536  if (!existing_boundary)
5537  {
5538  // ...if not, clone the input (and its inverse) and add them to the PeriodicBoundaries object
5539  // Pass the pairedboundary of the original as the boundary id of the inverse clone.
5540  _periodic_boundaries->emplace(periodic_boundary.myboundary, periodic_boundary.clone());
5541  _periodic_boundaries->emplace(periodic_boundary.pairedboundary, periodic_boundary.clone(PeriodicBoundaryBase::INVERSE));
5542  }
5543  else
5544  {
5545  // ...otherwise, merge this object's variable IDs with the existing boundary object's.
5546  existing_boundary->merge(periodic_boundary);
5547 
5548  // Do the same merging process for the inverse boundary. Note: the inverse better already exist!
5549  PeriodicBoundaryBase * inverse_boundary = _periodic_boundaries->boundary(periodic_boundary.pairedboundary);
5550  libmesh_assert(inverse_boundary);
5551  inverse_boundary->merge(periodic_boundary);
5552  }
5553 }
5554 
5555 
5556 
5557 
5559  const PeriodicBoundaryBase & inverse_boundary)
5560 {
5561  libmesh_assert_equal_to (boundary.myboundary, inverse_boundary.pairedboundary);
5562  libmesh_assert_equal_to (boundary.pairedboundary, inverse_boundary.myboundary);
5563 
5564  // Store clones of the passed-in objects. These will be cleaned up
5565  // automatically in the _periodic_boundaries destructor.
5566  _periodic_boundaries->emplace(boundary.myboundary, boundary.clone());
5567  _periodic_boundaries->emplace(inverse_boundary.myboundary, inverse_boundary.clone());
5568 }
5569 
5570 
5571 #endif
5572 
5573 
5574 } // namespace libMesh
void remove_adjoint_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary, unsigned int q)
Removes from the system the specified Dirichlet boundary for the adjoint equation defined by Quantity...
std::unique_ptr< FEMFunctionBase< Gradient > > g_fem
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:196
std::string get_local_constraints(bool print_nonlocal=false) const
Gets a string reporting all DoF and Node constraints local to this processor.
FEFamily family
The type of finite element.
Definition: fe_type.h:221
void constrain_element_dyad_matrix(DenseVector< Number > &v, DenseVector< Number > &w, std::vector< dof_id_type > &row_dofs, bool asymmetric_constraint_rows=true) const
Constrains a dyadic element matrix B = v w&#39;.
Definition: dof_map.h:2343
bool is_constrained_node(const Node *node) const
Definition: dof_map.h:2242
void add_adjoint_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary, unsigned int q)
Adds a copy of the specified Dirichlet boundary to the system, corresponding to the adjoint problem d...
bool empty() const
Definition: stored_range.h:300
void check_dirichlet_bcid_consistency(const MeshBase &mesh, const DirichletBoundary &boundary) const
Check that all the ids in dirichlet_bcids are actually present in the mesh.
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:1032
bool is_prepared() const
Definition: mesh_base.h:198
virtual Output component(unsigned int i, const Point &p, Real time=0.)
constraint_rows_type & get_constraint_rows()
Constraint rows accessors.
Definition: mesh_base.h:1703
A Node is like a Point, but with more information.
Definition: node.h:52
std::set< GhostingFunctor * >::const_iterator coupling_functors_begin() const
Beginning of range of coupling functors.
Definition: dof_map.h:359
void build_constraint_matrix_and_vector(DenseMatrix< Number > &C, DenseVector< Number > &H, std::vector< dof_id_type > &elem_dofs, int qoi_index=-1, const bool called_recursively=false) const
Build the constraint matrix C and the forcing vector H associated with the element degree of freedom ...
void get_side_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge/face (2D/3D) finite element object for variable var for the largest dimension in th...
Definition: fem_context.h:317
virtual void zero() override final
Set every element in the vector to 0.
Definition: dense_vector.h:420
virtual void pre_fe_reinit(const System &, const Elem *e)
Reinitializes local data vectors/matrices on the current geometric element.
Definition: fem_context.C:1683
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:1002
std::unique_ptr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:1638
std::unique_ptr< FunctionBase< Number > > f
void scatter_constraints(MeshBase &)
Sends constraint equations to constraining processors.
void add_periodic_boundary(const PeriodicBoundaryBase &periodic_boundary)
Adds a copy of the specified periodic boundary to the system.
static void compute_node_constraints(NodeConstraints &constraints, const Elem *elem)
Computes the nodal constraint contributions (for non-conforming adapted meshes), using Lagrange geome...
Definition: fe_abstract.C:886
bool _error_on_constraint_loop
This flag indicates whether or not we do an opt-mode check for the presence of constraint loops...
Definition: dof_map.h:1934
void add_adjoint_constraint_row(const unsigned int qoi_index, const dof_id_type dof_number, const DofConstraintRow &constraint_row, const Number constraint_rhs, const bool forbid_constraint_overwrite)
Adds a copy of the user-defined row to the constraint matrix, using an inhomogeneous right-hand-side ...
void get_edge_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for edge (3D only!) finite element object for variable var.
Definition: fem_context.h:1332
void process_mesh_constraint_rows(const MeshBase &mesh)
Adds any spline constraints from the Mesh to our DoF constraints.
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:53
void constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix and vector.
Definition: dof_map.h:2330
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
static std::unique_ptr< LinearSolver< T > > build(const libMesh::Parallel::Communicator &comm_in, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a LinearSolver using the linear solver package specified by solver_package.
Definition: linear_solver.C:59
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:2164
virtual void zero() override final
Sets all elements of the matrix to 0 and resets any decomposition flag which may have been previously...
Definition: dense_matrix.h:911
void parallel_for(const Range &range, const Body &body)
Execute the provided function object in parallel on the specified range.
Definition: threads_none.h:73
dof_id_type n_dofs() const
Definition: dof_map_base.h:105
static constexpr Real TOLERANCE
unsigned int dim
Real jacobian_tolerance
Defaults to zero, but can be set to a custom small negative value to try and avoid spurious zero (or ...
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
virtual bool is_edge_on_side(const unsigned int e, const unsigned int s) const =0
MessageTag get_unique_tag(int tagvalue=MessageTag::invalid_tag) const
dof_id_type n_local_constrained_dofs() const
virtual numeric_index_type size() const =0
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
void remove_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Removes the specified Dirichlet boundary from the system.
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, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
Definition: fe_interface.C:611
std::vector< dof_id_type > _send_list
A list containing all the global DOF indices that affect the solution on my processor.
Definition: dof_map.h:1989
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:2220
std::unique_ptr< SparsityPattern::Build > build_sparsity(const MeshBase &mesh, bool calculate_constrained=false, bool use_condensed_system=false) const
Builds a sparsity pattern for matrices using the current degree-of-freedom numbering and coupling...
Definition: dof_map.C:61
void sum(T &r) const
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
MeshBase & mesh
void gather_constraints(MeshBase &mesh, std::set< dof_id_type > &unexpanded_dofs, bool look_for_constrainees)
Helper function for querying about constraint equations on other processors.
static FEFieldType field_type(const FEType &fe_type)
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
std::unique_ptr< DirichletBoundaries > _dirichlet_boundaries
Data structure containing Dirichlet functions.
Definition: dof_map.h:2134
NumericVector< Number > * rhs
The system matrix.
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
const Parallel::Communicator & comm() const
The Node constraint storage format.
Definition: dof_map.h:156
unsigned int p_level() const
Definition: elem.h:3108
std::vector< unsigned int > variables
unsigned int m() const
void set_algebraic_type(const AlgebraicType atype)
Setting which determines whether to initialize algebraic structures (elem_*) on each element and set ...
Definition: fem_context.h:986
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
std::map< const Elem *, const CouplingMatrix *, CompareDofObjectsByPIDAndThenID > map_type
What elements do we care about and what variables do we care about on each element?
This class defines a vector in LIBMESH_DIM dimensional Real or Complex space.
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:54
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
bool libmesh_isnan(T x)
const MeshBase & get_mesh() const
Definition: system.h:2358
void send_packed_range(const unsigned int dest_processor_id, const Context *context, Iter range_begin, const Iter range_end, const MessageTag &tag=no_tag, std::size_t approx_buffer_size=1000000) const
Real distance(const Point &p)
Number exact_value(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
static std::unique_ptr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package(), const MatrixBuildType matrix_build_type=MatrixBuildType::AUTOMATIC)
Builds a SparseMatrix<T> using the linear solver package specified by solver_package.
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota was created back when std::iota was just an SGI STL extension.
Definition: utility.h:229
unsigned int sys_number() const
Definition: dof_map.h:2172
uint8_t processor_id_type
Definition: id_types.h:104
This is the MeshBase class.
Definition: mesh_base.h:75
const Variable & variable(const unsigned int c) const override
Definition: dof_map.h:2190
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value)=0
Set the element (i,j) to value.
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
Heterogeneously constrains the numeric vector v, which represents an adjoint solution defined on the ...
Definition: dof_map.h:2354
AdjointDofConstraintValues _adjoint_constraint_values
Definition: dof_map.h:2110
void vector_mult_transpose(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this)^T * arg.
dof_id_type end_dof() const
Definition: dof_map_base.h:83
static void compute_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...
void merge(const PeriodicBoundaryBase &pb)
const std::set< boundary_id_type > & get_node_boundary_ids() const
void enforce_constraints_on_residual(const NonlinearImplicitSystem &system, NumericVector< Number > *rhs, NumericVector< Number > const *solution, bool homogeneous=true) const
Definition: dof_map.h:2359
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
std::unique_ptr< FEMFunctionBase< Number > > f_fem
virtual bool is_node_on_edge(const unsigned int n, const unsigned int e) const =0
const std::vector< std::vector< OutputGradient > > & get_dphi() const
Definition: fe_base.h:230
unsigned int first_scalar_number() const
Definition: variable.h:138
processor_id_type n_processors() const
virtual bool is_serial() const
Definition: mesh_base.h:211
void set_jacobian_tolerance(Real tol)
Set the Jacobian tolerance used for determining when the mapping fails.
Definition: fe_map.h:641
boundary_id_type myboundary
The boundary ID of this boundary and its counterpart.
const dof_id_type n_nodes
Definition: tecplot_io.C:67
dof_id_type first_dof() const
Definition: dof_map_base.h:73
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
This class defines the notion of a variable in the system.
Definition: variable.h:50
static void merge_ghost_functor_outputs(GhostingFunctor::map_type &elements_to_ghost, CouplingMatricesSet &temporary_coupling_matrices, const std::set< GhostingFunctor *>::iterator &gf_begin, const std::set< GhostingFunctor *>::iterator &gf_end, const MeshBase::const_element_iterator &elems_begin, const MeshBase::const_element_iterator &elems_end, processor_id_type p)
Definition: dof_map.C:1614
int8_t boundary_id_type
Definition: id_types.h:51
void vector_mult_add(DenseVector< T > &dest, const T factor, const DenseVector< T > &arg) const
Performs the scaled matrix-vector multiplication, dest += factor * (*this) * arg. ...
bool has_adjoint_dirichlet_boundaries(unsigned int q) const
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2529
dof_id_type id() const
Definition: dof_object.h:828
void min(const T &r, T &o, Request &req) const
void add_constraint_row(const dof_id_type dof_number, const DofConstraintRow &constraint_row, const Number constraint_rhs, const bool forbid_constraint_overwrite)
Adds a copy of the user-defined row to the constraint matrix, using an inhomogeneous right-hand-side ...
unsigned int n_vars
const System & get_system() const
Accessor for associated system.
Definition: diff_context.h:106
virtual unsigned int n_nodes() const =0
std::set< std::unique_ptr< CouplingMatrix >, Utility::CompareUnderlying > CouplingMatricesSet
Definition: dof_map.h:1848
bool is_sorted(const std::vector< KeyType > &v)
void heterogeneously_constrain_element_vector(const DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
Constrains the element vector.
Definition: dof_map.h:2340
void receive_packed_range(const unsigned int dest_processor_id, Context *context, OutputIter out, const T *output_type, const MessageTag &tag=any_tag) const
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:967
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:96
unsigned int n_variables() const override
Definition: dof_map.h:628
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
virtual void init_context(const FEMContext &)
Prepares a context object for use.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
void print_dof_constraints(std::ostream &os=libMesh::out, bool print_nonlocal=false) const
Prints (from processor 0) all DoF and Node constraints.
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libmesh_assert(ctx)
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:2258
Manages consistently variables, degrees of freedom, coefficient vectors, matrices and non-linear solv...
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
virtual_for_inffe const std::vector< Real > & get_JxW() const
Definition: fe_abstract.h:303
void libmesh_assert_valid_boundary_ids(const MeshBase &mesh)
A function for verifying that boundary condition ids match across processors.
Definition: mesh_tools.C:1708
void allgather_recursive_constraints(MeshBase &)
Gathers constraint equation dependencies from other processors.
This class provides all data required for a physics package (e.g.
Definition: fem_context.h:62
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr)=0
This is at the core of this class.
const std::vector< dof_id_type > & get_dof_indices() const
Accessor for element dof indices.
Definition: diff_context.h:363
DofConstraints _dof_constraints
Data structure containing DOF constraints.
Definition: dof_map.h:2106
unsigned int n_points() const
Definition: quadrature.h:131
virtual_for_inffe const std::vector< Point > & get_xyz() const
Definition: fe_abstract.h:280
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:167
virtual unsigned int n_edges() const =0
virtual std::unique_ptr< PeriodicBoundaryBase > clone(TransformationType t=FORWARD) const =0
If we want the DofMap to be able to make copies of references and store them in the underlying map...
static unsigned int n_vec_dim(const MeshBase &mesh, const FEType &fe_type)
virtual void close()=0
Calls the NumericVector&#39;s internal assembly routines, ensuring that the values are consistent across ...
std::pair< Real, Real > max_constraint_error(const System &system, NumericVector< Number > *v=nullptr) const
Tests the constrained degrees of freedom on the numeric vector v, which represents a solution defined...
A do-nothing class for templated methods that expect output iterator arguments.
Storage for DofConstraint right hand sides for a particular problem.
Definition: dof_map.h:120
std::string enum_to_string(const T e)
virtual bool closed() const
static FEContinuity get_continuity(const FEType &fe_type)
Returns the input FEType&#39;s FEContinuity based on the underlying FEFamily and potentially the Order...
const std::set< boundary_id_type > & get_boundary_ids() const
std::vector< std::unique_ptr< DirichletBoundaries > > _adjoint_dirichlet_boundaries
Data structure containing Dirichlet functions.
Definition: dof_map.h:2140
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
virtual const Elem * elem_ptr(const dof_id_type i) const =0
virtual numeric_index_type first_local_index() const =0
virtual unsigned int n_sides() const =0
virtual void edge_reinit(const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *pts=nullptr, const std::vector< Real > *weights=nullptr)=0
Reinitializes all the physical element-dependent data based on the edge of the element elem...
ParallelType type() const
std::unique_ptr< PeriodicBoundaries > _periodic_boundaries
Data structure containing periodic boundaries.
Definition: dof_map.h:2126
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void left_multiply_transpose(const DenseMatrix< T > &A)
Left multiplies by the transpose of the matrix A.
void process_constraints(MeshBase &)
Postprocesses any constrained degrees of freedom to be constrained only in terms of unconstrained dof...
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
std::unique_ptr< FunctionBase< Gradient > > g
timpi_pure bool verify(const T &r) const
dof_id_type n_constrained_dofs() const
void constrain_element_vector(DenseVector< Number > &rhs, std::vector< dof_id_type > &dofs, bool asymmetric_constraint_rows=true) const
Constrains the element vector.
Definition: dof_map.h:2326
void max(const T &r, T &o, Request &req) const
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
System * system() const
Definition: variable.h:114
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:2507
const DirichletBoundaries * get_adjoint_dirichlet_boundaries(unsigned int q) const
virtual bool is_vertex(const unsigned int i) const =0
bool _verify_dirichlet_bc_consistency
Flag which determines whether we should do some additional checking of the consistency of the Dirichl...
Definition: dof_map.h:2162
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
void add_constraints_to_send_list()
Adds entries to the _send_list vector corresponding to DoFs which are dependencies for constraint equ...
virtual void right_multiply(const DenseMatrixBase< T > &M2) override final
Performs the operation: (*this) <- (*this) * M3.
SparseMatrix< Number > * matrix
The system matrix.
DofConstraintValueMap _primal_constraint_values
Definition: dof_map.h:2108
virtual numeric_index_type local_size() const =0
void heterogeneously_constrain_element_matrix_and_vector(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true, int qoi_index=-1) const
Constrains the element matrix and vector.
Definition: dof_map.h:2336
const std::set< boundary_id_type > & get_edge_boundary_ids() const
virtual Output component(const FEMContext &, unsigned int i, const Point &p, Real time=0.)
void resize(const unsigned int new_m, const unsigned int new_n)
Resizes the matrix to the specified size and calls zero().
Definition: dense_matrix.h:895
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:100
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
dof_id_type n_local_dofs() const
Definition: dof_map_base.h:115
virtual unsigned int size() const override final
Definition: dense_vector.h:104
std::unique_ptr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1605
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
void build_constraint_matrix(DenseMatrix< Number > &C, std::vector< dof_id_type > &elem_dofs, const bool called_recursively=false) const
Build the constraint matrix C associated with the element degree of freedom indices elem_dofs...
void get_element_fe(unsigned int var, FEGenericBase< OutputShape > *&fe) const
Accessor for interior finite element object for variable var for the largest dimension in the mesh...
Definition: fem_context.h:277
void heterogeneously_constrain_element_residual(DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, NumericVector< Number > &solution_local) const
Constrains the element residual.
The base class for defining periodic boundaries.
void create_dof_constraints(const MeshBase &, Real time=0)
Rebuilds the raw degree of freedom and DofObject constraints, based on attached DirichletBoundary obj...
std::set< GhostingFunctor * >::const_iterator coupling_functors_end() const
End of range of coupling functors.
Definition: dof_map.h:365
std::set< boundary_id_type > b
std::map< const Node *, Real, std::less< const Node * >, Threads::scalable_allocator< std::pair< const Node *const, Real > > > NodeConstraintRow
A row of the Node constraint mapping.
Definition: dof_map.h:148
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
void edge_boundary_ids(const Elem *const elem, const unsigned short int edge, std::vector< boundary_id_type > &vec_to_fill) const
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
const QBase & get_edge_qrule() const
Accessor for element edge quadrature rule.
Definition: fem_context.h:831
void heterogeneously_constrain_element_jacobian_and_residual(DenseMatrix< Number > &matrix, DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, NumericVector< Number > &solution_local) const
Constrains the element Jacobian and residual.
unsigned int number() const
Definition: variable.h:128
const FEMap & get_fe_map() const
Definition: fe_abstract.h:554
FEFamily
defines an enum for finite element families.
unsigned int n() const
virtual dof_id_type n_elem() const =0
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
void constrain_element_residual(DenseVector< Number > &rhs, std::vector< dof_id_type > &elem_dofs, NumericVector< Number > &solution_local) const
Constrains the element residual.
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 periodic boundary conditions) corresponding to vari...
bool active() const
Definition: elem.h:2941
const DofMap & get_dof_map() const
Definition: system.h:2374
processor_id_type processor_id() const
Definition: dof_object.h:905
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, const bool add_p_level=true)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
Definition: fe_interface.C:597
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:2475
const Point & point(const unsigned int i) const
Definition: elem.h:2453
virtual numeric_index_type last_local_index() const =0
The constraint matrix storage format.
Definition: dof_map.h:108
const std::vector< dof_id_type > & get_send_list() const
Definition: dof_map.h:526
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
static void compute_periodic_node_constraints(NodeConstraints &constraints, const PeriodicBoundaries &boundaries, const MeshBase &mesh, const PointLocatorBase *point_locator, const Elem *elem)
Computes the node position constraint equation contributions (for meshes with periodic boundary condi...
Definition: fe_abstract.C:1078
This class forms the foundation from which generic finite elements may be derived.
void constrain_nothing(std::vector< dof_id_type > &dofs) const
Does not actually constrain anything, but modifies dofs in the same way as any of the constrain funct...
Definition: dof_map.h:2348
const QBase & get_element_qrule() const
Accessor for element interior quadrature rule for the dimension of the current _elem.
Definition: fem_context.h:802
This class provides single index access to FieldType (i.e.
Definition: raw_accessor.h:93
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this) * arg.
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...
void check_for_cyclic_constraints()
Throw an error if we detect any constraint loops, i.e.
const QBase & get_side_qrule() const
Accessor for element side quadrature rule for the dimension of the current _elem. ...
Definition: fem_context.h:809
uint8_t dof_id_type
Definition: id_types.h:67
std::vector< dof_id_type > _end_df
Last DOF index (plus 1) on processor p.
Definition: dof_map_base.h:159
bool local_index(dof_id_type dof_index) const
Definition: dof_map.h:839
const std::vector< std::vector< OutputShape > > & get_phi() const
Definition: fe_base.h:207
const FEType & type() const
Definition: variable.h:144
NodeConstraints _node_constraints
Data structure containing DofObject constraints.
Definition: dof_map.h:2117
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=nullptr, bool homogeneous=false) const
Constrains the numeric vector v, which represents a solution defined on the mesh. ...
Definition: dof_map.h:2350
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:30
void constrain_element_matrix(DenseMatrix< Number > &matrix, std::vector< dof_id_type > &elem_dofs, bool asymmetric_constraint_rows=true) const
Constrains the element matrix.
Definition: dof_map.h:2317
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
void enforce_constraints_on_jacobian(const NonlinearImplicitSystem &system, SparseMatrix< Number > *jac) const
Definition: dof_map.h:2365