libMesh
dof_map_constraints.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // Local Includes
19 #include "libmesh/boundary_info.h" // needed for dirichlet constraints
20 #include "libmesh/dense_matrix.h"
21 #include "libmesh/dense_vector.h"
22 #include "libmesh/dirichlet_boundaries.h"
23 #include "libmesh/dof_map.h"
24 #include "libmesh/elem.h"
25 #include "libmesh/elem_range.h"
26 #include "libmesh/fe_base.h"
27 #include "libmesh/fe_interface.h"
28 #include "libmesh/fe_type.h"
29 #include "libmesh/function_base.h"
30 #include "libmesh/libmesh_logging.h"
31 #include "libmesh/mesh_base.h"
32 #include "libmesh/mesh_inserter_iterator.h"
33 #include "libmesh/mesh_tools.h" // for libmesh_assert_valid_boundary_ids()
34 #include "libmesh/numeric_vector.h" // for enforce_constraints_exactly()
35 #include "libmesh/parallel.h"
36 #include "libmesh/parallel_algebra.h"
37 #include "libmesh/parallel_elem.h"
38 #include "libmesh/parallel_node.h"
39 #include "libmesh/periodic_boundaries.h"
40 #include "libmesh/periodic_boundary.h"
41 #include "libmesh/periodic_boundary_base.h"
42 #include "libmesh/point_locator_base.h"
43 #include "libmesh/quadrature.h" // for dirichlet constraints
44 #include "libmesh/raw_accessor.h"
45 #include "libmesh/sparse_matrix.h" // needed to constrain adjoint rhs
46 #include "libmesh/system.h" // needed by enforce_constraints_exactly()
47 #include "libmesh/threads.h"
48 #include "libmesh/tensor_tools.h"
49 #include "libmesh/utility.h" // Utility::iota()
50 
51 // C++ Includes
52 #include <set>
53 #include <algorithm> // for std::count, std::fill
54 #include <sstream>
55 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
56 #include <cmath>
57 
58 
59 // Anonymous namespace to hold helper classes
60 namespace {
61 
62 using namespace libMesh;
63 
64 class ComputeConstraints
65 {
66 public:
67  ComputeConstraints (DofConstraints & constraints,
68  DofMap & dof_map,
69 #ifdef LIBMESH_ENABLE_PERIODIC
70  PeriodicBoundaries & periodic_boundaries,
71 #endif
72  const MeshBase & mesh,
73  const unsigned int variable_number) :
74  _constraints(constraints),
75  _dof_map(dof_map),
76 #ifdef LIBMESH_ENABLE_PERIODIC
77  _periodic_boundaries(periodic_boundaries),
78 #endif
79  _mesh(mesh),
80  _variable_number(variable_number)
81  {}
82 
83  void operator()(const ConstElemRange & range) const
84  {
85  const Variable & var_description = _dof_map.variable(_variable_number);
86 
87 #ifdef LIBMESH_ENABLE_PERIODIC
88  UniquePtr<PointLocatorBase> point_locator;
89  const bool have_periodic_boundaries =
90  !_periodic_boundaries.empty();
91  if (have_periodic_boundaries && !range.empty())
92  point_locator = _mesh.sub_point_locator();
93 #endif
94 
95  for (ConstElemRange::const_iterator it = range.begin(); it!=range.end(); ++it)
96  if (var_description.active_on_subdomain((*it)->subdomain_id()))
97  {
98 #ifdef LIBMESH_ENABLE_AMR
100  _dof_map,
101  _variable_number,
102  *it);
103 #endif
104 #ifdef LIBMESH_ENABLE_PERIODIC
105  // FIXME: periodic constraints won't work on a non-serial
106  // mesh unless it's kept ghost elements from opposing
107  // boundaries!
108  if (have_periodic_boundaries)
110  _dof_map,
111  _periodic_boundaries,
112  _mesh,
113  point_locator.get(),
114  _variable_number,
115  *it);
116 #endif
117  }
118  }
119 
120 private:
121  DofConstraints & _constraints;
122  DofMap & _dof_map;
123 #ifdef LIBMESH_ENABLE_PERIODIC
124  PeriodicBoundaries & _periodic_boundaries;
125 #endif
126  const MeshBase & _mesh;
127  const unsigned int _variable_number;
128 };
129 
130 
131 
132 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
133 class ComputeNodeConstraints
134 {
135 public:
136  ComputeNodeConstraints (NodeConstraints & node_constraints,
137  DofMap & dof_map,
138 #ifdef LIBMESH_ENABLE_PERIODIC
139  PeriodicBoundaries & periodic_boundaries,
140 #endif
141  const MeshBase & mesh) :
142  _node_constraints(node_constraints),
143  _dof_map(dof_map),
144 #ifdef LIBMESH_ENABLE_PERIODIC
145  _periodic_boundaries(periodic_boundaries),
146 #endif
147  _mesh(mesh)
148  {
149  // Clang detects that _dof_map is not used for anything (other
150  // than being initialized) so let's prevent that warning.
151  libmesh_ignore(_dof_map);
152  }
153 
154  void operator()(const ConstElemRange & range) const
155  {
156 #ifdef LIBMESH_ENABLE_PERIODIC
157  UniquePtr<PointLocatorBase> point_locator;
158  bool have_periodic_boundaries = !_periodic_boundaries.empty();
159  if (have_periodic_boundaries && !range.empty())
160  point_locator = _mesh.sub_point_locator();
161 #endif
162 
163  for (ConstElemRange::const_iterator it = range.begin(); it!=range.end(); ++it)
164  {
165 #ifdef LIBMESH_ENABLE_AMR
166  FEBase::compute_node_constraints (_node_constraints, *it);
167 #endif
168 #ifdef LIBMESH_ENABLE_PERIODIC
169  // FIXME: periodic constraints won't work on a non-serial
170  // mesh unless it's kept ghost elements from opposing
171  // boundaries!
172  if (have_periodic_boundaries)
174  _periodic_boundaries,
175  _mesh,
176  point_locator.get(),
177  *it);
178 #endif
179  }
180  }
181 
182 private:
183  NodeConstraints & _node_constraints;
184  DofMap & _dof_map;
185 #ifdef LIBMESH_ENABLE_PERIODIC
186  PeriodicBoundaries & _periodic_boundaries;
187 #endif
188  const MeshBase & _mesh;
189 };
190 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
191 
192 
193 #ifdef LIBMESH_ENABLE_DIRICHLET
194 
198 class AddConstraint
199 {
200 protected:
201  DofMap & dof_map;
202 
203 public:
204  AddConstraint(DofMap & dof_map_in) : dof_map(dof_map_in) {}
205 
206  virtual void operator()(dof_id_type dof_number,
207  const DofConstraintRow & constraint_row,
208  const Number constraint_rhs) const = 0;
209 };
210 
211 class AddPrimalConstraint : public AddConstraint
212 {
213 public:
214  AddPrimalConstraint(DofMap & dof_map_in) : AddConstraint(dof_map_in) {}
215 
216  virtual void operator()(dof_id_type dof_number,
217  const DofConstraintRow & constraint_row,
218  const Number constraint_rhs) const
219  {
220  if (!dof_map.is_constrained_dof(dof_number))
221  dof_map.add_constraint_row (dof_number, constraint_row,
222  constraint_rhs, true);
223  }
224 };
225 
226 class AddAdjointConstraint : public AddConstraint
227 {
228 private:
229  const unsigned int qoi_index;
230 
231 public:
232  AddAdjointConstraint(DofMap & dof_map_in, unsigned int qoi_index_in)
233  : AddConstraint(dof_map_in), qoi_index(qoi_index_in) {}
234 
235  virtual void operator()(dof_id_type dof_number,
236  const DofConstraintRow & constraint_row,
237  const Number constraint_rhs) const
238  {
239  dof_map.add_adjoint_constraint_row
240  (qoi_index, dof_number, constraint_row, constraint_rhs,
241  true);
242  }
243 };
244 
245 
246 
252 class ConstrainDirichlet
253 {
254 private:
255  DofMap & dof_map;
256  const MeshBase & mesh;
257  const Real time;
258  const DirichletBoundary dirichlet;
259 
260  const AddConstraint & add_fn;
261 
262  static Number f_component (FunctionBase<Number> * f,
263  FEMFunctionBase<Number> * f_fem,
264  const FEMContext * c,
265  unsigned int i,
266  const Point & p,
267  Real time)
268  {
269  if (f_fem)
270  {
271  if (c)
272  return f_fem->component(*c, i, p, time);
273  else
274  return std::numeric_limits<Real>::quiet_NaN();
275  }
276  return f->component(i, p, time);
277  }
278 
279  static Gradient g_component (FunctionBase<Gradient> * g,
281  const FEMContext * c,
282  unsigned int i,
283  const Point & p,
284  Real time)
285  {
286  if (g_fem)
287  {
288  if (c)
289  return g_fem->component(*c, i, p, time);
290  else
291  return std::numeric_limits<Number>::quiet_NaN();
292  }
293  return g->component(i, p, time);
294  }
295 
296  template<typename OutputType>
297  void apply_dirichlet_impl(const ConstElemRange & range,
298  const unsigned int var,
299  const Variable & variable,
300  const FEType & fe_type) const
301  {
302  typedef OutputType OutputShape;
303  typedef typename TensorTools::IncrementRank<OutputShape>::type OutputGradient;
304  //typedef typename TensorTools::IncrementRank<OutputGradient>::type OutputTensor;
305  typedef typename TensorTools::MakeNumber<OutputShape>::type OutputNumber;
306  typedef typename TensorTools::IncrementRank<OutputNumber>::type OutputNumberGradient;
307  //typedef typename TensorTools::IncrementRank<OutputNumberGradient>::type OutputNumberTensor;
308 
309  FunctionBase<Number> * f = dirichlet.f.get();
310  FunctionBase<Gradient> * g = dirichlet.g.get();
311 
312  FEMFunctionBase<Number> * f_fem = dirichlet.f_fem.get();
313  FEMFunctionBase<Gradient> * g_fem = dirichlet.g_fem.get();
314 
315  const System * f_system = dirichlet.f_system;
316 
317  const std::set<boundary_id_type> & b = dirichlet.b;
318 
319  // We need data to project
320  libmesh_assert(f || f_fem);
321  libmesh_assert(!(f && f_fem));
322 
323  // Iff our data depends on a system, we should have it.
324  libmesh_assert(!(f && f_system));
325  libmesh_assert(!(f_fem && !f_system));
326 
327  // The element matrix and RHS for projections.
328  // Note that Ke is always real-valued, whereas
329  // Fe may be complex valued if complex number
330  // support is enabled
333  // The new element coefficients
335 
336  // The dimensionality of the current mesh
337  const unsigned int dim = mesh.mesh_dimension();
338 
339  // Boundary info for the current mesh
340  const BoundaryInfo & boundary_info = mesh.get_boundary_info();
341 
342  unsigned int n_vec_dim = FEInterface::n_vec_dim(mesh, fe_type);
343 
344  const unsigned int var_component =
345  variable.first_scalar_number();
346 
347  // Get FE objects of the appropriate type
349 
350  // Prepare variables for projection
351  UniquePtr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
352  UniquePtr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));
353  UniquePtr<QBase> qrule (fe_type.default_quadrature_rule(dim));
354 
355  // The values of the shape functions at the quadrature
356  // points
357  const std::vector<std::vector<OutputShape>> & phi = fe->get_phi();
358 
359  // The gradients of the shape functions at the quadrature
360  // points on the child element.
361  const std::vector<std::vector<OutputGradient>> * dphi = libmesh_nullptr;
362 
363  const FEContinuity cont = fe->get_continuity();
364 
365  if ((cont == C_ONE) && (fe_type.family != SUBDIVISION))
366  {
367  // We'll need gradient data for a C1 projection
368  libmesh_assert(g || g_fem);
369 
370  // We currently demand that either neither nor both function
371  // object depend on current FEM data.
372  libmesh_assert(!(g && g_fem));
373  libmesh_assert(!(f && g_fem));
374  libmesh_assert(!(f_fem && g));
375 
376  const std::vector<std::vector<OutputGradient>> & ref_dphi = fe->get_dphi();
377  dphi = &ref_dphi;
378  }
379 
380  // The Jacobian * quadrature weight at the quadrature points
381  const std::vector<Real> & JxW = fe->get_JxW();
382 
383  // The XYZ locations of the quadrature points
384  const std::vector<Point> & xyz_values = fe->get_xyz();
385 
386  // The global DOF indices
387  std::vector<dof_id_type> dof_indices;
388  // Side/edge local DOF indices
389  std::vector<unsigned int> side_dofs;
390 
391  // If our supplied functions require a FEMContext, and if we have
392  // an initialized solution to use with that FEMContext, then
393  // create one
394  UniquePtr<FEMContext> context;
395  if (f_fem)
396  {
397  libmesh_assert(f_system);
398  if (f_system->current_local_solution->initialized())
399  {
400  context = UniquePtr<FEMContext>(new FEMContext(*f_system));
401  f_fem->init_context(*context);
402  if (g_fem)
403  g_fem->init_context(*context);
404  }
405  }
406 
407  // Iterate over all the elements in the range
408  for (ConstElemRange::const_iterator elem_it=range.begin(); elem_it != range.end(); ++elem_it)
409  {
410  const Elem * elem = *elem_it;
411 
412  // We only calculate Dirichlet constraints on active
413  // elements
414  if (!elem->active())
415  continue;
416 
417  // Per-subdomain variables don't need to be projected on
418  // elements where they're not active
419  if (!variable.active_on_subdomain(elem->subdomain_id()))
420  continue;
421 
422  // There's a chicken-and-egg problem with FEMFunction-based
423  // Dirichlet constraints: we can't evaluate the FEMFunction
424  // until we have an initialized local solution vector, we
425  // can't initialize local solution vectors until we have a
426  // send list, and we can't generate a send list until we know
427  // all our constraints
428  //
429  // We don't generate constraints on uninitialized systems;
430  // currently user code will have to reinit() before any
431  // FEMFunction-based constraints will be correct. This should
432  // be fine, since user code would want to reinit() after
433  // setting initial conditions anyway.
434  if (f_system && context.get())
435  context->pre_fe_reinit(*f_system, elem);
436 
437  const unsigned short n_sides = elem->n_sides();
438  const unsigned short n_edges = elem->n_edges();
439  const unsigned short n_nodes = elem->n_nodes();
440 
441  // Find out which nodes, edges, sides and shellfaces are on a requested
442  // boundary:
443  std::vector<bool> is_boundary_node(n_nodes, false),
444  is_boundary_edge(n_edges, false),
445  is_boundary_side(n_sides, false),
446  is_boundary_shellface(2, false);
447 
448  // We also maintain a separate list of nodeset-based boundary nodes
449  std::vector<bool> is_boundary_nodeset(n_nodes, false);
450 
451  // Container to catch boundary ids handed back for sides,
452  // nodes, and edges in the loops below.
453  std::vector<boundary_id_type> ids_vec;
454 
455  for (unsigned char s = 0; s != n_sides; ++s)
456  {
457  // First see if this side has been requested
458  boundary_info.boundary_ids (elem, s, ids_vec);
459 
460  bool do_this_side = false;
461  for (std::vector<boundary_id_type>::iterator bc_it = ids_vec.begin();
462  bc_it != ids_vec.end(); ++bc_it)
463  if (b.count(*bc_it))
464  {
465  do_this_side = true;
466  break;
467  }
468  if (!do_this_side)
469  continue;
470 
471  is_boundary_side[s] = true;
472 
473  // Then see what nodes and what edges are on it
474  for (unsigned int n = 0; n != n_nodes; ++n)
475  if (elem->is_node_on_side(n,s))
476  is_boundary_node[n] = true;
477  for (unsigned int e = 0; e != n_edges; ++e)
478  if (elem->is_edge_on_side(e,s))
479  is_boundary_edge[e] = true;
480  }
481 
482  // We can also impose Dirichlet boundary conditions on nodes, so we should
483  // also independently check whether the nodes have been requested
484  for (unsigned int n=0; n != n_nodes; ++n)
485  {
486  boundary_info.boundary_ids (elem->node_ptr(n), ids_vec);
487 
488  for (std::vector<boundary_id_type>::iterator bc_it = ids_vec.begin();
489  bc_it != ids_vec.end(); ++bc_it)
490  if (b.count(*bc_it))
491  {
492  is_boundary_node[n] = true;
493  is_boundary_nodeset[n] = true;
494  }
495  }
496 
497  // We can also impose Dirichlet boundary conditions on edges, so we should
498  // also independently check whether the edges have been requested
499  for (unsigned short e=0; e != n_edges; ++e)
500  {
501  boundary_info.edge_boundary_ids (elem, e, ids_vec);
502 
503  for (std::vector<boundary_id_type>::iterator bc_it = ids_vec.begin();
504  bc_it != ids_vec.end(); ++bc_it)
505  if (b.count(*bc_it))
506  is_boundary_edge[e] = true;
507  }
508 
509  // We can also impose Dirichlet boundary conditions on shellfaces, so we should
510  // also independently check whether the shellfaces have been requested
511  for (unsigned short shellface=0; shellface != 2; ++shellface)
512  {
513  boundary_info.shellface_boundary_ids (elem, shellface, ids_vec);
514 
515  for (std::vector<boundary_id_type>::iterator bc_it = ids_vec.begin();
516  bc_it != ids_vec.end(); ++bc_it)
517  if (b.count(*bc_it))
518  is_boundary_shellface[shellface] = true;
519  }
520 
521  // Update the DOF indices for this element based on
522  // the current mesh
523  dof_map.dof_indices (elem, dof_indices, var);
524 
525  // The number of DOFs on the element
526  const unsigned int n_dofs =
527  cast_int<unsigned int>(dof_indices.size());
528 
529  // Fixed vs. free DoFs on edge/face projections
530  std::vector<char> dof_is_fixed(n_dofs, false); // bools
531  std::vector<int> free_dof(n_dofs, 0);
532 
533  // The element type
534  const ElemType elem_type = elem->type();
535 
536  // Zero the interpolated values
537  Ue.resize (n_dofs); Ue.zero();
538 
539  // In general, we need a series of
540  // projections to ensure a unique and continuous
541  // solution. We start by interpolating boundary nodes, then
542  // hold those fixed and project boundary edges, then hold
543  // those fixed and project boundary faces,
544 
545  // Interpolate node values first. Note that we have a special
546  // case for nodes that have a boundary nodeset, since we do
547  // need to interpolate them directly, even if they're non-vertex
548  // nodes.
549  unsigned int current_dof = 0;
550  for (unsigned int n=0; n!= n_nodes; ++n)
551  {
552  // FIXME: this should go through the DofMap,
553  // not duplicate dof_indices code badly!
554  const unsigned int nc =
555  FEInterface::n_dofs_at_node (dim, fe_type, elem_type,
556  n);
557  if ((!elem->is_vertex(n) || !is_boundary_node[n]) &&
558  !is_boundary_nodeset[n])
559  {
560  current_dof += nc;
561  continue;
562  }
563  if (cont == DISCONTINUOUS)
564  {
565  libmesh_assert_equal_to (nc, 0);
566  }
567  // Assume that C_ZERO elements have a single nodal
568  // value shape function
569  else if ((cont == C_ZERO) || (fe_type.family == SUBDIVISION))
570  {
571  libmesh_assert_equal_to (nc, n_vec_dim);
572  for (unsigned int c = 0; c < n_vec_dim; c++)
573  {
574  Ue(current_dof+c) =
575  f_component(f, f_fem, context.get(), var_component+c,
576  elem->point(n), time);
577  dof_is_fixed[current_dof+c] = true;
578  }
579  current_dof += n_vec_dim;
580  }
581  // The hermite element vertex shape functions are weird
582  else if (fe_type.family == HERMITE)
583  {
584  Ue(current_dof) =
585  f_component(f, f_fem, context.get(), var_component,
586  elem->point(n), time);
587  dof_is_fixed[current_dof] = true;
588  current_dof++;
589  Gradient grad =
590  g_component(g, g_fem, context.get(), var_component,
591  elem->point(n), time);
592  // x derivative
593  Ue(current_dof) = grad(0);
594  dof_is_fixed[current_dof] = true;
595  current_dof++;
596  if (dim > 1)
597  {
598  // We'll finite difference mixed derivatives
599  Point nxminus = elem->point(n),
600  nxplus = elem->point(n);
601  nxminus(0) -= TOLERANCE;
602  nxplus(0) += TOLERANCE;
603  Gradient gxminus =
604  g_component(g, g_fem, context.get(), var_component,
605  nxminus, time);
606  Gradient gxplus =
607  g_component(g, g_fem, context.get(), var_component,
608  nxplus, time);
609  // y derivative
610  Ue(current_dof) = grad(1);
611  dof_is_fixed[current_dof] = true;
612  current_dof++;
613  // xy derivative
614  Ue(current_dof) = (gxplus(1) - gxminus(1))
615  / 2. / TOLERANCE;
616  dof_is_fixed[current_dof] = true;
617  current_dof++;
618 
619  if (dim > 2)
620  {
621  // z derivative
622  Ue(current_dof) = grad(2);
623  dof_is_fixed[current_dof] = true;
624  current_dof++;
625  // xz derivative
626  Ue(current_dof) = (gxplus(2) - gxminus(2))
627  / 2. / TOLERANCE;
628  dof_is_fixed[current_dof] = true;
629  current_dof++;
630  // We need new points for yz
631  Point nyminus = elem->point(n),
632  nyplus = elem->point(n);
633  nyminus(1) -= TOLERANCE;
634  nyplus(1) += TOLERANCE;
635  Gradient gyminus =
636  g_component(g, g_fem, context.get(), var_component,
637  nyminus, time);
638  Gradient gyplus =
639  g_component(g, g_fem, context.get(), var_component,
640  nyplus, time);
641  // xz derivative
642  Ue(current_dof) = (gyplus(2) - gyminus(2))
643  / 2. / TOLERANCE;
644  dof_is_fixed[current_dof] = true;
645  current_dof++;
646  // Getting a 2nd order xyz is more tedious
647  Point nxmym = elem->point(n),
648  nxmyp = elem->point(n),
649  nxpym = elem->point(n),
650  nxpyp = elem->point(n);
651  nxmym(0) -= TOLERANCE;
652  nxmym(1) -= TOLERANCE;
653  nxmyp(0) -= TOLERANCE;
654  nxmyp(1) += TOLERANCE;
655  nxpym(0) += TOLERANCE;
656  nxpym(1) -= TOLERANCE;
657  nxpyp(0) += TOLERANCE;
658  nxpyp(1) += TOLERANCE;
659  Gradient gxmym =
660  g_component(g, g_fem, context.get(), var_component,
661  nxmym, time);
662  Gradient gxmyp =
663  g_component(g, g_fem, context.get(), var_component,
664  nxmyp, time);
665  Gradient gxpym =
666  g_component(g, g_fem, context.get(), var_component,
667  nxpym, time);
668  Gradient gxpyp =
669  g_component(g, g_fem, context.get(), var_component,
670  nxpyp, time);
671  Number gxzplus = (gxpyp(2) - gxmyp(2))
672  / 2. / TOLERANCE;
673  Number gxzminus = (gxpym(2) - gxmym(2))
674  / 2. / TOLERANCE;
675  // xyz derivative
676  Ue(current_dof) = (gxzplus - gxzminus)
677  / 2. / TOLERANCE;
678  dof_is_fixed[current_dof] = true;
679  current_dof++;
680  }
681  }
682  }
683  // Assume that other C_ONE elements have a single nodal
684  // value shape function and nodal gradient component
685  // shape functions
686  else if (cont == C_ONE)
687  {
688  libmesh_assert_equal_to (nc, 1 + dim);
689  Ue(current_dof) =
690  f_component(f, f_fem, context.get(), var_component,
691  elem->point(n), time);
692  dof_is_fixed[current_dof] = true;
693  current_dof++;
694  Gradient grad =
695  g_component(g, g_fem, context.get(), var_component,
696  elem->point(n), time);
697  for (unsigned int i=0; i!= dim; ++i)
698  {
699  Ue(current_dof) = grad(i);
700  dof_is_fixed[current_dof] = true;
701  current_dof++;
702  }
703  }
704  else
705  libmesh_error_msg("Unknown continuity cont = " << cont);
706  }
707 
708  // In 3D, project any edge values next
709  if (dim > 2 && cont != DISCONTINUOUS)
710  for (unsigned int e=0; e != n_edges; ++e)
711  {
712  if (!is_boundary_edge[e])
713  continue;
714 
715  FEInterface::dofs_on_edge(elem, dim, fe_type, e,
716  side_dofs);
717 
718  // Some edge dofs are on nodes and already
719  // fixed, others are free to calculate
720  unsigned int free_dofs = 0;
721  for (std::size_t i=0; i != side_dofs.size(); ++i)
722  if (!dof_is_fixed[side_dofs[i]])
723  free_dof[free_dofs++] = i;
724 
725  // There may be nothing to project
726  if (!free_dofs)
727  continue;
728 
729  Ke.resize (free_dofs, free_dofs); Ke.zero();
730  Fe.resize (free_dofs); Fe.zero();
731  // The new edge coefficients
732  DenseVector<Number> Uedge(free_dofs);
733 
734  // Initialize FE data on the edge
735  fe->attach_quadrature_rule (qedgerule.get());
736  fe->edge_reinit (elem, e);
737  const unsigned int n_qp = qedgerule->n_points();
738 
739  // Loop over the quadrature points
740  for (unsigned int qp=0; qp<n_qp; qp++)
741  {
742  // solution at the quadrature point
743  OutputNumber fineval(0);
744  libMesh::RawAccessor<OutputNumber> f_accessor( fineval, dim );
745 
746  for (unsigned int c = 0; c < n_vec_dim; c++)
747  f_accessor(c) =
748  f_component(f, f_fem, context.get(), var_component+c,
749  xyz_values[qp], time);
750 
751  // solution grad at the quadrature point
752  OutputNumberGradient finegrad;
753  libMesh::RawAccessor<OutputNumberGradient> g_accessor( finegrad, dim );
754 
755  unsigned int g_rank;
756  switch( FEInterface::field_type( fe_type ) )
757  {
758  case TYPE_SCALAR:
759  {
760  g_rank = 1;
761  break;
762  }
763  case TYPE_VECTOR:
764  {
765  g_rank = 2;
766  break;
767  }
768  default:
769  libmesh_error_msg("Unknown field type!");
770  }
771 
772  if (cont == C_ONE)
773  for (unsigned int c = 0; c < n_vec_dim; c++)
774  for (unsigned int d = 0; d < g_rank; d++)
775  g_accessor(c + d*dim ) =
776  g_component(g, g_fem, context.get(), var_component,
777  xyz_values[qp], time)(c);
778 
779  // Form edge projection matrix
780  for (std::size_t sidei=0, freei=0; sidei != side_dofs.size(); ++sidei)
781  {
782  unsigned int i = side_dofs[sidei];
783  // fixed DoFs aren't test functions
784  if (dof_is_fixed[i])
785  continue;
786  for (std::size_t sidej=0, freej=0; sidej != side_dofs.size(); ++sidej)
787  {
788  unsigned int j = side_dofs[sidej];
789  if (dof_is_fixed[j])
790  Fe(freei) -= phi[i][qp] * phi[j][qp] *
791  JxW[qp] * Ue(j);
792  else
793  Ke(freei,freej) += phi[i][qp] *
794  phi[j][qp] * JxW[qp];
795  if (cont == C_ONE)
796  {
797  if (dof_is_fixed[j])
798  Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp]) ) *
799  JxW[qp] * Ue(j);
800  else
801  Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
802  * JxW[qp];
803  }
804  if (!dof_is_fixed[j])
805  freej++;
806  }
807  Fe(freei) += phi[i][qp] * fineval * JxW[qp];
808  if (cont == C_ONE)
809  Fe(freei) += (finegrad.contract( (*dphi)[i][qp]) ) *
810  JxW[qp];
811  freei++;
812  }
813  }
814 
815  Ke.cholesky_solve(Fe, Uedge);
816 
817  // Transfer new edge solutions to element
818  for (unsigned int i=0; i != free_dofs; ++i)
819  {
820  Number & ui = Ue(side_dofs[free_dof[i]]);
822  std::abs(ui - Uedge(i)) < TOLERANCE);
823  ui = Uedge(i);
824  dof_is_fixed[side_dofs[free_dof[i]]] = true;
825  }
826  }
827 
828  // Project any side values (edges in 2D, faces in 3D)
829  if (dim > 1 && cont != DISCONTINUOUS)
830  for (unsigned int s=0; s != n_sides; ++s)
831  {
832  if (!is_boundary_side[s])
833  continue;
834 
835  FEInterface::dofs_on_side(elem, dim, fe_type, s,
836  side_dofs);
837 
838  // Some side dofs are on nodes/edges and already
839  // fixed, others are free to calculate
840  unsigned int free_dofs = 0;
841  for (std::size_t i=0; i != side_dofs.size(); ++i)
842  if (!dof_is_fixed[side_dofs[i]])
843  free_dof[free_dofs++] = i;
844 
845  // There may be nothing to project
846  if (!free_dofs)
847  continue;
848 
849  Ke.resize (free_dofs, free_dofs); Ke.zero();
850  Fe.resize (free_dofs); Fe.zero();
851  // The new side coefficients
852  DenseVector<Number> Uside(free_dofs);
853 
854  // Initialize FE data on the side
855  fe->attach_quadrature_rule (qsiderule.get());
856  fe->reinit (elem, s);
857  const unsigned int n_qp = qsiderule->n_points();
858 
859  // Loop over the quadrature points
860  for (unsigned int qp=0; qp<n_qp; qp++)
861  {
862  // solution at the quadrature point
863  OutputNumber fineval(0);
864  libMesh::RawAccessor<OutputNumber> f_accessor( fineval, dim );
865 
866  for (unsigned int c = 0; c < n_vec_dim; c++)
867  f_accessor(c) =
868  f_component(f, f_fem, context.get(), var_component+c,
869  xyz_values[qp], time);
870 
871  // solution grad at the quadrature point
872  OutputNumberGradient finegrad;
873  libMesh::RawAccessor<OutputNumberGradient> g_accessor( finegrad, dim );
874 
875  unsigned int g_rank;
876  switch( FEInterface::field_type( fe_type ) )
877  {
878  case TYPE_SCALAR:
879  {
880  g_rank = 1;
881  break;
882  }
883  case TYPE_VECTOR:
884  {
885  g_rank = 2;
886  break;
887  }
888  default:
889  libmesh_error_msg("Unknown field type!");
890  }
891 
892  if (cont == C_ONE)
893  for (unsigned int c = 0; c < n_vec_dim; c++)
894  for (unsigned int d = 0; d < g_rank; d++)
895  g_accessor(c + d*dim ) =
896  g_component(g, g_fem, context.get(), var_component,
897  xyz_values[qp], time)(c);
898 
899  // Form side projection matrix
900  for (std::size_t sidei=0, freei=0; sidei != side_dofs.size(); ++sidei)
901  {
902  unsigned int i = side_dofs[sidei];
903  // fixed DoFs aren't test functions
904  if (dof_is_fixed[i])
905  continue;
906  for (std::size_t sidej=0, freej=0; sidej != side_dofs.size(); ++sidej)
907  {
908  unsigned int j = side_dofs[sidej];
909  if (dof_is_fixed[j])
910  Fe(freei) -= phi[i][qp] * phi[j][qp] *
911  JxW[qp] * Ue(j);
912  else
913  Ke(freei,freej) += phi[i][qp] *
914  phi[j][qp] * JxW[qp];
915  if (cont == C_ONE)
916  {
917  if (dof_is_fixed[j])
918  Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp])) *
919  JxW[qp] * Ue(j);
920  else
921  Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
922  * JxW[qp];
923  }
924  if (!dof_is_fixed[j])
925  freej++;
926  }
927  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
928  if (cont == C_ONE)
929  Fe(freei) += (finegrad.contract((*dphi)[i][qp])) *
930  JxW[qp];
931  freei++;
932  }
933  }
934 
935  Ke.cholesky_solve(Fe, Uside);
936 
937  // Transfer new side solutions to element
938  for (unsigned int i=0; i != free_dofs; ++i)
939  {
940  Number & ui = Ue(side_dofs[free_dof[i]]);
942  std::abs(ui - Uside(i)) < TOLERANCE);
943  ui = Uside(i);
944  dof_is_fixed[side_dofs[free_dof[i]]] = true;
945  }
946  }
947 
948  // Project any shellface values
949  if (dim == 2 && cont != DISCONTINUOUS)
950  for (unsigned int shellface=0; shellface != 2; ++shellface)
951  {
952  if (!is_boundary_shellface[shellface])
953  continue;
954 
955  // A shellface has the same dof indices as the element itself
956  std::vector<unsigned int> shellface_dofs(dof_indices.size());
957  Utility::iota(shellface_dofs.begin(), shellface_dofs.end(), 0);
958 
959  // Some shellface dofs are on nodes/edges and already
960  // fixed, others are free to calculate
961  unsigned int free_dofs = 0;
962  for (std::size_t i=0; i != shellface_dofs.size(); ++i)
963  if (!dof_is_fixed[shellface_dofs[i]])
964  free_dof[free_dofs++] = i;
965 
966  // There may be nothing to project
967  if (!free_dofs)
968  continue;
969 
970  Ke.resize (free_dofs, free_dofs); Ke.zero();
971  Fe.resize (free_dofs); Fe.zero();
972  // The new shellface coefficients
973  DenseVector<Number> Ushellface(free_dofs);
974 
975  // Initialize FE data on the element
976  fe->attach_quadrature_rule (qrule.get());
977  fe->reinit (elem);
978  const unsigned int n_qp = qrule->n_points();
979 
980  // Loop over the quadrature points
981  for (unsigned int qp=0; qp<n_qp; qp++)
982  {
983  // solution at the quadrature point
984  OutputNumber fineval(0);
985  libMesh::RawAccessor<OutputNumber> f_accessor( fineval, dim );
986 
987  for (unsigned int c = 0; c < n_vec_dim; c++)
988  f_accessor(c) =
989  f_component(f, f_fem, context.get(), var_component+c,
990  xyz_values[qp], time);
991 
992  // solution grad at the quadrature point
993  OutputNumberGradient finegrad;
994  libMesh::RawAccessor<OutputNumberGradient> g_accessor( finegrad, dim );
995 
996  unsigned int g_rank;
997  switch( FEInterface::field_type( fe_type ) )
998  {
999  case TYPE_SCALAR:
1000  {
1001  g_rank = 1;
1002  break;
1003  }
1004  case TYPE_VECTOR:
1005  {
1006  g_rank = 2;
1007  break;
1008  }
1009  default:
1010  libmesh_error_msg("Unknown field type!");
1011  }
1012 
1013  if (cont == C_ONE)
1014  for (unsigned int c = 0; c < n_vec_dim; c++)
1015  for (unsigned int d = 0; d < g_rank; d++)
1016  g_accessor(c + d*dim ) =
1017  g_component(g, g_fem, context.get(), var_component,
1018  xyz_values[qp], time)(c);
1019 
1020  // Form shellface projection matrix
1021  for (std::size_t shellfacei=0, freei=0;
1022  shellfacei != shellface_dofs.size(); ++shellfacei)
1023  {
1024  unsigned int i = shellface_dofs[shellfacei];
1025  // fixed DoFs aren't test functions
1026  if (dof_is_fixed[i])
1027  continue;
1028  for (std::size_t shellfacej=0, freej=0;
1029  shellfacej != shellface_dofs.size(); ++shellfacej)
1030  {
1031  unsigned int j = shellface_dofs[shellfacej];
1032  if (dof_is_fixed[j])
1033  Fe(freei) -= phi[i][qp] * phi[j][qp] *
1034  JxW[qp] * Ue(j);
1035  else
1036  Ke(freei,freej) += phi[i][qp] *
1037  phi[j][qp] * JxW[qp];
1038  if (cont == C_ONE)
1039  {
1040  if (dof_is_fixed[j])
1041  Fe(freei) -= ((*dphi)[i][qp].contract((*dphi)[j][qp])) *
1042  JxW[qp] * Ue(j);
1043  else
1044  Ke(freei,freej) += ((*dphi)[i][qp].contract((*dphi)[j][qp]))
1045  * JxW[qp];
1046  }
1047  if (!dof_is_fixed[j])
1048  freej++;
1049  }
1050  Fe(freei) += (fineval * phi[i][qp]) * JxW[qp];
1051  if (cont == C_ONE)
1052  Fe(freei) += (finegrad.contract((*dphi)[i][qp])) *
1053  JxW[qp];
1054  freei++;
1055  }
1056  }
1057 
1058  Ke.cholesky_solve(Fe, Ushellface);
1059 
1060  // Transfer new shellface solutions to element
1061  for (unsigned int i=0; i != free_dofs; ++i)
1062  {
1063  Number & ui = Ue(shellface_dofs[free_dof[i]]);
1065  std::abs(ui - Ushellface(i)) < TOLERANCE);
1066  ui = Ushellface(i);
1067  dof_is_fixed[shellface_dofs[free_dof[i]]] = true;
1068  }
1069  }
1070 
1071  // Lock the DofConstraints since it is shared among threads.
1072  {
1073  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
1074 
1075  for (unsigned int i = 0; i < n_dofs; i++)
1076  {
1077  DofConstraintRow empty_row;
1078  if (dof_is_fixed[i] && !libmesh_isnan(Ue(i)))
1079  add_fn (dof_indices[i], empty_row, Ue(i));
1080  }
1081  }
1082  }
1083 
1084  } // apply_dirichlet_impl
1085 
1086 public:
1087  ConstrainDirichlet (DofMap & dof_map_in,
1088  const MeshBase & mesh_in,
1089  const Real time_in,
1090  const DirichletBoundary & dirichlet_in,
1091  const AddConstraint & add_in) :
1092  dof_map(dof_map_in),
1093  mesh(mesh_in),
1094  time(time_in),
1095  dirichlet(dirichlet_in),
1096  add_fn(add_in) { }
1097 
1098  ConstrainDirichlet (const ConstrainDirichlet & in) :
1099  dof_map(in.dof_map),
1100  mesh(in.mesh),
1101  time(in.time),
1102  dirichlet(in.dirichlet),
1103  add_fn(in.add_fn) { }
1104 
1105  void operator()(const ConstElemRange & range) const
1106  {
1113  // Loop over all the variables we've been requested to project
1114  for (std::size_t v=0; v!=dirichlet.variables.size(); v++)
1115  {
1116  const unsigned int var = dirichlet.variables[v];
1117 
1118  const Variable & variable = dof_map.variable(var);
1119 
1120  const FEType & fe_type = variable.type();
1121 
1122  if (fe_type.family == SCALAR)
1123  continue;
1124 
1125  switch( FEInterface::field_type( fe_type ) )
1126  {
1127  case TYPE_SCALAR:
1128  {
1129  this->apply_dirichlet_impl<Real>( range, var, variable, fe_type );
1130  break;
1131  }
1132  case TYPE_VECTOR:
1133  {
1134  this->apply_dirichlet_impl<RealGradient>( range, var, variable, fe_type );
1135  break;
1136  }
1137  default:
1138  libmesh_error_msg("Unknown field type!");
1139  }
1140 
1141  }
1142  }
1143 
1144 }; // class ConstrainDirichlet
1145 
1146 
1147 #endif // LIBMESH_ENABLE_DIRICHLET
1148 
1149 
1150 } // anonymous namespace
1151 
1152 
1153 
1154 namespace libMesh
1155 {
1156 
1157 // ------------------------------------------------------------
1158 // DofMap member functions
1159 
1160 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1161 
1162 
1164 {
1165  parallel_object_only();
1166 
1167  dof_id_type nc_dofs = this->n_local_constrained_dofs();
1168  this->comm().sum(nc_dofs);
1169  return nc_dofs;
1170 }
1171 
1172 
1174 {
1175  const DofConstraints::const_iterator lower =
1176  _dof_constraints.lower_bound(this->first_dof()),
1177  upper =
1178  _dof_constraints.upper_bound(this->end_dof()-1);
1179 
1180  return cast_int<dof_id_type>(std::distance(lower, upper));
1181 }
1182 
1183 
1184 
1186 {
1187  parallel_object_only();
1188 
1189  LOG_SCOPE("create_dof_constraints()", "DofMap");
1190 
1191  libmesh_assert (mesh.is_prepared());
1192 
1193  // The user might have set boundary conditions after the mesh was
1194  // prepared; we should double-check that those boundary conditions
1195  // are still consistent.
1196 #ifdef DEBUG
1198 #endif
1199 
1200  // We might get constraint equations from AMR hanging nodes in 2D/3D
1201  // or from boundary conditions in any dimension
1202  const bool possible_local_constraints = false
1203  || !mesh.n_elem()
1204 #ifdef LIBMESH_ENABLE_AMR
1205  || mesh.mesh_dimension() > 1
1206 #endif
1207 #ifdef LIBMESH_ENABLE_PERIODIC
1208  || !_periodic_boundaries->empty()
1209 #endif
1210 #ifdef LIBMESH_ENABLE_DIRICHLET
1211  || !_dirichlet_boundaries->empty()
1212 #endif
1213  ;
1214 
1215  // Even if we don't have constraints, another processor might.
1216  bool possible_global_constraints = possible_local_constraints;
1217 #if defined(LIBMESH_ENABLE_PERIODIC) || defined(LIBMESH_ENABLE_DIRICHLET) || defined(LIBMESH_ENABLE_AMR)
1218  libmesh_assert(this->comm().verify(mesh.is_serial()));
1219 
1220  this->comm().max(possible_global_constraints);
1221 #endif
1222 
1223  if (!possible_global_constraints)
1224  {
1225  // Clear out any old constraints; maybe the user just deleted
1226  // their last remaining dirichlet/periodic/user constraint?
1227 #ifdef LIBMESH_ENABLE_CONSTRAINTS
1228  _dof_constraints.clear();
1229  _stashed_dof_constraints.clear();
1230  _primal_constraint_values.clear();
1231  _adjoint_constraint_values.clear();
1232 #endif
1233 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1234  _node_constraints.clear();
1235 #endif
1236 
1237  return;
1238  }
1239 
1240  // Here we build the hanging node constraints. This is done
1241  // by enforcing the condition u_a = u_b along hanging sides.
1242  // u_a = u_b is collocated at the nodes of side a, which gives
1243  // one row of the constraint matrix.
1244 
1245  // Processors only compute their local constraints
1246  ConstElemRange range (mesh.local_elements_begin(),
1247  mesh.local_elements_end());
1248 
1249  // Global computation fails if we're using a FEMFunctionBase BC on a
1250  // ReplicatedMesh in parallel
1251  // ConstElemRange range (mesh.elements_begin(),
1252  // mesh.elements_end());
1253 
1254  // compute_periodic_constraints requires a point_locator() from our
1255  // Mesh, but point_locator() construction is parallel and threaded.
1256  // Rather than nest threads within threads we'll make sure it's
1257  // preconstructed.
1258 #ifdef LIBMESH_ENABLE_PERIODIC
1259  bool need_point_locator = !_periodic_boundaries->empty() && !range.empty();
1260 
1261  this->comm().max(need_point_locator);
1262 
1263  if (need_point_locator)
1264  mesh.sub_point_locator();
1265 #endif
1266 
1267 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1268  // recalculate node constraints from scratch
1269  _node_constraints.clear();
1270 
1271  Threads::parallel_for (range,
1272  ComputeNodeConstraints (_node_constraints,
1273  *this,
1274 #ifdef LIBMESH_ENABLE_PERIODIC
1275  *_periodic_boundaries,
1276 #endif
1277  mesh));
1278 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1279 
1280 
1281  // recalculate dof constraints from scratch
1282  _dof_constraints.clear();
1283  _stashed_dof_constraints.clear();
1284  _primal_constraint_values.clear();
1285  _adjoint_constraint_values.clear();
1286 
1287  // Look at all the variables in the system. Reset the element
1288  // range at each iteration -- there is no need to reconstruct it.
1289  for (unsigned int variable_number=0; variable_number<this->n_variables();
1290  ++variable_number, range.reset())
1291  Threads::parallel_for (range,
1292  ComputeConstraints (_dof_constraints,
1293  *this,
1294 #ifdef LIBMESH_ENABLE_PERIODIC
1295  *_periodic_boundaries,
1296 #endif
1297  mesh,
1298  variable_number));
1299 
1300 #ifdef LIBMESH_ENABLE_DIRICHLET
1301  for (DirichletBoundaries::iterator
1302  i = _dirichlet_boundaries->begin();
1303  i != _dirichlet_boundaries->end(); ++i, range.reset())
1304  {
1305  // Sanity check that the boundary ids associated with the DirichletBoundary
1306  // objects are actually present in the mesh
1307  this->check_dirichlet_bcid_consistency(mesh,**i);
1308 
1310  (range,
1311  ConstrainDirichlet(*this, mesh, time, **i,
1312  AddPrimalConstraint(*this))
1313  );
1314  }
1315 
1316  for (std::size_t qoi_index = 0;
1317  qoi_index != _adjoint_dirichlet_boundaries.size();
1318  ++qoi_index)
1319  {
1320  for (DirichletBoundaries::iterator
1321  i = _adjoint_dirichlet_boundaries[qoi_index]->begin();
1322  i != _adjoint_dirichlet_boundaries[qoi_index]->end();
1323  ++i, range.reset())
1324  {
1325  // Sanity check that the boundary ids associated with the DirichletBoundary
1326  // objects are actually present in the mesh
1327  this->check_dirichlet_bcid_consistency(mesh,**i);
1328 
1330  (range,
1331  ConstrainDirichlet(*this, mesh, time, **i,
1332  AddAdjointConstraint(*this, qoi_index))
1333  );
1334  }
1335  }
1336 
1337 #endif // LIBMESH_ENABLE_DIRICHLET
1338 }
1339 
1340 
1341 
1343  const DofConstraintRow & constraint_row,
1344  const Number constraint_rhs,
1345  const bool forbid_constraint_overwrite)
1346 {
1347  // Optionally allow the user to overwrite constraints. Defaults to false.
1348  if (forbid_constraint_overwrite)
1349  if (this->is_constrained_dof(dof_number))
1350  libmesh_error_msg("ERROR: DOF " << dof_number << " was already constrained!");
1351 
1352  // We don't get insert_or_assign until C++17 so we make do.
1353  std::pair<DofConstraints::iterator, bool> it =
1354  _dof_constraints.insert(std::make_pair(dof_number, constraint_row));
1355  if (!it.second)
1356  it.first->second = constraint_row;
1357 
1358  std::pair<DofConstraintValueMap::iterator, bool> rhs_it =
1359  _primal_constraint_values.insert(std::make_pair(dof_number, constraint_rhs));
1360  if (!rhs_it.second)
1361  rhs_it.first->second = constraint_rhs;
1362 }
1363 
1364 
1365 void DofMap::add_adjoint_constraint_row (const unsigned int qoi_index,
1366  const dof_id_type dof_number,
1367  const DofConstraintRow & /*constraint_row*/,
1368  const Number constraint_rhs,
1369  const bool forbid_constraint_overwrite)
1370 {
1371  // Optionally allow the user to overwrite constraints. Defaults to false.
1372  if (forbid_constraint_overwrite)
1373  {
1374  if (!this->is_constrained_dof(dof_number))
1375  libmesh_error_msg("ERROR: DOF " << dof_number << " has no corresponding primal constraint!");
1376 #ifndef NDEBUG
1377  // No way to do this without a non-normalized tolerance?
1378  /*
1379  // If the user passed in more than just the rhs, let's check the
1380  // coefficients for consistency
1381  if (!constraint_row.empty())
1382  {
1383  DofConstraintRow row = _dof_constraints[dof_number];
1384  for (DofConstraintRow::const_iterator pos=row.begin();
1385  pos != row.end(); ++pos)
1386  {
1387  libmesh_assert(constraint_row.find(pos->first)->second
1388  == pos->second);
1389  }
1390  }
1391 
1392  if (_adjoint_constraint_values[qoi_index].find(dof_number) !=
1393  _adjoint_constraint_values[qoi_index].end())
1394  libmesh_assert_equal_to
1395  (_adjoint_constraint_values[qoi_index][dof_number],
1396  constraint_rhs);
1397  */
1398 #endif
1399  }
1400 
1401  // Creates the map of rhs values if it doesn't already exist; then
1402  // adds the current value to that map
1403 
1404  // We don't get insert_or_assign until C++17 so we make do.
1405  std::pair<DofConstraintValueMap::iterator, bool> rhs_it =
1406  _adjoint_constraint_values[qoi_index].insert(std::make_pair(dof_number,
1407  constraint_rhs));
1408  if (!rhs_it.second)
1409  rhs_it.first->second = constraint_rhs;
1410 }
1411 
1412 
1413 
1414 
1415 void DofMap::print_dof_constraints(std::ostream & os,
1416  bool print_nonlocal) const
1417 {
1418  parallel_object_only();
1419 
1420  std::string local_constraints =
1421  this->get_local_constraints(print_nonlocal);
1422 
1423  if (this->processor_id())
1424  {
1425  this->comm().send(0, local_constraints);
1426  }
1427  else
1428  {
1429  os << "Processor 0:\n";
1430  os << local_constraints;
1431 
1432  for (processor_id_type i=1; i<this->n_processors(); ++i)
1433  {
1434  this->comm().receive(i, local_constraints);
1435  os << "Processor " << i << ":\n";
1436  os << local_constraints;
1437  }
1438  }
1439 }
1440 
1441 std::string DofMap::get_local_constraints(bool print_nonlocal) const
1442 {
1443  std::ostringstream os;
1444 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
1445  if (print_nonlocal)
1446  os << "All ";
1447  else
1448  os << "Local ";
1449 
1450  os << "Node Constraints:"
1451  << std::endl;
1452 
1453  for (NodeConstraints::const_iterator it=_node_constraints.begin();
1454  it != _node_constraints.end(); ++it)
1455  {
1456  const Node * node = it->first;
1457 
1458  // Skip non-local nodes if requested
1459  if (!print_nonlocal &&
1460  node->processor_id() != this->processor_id())
1461  continue;
1462 
1463  const NodeConstraintRow & row = it->second.first;
1464  const Point & offset = it->second.second;
1465 
1466  os << "Constraints for Node id " << node->id()
1467  << ": \t";
1468 
1469  for (NodeConstraintRow::const_iterator pos=row.begin();
1470  pos != row.end(); ++pos)
1471  os << " (" << pos->first->id() << ","
1472  << pos->second << ")\t";
1473 
1474  os << "rhs: " << offset;
1475 
1476  os << std::endl;
1477  }
1478 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
1479 
1480  if (print_nonlocal)
1481  os << "All ";
1482  else
1483  os << "Local ";
1484 
1485  os << "DoF Constraints:"
1486  << std::endl;
1487 
1488  for (DofConstraints::const_iterator it=_dof_constraints.begin();
1489  it != _dof_constraints.end(); ++it)
1490  {
1491  const dof_id_type i = it->first;
1492 
1493  // Skip non-local dofs if requested
1494  if (!print_nonlocal &&
1495  ((i < this->first_dof()) ||
1496  (i >= this->end_dof())))
1497  continue;
1498 
1499  const DofConstraintRow & row = it->second;
1500  DofConstraintValueMap::const_iterator rhsit =
1501  _primal_constraint_values.find(i);
1502  const Number rhs = (rhsit == _primal_constraint_values.end()) ?
1503  0 : rhsit->second;
1504 
1505  os << "Constraints for DoF " << i
1506  << ": \t";
1507 
1508  for (DofConstraintRow::const_iterator pos=row.begin();
1509  pos != row.end(); ++pos)
1510  os << " (" << pos->first << ","
1511  << pos->second << ")\t";
1512 
1513  os << "rhs: " << rhs;
1514 
1515  os << std::endl;
1516  }
1517 
1518  for (std::size_t qoi_index = 0;
1519  qoi_index != _adjoint_dirichlet_boundaries.size();
1520  ++qoi_index)
1521  {
1522  os << "Adjoint " << qoi_index << " DoF rhs values:"
1523  << std::endl;
1524 
1525  AdjointDofConstraintValues::const_iterator adjoint_map_it =
1526  _adjoint_constraint_values.find(qoi_index);
1527 
1528  if (adjoint_map_it != _adjoint_constraint_values.end())
1529  for (DofConstraintValueMap::const_iterator
1530  it=adjoint_map_it->second.begin();
1531  it != adjoint_map_it->second.end(); ++it)
1532  {
1533  const dof_id_type i = it->first;
1534 
1535  // Skip non-local dofs if requested
1536  if (!print_nonlocal &&
1537  ((i < this->first_dof()) ||
1538  (i >= this->end_dof())))
1539  continue;
1540 
1541  const Number rhs = it->second;
1542 
1543  os << "RHS for DoF " << i
1544  << ": " << rhs;
1545 
1546  os << std::endl;
1547  }
1548  }
1549 
1550  return os.str();
1551 }
1552 
1553 
1554 
1556  std::vector<dof_id_type> & elem_dofs,
1557  bool asymmetric_constraint_rows) const
1558 {
1559  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
1560  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
1561 
1562  // check for easy return
1563  if (this->_dof_constraints.empty())
1564  return;
1565 
1566  // The constrained matrix is built up as C^T K C.
1568 
1569 
1570  this->build_constraint_matrix (C, elem_dofs);
1571 
1572  LOG_SCOPE("constrain_elem_matrix()", "DofMap");
1573 
1574  // It is possible that the matrix is not constrained at all.
1575  if ((C.m() == matrix.m()) &&
1576  (C.n() == elem_dofs.size())) // It the matrix is constrained
1577  {
1578  // Compute the matrix-matrix-matrix product C^T K C
1579  matrix.left_multiply_transpose (C);
1580  matrix.right_multiply (C);
1581 
1582 
1583  libmesh_assert_equal_to (matrix.m(), matrix.n());
1584  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
1585  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
1586 
1587 
1588  for (std::size_t i=0; i<elem_dofs.size(); i++)
1589  // If the DOF is constrained
1590  if (this->is_constrained_dof(elem_dofs[i]))
1591  {
1592  for (unsigned int j=0; j<matrix.n(); j++)
1593  matrix(i,j) = 0.;
1594 
1595  matrix(i,i) = 1.;
1596 
1597  if (asymmetric_constraint_rows)
1598  {
1599  DofConstraints::const_iterator
1600  pos = _dof_constraints.find(elem_dofs[i]);
1601 
1602  libmesh_assert (pos != _dof_constraints.end());
1603 
1604  const DofConstraintRow & constraint_row = pos->second;
1605 
1606  // This is an overzealous assertion in the presence of
1607  // heterogenous constraints: we now can constrain "u_i = c"
1608  // with no other u_j terms involved.
1609  //
1610  // libmesh_assert (!constraint_row.empty());
1611 
1612  for (DofConstraintRow::const_iterator
1613  it=constraint_row.begin(); it != constraint_row.end();
1614  ++it)
1615  for (std::size_t j=0; j<elem_dofs.size(); j++)
1616  if (elem_dofs[j] == it->first)
1617  matrix(i,j) = -it->second;
1618  }
1619  }
1620  } // end if is constrained...
1621 }
1622 
1623 
1624 
1626  DenseVector<Number> & rhs,
1627  std::vector<dof_id_type> & elem_dofs,
1628  bool asymmetric_constraint_rows) const
1629 {
1630  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
1631  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
1632  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
1633 
1634  // check for easy return
1635  if (this->_dof_constraints.empty())
1636  return;
1637 
1638  // The constrained matrix is built up as C^T K C.
1639  // The constrained RHS is built up as C^T F
1641 
1642  this->build_constraint_matrix (C, elem_dofs);
1643 
1644  LOG_SCOPE("cnstrn_elem_mat_vec()", "DofMap");
1645 
1646  // It is possible that the matrix is not constrained at all.
1647  if ((C.m() == matrix.m()) &&
1648  (C.n() == elem_dofs.size())) // It the matrix is constrained
1649  {
1650  // Compute the matrix-matrix-matrix product C^T K C
1651  matrix.left_multiply_transpose (C);
1652  matrix.right_multiply (C);
1653 
1654 
1655  libmesh_assert_equal_to (matrix.m(), matrix.n());
1656  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
1657  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
1658 
1659 
1660  for (std::size_t i=0; i<elem_dofs.size(); i++)
1661  if (this->is_constrained_dof(elem_dofs[i]))
1662  {
1663  for (unsigned int j=0; j<matrix.n(); j++)
1664  matrix(i,j) = 0.;
1665 
1666  // If the DOF is constrained
1667  matrix(i,i) = 1.;
1668 
1669  // This will put a nonsymmetric entry in the constraint
1670  // row to ensure that the linear system produces the
1671  // correct value for the constrained DOF.
1672  if (asymmetric_constraint_rows)
1673  {
1674  DofConstraints::const_iterator
1675  pos = _dof_constraints.find(elem_dofs[i]);
1676 
1677  libmesh_assert (pos != _dof_constraints.end());
1678 
1679  const DofConstraintRow & constraint_row = pos->second;
1680 
1681  // p refinement creates empty constraint rows
1682  // libmesh_assert (!constraint_row.empty());
1683 
1684  for (DofConstraintRow::const_iterator
1685  it=constraint_row.begin(); it != constraint_row.end();
1686  ++it)
1687  for (std::size_t j=0; j<elem_dofs.size(); j++)
1688  if (elem_dofs[j] == it->first)
1689  matrix(i,j) = -it->second;
1690  }
1691  }
1692 
1693 
1694  // Compute the matrix-vector product C^T F
1695  DenseVector<Number> old_rhs(rhs);
1696 
1697  // compute matrix/vector product
1698  C.vector_mult_transpose(rhs, old_rhs);
1699  } // end if is constrained...
1700 }
1701 
1702 
1703 
1705  DenseVector<Number> & rhs,
1706  std::vector<dof_id_type> & elem_dofs,
1707  bool asymmetric_constraint_rows,
1708  int qoi_index) const
1709 {
1710  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
1711  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
1712  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
1713 
1714  // check for easy return
1715  if (this->_dof_constraints.empty())
1716  return;
1717 
1718  // The constrained matrix is built up as C^T K C.
1719  // The constrained RHS is built up as C^T (F - K H)
1722 
1723  this->build_constraint_matrix_and_vector (C, H, elem_dofs, qoi_index);
1724 
1725  LOG_SCOPE("hetero_cnstrn_elem_mat_vec()", "DofMap");
1726 
1727  // It is possible that the matrix is not constrained at all.
1728  if ((C.m() == matrix.m()) &&
1729  (C.n() == elem_dofs.size())) // It the matrix is constrained
1730  {
1731  // We may have rhs values to use later
1732  const DofConstraintValueMap * rhs_values = libmesh_nullptr;
1733  if (qoi_index < 0)
1734  rhs_values = &_primal_constraint_values;
1735  else
1736  {
1737  const AdjointDofConstraintValues::const_iterator
1738  it = _adjoint_constraint_values.find(qoi_index);
1739  if (it != _adjoint_constraint_values.end())
1740  rhs_values = &it->second;
1741  }
1742 
1743  // Compute matrix/vector product K H
1745  matrix.vector_mult(KH, H);
1746 
1747  // Compute the matrix-vector product C^T (F - KH)
1748  DenseVector<Number> F_minus_KH(rhs);
1749  F_minus_KH -= KH;
1750  C.vector_mult_transpose(rhs, F_minus_KH);
1751 
1752  // Compute the matrix-matrix-matrix product C^T K C
1753  matrix.left_multiply_transpose (C);
1754  matrix.right_multiply (C);
1755 
1756  libmesh_assert_equal_to (matrix.m(), matrix.n());
1757  libmesh_assert_equal_to (matrix.m(), elem_dofs.size());
1758  libmesh_assert_equal_to (matrix.n(), elem_dofs.size());
1759 
1760  for (std::size_t i=0; i<elem_dofs.size(); i++)
1761  {
1762  const dof_id_type dof_id = elem_dofs[i];
1763 
1764  if (this->is_constrained_dof(dof_id))
1765  {
1766  for (unsigned int j=0; j<matrix.n(); j++)
1767  matrix(i,j) = 0.;
1768 
1769  // If the DOF is constrained
1770  matrix(i,i) = 1.;
1771 
1772  // This will put a nonsymmetric entry in the constraint
1773  // row to ensure that the linear system produces the
1774  // correct value for the constrained DOF.
1775  if (asymmetric_constraint_rows)
1776  {
1777  DofConstraints::const_iterator
1778  pos = _dof_constraints.find(dof_id);
1779 
1780  libmesh_assert (pos != _dof_constraints.end());
1781 
1782  const DofConstraintRow & constraint_row = pos->second;
1783 
1784  for (DofConstraintRow::const_iterator
1785  it=constraint_row.begin(); it != constraint_row.end();
1786  ++it)
1787  for (std::size_t j=0; j<elem_dofs.size(); j++)
1788  if (elem_dofs[j] == it->first)
1789  matrix(i,j) = -it->second;
1790 
1791  if (rhs_values)
1792  {
1793  const DofConstraintValueMap::const_iterator valpos =
1794  rhs_values->find(dof_id);
1795 
1796  rhs(i) = (valpos == rhs_values->end()) ?
1797  0 : valpos->second;
1798  }
1799  }
1800  else
1801  rhs(i) = 0.;
1802  }
1803  }
1804 
1805  } // end if is constrained...
1806 }
1807 
1808 
1809 
1811  DenseVector<Number> & rhs,
1812  std::vector<dof_id_type> & elem_dofs,
1813  bool asymmetric_constraint_rows,
1814  int qoi_index) const
1815 {
1816  libmesh_assert_equal_to (elem_dofs.size(), matrix.m());
1817  libmesh_assert_equal_to (elem_dofs.size(), matrix.n());
1818  libmesh_assert_equal_to (elem_dofs.size(), rhs.size());
1819 
1820  // check for easy return
1821  if (this->_dof_constraints.empty())
1822  return;
1823 
1824  // The constrained matrix is built up as C^T K C.
1825  // The constrained RHS is built up as C^T (F - K H)
1828 
1829  this->build_constraint_matrix_and_vector (C, H, elem_dofs, qoi_index);
1830 
1831  LOG_SCOPE("hetero_cnstrn_elem_vec()", "DofMap");
1832 
1833  // It is possible that the matrix is not constrained at all.
1834  if ((C.m() == matrix.m()) &&
1835  (C.n() == elem_dofs.size())) // It the matrix is constrained
1836  {
1837  // We may have rhs values to use later
1838  const DofConstraintValueMap * rhs_values = libmesh_nullptr;
1839  if (qoi_index < 0)
1840  rhs_values = &_primal_constraint_values;
1841  else
1842  {
1843  const AdjointDofConstraintValues::const_iterator
1844  it = _adjoint_constraint_values.find(qoi_index);
1845  if (it != _adjoint_constraint_values.end())
1846  rhs_values = &it->second;
1847  }
1848 
1849  // Compute matrix/vector product K H
1851  matrix.vector_mult(KH, H);
1852 
1853  // Compute the matrix-vector product C^T (F - KH)
1854  DenseVector<Number> F_minus_KH(rhs);
1855  F_minus_KH -= KH;
1856  C.vector_mult_transpose(rhs, F_minus_KH);
1857 
1858  for (std::size_t i=0; i<elem_dofs.size(); i++)
1859  {
1860  const dof_id_type dof_id = elem_dofs[i];
1861 
1862  if (this->is_constrained_dof(dof_id))
1863  {
1864  // This will put a nonsymmetric entry in the constraint
1865  // row to ensure that the linear system produces the
1866  // correct value for the constrained DOF.
1867  if (asymmetric_constraint_rows && rhs_values)
1868  {
1869  const DofConstraintValueMap::const_iterator valpos =
1870  rhs_values->find(dof_id);
1871 
1872  rhs(i) = (valpos == rhs_values->end()) ?
1873  0 : valpos->second;
1874  }
1875  else
1876  rhs(i) = 0.;
1877  }
1878  }
1879 
1880  } // end if is constrained...
1881 }
1882 
1883 
1884 
1885 
1887  std::vector<dof_id_type> & row_dofs,
1888  std::vector<dof_id_type> & col_dofs,
1889  bool asymmetric_constraint_rows) const
1890 {
1891  libmesh_assert_equal_to (row_dofs.size(), matrix.m());
1892  libmesh_assert_equal_to (col_dofs.size(), matrix.n());
1893 
1894  // check for easy return
1895  if (this->_dof_constraints.empty())
1896  return;
1897 
1898  // The constrained matrix is built up as R^T K C.
1901 
1902  // Safeguard against the user passing us the same
1903  // object for row_dofs and col_dofs. If that is done
1904  // the calls to build_matrix would fail
1905  std::vector<dof_id_type> orig_row_dofs(row_dofs);
1906  std::vector<dof_id_type> orig_col_dofs(col_dofs);
1907 
1908  this->build_constraint_matrix (R, orig_row_dofs);
1909  this->build_constraint_matrix (C, orig_col_dofs);
1910 
1911  LOG_SCOPE("constrain_elem_matrix()", "DofMap");
1912 
1913  row_dofs = orig_row_dofs;
1914  col_dofs = orig_col_dofs;
1915 
1916  bool constraint_found = false;
1917 
1918  // K_constrained = R^T K C
1919 
1920  if ((R.m() == matrix.m()) &&
1921  (R.n() == row_dofs.size()))
1922  {
1923  matrix.left_multiply_transpose (R);
1924  constraint_found = true;
1925  }
1926 
1927  if ((C.m() == matrix.n()) &&
1928  (C.n() == col_dofs.size()))
1929  {
1930  matrix.right_multiply (C);
1931  constraint_found = true;
1932  }
1933 
1934  // It is possible that the matrix is not constrained at all.
1935  if (constraint_found)
1936  {
1937  libmesh_assert_equal_to (matrix.m(), row_dofs.size());
1938  libmesh_assert_equal_to (matrix.n(), col_dofs.size());
1939 
1940 
1941  for (std::size_t i=0; i<row_dofs.size(); i++)
1942  if (this->is_constrained_dof(row_dofs[i]))
1943  {
1944  for (unsigned int j=0; j<matrix.n(); j++)
1945  {
1946  if (row_dofs[i] != col_dofs[j])
1947  matrix(i,j) = 0.;
1948  else // If the DOF is constrained
1949  matrix(i,j) = 1.;
1950  }
1951 
1952  if (asymmetric_constraint_rows)
1953  {
1954  DofConstraints::const_iterator
1955  pos = _dof_constraints.find(row_dofs[i]);
1956 
1957  libmesh_assert (pos != _dof_constraints.end());
1958 
1959  const DofConstraintRow & constraint_row = pos->second;
1960 
1961  libmesh_assert (!constraint_row.empty());
1962 
1963  for (DofConstraintRow::const_iterator
1964  it=constraint_row.begin(); it != constraint_row.end();
1965  ++it)
1966  for (std::size_t j=0; j<col_dofs.size(); j++)
1967  if (col_dofs[j] == it->first)
1968  matrix(i,j) = -it->second;
1969  }
1970  }
1971  } // end if is constrained...
1972 }
1973 
1974 
1975 
1977  std::vector<dof_id_type> & row_dofs,
1978  bool) const
1979 {
1980  libmesh_assert_equal_to (rhs.size(), row_dofs.size());
1981 
1982  // check for easy return
1983  if (this->_dof_constraints.empty())
1984  return;
1985 
1986  // The constrained RHS is built up as R^T F.
1988 
1989  this->build_constraint_matrix (R, row_dofs);
1990 
1991  LOG_SCOPE("constrain_elem_vector()", "DofMap");
1992 
1993  // It is possible that the vector is not constrained at all.
1994  if ((R.m() == rhs.size()) &&
1995  (R.n() == row_dofs.size())) // if the RHS is constrained
1996  {
1997  // Compute the matrix-vector product
1998  DenseVector<Number> old_rhs(rhs);
1999  R.vector_mult_transpose(rhs, old_rhs);
2000 
2001  libmesh_assert_equal_to (row_dofs.size(), rhs.size());
2002 
2003  for (std::size_t i=0; i<row_dofs.size(); i++)
2004  if (this->is_constrained_dof(row_dofs[i]))
2005  {
2006  // If the DOF is constrained
2007  libmesh_assert (_dof_constraints.find(row_dofs[i]) != _dof_constraints.end());
2008 
2009  rhs(i) = 0;
2010  }
2011  } // end if the RHS is constrained.
2012 }
2013 
2014 
2015 
2017  DenseVector<Number> & w,
2018  std::vector<dof_id_type> & row_dofs,
2019  bool) const
2020 {
2021  libmesh_assert_equal_to (v.size(), row_dofs.size());
2022  libmesh_assert_equal_to (w.size(), row_dofs.size());
2023 
2024  // check for easy return
2025  if (this->_dof_constraints.empty())
2026  return;
2027 
2028  // The constrained RHS is built up as R^T F.
2030 
2031  this->build_constraint_matrix (R, row_dofs);
2032 
2033  LOG_SCOPE("cnstrn_elem_dyad_mat()", "DofMap");
2034 
2035  // It is possible that the vector is not constrained at all.
2036  if ((R.m() == v.size()) &&
2037  (R.n() == row_dofs.size())) // if the RHS is constrained
2038  {
2039  // Compute the matrix-vector products
2040  DenseVector<Number> old_v(v);
2041  DenseVector<Number> old_w(w);
2042 
2043  // compute matrix/vector product
2044  R.vector_mult_transpose(v, old_v);
2045  R.vector_mult_transpose(w, old_w);
2046 
2047  libmesh_assert_equal_to (row_dofs.size(), v.size());
2048  libmesh_assert_equal_to (row_dofs.size(), w.size());
2049 
2050  /* Constrain only v, not w. */
2051 
2052  for (std::size_t i=0; i<row_dofs.size(); i++)
2053  if (this->is_constrained_dof(row_dofs[i]))
2054  {
2055  // If the DOF is constrained
2056  libmesh_assert (_dof_constraints.find(row_dofs[i]) != _dof_constraints.end());
2057 
2058  v(i) = 0;
2059  }
2060  } // end if the RHS is constrained.
2061 }
2062 
2063 
2064 
2065 void DofMap::constrain_nothing (std::vector<dof_id_type> & dofs) const
2066 {
2067  // check for easy return
2068  if (this->_dof_constraints.empty())
2069  return;
2070 
2071  // All the work is done by \p build_constraint_matrix. We just need
2072  // a dummy matrix.
2074  this->build_constraint_matrix (R, dofs);
2075 }
2076 
2077 
2078 
2079 void DofMap::enforce_constraints_exactly (const System & system,
2081  bool homogeneous) const
2082 {
2083  parallel_object_only();
2084 
2085  if (!this->n_constrained_dofs())
2086  return;
2087 
2088  LOG_SCOPE("enforce_constraints_exactly()","DofMap");
2089 
2090  if (!v)
2091  v = system.solution.get();
2092 
2093  NumericVector<Number> * v_local = libmesh_nullptr; // will be initialized below
2094  NumericVector<Number> * v_global = libmesh_nullptr; // will be initialized below
2096  if (v->type() == SERIAL)
2097  {
2098  v_built = NumericVector<Number>::build(this->comm());
2099  v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
2100  v_built->close();
2101 
2102  for (dof_id_type i=v_built->first_local_index();
2103  i<v_built->last_local_index(); i++)
2104  v_built->set(i, (*v)(i));
2105  v_built->close();
2106  v_global = v_built.get();
2107 
2108  v_local = v;
2109  libmesh_assert (v_local->closed());
2110  }
2111  else if (v->type() == PARALLEL)
2112  {
2113  v_built = NumericVector<Number>::build(this->comm());
2114  v_built->init (v->size(), v->size(), true, SERIAL);
2115  v->localize(*v_built);
2116  v_built->close();
2117  v_local = v_built.get();
2118 
2119  v_global = v;
2120  }
2121  else if (v->type() == GHOSTED)
2122  {
2123  v_local = v;
2124  v_global = v;
2125  }
2126  else // unknown v->type()
2127  libmesh_error_msg("ERROR: Unknown v->type() == " << v->type());
2128 
2129  // We should never hit these asserts because we should error-out in
2130  // else clause above. Just to be sure we don't try to use v_local
2131  // and v_global uninitialized...
2132  libmesh_assert(v_local);
2133  libmesh_assert(v_global);
2134  libmesh_assert_equal_to (this, &(system.get_dof_map()));
2135 
2136  DofConstraints::const_iterator c_it = _dof_constraints.begin();
2137  const DofConstraints::const_iterator c_end = _dof_constraints.end();
2138 
2139  for ( ; c_it != c_end; ++c_it)
2140  {
2141  dof_id_type constrained_dof = c_it->first;
2142  if (constrained_dof < this->first_dof() ||
2143  constrained_dof >= this->end_dof())
2144  continue;
2145 
2146  const DofConstraintRow constraint_row = c_it->second;
2147 
2148  Number exact_value = 0;
2149  if (!homogeneous)
2150  {
2151  DofConstraintValueMap::const_iterator rhsit =
2152  _primal_constraint_values.find(constrained_dof);
2153  if (rhsit != _primal_constraint_values.end())
2154  exact_value = rhsit->second;
2155  }
2156  for (DofConstraintRow::const_iterator
2157  j=constraint_row.begin(); j != constraint_row.end();
2158  ++j)
2159  exact_value += j->second * (*v_local)(j->first);
2160 
2161  v_global->set(constrained_dof, exact_value);
2162  }
2163 
2164  // If the old vector was serial, we probably need to send our values
2165  // to other processors
2166  if (v->type() == SERIAL)
2167  {
2168 #ifndef NDEBUG
2169  v_global->close();
2170 #endif
2171  v_global->localize (*v);
2172  }
2173  v->close();
2174 }
2175 
2176 
2177 
2179  unsigned int q) const
2180 {
2181  parallel_object_only();
2182 
2183  if (!this->n_constrained_dofs())
2184  return;
2185 
2186  LOG_SCOPE("enforce_adjoint_constraints_exactly()", "DofMap");
2187 
2188  NumericVector<Number> * v_local = libmesh_nullptr; // will be initialized below
2189  NumericVector<Number> * v_global = libmesh_nullptr; // will be initialized below
2191  if (v.type() == SERIAL)
2192  {
2193  v_built = NumericVector<Number>::build(this->comm());
2194  v_built->init(this->n_dofs(), this->n_local_dofs(), true, PARALLEL);
2195  v_built->close();
2196 
2197  for (dof_id_type i=v_built->first_local_index();
2198  i<v_built->last_local_index(); i++)
2199  v_built->set(i, v(i));
2200  v_built->close();
2201  v_global = v_built.get();
2202 
2203  v_local = &v;
2204  libmesh_assert (v_local->closed());
2205  }
2206  else if (v.type() == PARALLEL)
2207  {
2208  v_built = NumericVector<Number>::build(this->comm());
2209  v_built->init (v.size(), v.size(), true, SERIAL);
2210  v.localize(*v_built);
2211  v_built->close();
2212  v_local = v_built.get();
2213 
2214  v_global = &v;
2215  }
2216  else if (v.type() == GHOSTED)
2217  {
2218  v_local = &v;
2219  v_global = &v;
2220  }
2221  else // unknown v.type()
2222  libmesh_error_msg("ERROR: Unknown v.type() == " << v.type());
2223 
2224  // We should never hit these asserts because we should error-out in
2225  // else clause above. Just to be sure we don't try to use v_local
2226  // and v_global uninitialized...
2227  libmesh_assert(v_local);
2228  libmesh_assert(v_global);
2229 
2230  // Do we have any non_homogeneous constraints?
2231  const AdjointDofConstraintValues::const_iterator
2232  adjoint_constraint_map_it = _adjoint_constraint_values.find(q);
2233  const DofConstraintValueMap * constraint_map =
2234  (adjoint_constraint_map_it == _adjoint_constraint_values.end()) ?
2235  libmesh_nullptr : &adjoint_constraint_map_it->second;
2236 
2237  DofConstraints::const_iterator c_it = _dof_constraints.begin();
2238  const DofConstraints::const_iterator c_end = _dof_constraints.end();
2239 
2240  for ( ; c_it != c_end; ++c_it)
2241  {
2242  dof_id_type constrained_dof = c_it->first;
2243  if (constrained_dof < this->first_dof() ||
2244  constrained_dof >= this->end_dof())
2245  continue;
2246 
2247  const DofConstraintRow constraint_row = c_it->second;
2248 
2249  Number exact_value = 0;
2250  if (constraint_map)
2251  {
2252  const DofConstraintValueMap::const_iterator
2253  adjoint_constraint_it =
2254  constraint_map->find(constrained_dof);
2255  if (adjoint_constraint_it != constraint_map->end())
2256  exact_value = adjoint_constraint_it->second;
2257  }
2258 
2259  for (DofConstraintRow::const_iterator
2260  j=constraint_row.begin(); j != constraint_row.end();
2261  ++j)
2262  exact_value += j->second * (*v_local)(j->first);
2263 
2264  v_global->set(constrained_dof, exact_value);
2265  }
2266 
2267  // If the old vector was serial, we probably need to send our values
2268  // to other processors
2269  if (v.type() == SERIAL)
2270  {
2271 #ifndef NDEBUG
2272  v_global->close();
2273 #endif
2274  v_global->localize (v);
2275  }
2276  v.close();
2277 }
2278 
2279 
2280 
2281 std::pair<Real, Real>
2283  NumericVector<Number> * v) const
2284 {
2285  if (!v)
2286  v = system.solution.get();
2287  NumericVector<Number> & vec = *v;
2288 
2289  // We'll assume the vector is closed
2290  libmesh_assert (vec.closed());
2291 
2292  Real max_absolute_error = 0., max_relative_error = 0.;
2293 
2294  const MeshBase & mesh = system.get_mesh();
2295 
2296  libmesh_assert_equal_to (this, &(system.get_dof_map()));
2297 
2298  // indices on each element
2299  std::vector<dof_id_type> local_dof_indices;
2300 
2301  for (const auto & elem : mesh.active_local_element_ptr_range())
2302  {
2303  this->dof_indices(elem, local_dof_indices);
2304  std::vector<dof_id_type> raw_dof_indices = local_dof_indices;
2305 
2306  // Constraint matrix for each element
2308 
2309  this->build_constraint_matrix (C, local_dof_indices);
2310 
2311  // Continue if the element is unconstrained
2312  if (!C.m())
2313  continue;
2314 
2315  libmesh_assert_equal_to (C.m(), raw_dof_indices.size());
2316  libmesh_assert_equal_to (C.n(), local_dof_indices.size());
2317 
2318  for (unsigned int i=0; i!=C.m(); ++i)
2319  {
2320  // Recalculate any constrained dof owned by this processor
2321  dof_id_type global_dof = raw_dof_indices[i];
2322  if (this->is_constrained_dof(global_dof) &&
2323  global_dof >= vec.first_local_index() &&
2324  global_dof < vec.last_local_index())
2325  {
2326 #ifndef NDEBUG
2327  DofConstraints::const_iterator
2328  pos = _dof_constraints.find(global_dof);
2329 
2330  libmesh_assert (pos != _dof_constraints.end());
2331 #endif
2332 
2333  Number exact_value = 0;
2334  DofConstraintValueMap::const_iterator rhsit =
2335  _primal_constraint_values.find(global_dof);
2336  if (rhsit != _primal_constraint_values.end())
2337  exact_value = rhsit->second;
2338 
2339  for (unsigned int j=0; j!=C.n(); ++j)
2340  {
2341  if (local_dof_indices[j] != global_dof)
2342  exact_value += C(i,j) *
2343  vec(local_dof_indices[j]);
2344  }
2345 
2346  max_absolute_error = std::max(max_absolute_error,
2347  std::abs(vec(global_dof) - exact_value));
2348  max_relative_error = std::max(max_relative_error,
2349  std::abs(vec(global_dof) - exact_value)
2350  / std::abs(exact_value));
2351  }
2352  }
2353  }
2354 
2355  return std::pair<Real, Real>(max_absolute_error, max_relative_error);
2356 }
2357 
2358 
2359 
2361  std::vector<dof_id_type> & elem_dofs,
2362  const bool called_recursively) const
2363 {
2364  LOG_SCOPE_IF("build_constraint_matrix()", "DofMap", !called_recursively);
2365 
2366  // Create a set containing the DOFs we already depend on
2367  typedef std::set<dof_id_type> RCSet;
2368  RCSet dof_set;
2369 
2370  bool we_have_constraints = false;
2371 
2372  // Next insert any other dofs the current dofs might be constrained
2373  // in terms of. Note that in this case we may not be done: Those
2374  // may in turn depend on others. So, we need to repeat this process
2375  // in that case until the system depends only on unconstrained
2376  // degrees of freedom.
2377  for (std::size_t i=0; i<elem_dofs.size(); i++)
2378  if (this->is_constrained_dof(elem_dofs[i]))
2379  {
2380  we_have_constraints = true;
2381 
2382  // If the DOF is constrained
2383  DofConstraints::const_iterator
2384  pos = _dof_constraints.find(elem_dofs[i]);
2385 
2386  libmesh_assert (pos != _dof_constraints.end());
2387 
2388  const DofConstraintRow & constraint_row = pos->second;
2389 
2390  // Constraint rows in p refinement may be empty
2391  //libmesh_assert (!constraint_row.empty());
2392 
2393  for (DofConstraintRow::const_iterator
2394  it=constraint_row.begin(); it != constraint_row.end();
2395  ++it)
2396  dof_set.insert (it->first);
2397  }
2398 
2399  // May be safe to return at this point
2400  // (but remember to stop the perflog)
2401  if (!we_have_constraints)
2402  return;
2403 
2404  for (std::size_t i=0; i != elem_dofs.size(); ++i)
2405  dof_set.erase (elem_dofs[i]);
2406 
2407  // If we added any DOFS then we need to do this recursively.
2408  // It is possible that we just added a DOF that is also
2409  // constrained!
2410  //
2411  // Also, we need to handle the special case of an element having DOFs
2412  // constrained in terms of other, local DOFs
2413  if (!dof_set.empty() || // case 1: constrained in terms of other DOFs
2414  !called_recursively) // case 2: constrained in terms of our own DOFs
2415  {
2416  const unsigned int old_size =
2417  cast_int<unsigned int>(elem_dofs.size());
2418 
2419  // Add new dependency dofs to the end of the current dof set
2420  elem_dofs.insert(elem_dofs.end(),
2421  dof_set.begin(), dof_set.end());
2422 
2423  // Now we can build the constraint matrix.
2424  // Note that resize also zeros for a DenseMatrix<Number>.
2425  C.resize (old_size,
2426  cast_int<unsigned int>(elem_dofs.size()));
2427 
2428  // Create the C constraint matrix.
2429  for (unsigned int i=0; i != old_size; i++)
2430  if (this->is_constrained_dof(elem_dofs[i]))
2431  {
2432  // If the DOF is constrained
2433  DofConstraints::const_iterator
2434  pos = _dof_constraints.find(elem_dofs[i]);
2435 
2436  libmesh_assert (pos != _dof_constraints.end());
2437 
2438  const DofConstraintRow & constraint_row = pos->second;
2439 
2440  // p refinement creates empty constraint rows
2441  // libmesh_assert (!constraint_row.empty());
2442 
2443  for (DofConstraintRow::const_iterator
2444  it=constraint_row.begin(); it != constraint_row.end();
2445  ++it)
2446  for (std::size_t j=0; j != elem_dofs.size(); j++)
2447  if (elem_dofs[j] == it->first)
2448  C(i,j) = it->second;
2449  }
2450  else
2451  {
2452  C(i,i) = 1.;
2453  }
2454 
2455  // May need to do this recursively. It is possible
2456  // that we just replaced a constrained DOF with another
2457  // constrained DOF.
2458  DenseMatrix<Number> Cnew;
2459 
2460  this->build_constraint_matrix (Cnew, elem_dofs, true);
2461 
2462  if ((C.n() == Cnew.m()) &&
2463  (Cnew.n() == elem_dofs.size())) // If the constraint matrix
2464  C.right_multiply(Cnew); // is constrained...
2465 
2466  libmesh_assert_equal_to (C.n(), elem_dofs.size());
2467  }
2468 }
2469 
2470 
2471 
2473  DenseVector<Number> & H,
2474  std::vector<dof_id_type> & elem_dofs,
2475  int qoi_index,
2476  const bool called_recursively) const
2477 {
2478  LOG_SCOPE_IF("build_constraint_matrix_and_vector()", "DofMap", !called_recursively);
2479 
2480  // Create a set containing the DOFs we already depend on
2481  typedef std::set<dof_id_type> RCSet;
2482  RCSet dof_set;
2483 
2484  bool we_have_constraints = false;
2485 
2486  // Next insert any other dofs the current dofs might be constrained
2487  // in terms of. Note that in this case we may not be done: Those
2488  // may in turn depend on others. So, we need to repeat this process
2489  // in that case until the system depends only on unconstrained
2490  // degrees of freedom.
2491  for (std::size_t i=0; i<elem_dofs.size(); i++)
2492  if (this->is_constrained_dof(elem_dofs[i]))
2493  {
2494  we_have_constraints = true;
2495 
2496  // If the DOF is constrained
2497  DofConstraints::const_iterator
2498  pos = _dof_constraints.find(elem_dofs[i]);
2499 
2500  libmesh_assert (pos != _dof_constraints.end());
2501 
2502  const DofConstraintRow & constraint_row = pos->second;
2503 
2504  // Constraint rows in p refinement may be empty
2505  //libmesh_assert (!constraint_row.empty());
2506 
2507  for (DofConstraintRow::const_iterator
2508  it=constraint_row.begin(); it != constraint_row.end();
2509  ++it)
2510  dof_set.insert (it->first);
2511  }
2512 
2513  // May be safe to return at this point
2514  // (but remember to stop the perflog)
2515  if (!we_have_constraints)
2516  return;
2517 
2518  for (std::size_t i=0; i != elem_dofs.size(); ++i)
2519  dof_set.erase (elem_dofs[i]);
2520 
2521  // If we added any DOFS then we need to do this recursively.
2522  // It is possible that we just added a DOF that is also
2523  // constrained!
2524  //
2525  // Also, we need to handle the special case of an element having DOFs
2526  // constrained in terms of other, local DOFs
2527  if (!dof_set.empty() || // case 1: constrained in terms of other DOFs
2528  !called_recursively) // case 2: constrained in terms of our own DOFs
2529  {
2530  const DofConstraintValueMap * rhs_values = libmesh_nullptr;
2531  if (qoi_index < 0)
2532  rhs_values = &_primal_constraint_values;
2533  else
2534  {
2535  const AdjointDofConstraintValues::const_iterator
2536  it = _adjoint_constraint_values.find(qoi_index);
2537  if (it != _adjoint_constraint_values.end())
2538  rhs_values = &it->second;
2539  }
2540 
2541  const unsigned int old_size =
2542  cast_int<unsigned int>(elem_dofs.size());
2543 
2544  // Add new dependency dofs to the end of the current dof set
2545  elem_dofs.insert(elem_dofs.end(),
2546  dof_set.begin(), dof_set.end());
2547 
2548  // Now we can build the constraint matrix and vector.
2549  // Note that resize also zeros for a DenseMatrix and DenseVector
2550  C.resize (old_size,
2551  cast_int<unsigned int>(elem_dofs.size()));
2552  H.resize (old_size);
2553 
2554  // Create the C constraint matrix.
2555  for (unsigned int i=0; i != old_size; i++)
2556  if (this->is_constrained_dof(elem_dofs[i]))
2557  {
2558  // If the DOF is constrained
2559  DofConstraints::const_iterator
2560  pos = _dof_constraints.find(elem_dofs[i]);
2561 
2562  libmesh_assert (pos != _dof_constraints.end());
2563 
2564  const DofConstraintRow & constraint_row = pos->second;
2565 
2566  // p refinement creates empty constraint rows
2567  // libmesh_assert (!constraint_row.empty());
2568 
2569  for (DofConstraintRow::const_iterator
2570  it=constraint_row.begin(); it != constraint_row.end();
2571  ++it)
2572  for (std::size_t j=0; j != elem_dofs.size(); j++)
2573  if (elem_dofs[j] == it->first)
2574  C(i,j) = it->second;
2575 
2576  if (rhs_values)
2577  {
2578  DofConstraintValueMap::const_iterator rhsit =
2579  rhs_values->find(elem_dofs[i]);
2580  if (rhsit != rhs_values->end())
2581  H(i) = rhsit->second;
2582  }
2583  }
2584  else
2585  {
2586  C(i,i) = 1.;
2587  }
2588 
2589  // May need to do this recursively. It is possible
2590  // that we just replaced a constrained DOF with another
2591  // constrained DOF.
2592  DenseMatrix<Number> Cnew;
2593  DenseVector<Number> Hnew;
2594 
2595  this->build_constraint_matrix_and_vector (Cnew, Hnew, elem_dofs,
2596  qoi_index, true);
2597 
2598  if ((C.n() == Cnew.m()) && // If the constraint matrix
2599  (Cnew.n() == elem_dofs.size())) // is constrained...
2600  {
2601  // If x = Cy + h and y = Dz + g
2602  // Then x = (CD)z + (Cg + h)
2603  C.vector_mult_add(H, 1, Hnew);
2604 
2605  C.right_multiply(Cnew);
2606  }
2607 
2608  libmesh_assert_equal_to (C.n(), elem_dofs.size());
2609  }
2610 }
2611 
2612 
2614 {
2615  // This function must be run on all processors at once
2616  parallel_object_only();
2617 
2618  // Return immediately if there's nothing to gather
2619  if (this->n_processors() == 1)
2620  return;
2621 
2622  // We might get to return immediately if none of the processors
2623  // found any constraints
2624  unsigned int has_constraints = !_dof_constraints.empty()
2625 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2626  || !_node_constraints.empty()
2627 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2628  ;
2629  this->comm().max(has_constraints);
2630  if (!has_constraints)
2631  return;
2632 
2633  // If we have heterogenous adjoint constraints we need to
2634  // communicate those too.
2635  const unsigned int max_qoi_num =
2636  _adjoint_constraint_values.empty() ?
2637  0 : _adjoint_constraint_values.rbegin()->first;
2638 
2639  // We might have calculated constraints for constrained dofs
2640  // which have support on other processors.
2641  // Push these out first.
2642  {
2643  std::vector<std::set<dof_id_type>> pushed_ids(this->n_processors());
2644 
2645 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2646  std::vector<std::set<dof_id_type>> pushed_node_ids(this->n_processors());
2647 #endif
2648 
2649  const unsigned int sys_num = this->sys_number();
2650 
2652  foreign_elem_it = mesh.active_not_local_elements_begin(),
2653  foreign_elem_end = mesh.active_not_local_elements_end();
2654 
2655  // Collect the constraints to push to each processor
2656  for (; foreign_elem_it != foreign_elem_end; ++foreign_elem_it)
2657  {
2658  const Elem * elem = *foreign_elem_it;
2659 
2660  const unsigned short n_nodes = elem->n_nodes();
2661 
2662  // Just checking dof_indices on the foreign element isn't
2663  // enough. Consider a central hanging node between a coarse
2664  // Q2/Q1 element and its finer neighbors on a higher-ranked
2665  // processor. The coarse element's processor will own the node,
2666  // and will thereby own the pressure dof on that node, despite
2667  // the fact that that pressure dof doesn't directly exist on the
2668  // coarse element!
2669  //
2670  // So, we loop through dofs manually.
2671 
2672  {
2673  const unsigned int n_vars = elem->n_vars(sys_num);
2674  for (unsigned int v=0; v != n_vars; ++v)
2675  {
2676  const unsigned int n_comp = elem->n_comp(sys_num,v);
2677  for (unsigned int c=0; c != n_comp; ++c)
2678  {
2679  const unsigned int id =
2680  elem->dof_number(sys_num,v,c);
2681  if (this->is_constrained_dof(id))
2682  pushed_ids[elem->processor_id()].insert(id);
2683  }
2684  }
2685  }
2686 
2687  for (unsigned short n = 0; n != n_nodes; ++n)
2688  {
2689  const Node & node = elem->node_ref(n);
2690  const unsigned int n_vars = node.n_vars(sys_num);
2691  for (unsigned int v=0; v != n_vars; ++v)
2692  {
2693  const unsigned int n_comp = node.n_comp(sys_num,v);
2694  for (unsigned int c=0; c != n_comp; ++c)
2695  {
2696  const unsigned int id =
2697  node.dof_number(sys_num,v,c);
2698  if (this->is_constrained_dof(id))
2699  pushed_ids[elem->processor_id()].insert(id);
2700  }
2701  }
2702  }
2703 
2704 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2705  for (unsigned short n = 0; n != n_nodes; ++n)
2706  if (this->is_constrained_node(elem->node_ptr(n)))
2707  pushed_node_ids[elem->processor_id()].insert(elem->node_id(n));
2708 #endif
2709  }
2710 
2711  // Now trade constraint rows
2712  for (processor_id_type p = 0; p != this->n_processors(); ++p)
2713  {
2714  // Push to processor procup while receiving from procdown
2715  processor_id_type procup =
2716  cast_int<processor_id_type>((this->processor_id() + p) %
2717  this->n_processors());
2718  processor_id_type procdown =
2719  cast_int<processor_id_type>((this->n_processors() +
2720  this->processor_id() - p) %
2721  this->n_processors());
2722 
2723  // Pack the dof constraint rows and rhs's to push to procup
2724  const std::size_t pushed_ids_size = pushed_ids[procup].size();
2725  std::vector<std::vector<dof_id_type>> pushed_keys(pushed_ids_size);
2726  std::vector<std::vector<Real>> pushed_vals(pushed_ids_size);
2727  std::vector<Number> pushed_rhss(pushed_ids_size);
2728 
2729  std::vector<std::vector<Number>>
2730  pushed_adj_rhss(max_qoi_num,
2731  std::vector<Number>(pushed_ids_size));
2732  std::set<dof_id_type>::const_iterator it = pushed_ids[procup].begin();
2733  for (std::size_t i = 0; it != pushed_ids[procup].end();
2734  ++i, ++it)
2735  {
2736  const dof_id_type pushed_id = *it;
2737  DofConstraintRow & row = _dof_constraints[pushed_id];
2738  std::size_t row_size = row.size();
2739  pushed_keys[i].reserve(row_size);
2740  pushed_vals[i].reserve(row_size);
2741  for (DofConstraintRow::iterator j = row.begin();
2742  j != row.end(); ++j)
2743  {
2744  pushed_keys[i].push_back(j->first);
2745  pushed_vals[i].push_back(j->second);
2746  }
2747 
2748  DofConstraintValueMap::const_iterator rhsit =
2749  _primal_constraint_values.find(pushed_id);
2750  pushed_rhss[i] =
2751  (rhsit == _primal_constraint_values.end()) ?
2752  0 : rhsit->second;
2753 
2754  for (unsigned int q = 0; q != max_qoi_num; ++q)
2755  {
2756  AdjointDofConstraintValues::const_iterator adjoint_map_it =
2757  _adjoint_constraint_values.find(q);
2758 
2759  if (adjoint_map_it == _adjoint_constraint_values.end())
2760  continue;
2761 
2762  const DofConstraintValueMap & constraint_map =
2763  adjoint_map_it->second;
2764 
2765  DofConstraintValueMap::const_iterator rhsit =
2766  constraint_map.find(pushed_id);
2767 
2768  pushed_adj_rhss[q][i] =
2769  (rhsit == constraint_map.end()) ?
2770  0 : rhsit->second;
2771  }
2772  }
2773 
2774 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2775  // Pack the node constraint rows to push to procup
2776  const std::size_t pushed_nodes_size = pushed_node_ids[procup].size();
2777  std::vector<std::vector<dof_id_type>> pushed_node_keys(pushed_nodes_size);
2778  std::vector<std::vector<Real>> pushed_node_vals(pushed_nodes_size);
2779  std::vector<Point> pushed_node_offsets(pushed_nodes_size);
2780  std::set<dof_id_type>::const_iterator node_it = pushed_node_ids[procup].begin();
2781  for (std::size_t i = 0; node_it != pushed_node_ids[procup].end();
2782  ++i, ++node_it)
2783  {
2784  const Node * node = mesh.node_ptr(*node_it);
2785  NodeConstraintRow & row = _node_constraints[node].first;
2786  std::size_t row_size = row.size();
2787  pushed_node_keys[i].reserve(row_size);
2788  pushed_node_vals[i].reserve(row_size);
2789  for (NodeConstraintRow::iterator j = row.begin();
2790  j != row.end(); ++j)
2791  {
2792  pushed_node_keys[i].push_back(j->first->id());
2793  pushed_node_vals[i].push_back(j->second);
2794  }
2795  pushed_node_offsets[i] = _node_constraints[node].second;
2796  }
2797 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2798 
2799  // Trade pushed dof constraint rows
2800  std::vector<dof_id_type> pushed_ids_from_me
2801  (pushed_ids[procup].begin(), pushed_ids[procup].end());
2802  std::vector<dof_id_type> pushed_ids_to_me;
2803  std::vector<std::vector<dof_id_type>> pushed_keys_to_me;
2804  std::vector<std::vector<Real>> pushed_vals_to_me;
2805  std::vector<Number> pushed_rhss_to_me;
2806  std::vector<std::vector<Number>> pushed_adj_rhss_to_me;
2807  this->comm().send_receive(procup, pushed_ids_from_me,
2808  procdown, pushed_ids_to_me);
2809  this->comm().send_receive(procup, pushed_keys,
2810  procdown, pushed_keys_to_me);
2811  this->comm().send_receive(procup, pushed_vals,
2812  procdown, pushed_vals_to_me);
2813  this->comm().send_receive(procup, pushed_rhss,
2814  procdown, pushed_rhss_to_me);
2815  this->comm().send_receive(procup, pushed_adj_rhss,
2816  procdown, pushed_adj_rhss_to_me);
2817  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_keys_to_me.size());
2818  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_vals_to_me.size());
2819  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_rhss_to_me.size());
2820 
2821 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2822  // Trade pushed node constraint rows
2823  std::vector<dof_id_type> pushed_node_ids_from_me
2824  (pushed_node_ids[procup].begin(), pushed_node_ids[procup].end());
2825  std::vector<dof_id_type> pushed_node_ids_to_me;
2826  std::vector<std::vector<dof_id_type>> pushed_node_keys_to_me;
2827  std::vector<std::vector<Real>> pushed_node_vals_to_me;
2828  std::vector<Point> pushed_node_offsets_to_me;
2829  this->comm().send_receive(procup, pushed_node_ids_from_me,
2830  procdown, pushed_node_ids_to_me);
2831  this->comm().send_receive(procup, pushed_node_keys,
2832  procdown, pushed_node_keys_to_me);
2833  this->comm().send_receive(procup, pushed_node_vals,
2834  procdown, pushed_node_vals_to_me);
2835  this->comm().send_receive(procup, pushed_node_offsets,
2836  procdown, pushed_node_offsets_to_me);
2837 
2838  // Note that we aren't doing a send_receive on the Nodes
2839  // themselves. At this point we should only be pushing out
2840  // "raw" constraints, and there should be no
2841  // constrained-by-constrained-by-etc. situations that could
2842  // involve non-semilocal nodes.
2843 
2844  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_keys_to_me.size());
2845  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_vals_to_me.size());
2846  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_offsets_to_me.size());
2847 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2848 
2849  // Add the dof constraints that I've been sent
2850  for (std::size_t i = 0; i != pushed_ids_to_me.size(); ++i)
2851  {
2852  libmesh_assert_equal_to (pushed_keys_to_me[i].size(), pushed_vals_to_me[i].size());
2853 
2854  dof_id_type constrained = pushed_ids_to_me[i];
2855 
2856  // If we don't already have a constraint for this dof,
2857  // add the one we were sent
2858  if (!this->is_constrained_dof(constrained))
2859  {
2860  DofConstraintRow & row = _dof_constraints[constrained];
2861  for (std::size_t j = 0; j != pushed_keys_to_me[i].size(); ++j)
2862  {
2863  row[pushed_keys_to_me[i][j]] = pushed_vals_to_me[i][j];
2864  }
2865  if (libmesh_isnan(pushed_rhss_to_me[i]))
2866  libmesh_assert(pushed_keys_to_me[i].empty());
2867  if (pushed_rhss_to_me[i] != Number(0))
2868  _primal_constraint_values[constrained] = pushed_rhss_to_me[i];
2869  else
2870  _primal_constraint_values.erase(constrained);
2871 
2872  for (unsigned int q = 0; q != max_qoi_num; ++q)
2873  {
2874  AdjointDofConstraintValues::iterator adjoint_map_it =
2875  _adjoint_constraint_values.find(q);
2876 
2877  if ((adjoint_map_it == _adjoint_constraint_values.end()) &&
2878  pushed_adj_rhss_to_me[q][constrained] == Number(0))
2879  continue;
2880 
2881  if (adjoint_map_it == _adjoint_constraint_values.end())
2882  adjoint_map_it = _adjoint_constraint_values.insert
2883  (std::make_pair(q,DofConstraintValueMap())).first;
2884 
2885  DofConstraintValueMap & constraint_map =
2886  adjoint_map_it->second;
2887 
2888  if (pushed_adj_rhss_to_me[q][i] != Number(0))
2889  constraint_map[constrained] =
2890  pushed_adj_rhss_to_me[q][i];
2891  else
2892  constraint_map.erase(constrained);
2893  }
2894  }
2895  }
2896 
2897 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2898  // Add the node constraints that I've been sent
2899  for (std::size_t i = 0; i != pushed_node_ids_to_me.size(); ++i)
2900  {
2901  libmesh_assert_equal_to (pushed_node_keys_to_me[i].size(), pushed_node_vals_to_me[i].size());
2902 
2903  dof_id_type constrained_id = pushed_node_ids_to_me[i];
2904 
2905  // If we don't already have a constraint for this node,
2906  // add the one we were sent
2907  const Node * constrained = mesh.node_ptr(constrained_id);
2908  if (!this->is_constrained_node(constrained))
2909  {
2910  NodeConstraintRow & row = _node_constraints[constrained].first;
2911  for (std::size_t j = 0; j != pushed_node_keys_to_me[i].size(); ++j)
2912  {
2913  const Node * key_node = mesh.node_ptr(pushed_node_keys_to_me[i][j]);
2914  libmesh_assert(key_node);
2915  row[key_node] = pushed_node_vals_to_me[i][j];
2916  }
2917  _node_constraints[constrained].second = pushed_node_offsets_to_me[i];
2918  }
2919  }
2920 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
2921  }
2922  }
2923 
2924  // Now start checking for any other constraints we need
2925  // to know about, requesting them recursively.
2926 
2927  // Create sets containing the DOFs and nodes we already depend on
2928  typedef std::set<dof_id_type> DoF_RCSet;
2929  DoF_RCSet unexpanded_dofs;
2930 
2931  for (DofConstraints::iterator i = _dof_constraints.begin();
2932  i != _dof_constraints.end(); ++i)
2933  {
2934  unexpanded_dofs.insert(i->first);
2935  }
2936 
2937  // Gather all the dof constraints we need
2938  this->gather_constraints(mesh, unexpanded_dofs, false);
2939 
2940  // Gather all the node constraints we need
2941 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
2942  typedef std::set<const Node *> Node_RCSet;
2943  Node_RCSet unexpanded_nodes;
2944 
2945  for (NodeConstraints::iterator i = _node_constraints.begin();
2946  i != _node_constraints.end(); ++i)
2947  {
2948  unexpanded_nodes.insert(i->first);
2949  }
2950 
2951  // We have to keep recursing while the unexpanded set is
2952  // nonempty on *any* processor
2953  bool unexpanded_set_nonempty = !unexpanded_nodes.empty();
2954  this->comm().max(unexpanded_set_nonempty);
2955 
2956  while (unexpanded_set_nonempty)
2957  {
2958  // Let's make sure we don't lose sync in this loop.
2959  parallel_object_only();
2960 
2961  // Request sets
2962  Node_RCSet node_request_set;
2963 
2964  // Request sets to send to each processor
2965  std::vector<std::vector<dof_id_type>>
2966  requested_node_ids(this->n_processors());
2967 
2968  // And the sizes of each
2969  std::vector<dof_id_type>
2970  node_ids_on_proc(this->n_processors(), 0);
2971 
2972  // Fill (and thereby sort and uniq!) the main request sets
2973  for (Node_RCSet::iterator i = unexpanded_nodes.begin();
2974  i != unexpanded_nodes.end(); ++i)
2975  {
2976  NodeConstraintRow & row = _node_constraints[*i].first;
2977  for (NodeConstraintRow::iterator j = row.begin();
2978  j != row.end(); ++j)
2979  {
2980  const Node * const node = j->first;
2981  libmesh_assert(node);
2982 
2983  // If it's non-local and we haven't already got a
2984  // constraint for it, we might need to ask for one
2985  if ((node->processor_id() != this->processor_id()) &&
2986  !_node_constraints.count(node))
2987  node_request_set.insert(node);
2988  }
2989  }
2990 
2991  // Clear the unexpanded constraint sets; we're about to expand
2992  // them
2993  unexpanded_nodes.clear();
2994 
2995  // Count requests by processor
2996  for (Node_RCSet::iterator i = node_request_set.begin();
2997  i != node_request_set.end(); ++i)
2998  {
2999  libmesh_assert(*i);
3000  libmesh_assert_less ((*i)->processor_id(), this->n_processors());
3001  node_ids_on_proc[(*i)->processor_id()]++;
3002  }
3003 
3004  for (processor_id_type p = 0; p != this->n_processors(); ++p)
3005  {
3006  requested_node_ids[p].reserve(node_ids_on_proc[p]);
3007  }
3008 
3009  // Prepare each processor's request set
3010  for (Node_RCSet::iterator i = node_request_set.begin();
3011  i != node_request_set.end(); ++i)
3012  {
3013  requested_node_ids[(*i)->processor_id()].push_back((*i)->id());
3014  }
3015 
3016  // Now request constraint rows from other processors
3017  for (processor_id_type p=1; p != this->n_processors(); ++p)
3018  {
3019  // Trade my requests with processor procup and procdown
3020  processor_id_type procup =
3021  cast_int<processor_id_type>((this->processor_id() + p) %
3022  this->n_processors());
3023  processor_id_type procdown =
3024  cast_int<processor_id_type>((this->n_processors() +
3025  this->processor_id() - p) %
3026  this->n_processors());
3027  std::vector<dof_id_type> node_request_to_fill;
3028 
3029  this->comm().send_receive(procup, requested_node_ids[procup],
3030  procdown, node_request_to_fill);
3031 
3032  // Fill those requests
3033  std::vector<std::vector<dof_id_type>>
3034  node_row_keys(node_request_to_fill.size());
3035  std::vector<std::vector<Real>>
3036  node_row_vals(node_request_to_fill.size());
3037  std::vector<Point>
3038  node_row_rhss(node_request_to_fill.size());
3039 
3040  // FIXME - this could be an unordered set, given a
3041  // hash<pointers> specialization
3042  std::set<const Node *> nodes_requested;
3043 
3044  for (std::size_t i=0; i != node_request_to_fill.size(); ++i)
3045  {
3046  dof_id_type constrained_id = node_request_to_fill[i];
3047  const Node * constrained_node = mesh.node_ptr(constrained_id);
3048  if (_node_constraints.count(constrained_node))
3049  {
3050  const NodeConstraintRow & row = _node_constraints[constrained_node].first;
3051  std::size_t row_size = row.size();
3052  node_row_keys[i].reserve(row_size);
3053  node_row_vals[i].reserve(row_size);
3054  for (NodeConstraintRow::const_iterator j = row.begin();
3055  j != row.end(); ++j)
3056  {
3057  const Node * node = j->first;
3058  node_row_keys[i].push_back(node->id());
3059  node_row_vals[i].push_back(j->second);
3060 
3061  // If we're not sure whether our send
3062  // destination already has this node, let's give
3063  // it a copy.
3064  if (node->processor_id() != procdown)
3065  nodes_requested.insert(node);
3066 
3067  // We can have 0 nodal constraint
3068  // coefficients, where no Lagrange constraint
3069  // exists but non-Lagrange basis constraints
3070  // might.
3071  // libmesh_assert(j->second);
3072  }
3073  node_row_rhss[i] = _node_constraints[constrained_node].second;
3074  }
3075  }
3076 
3077  // Trade back the results
3078  std::vector<std::vector<dof_id_type>> node_filled_keys;
3079  std::vector<std::vector<Real>> node_filled_vals;
3080  std::vector<Point> node_filled_rhss;
3081 
3082  this->comm().send_receive(procdown, node_row_keys,
3083  procup, node_filled_keys);
3084  this->comm().send_receive(procdown, node_row_vals,
3085  procup, node_filled_vals);
3086  this->comm().send_receive(procdown, node_row_rhss,
3087  procup, node_filled_rhss);
3088 
3089  // Constraining nodes might not even exist on our subset of
3090  // a distributed mesh, so let's make them exist.
3091  if (!mesh.is_serial())
3092  this->comm().send_receive_packed_range
3093  (procdown, &mesh, nodes_requested.begin(), nodes_requested.end(),
3095  (Node**)libmesh_nullptr);
3096 
3097  libmesh_assert_equal_to (node_filled_keys.size(), requested_node_ids[procup].size());
3098  libmesh_assert_equal_to (node_filled_vals.size(), requested_node_ids[procup].size());
3099  libmesh_assert_equal_to (node_filled_rhss.size(), requested_node_ids[procup].size());
3100 
3101  for (std::size_t i=0; i != requested_node_ids[procup].size(); ++i)
3102  {
3103  libmesh_assert_equal_to (node_filled_keys[i].size(), node_filled_vals[i].size());
3104  if (!node_filled_keys[i].empty())
3105  {
3106  dof_id_type constrained_id = requested_node_ids[procup][i];
3107  const Node * constrained_node = mesh.node_ptr(constrained_id);
3108  NodeConstraintRow & row = _node_constraints[constrained_node].first;
3109  for (std::size_t j = 0; j != node_filled_keys[i].size(); ++j)
3110  {
3111  const Node * key_node =
3112  mesh.node_ptr(node_filled_keys[i][j]);
3113  libmesh_assert(key_node);
3114  row[key_node] = node_filled_vals[i][j];
3115  }
3116  _node_constraints[constrained_node].second = node_filled_rhss[i];
3117 
3118  // And prepare to check for more recursive constraints
3119  unexpanded_nodes.insert(constrained_node);
3120  }
3121  }
3122  }
3123 
3124  // We have to keep recursing while the unexpanded set is
3125  // nonempty on *any* processor
3126  unexpanded_set_nonempty = !unexpanded_nodes.empty();
3127  this->comm().max(unexpanded_set_nonempty);
3128  }
3129 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3130 }
3131 
3132 
3133 
3135 {
3136  // We've computed our local constraints, but they may depend on
3137  // non-local constraints that we'll need to take into account.
3138  this->allgather_recursive_constraints(mesh);
3139 
3140  // Create a set containing the DOFs we already depend on
3141  typedef std::set<dof_id_type> RCSet;
3142  RCSet unexpanded_set;
3143 
3144  for (DofConstraints::iterator i = _dof_constraints.begin();
3145  i != _dof_constraints.end(); ++i)
3146  unexpanded_set.insert(i->first);
3147 
3148  while (!unexpanded_set.empty())
3149  for (RCSet::iterator i = unexpanded_set.begin();
3150  i != unexpanded_set.end(); /* nothing */)
3151  {
3152  // If the DOF is constrained
3153  DofConstraints::iterator
3154  pos = _dof_constraints.find(*i);
3155 
3156  libmesh_assert (pos != _dof_constraints.end());
3157 
3158  DofConstraintRow & constraint_row = pos->second;
3159 
3160  DofConstraintValueMap::iterator rhsit =
3161  _primal_constraint_values.find(*i);
3162  Number constraint_rhs = (rhsit == _primal_constraint_values.end()) ?
3163  0 : rhsit->second;
3164 
3165  std::vector<dof_id_type> constraints_to_expand;
3166 
3167  for (DofConstraintRow::const_iterator
3168  it=constraint_row.begin(); it != constraint_row.end();
3169  ++it)
3170  if (it->first != *i && this->is_constrained_dof(it->first))
3171  {
3172  unexpanded_set.insert(it->first);
3173  constraints_to_expand.push_back(it->first);
3174  }
3175 
3176  for (std::size_t j = 0; j != constraints_to_expand.size();
3177  ++j)
3178  {
3179  dof_id_type expandable = constraints_to_expand[j];
3180 
3181  const Real this_coef = constraint_row[expandable];
3182 
3183  DofConstraints::const_iterator
3184  subpos = _dof_constraints.find(expandable);
3185 
3186  libmesh_assert (subpos != _dof_constraints.end());
3187 
3188  const DofConstraintRow & subconstraint_row = subpos->second;
3189 
3190  for (DofConstraintRow::const_iterator
3191  it=subconstraint_row.begin();
3192  it != subconstraint_row.end(); ++it)
3193  {
3194  constraint_row[it->first] += it->second * this_coef;
3195  }
3196  DofConstraintValueMap::const_iterator subrhsit =
3197  _primal_constraint_values.find(expandable);
3198  if (subrhsit != _primal_constraint_values.end())
3199  constraint_rhs += subrhsit->second * this_coef;
3200 
3201  constraint_row.erase(expandable);
3202  }
3203 
3204  if (rhsit == _primal_constraint_values.end())
3205  {
3206  if (constraint_rhs != Number(0))
3207  _primal_constraint_values[*i] = constraint_rhs;
3208  else
3209  _primal_constraint_values.erase(*i);
3210  }
3211  else
3212  {
3213  if (constraint_rhs != Number(0))
3214  rhsit->second = constraint_rhs;
3215  else
3216  _primal_constraint_values.erase(rhsit);
3217  }
3218 
3219  if (constraints_to_expand.empty())
3220  unexpanded_set.erase(i++);
3221  else
3222  ++i;
3223  }
3224 
3225  // In parallel we can't guarantee that nodes/dofs which constrain
3226  // others are on processors which are aware of that constraint, yet
3227  // we need such awareness for sparsity pattern generation. So send
3228  // other processors any constraints they might need to know about.
3229  this->scatter_constraints(mesh);
3230 
3231  // Now that we have our root constraint dependencies sorted out, add
3232  // them to the send_list
3233  this->add_constraints_to_send_list();
3234 }
3235 
3236 
3238 {
3239  // At this point each processor with a constrained node knows
3240  // the corresponding constraint row, but we also need each processor
3241  // with a constrainer node to know the corresponding row(s).
3242 
3243  // This function must be run on all processors at once
3244  parallel_object_only();
3245 
3246  // Return immediately if there's nothing to gather
3247  if (this->n_processors() == 1)
3248  return;
3249 
3250  // We might get to return immediately if none of the processors
3251  // found any constraints
3252  unsigned int has_constraints = !_dof_constraints.empty()
3253 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3254  || !_node_constraints.empty()
3255 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3256  ;
3257  this->comm().max(has_constraints);
3258  if (!has_constraints)
3259  return;
3260 
3261 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3262  std::vector<std::set<dof_id_type>> pushed_node_ids(this->n_processors());
3263 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3264 
3265  std::vector<std::set<dof_id_type>> pushed_ids(this->n_processors());
3266 
3267  // Collect the dof constraints I need to push to each processor
3268  dof_id_type constrained_proc_id = 0;
3269  for (DofConstraints::iterator i = _dof_constraints.begin();
3270  i != _dof_constraints.end(); ++i)
3271  {
3272  const dof_id_type constrained = i->first;
3273  while (constrained >= _end_df[constrained_proc_id])
3274  constrained_proc_id++;
3275 
3276  if (constrained_proc_id != this->processor_id())
3277  continue;
3278 
3279  DofConstraintRow & row = i->second;
3280  for (DofConstraintRow::iterator j = row.begin();
3281  j != row.end(); ++j)
3282  {
3283  const dof_id_type constraining = j->first;
3284 
3285  processor_id_type constraining_proc_id = 0;
3286  while (constraining >= _end_df[constraining_proc_id])
3287  constraining_proc_id++;
3288 
3289  if (constraining_proc_id != this->processor_id() &&
3290  constraining_proc_id != constrained_proc_id)
3291  pushed_ids[constraining_proc_id].insert(constrained);
3292  }
3293  }
3294 
3295 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3296  // Collect the node constraints to push to each processor
3297  for (NodeConstraints::iterator i = _node_constraints.begin();
3298  i != _node_constraints.end(); ++i)
3299  {
3300  const Node * constrained = i->first;
3301 
3302  if (constrained->processor_id() != this->processor_id())
3303  continue;
3304 
3305  NodeConstraintRow & row = i->second.first;
3306  for (NodeConstraintRow::iterator j = row.begin();
3307  j != row.end(); ++j)
3308  {
3309  const Node * constraining = j->first;
3310 
3311  if (constraining->processor_id() != this->processor_id() &&
3312  constraining->processor_id() != constrained->processor_id())
3313  pushed_node_ids[constraining->processor_id()].insert(constrained->id());
3314  }
3315  }
3316 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3317 
3318  // Now trade constraint rows
3319  for (processor_id_type p = 0; p != this->n_processors(); ++p)
3320  {
3321  // Push to processor procup while receiving from procdown
3322  processor_id_type procup =
3323  cast_int<processor_id_type>((this->processor_id() + p) %
3324  this->n_processors());
3325  processor_id_type procdown =
3326  cast_int<processor_id_type>((this->n_processors() +
3327  this->processor_id() - p) %
3328  this->n_processors());
3329 
3330  // Pack the dof constraint rows and rhs's to push to procup
3331  const std::size_t pushed_ids_size = pushed_ids[procup].size();
3332  std::vector<std::vector<dof_id_type>> pushed_keys(pushed_ids_size);
3333  std::vector<std::vector<Real>> pushed_vals(pushed_ids_size);
3334  std::vector<Number> pushed_rhss(pushed_ids_size);
3335 
3336  std::set<dof_id_type>::const_iterator it;
3337  std::size_t push_i;
3338  for (push_i = 0, it = pushed_ids[procup].begin();
3339  it != pushed_ids[procup].end(); ++push_i, ++it)
3340  {
3341  const dof_id_type constrained = *it;
3342  DofConstraintRow & row = _dof_constraints[constrained];
3343  std::size_t row_size = row.size();
3344  pushed_keys[push_i].reserve(row_size);
3345  pushed_vals[push_i].reserve(row_size);
3346  for (DofConstraintRow::iterator j = row.begin();
3347  j != row.end(); ++j)
3348  {
3349  pushed_keys[push_i].push_back(j->first);
3350  pushed_vals[push_i].push_back(j->second);
3351  }
3352 
3353  DofConstraintValueMap::const_iterator rhsit =
3354  _primal_constraint_values.find(constrained);
3355  pushed_rhss[push_i] = (rhsit == _primal_constraint_values.end()) ?
3356  0 : rhsit->second;
3357  }
3358 
3359 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3360  // Pack the node constraint rows to push to procup
3361  const std::size_t pushed_node_ids_size = pushed_node_ids[procup].size();
3362  std::vector<std::vector<dof_id_type>> pushed_node_keys(pushed_node_ids_size);
3363  std::vector<std::vector<Real>> pushed_node_vals(pushed_node_ids_size);
3364  std::vector<Point> pushed_node_offsets(pushed_node_ids_size);
3365  std::set<const Node *> pushed_nodes;
3366 
3367  for (push_i = 0, it = pushed_node_ids[procup].begin();
3368  it != pushed_node_ids[procup].end(); ++push_i, ++it)
3369  {
3370  const Node * constrained = mesh.node_ptr(*it);
3371 
3372  if (constrained->processor_id() != procdown)
3373  pushed_nodes.insert(constrained);
3374 
3375  NodeConstraintRow & row = _node_constraints[constrained].first;
3376  std::size_t row_size = row.size();
3377  pushed_node_keys[push_i].reserve(row_size);
3378  pushed_node_vals[push_i].reserve(row_size);
3379  for (NodeConstraintRow::iterator j = row.begin();
3380  j != row.end(); ++j)
3381  {
3382  const Node * constraining = j->first;
3383 
3384  pushed_node_keys[push_i].push_back(constraining->id());
3385  pushed_node_vals[push_i].push_back(j->second);
3386 
3387  if (constraining->processor_id() != procup)
3388  pushed_nodes.insert(constraining);
3389  }
3390  pushed_node_offsets[push_i] = _node_constraints[constrained].second;
3391  }
3392 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3393 
3394  // Trade pushed dof constraint rows
3395  std::vector<dof_id_type> pushed_ids_from_me
3396  (pushed_ids[procup].begin(), pushed_ids[procup].end());
3397  std::vector<dof_id_type> pushed_ids_to_me;
3398  std::vector<std::vector<dof_id_type>> pushed_keys_to_me;
3399  std::vector<std::vector<Real>> pushed_vals_to_me;
3400  std::vector<Number> pushed_rhss_to_me;
3401  this->comm().send_receive(procup, pushed_ids_from_me,
3402  procdown, pushed_ids_to_me);
3403  this->comm().send_receive(procup, pushed_keys,
3404  procdown, pushed_keys_to_me);
3405  this->comm().send_receive(procup, pushed_vals,
3406  procdown, pushed_vals_to_me);
3407  this->comm().send_receive(procup, pushed_rhss,
3408  procdown, pushed_rhss_to_me);
3409  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_keys_to_me.size());
3410  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_vals_to_me.size());
3411  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_rhss_to_me.size());
3412 
3413 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3414  // Trade pushed node constraint rows
3415  std::vector<dof_id_type> pushed_node_ids_from_me
3416  (pushed_node_ids[procup].begin(), pushed_node_ids[procup].end());
3417  std::vector<dof_id_type> pushed_node_ids_to_me;
3418  std::vector<std::vector<dof_id_type>> pushed_node_keys_to_me;
3419  std::vector<std::vector<Real>> pushed_node_vals_to_me;
3420  std::vector<Point> pushed_node_offsets_to_me;
3421  this->comm().send_receive(procup, pushed_node_ids_from_me,
3422  procdown, pushed_node_ids_to_me);
3423  this->comm().send_receive(procup, pushed_node_keys,
3424  procdown, pushed_node_keys_to_me);
3425  this->comm().send_receive(procup, pushed_node_vals,
3426  procdown, pushed_node_vals_to_me);
3427  this->comm().send_receive(procup, pushed_node_offsets,
3428  procdown, pushed_node_offsets_to_me);
3429 
3430  // Constraining nodes might not even exist on our subset of
3431  // a distributed mesh, so let's make them exist.
3432  if (!mesh.is_serial())
3433  this->comm().send_receive_packed_range
3434  (procup, &mesh, pushed_nodes.begin(), pushed_nodes.end(),
3435  procdown, &mesh, mesh_inserter_iterator<Node>(mesh),
3436  (Node**)libmesh_nullptr);
3437 
3438  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_keys_to_me.size());
3439  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_vals_to_me.size());
3440  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_offsets_to_me.size());
3441 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3442 
3443  // Add the dof constraints that I've been sent
3444  for (std::size_t i = 0; i != pushed_ids_to_me.size(); ++i)
3445  {
3446  libmesh_assert_equal_to (pushed_keys_to_me[i].size(), pushed_vals_to_me[i].size());
3447 
3448  dof_id_type constrained = pushed_ids_to_me[i];
3449 
3450  // If we don't already have a constraint for this dof,
3451  // add the one we were sent
3452  if (!this->is_constrained_dof(constrained))
3453  {
3454  DofConstraintRow & row = _dof_constraints[constrained];
3455  for (std::size_t j = 0; j != pushed_keys_to_me[i].size(); ++j)
3456  {
3457  row[pushed_keys_to_me[i][j]] = pushed_vals_to_me[i][j];
3458  }
3459  if (pushed_rhss_to_me[i] != Number(0))
3460  _primal_constraint_values[constrained] = pushed_rhss_to_me[i];
3461  else
3462  _primal_constraint_values.erase(constrained);
3463  }
3464  }
3465 
3466 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3467  // Add the node constraints that I've been sent
3468  for (std::size_t i = 0; i != pushed_node_ids_to_me.size(); ++i)
3469  {
3470  libmesh_assert_equal_to (pushed_node_keys_to_me[i].size(), pushed_node_vals_to_me[i].size());
3471 
3472  dof_id_type constrained_id = pushed_node_ids_to_me[i];
3473 
3474  // If we don't already have a constraint for this node,
3475  // add the one we were sent
3476  const Node * constrained = mesh.node_ptr(constrained_id);
3477  if (!this->is_constrained_node(constrained))
3478  {
3479  NodeConstraintRow & row = _node_constraints[constrained].first;
3480  for (std::size_t j = 0; j != pushed_node_keys_to_me[i].size(); ++j)
3481  {
3482  const Node * key_node = mesh.node_ptr(pushed_node_keys_to_me[i][j]);
3483  row[key_node] = pushed_node_vals_to_me[i][j];
3484  }
3485  _node_constraints[constrained].second = pushed_node_offsets_to_me[i];
3486  }
3487  }
3488 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3489  }
3490 
3491  // Next we need to push constraints to processors which don't own
3492  // the constrained dof, don't own the constraining dof, but own an
3493  // element supporting the constraining dof.
3494  //
3495  // We need to be able to quickly look up constrained dof ids by what
3496  // constrains them, so that we can handle the case where we see a
3497  // foreign element containing one of our constraining DoF ids and we
3498  // need to push that constraint.
3499  //
3500  // Getting distributed adaptive sparsity patterns right is hard.
3501 
3502  typedef std::map<dof_id_type, std::set<dof_id_type>> DofConstrainsMap;
3503  DofConstrainsMap dof_id_constrains;
3504 
3505  for (DofConstraints::iterator i = _dof_constraints.begin();
3506  i != _dof_constraints.end(); ++i)
3507  {
3508  const dof_id_type constrained = i->first;
3509  DofConstraintRow & row = i->second;
3510  for (DofConstraintRow::iterator j = row.begin();
3511  j != row.end(); ++j)
3512  {
3513  const dof_id_type constraining = j->first;
3514 
3515  dof_id_type constraining_proc_id = 0;
3516  while (constraining >= _end_df[constraining_proc_id])
3517  constraining_proc_id++;
3518 
3519  if (constraining_proc_id == this->processor_id())
3520  dof_id_constrains[constraining].insert(constrained);
3521  }
3522  }
3523 
3524  // Loop over all foreign elements, find any supporting our
3525  // constrained dof indices.
3526  pushed_ids.clear();
3527  pushed_ids.resize(this->n_processors());
3528 
3531  for (; it != end; ++it)
3532  {
3533  const Elem * elem = *it;
3534 
3535  std::vector<dof_id_type> my_dof_indices;
3536  this->dof_indices (elem, my_dof_indices);
3537 
3538  for (std::size_t i=0; i != my_dof_indices.size(); ++i)
3539  {
3540  DofConstrainsMap::const_iterator dcmi =
3541  dof_id_constrains.find(my_dof_indices[i]);
3542  if (dcmi != dof_id_constrains.end())
3543  {
3544  for (DofConstrainsMap::mapped_type::const_iterator mti =
3545  dcmi->second.begin();
3546  mti != dcmi->second.end(); ++mti)
3547  {
3548  const dof_id_type constrained = *mti;
3549 
3550  dof_id_type the_constrained_proc_id = 0;
3551  while (constrained >= _end_df[the_constrained_proc_id])
3552  the_constrained_proc_id++;
3553 
3554  const dof_id_type elemproc = elem->processor_id();
3555  if (elemproc != the_constrained_proc_id)
3556  pushed_ids[elemproc].insert(constrained);
3557  }
3558  }
3559  }
3560  }
3561 
3562  // One last trade of constraint rows
3563  for (processor_id_type p = 0; p != this->n_processors(); ++p)
3564  {
3565  // Push to processor procup while receiving from procdown
3566  processor_id_type procup =
3567  cast_int<processor_id_type>((this->processor_id() + p) %
3568  this->n_processors());
3569  processor_id_type procdown =
3570  cast_int<processor_id_type>((this->n_processors() +
3571  this->processor_id() - p) %
3572  this->n_processors());
3573 
3574  // Pack the dof constraint rows and rhs's to push to procup
3575  const std::size_t pushed_ids_size = pushed_ids[procup].size();
3576  std::vector<std::vector<dof_id_type>> pushed_keys(pushed_ids_size);
3577  std::vector<std::vector<Real>> pushed_vals(pushed_ids_size);
3578  std::vector<Number> pushed_rhss(pushed_ids_size);
3579 
3580  // As long as we're declaring them outside the loop, let's initialize them too!
3581  std::set<dof_id_type>::const_iterator pushed_ids_iter = pushed_ids[procup].begin();
3582  std::size_t push_i = 0;
3583  for ( ; pushed_ids_iter != pushed_ids[procup].end(); ++push_i, ++pushed_ids_iter)
3584  {
3585  const dof_id_type constrained = *pushed_ids_iter;
3586  DofConstraintRow & row = _dof_constraints[constrained];
3587  std::size_t row_size = row.size();
3588  pushed_keys[push_i].reserve(row_size);
3589  pushed_vals[push_i].reserve(row_size);
3590  for (DofConstraintRow::iterator j = row.begin();
3591  j != row.end(); ++j)
3592  {
3593  pushed_keys[push_i].push_back(j->first);
3594  pushed_vals[push_i].push_back(j->second);
3595  }
3596 
3597  DofConstraintValueMap::const_iterator rhsit =
3598  _primal_constraint_values.find(constrained);
3599  pushed_rhss[push_i] = (rhsit == _primal_constraint_values.end()) ?
3600  0 : rhsit->second;
3601  }
3602 
3603  // Trade pushed dof constraint rows
3604  std::vector<dof_id_type> pushed_ids_from_me
3605  (pushed_ids[procup].begin(), pushed_ids[procup].end());
3606  std::vector<dof_id_type> pushed_ids_to_me;
3607  std::vector<std::vector<dof_id_type>> pushed_keys_to_me;
3608  std::vector<std::vector<Real>> pushed_vals_to_me;
3609  std::vector<Number> pushed_rhss_to_me;
3610  this->comm().send_receive(procup, pushed_ids_from_me,
3611  procdown, pushed_ids_to_me);
3612  this->comm().send_receive(procup, pushed_keys,
3613  procdown, pushed_keys_to_me);
3614  this->comm().send_receive(procup, pushed_vals,
3615  procdown, pushed_vals_to_me);
3616  this->comm().send_receive(procup, pushed_rhss,
3617  procdown, pushed_rhss_to_me);
3618  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_keys_to_me.size());
3619  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_vals_to_me.size());
3620  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_rhss_to_me.size());
3621 
3622  // Add the dof constraints that I've been sent
3623  for (std::size_t i = 0; i != pushed_ids_to_me.size(); ++i)
3624  {
3625  libmesh_assert_equal_to (pushed_keys_to_me[i].size(), pushed_vals_to_me[i].size());
3626 
3627  dof_id_type constrained = pushed_ids_to_me[i];
3628 
3629  // If we don't already have a constraint for this dof,
3630  // add the one we were sent
3631  if (!this->is_constrained_dof(constrained))
3632  {
3633  DofConstraintRow & row = _dof_constraints[constrained];
3634  for (std::size_t j = 0; j != pushed_keys_to_me[i].size(); ++j)
3635  {
3636  row[pushed_keys_to_me[i][j]] = pushed_vals_to_me[i][j];
3637  }
3638 
3639  if (pushed_rhss_to_me[i] != Number(0))
3640  _primal_constraint_values[constrained] = pushed_rhss_to_me[i];
3641  else
3642  _primal_constraint_values.erase(constrained);
3643  }
3644  }
3645  }
3646 
3647  // Finally, we need to handle the case of remote dof coupling. If a
3648  // processor's element is coupled to a ghost element, then the
3649  // processor needs to know about all constraints which affect the
3650  // dofs on that ghost element, so we'll have to query the ghost
3651  // element's owner.
3652 
3653  GhostingFunctor::map_type elements_to_couple;
3654 
3655  // Man, I wish we had guaranteed unique_ptr availability...
3656  std::set<CouplingMatrix*> temporary_coupling_matrices;
3657 
3658  this->merge_ghost_functor_outputs
3659  (elements_to_couple,
3660  temporary_coupling_matrices,
3661  this->coupling_functors_begin(),
3662  this->coupling_functors_end(),
3665  this->processor_id());
3666 
3667  // Each ghost-coupled element's owner should get a request for its dofs
3668  std::set<dof_id_type> requested_dofs;
3669 
3670  GhostingFunctor::map_type::iterator etg_it = elements_to_couple.begin();
3671  const GhostingFunctor::map_type::iterator etg_end = elements_to_couple.end();
3672  for (; etg_it != etg_end; ++etg_it)
3673  {
3674  const Elem *elem = etg_it->first;
3675 
3676  // FIXME - optimize for the non-fully-coupled case?
3677  std::vector<dof_id_type> element_dofs;
3678  this->dof_indices(elem, element_dofs);
3679 
3680  for (std::size_t i=0; i != element_dofs.size(); ++i)
3681  requested_dofs.insert(element_dofs[i]);
3682  }
3683 
3684  this->gather_constraints(mesh, requested_dofs, false);
3685 }
3686 
3687 
3689  std::set<dof_id_type> & unexpanded_dofs,
3690  bool /*look_for_constrainees*/)
3691 {
3692  typedef std::set<dof_id_type> DoF_RCSet;
3693 
3694  // If we have heterogenous adjoint constraints we need to
3695  // communicate those too.
3696  const unsigned int max_qoi_num =
3697  _adjoint_constraint_values.empty() ?
3698  0 : _adjoint_constraint_values.rbegin()->first;
3699 
3700  // We have to keep recursing while the unexpanded set is
3701  // nonempty on *any* processor
3702  bool unexpanded_set_nonempty = !unexpanded_dofs.empty();
3703  this->comm().max(unexpanded_set_nonempty);
3704 
3705  while (unexpanded_set_nonempty)
3706  {
3707  // Let's make sure we don't lose sync in this loop.
3708  parallel_object_only();
3709 
3710  // Request sets
3711  DoF_RCSet dof_request_set;
3712 
3713  // Request sets to send to each processor
3714  std::vector<std::vector<dof_id_type>>
3715  requested_dof_ids(this->n_processors());
3716 
3717  // And the sizes of each
3718  std::vector<dof_id_type>
3719  dof_ids_on_proc(this->n_processors(), 0);
3720 
3721  // Fill (and thereby sort and uniq!) the main request sets
3722  for (DoF_RCSet::iterator i = unexpanded_dofs.begin();
3723  i != unexpanded_dofs.end(); ++i)
3724  {
3725  const dof_id_type unexpanded_dof = *i;
3726  DofConstraints::const_iterator
3727  pos = _dof_constraints.find(unexpanded_dof);
3728 
3729  // If we were asked for a DoF and we don't already have a
3730  // constraint for it, then we need to check for one.
3731  if (pos == _dof_constraints.end())
3732  {
3733  if (((unexpanded_dof < this->first_dof()) ||
3734  (unexpanded_dof >= this->end_dof())) &&
3735  !_dof_constraints.count(unexpanded_dof))
3736  dof_request_set.insert(unexpanded_dof);
3737  }
3738  // If we were asked for a DoF and we already have a
3739  // constraint for it, then we need to check if the
3740  // constraint is recursive.
3741  else
3742  {
3743  const DofConstraintRow & row = pos->second;
3744  for (DofConstraintRow::const_iterator j = row.begin();
3745  j != row.end(); ++j)
3746  {
3747  const dof_id_type constraining_dof = j->first;
3748 
3749  // If it's non-local and we haven't already got a
3750  // constraint for it, we might need to ask for one
3751  if (((constraining_dof < this->first_dof()) ||
3752  (constraining_dof >= this->end_dof())) &&
3753  !_dof_constraints.count(constraining_dof))
3754  dof_request_set.insert(constraining_dof);
3755  }
3756  }
3757  }
3758 
3759  // Clear the unexpanded constraint set; we're about to expand it
3760  unexpanded_dofs.clear();
3761 
3762  // Count requests by processor
3763  processor_id_type proc_id = 0;
3764  for (DoF_RCSet::iterator i = dof_request_set.begin();
3765  i != dof_request_set.end(); ++i)
3766  {
3767  while (*i >= _end_df[proc_id])
3768  proc_id++;
3769  dof_ids_on_proc[proc_id]++;
3770  }
3771 
3772  for (processor_id_type p = 0; p != this->n_processors(); ++p)
3773  {
3774  requested_dof_ids[p].reserve(dof_ids_on_proc[p]);
3775  }
3776 
3777  // Prepare each processor's request set
3778  proc_id = 0;
3779  for (DoF_RCSet::iterator i = dof_request_set.begin();
3780  i != dof_request_set.end(); ++i)
3781  {
3782  while (*i >= _end_df[proc_id])
3783  proc_id++;
3784  requested_dof_ids[proc_id].push_back(*i);
3785  }
3786 
3787  // Now request constraint rows from other processors
3788  for (processor_id_type p=1; p != this->n_processors(); ++p)
3789  {
3790  // Trade my requests with processor procup and procdown
3791  processor_id_type procup =
3792  cast_int<processor_id_type>((this->processor_id() + p) %
3793  this->n_processors());
3794  processor_id_type procdown =
3795  cast_int<processor_id_type>((this->n_processors() +
3796  this->processor_id() - p) %
3797  this->n_processors());
3798  std::vector<dof_id_type> dof_request_to_fill;
3799 
3800  this->comm().send_receive(procup, requested_dof_ids[procup],
3801  procdown, dof_request_to_fill);
3802 
3803  // Fill those requests
3804  std::vector<std::vector<dof_id_type>> dof_row_keys(dof_request_to_fill.size());
3805 
3806  std::vector<std::vector<Real>> dof_row_vals(dof_request_to_fill.size());
3807  std::vector<Number> dof_row_rhss(dof_request_to_fill.size());
3808  std::vector<std::vector<Number>>
3809  dof_adj_rhss(max_qoi_num,
3810  std::vector<Number>(dof_request_to_fill.size()));
3811  for (std::size_t i=0; i != dof_request_to_fill.size(); ++i)
3812  {
3813  dof_id_type constrained = dof_request_to_fill[i];
3814  if (_dof_constraints.count(constrained))
3815  {
3816  DofConstraintRow & row = _dof_constraints[constrained];
3817  std::size_t row_size = row.size();
3818  dof_row_keys[i].reserve(row_size);
3819  dof_row_vals[i].reserve(row_size);
3820  for (DofConstraintRow::iterator j = row.begin();
3821  j != row.end(); ++j)
3822  {
3823  dof_row_keys[i].push_back(j->first);
3824  dof_row_vals[i].push_back(j->second);
3825 
3826  // We should never have a 0 constraint
3827  // coefficient; that's implicit via sparse
3828  // constraint storage
3829  libmesh_assert(j->second);
3830  }
3831  DofConstraintValueMap::const_iterator rhsit =
3832  _primal_constraint_values.find(constrained);
3833  dof_row_rhss[i] = (rhsit == _primal_constraint_values.end()) ?
3834  0 : rhsit->second;
3835 
3836  for (unsigned int q = 0; q != max_qoi_num; ++q)
3837  {
3838  AdjointDofConstraintValues::const_iterator adjoint_map_it =
3839  _adjoint_constraint_values.find(q);
3840 
3841  if (adjoint_map_it == _adjoint_constraint_values.end())
3842  continue;
3843 
3844  const DofConstraintValueMap & constraint_map =
3845  adjoint_map_it->second;
3846 
3847  DofConstraintValueMap::const_iterator rhsit =
3848  constraint_map.find(constrained);
3849  dof_adj_rhss[q][i] = (rhsit == constraint_map.end()) ?
3850  0 : rhsit->second;
3851  }
3852  }
3853  else
3854  {
3855  // Get NaN from Real, where it should exist, not
3856  // from Number, which may be std::complex, in which
3857  // case quiet_NaN() silently returns zero, rather
3858  // than sanely returning NaN or throwing an
3859  // exception or sending Stroustrup hate mail.
3860  dof_row_rhss[i] =
3861  std::numeric_limits<Real>::quiet_NaN();
3862 
3863  // Make sure we don't get caught by "!isnan(NaN)"
3864  // bugs again.
3865  libmesh_assert(libmesh_isnan(dof_row_rhss[i]));
3866  }
3867  }
3868 
3869  // Trade back the results
3870  std::vector<std::vector<dof_id_type>> dof_filled_keys;
3871  std::vector<std::vector<Real>> dof_filled_vals;
3872  std::vector<Number> dof_filled_rhss;
3873  std::vector<std::vector<Number>> adj_filled_rhss;
3874  this->comm().send_receive(procdown, dof_row_keys,
3875  procup, dof_filled_keys);
3876  this->comm().send_receive(procdown, dof_row_vals,
3877  procup, dof_filled_vals);
3878  this->comm().send_receive(procdown, dof_row_rhss,
3879  procup, dof_filled_rhss);
3880  this->comm().send_receive(procdown, dof_adj_rhss,
3881  procup, adj_filled_rhss);
3882 
3883  libmesh_assert_equal_to (dof_filled_keys.size(), requested_dof_ids[procup].size());
3884  libmesh_assert_equal_to (dof_filled_vals.size(), requested_dof_ids[procup].size());
3885  libmesh_assert_equal_to (dof_filled_rhss.size(), requested_dof_ids[procup].size());
3886 #ifndef NDEBUG
3887  for (std::size_t q=0; q != adj_filled_rhss.size(); ++q)
3888  libmesh_assert_equal_to (adj_filled_rhss[q].size(), requested_dof_ids[procup].size());
3889 #endif
3890 
3891  // Add any new constraint rows we've found
3892  for (std::size_t i=0; i != requested_dof_ids[procup].size(); ++i)
3893  {
3894  libmesh_assert_equal_to (dof_filled_keys[i].size(), dof_filled_vals[i].size());
3895  if (!libmesh_isnan(dof_filled_rhss[i]))
3896  {
3897  dof_id_type constrained = requested_dof_ids[procup][i];
3898  DofConstraintRow & row = _dof_constraints[constrained];
3899  for (std::size_t j = 0; j != dof_filled_keys[i].size(); ++j)
3900  row[dof_filled_keys[i][j]] = dof_filled_vals[i][j];
3901  if (dof_filled_rhss[i] != Number(0))
3902  _primal_constraint_values[constrained] = dof_filled_rhss[i];
3903  else
3904  _primal_constraint_values.erase(constrained);
3905 
3906  for (unsigned int q = 0; q != max_qoi_num; ++q)
3907  {
3908  AdjointDofConstraintValues::iterator adjoint_map_it =
3909  _adjoint_constraint_values.find(q);
3910 
3911  if ((adjoint_map_it == _adjoint_constraint_values.end()) &&
3912  adj_filled_rhss[q][constrained] == Number(0))
3913  continue;
3914 
3915  if (adjoint_map_it == _adjoint_constraint_values.end())
3916  adjoint_map_it = _adjoint_constraint_values.insert
3917  (std::make_pair(q,DofConstraintValueMap())).first;
3918 
3919  DofConstraintValueMap & constraint_map =
3920  adjoint_map_it->second;
3921 
3922  if (adj_filled_rhss[q][i] != Number(0))
3923  constraint_map[constrained] =
3924  adj_filled_rhss[q][i];
3925  else
3926  constraint_map.erase(constrained);
3927  }
3928 
3929  // And prepare to check for more recursive constraints
3930  if (!dof_filled_keys[i].empty())
3931  unexpanded_dofs.insert(constrained);
3932  }
3933  }
3934  }
3935 
3936  // We have to keep recursing while the unexpanded set is
3937  // nonempty on *any* processor
3938  unexpanded_set_nonempty = !unexpanded_dofs.empty();
3939  this->comm().max(unexpanded_set_nonempty);
3940  }
3941 }
3942 
3944 {
3945  // This function must be run on all processors at once
3946  parallel_object_only();
3947 
3948  // Return immediately if there's nothing to gather
3949  if (this->n_processors() == 1)
3950  return;
3951 
3952  // We might get to return immediately if none of the processors
3953  // found any constraints
3954  unsigned int has_constraints = !_dof_constraints.empty();
3955  this->comm().max(has_constraints);
3956  if (!has_constraints)
3957  return;
3958 
3959  for (DofConstraints::iterator i = _dof_constraints.begin();
3960  i != _dof_constraints.end(); ++i)
3961  {
3962  dof_id_type constrained_dof = i->first;
3963 
3964  // We only need the dependencies of our own constrained dofs
3965  if (constrained_dof < this->first_dof() ||
3966  constrained_dof >= this->end_dof())
3967  continue;
3968 
3969  DofConstraintRow & constraint_row = i->second;
3970  for (DofConstraintRow::const_iterator
3971  j=constraint_row.begin(); j != constraint_row.end();
3972  ++j)
3973  {
3974  dof_id_type constraint_dependency = j->first;
3975 
3976  // No point in adding one of our own dofs to the send_list
3977  if (constraint_dependency >= this->first_dof() &&
3978  constraint_dependency < this->end_dof())
3979  continue;
3980 
3981  _send_list.push_back(constraint_dependency);
3982  }
3983  }
3984 }
3985 
3986 
3987 
3988 #endif // LIBMESH_ENABLE_CONSTRAINTS
3989 
3990 
3991 #ifdef LIBMESH_ENABLE_AMR
3992 
3993 void DofMap::constrain_p_dofs (unsigned int var,
3994  const Elem * elem,
3995  unsigned int s,
3996  unsigned int p)
3997 {
3998  // We're constraining dofs on elem which correspond to p refinement
3999  // levels above p - this only makes sense if elem's p refinement
4000  // level is above p.
4001  libmesh_assert_greater (elem->p_level(), p);
4002  libmesh_assert_less (s, elem->n_sides());
4003 
4004  const unsigned int sys_num = this->sys_number();
4005  const unsigned int dim = elem->dim();
4006  ElemType type = elem->type();
4007  FEType low_p_fe_type = this->variable_type(var);
4008  FEType high_p_fe_type = this->variable_type(var);
4009  low_p_fe_type.order = static_cast<Order>(low_p_fe_type.order + p);
4010  high_p_fe_type.order = static_cast<Order>(high_p_fe_type.order +
4011  elem->p_level());
4012 
4013  const unsigned int n_nodes = elem->n_nodes();
4014  for (unsigned int n = 0; n != n_nodes; ++n)
4015  if (elem->is_node_on_side(n, s))
4016  {
4017  const Node & node = elem->node_ref(n);
4018  const unsigned int low_nc =
4019  FEInterface::n_dofs_at_node (dim, low_p_fe_type, type, n);
4020  const unsigned int high_nc =
4021  FEInterface::n_dofs_at_node (dim, high_p_fe_type, type, n);
4022 
4023  // since we may be running this method concurrently
4024  // on multiple threads we need to acquire a lock
4025  // before modifying the _dof_constraints object.
4026  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
4027 
4028  if (elem->is_vertex(n))
4029  {
4030  // Add "this is zero" constraint rows for high p vertex
4031  // dofs
4032  for (unsigned int i = low_nc; i != high_nc; ++i)
4033  {
4034  _dof_constraints[node.dof_number(sys_num,var,i)].clear();
4035  _primal_constraint_values.erase(node.dof_number(sys_num,var,i));
4036  }
4037  }
4038  else
4039  {
4040  const unsigned int total_dofs = node.n_comp(sys_num, var);
4041  libmesh_assert_greater_equal (total_dofs, high_nc);
4042  // Add "this is zero" constraint rows for high p
4043  // non-vertex dofs, which are numbered in reverse
4044  for (unsigned int j = low_nc; j != high_nc; ++j)
4045  {
4046  const unsigned int i = total_dofs - j - 1;
4047  _dof_constraints[node.dof_number(sys_num,var,i)].clear();
4048  _primal_constraint_values.erase(node.dof_number(sys_num,var,i));
4049  }
4050  }
4051  }
4052 }
4053 
4054 #endif // LIBMESH_ENABLE_AMR
4055 
4056 
4057 #ifdef LIBMESH_ENABLE_DIRICHLET
4058 void DofMap::add_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary)
4059 {
4060  _dirichlet_boundaries->push_back(new DirichletBoundary(dirichlet_boundary));
4061 }
4062 
4063 
4065  unsigned int qoi_index)
4066 {
4067  unsigned int old_size = cast_int<unsigned int>
4068  (_adjoint_dirichlet_boundaries.size());
4069  for (unsigned int i = old_size; i <= qoi_index; ++i)
4070  _adjoint_dirichlet_boundaries.push_back(new DirichletBoundaries());
4071 
4072  _adjoint_dirichlet_boundaries[qoi_index]->push_back
4073  (new DirichletBoundary(dirichlet_boundary));
4074 }
4075 
4076 
4078 {
4079  if (_adjoint_dirichlet_boundaries.size() > q)
4080  return true;
4081 
4082  return false;
4083 }
4084 
4085 
4086 const DirichletBoundaries *
4088 {
4089  libmesh_assert_greater(_adjoint_dirichlet_boundaries.size(),q);
4090  return _adjoint_dirichlet_boundaries[q];
4091 }
4092 
4093 
4096 {
4097  unsigned int old_size = cast_int<unsigned int>
4098  (_adjoint_dirichlet_boundaries.size());
4099  for (unsigned int i = old_size; i <= q; ++i)
4100  _adjoint_dirichlet_boundaries.push_back(new DirichletBoundaries());
4101 
4102  return _adjoint_dirichlet_boundaries[q];
4103 }
4104 
4105 
4106 void DofMap::remove_dirichlet_boundary (const DirichletBoundary & boundary_to_remove)
4107 {
4108  // Find a boundary condition matching the one to be removed
4109  std::vector<DirichletBoundary *>::iterator it = _dirichlet_boundaries->begin();
4110  std::vector<DirichletBoundary *>::iterator end = _dirichlet_boundaries->end();
4111  for (; it != end; ++it)
4112  {
4113  DirichletBoundary * bdy = *it;
4114 
4115  if ((bdy->b == boundary_to_remove.b) &&
4116  bdy->variables == boundary_to_remove.variables)
4117  break;
4118  }
4119 
4120  // Delete it and remove it
4121  libmesh_assert (it != end);
4122  delete *it;
4123  _dirichlet_boundaries->erase(it);
4124 }
4125 
4126 
4128  unsigned int qoi_index)
4129 {
4130  libmesh_assert_greater(_adjoint_dirichlet_boundaries.size(),
4131  qoi_index);
4132 
4133  // Find a boundary condition matching the one to be removed
4134  std::vector<DirichletBoundary *>::iterator it =
4135  _adjoint_dirichlet_boundaries[qoi_index]->begin();
4136  std::vector<DirichletBoundary *>::iterator end =
4137  _adjoint_dirichlet_boundaries[qoi_index]->end();
4138  for (; it != end; ++it)
4139  {
4140  DirichletBoundary * bdy = *it;
4141 
4142  if ((bdy->b == boundary_to_remove.b) &&
4143  bdy->variables == boundary_to_remove.variables)
4144  break;
4145  }
4146 
4147  // Delete it and remove it
4148  libmesh_assert (it != end);
4149  delete *it;
4150  _adjoint_dirichlet_boundaries[qoi_index]->erase(it);
4151 }
4152 
4153 
4155 {
4156  for (std::vector<DirichletBoundary *>::iterator it = begin(); it != end(); ++it)
4157  delete *it;
4158 }
4159 
4161  const DirichletBoundary & boundary) const
4162 {
4163  const std::set<boundary_id_type>& mesh_bcids = mesh.get_boundary_info().get_boundary_ids();
4164  const std::set<boundary_id_type>& dbc_bcids = boundary.b;
4165 
4166  // DirichletBoundary id sets should be consistent across all ranks
4167  libmesh_assert(mesh.comm().verify(dbc_bcids.size()));
4168 
4169  for (std::set<boundary_id_type>::const_iterator bid = dbc_bcids.begin();
4170  bid != dbc_bcids.end(); ++bid)
4171  {
4172  // DirichletBoundary id sets should be consistent across all ranks
4173  libmesh_assert(mesh.comm().verify(*bid));
4174 
4175  bool found_bcid = (mesh_bcids.find(*bid) != mesh_bcids.end());
4176 
4177  // On a distributed mesh, boundary id sets may *not* be
4178  // consistent across all ranks, since not all ranks see all
4179  // boundaries
4180  mesh.comm().max(found_bcid);
4181 
4182  if (!found_bcid)
4183  libmesh_error_msg("Could not find Dirichlet boundary id " << *bid << " in mesh!");
4184  }
4185 }
4186 
4187 #endif // LIBMESH_ENABLE_DIRICHLET
4188 
4189 
4190 #ifdef LIBMESH_ENABLE_PERIODIC
4191 
4192 void DofMap::add_periodic_boundary (const PeriodicBoundaryBase & periodic_boundary)
4193 {
4194  // See if we already have a periodic boundary associated myboundary...
4195  PeriodicBoundaryBase * existing_boundary = _periodic_boundaries->boundary(periodic_boundary.myboundary);
4196 
4197  if (existing_boundary == libmesh_nullptr)
4198  {
4199  // ...if not, clone the input (and its inverse) and add them to the PeriodicBoundaries object
4200  PeriodicBoundaryBase * boundary = periodic_boundary.clone().release();
4201  PeriodicBoundaryBase * inverse_boundary = periodic_boundary.clone(PeriodicBoundaryBase::INVERSE).release();
4202 
4203  // _periodic_boundaries takes ownership of the pointers
4204  _periodic_boundaries->insert(std::make_pair(boundary->myboundary, boundary));
4205  _periodic_boundaries->insert(std::make_pair(inverse_boundary->myboundary, inverse_boundary));
4206  }
4207  else
4208  {
4209  // ...otherwise, merge this object's variable IDs with the existing boundary object's.
4210  existing_boundary->merge(periodic_boundary);
4211 
4212  // Do the same merging process for the inverse boundary. Note: the inverse better already exist!
4213  PeriodicBoundaryBase * inverse_boundary = _periodic_boundaries->boundary(periodic_boundary.pairedboundary);
4214  libmesh_assert(inverse_boundary);
4215  inverse_boundary->merge(periodic_boundary);
4216  }
4217 }
4218 
4219 
4220 
4221 
4223  const PeriodicBoundaryBase & inverse_boundary)
4224 {
4225  libmesh_assert_equal_to (boundary.myboundary, inverse_boundary.pairedboundary);
4226  libmesh_assert_equal_to (boundary.pairedboundary, inverse_boundary.myboundary);
4227 
4228  // Allocate copies on the heap. The _periodic_boundaries object will manage this memory.
4229  // Note: this also means that the copy constructor for the PeriodicBoundary (or user class
4230  // derived therefrom) must be implemented!
4231  // PeriodicBoundary * p_boundary = new PeriodicBoundary(boundary);
4232  // PeriodicBoundary * p_inverse_boundary = new PeriodicBoundary(inverse_boundary);
4233 
4234  // We can't use normal copy construction since this leads to slicing with derived classes.
4235  // Use clone() "virtual constructor" instead. But, this *requires* user to override the clone()
4236  // method. Note also that clone() allocates memory. In this case, the _periodic_boundaries object
4237  // takes responsibility for cleanup.
4238  PeriodicBoundaryBase * p_boundary = boundary.clone().release();
4239  PeriodicBoundaryBase * p_inverse_boundary = inverse_boundary.clone().release();
4240 
4241  // Add the periodic boundary and its inverse to the PeriodicBoundaries data structure. The
4242  // PeriodicBoundaries data structure takes ownership of the pointers.
4243  _periodic_boundaries->insert(std::make_pair(p_boundary->myboundary, p_boundary));
4244  _periodic_boundaries->insert(std::make_pair(p_inverse_boundary->myboundary, p_inverse_boundary));
4245 }
4246 
4247 
4248 #endif
4249 
4250 
4251 } // 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...
The definition of the element_iterator struct.
Definition: mesh_base.h:1476
static void dofs_on_edge(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int e, std::vector< unsigned int > &di)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem A...
Definition: fe_interface.C:510
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
FEFamily family
The type of finite element.
Definition: fe_type.h:203
virtual bool closed() const
bool verify(const T &r) const
Verify that a local variable has the same value on all processors.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
double abs(double a)
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...
const FEType & type() const
Definition: variable.h:119
virtual Output component(unsigned int i, const Point &p, Real time=0.)
A Node is like a Point, but with more information.
Definition: node.h:52
virtual bool is_serial() const
Definition: mesh_base.h:140
bool active() const
Definition: elem.h:2257
unsigned int n() const
virtual numeric_index_type size() const =0
virtual numeric_index_type last_local_index() const =0
void scatter_constraints(MeshBase &)
Sends constraint equations to constraining processors.
virtual ElemType type() const =0
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:796
virtual void zero() libmesh_override
Set every element in the matrix to 0.
Definition: dense_matrix.h:792
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 ...
virtual bool is_edge_on_side(const unsigned int e, const unsigned int s) const =0
unsigned int n_vars(const unsigned int s, const unsigned int vg) const
Definition: dof_object.h:745
unsigned int p_level() const
Definition: elem.h:2422
void max(T &r) const
Take a local variable and replace it with the maximum of it&#39;s values on all processors.
dof_id_type n_constrained_dofs() const
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
The definition of the const_element_iterator struct.
Definition: mesh_base.h:1494
static void dofs_on_side(const Elem *const elem, const unsigned int dim, const FEType &fe_t, unsigned int s, std::vector< unsigned int > &di)
Fills the vector di with the local degree of freedom indices associated with side s of element elem A...
Definition: fe_interface.C:495
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 ...
unsigned int dim
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
virtual unsigned int n_edges() const =0
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.
void vector_mult_transpose(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this)^T * arg.
Definition: dense_matrix.C:438
virtual element_iterator local_elements_begin()=0
ElemType
Defines an enum for geometric element types.
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:350
void remove_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Removes the specified Dirichlet boundary from the system.
virtual void zero() libmesh_override
Set every element in the vector to 0.
Definition: dense_vector.h:374
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:1802
UniquePtr< FunctionBase< Gradient > > g
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
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.
uint8_t processor_id_type
Definition: id_types.h:99
UniquePtr< NumericVector< Number > > current_local_solution
All the values I need to compute my contribution to the simulation at hand.
Definition: system.h:1535
static FEFieldType field_type(const FEType &fe_type)
const class libmesh_nullptr_t libmesh_nullptr
dof_id_type dof_number(const unsigned int s, const unsigned int var, const unsigned int comp) const
Definition: dof_object.h:810
This class allows one to associate Dirichlet boundary values with a given set of mesh boundary ids an...
virtual const Node * node_ptr(const dof_id_type i) const =0
The Node constraint storage format.
Definition: dof_map.h:144
std::vector< unsigned int > variables
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:197
static const Real TOLERANCE
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
std::pair< Real, Real > max_constraint_error(const System &system, NumericVector< Number > *v=libmesh_nullptr) const
Tests the constrained degrees of freedom on the numeric vector v, which represents a solution defined...
The StoredRange class defines a contiguous, divisible set of objects.
Definition: stored_range.h:52
const Variable & variable(const unsigned int c) const
Definition: dof_map.h:1667
const unsigned int n_vars
Definition: tecplot_io.C:68
The libMesh namespace provides an interface to certain functionality in the library.
dof_id_type dof_id
Definition: xdr_io.C:48
const_iterator begin() const
Beginning of the range.
Definition: stored_range.h:261
virtual unsigned int size() const libmesh_override
Definition: dense_vector.h:87
long double max(long double a, double b)
Real distance(const Point &p)
Number exact_value(const Point &p, const Parameters &parameters, const std::string &, const std::string &)
unsigned int first_scalar_number() const
Definition: variable.h:113
void iota(ForwardIter first, ForwardIter last, T value)
Utility::iota is a duplication of the SGI STL extension std::iota.
Definition: utility.h:57
const std::set< boundary_id_type > & get_boundary_ids() const
This is the MeshBase class.
Definition: mesh_base.h:68
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
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...
Definition: fe_interface.C:863
bool has_adjoint_dirichlet_boundaries(unsigned int q) const
void merge(const PeriodicBoundaryBase &pb)
virtual unsigned int n_nodes() const =0
std::vector< boundary_id_type > boundary_ids(const Node *node) const
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
const DirichletBoundaries * get_adjoint_dirichlet_boundaries(unsigned int q) const
boundary_id_type myboundary
The boundary ID of this boundary and its counterpart.
const dof_id_type n_nodes
Definition: tecplot_io.C:67
const MeshBase & get_mesh() const
Definition: system.h:2014
void print_dof_constraints(std::ostream &os=libMesh::out, bool print_nonlocal=false) const
Prints (from processor 0) all DoF and Node constraints.
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:29
This class defines the notion of a variable in the system.
Definition: variable.h:49
const DofMap & get_dof_map() const
Definition: system.h:2030
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
virtual void right_multiply(const DenseMatrixBase< T > &M2) libmesh_override
Performs the operation: (*this) <- (*this) * M3.
Definition: dense_matrix.C:210
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1874
std::string get_local_constraints(bool print_nonlocal=false) const
Gets a string reporting all DoF and Node constraints local to this processor.
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 m() const
virtual element_iterator active_local_elements_begin()=0
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. ...
Definition: dense_matrix.C:508
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
dof_id_type n_local_constrained_dofs() const
virtual void init_context(const FEMContext &)
Prepares a context object for use.
We&#39;re using a class instead of a typedef to allow forward declarations and future flexibility...
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:56
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
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:1806
std::unordered_map< const Elem *, const CouplingMatrix * > map_type
What elements do we care about and what variables do we care about on each element?
virtual element_iterator local_elements_end()=0
const Node & node_ref(const unsigned int i) const
Definition: elem.h:1896
void libmesh_assert_valid_boundary_ids(const MeshBase &mesh)
A function for verifying that boundary condition ids match across processors.
Definition: mesh_tools.C:1194
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:61
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...
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
subdomain_id_type subdomain_id() const
Definition: elem.h:1951
virtual unsigned int n_sides() const =0
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 enforce_constraints_exactly(const System &system, NumericVector< Number > *v=libmesh_nullptr, bool homogeneous=false) const
Constrains the numeric vector v, which represents a solution defined on the mesh. ...
Definition: dof_map.h:1816
bool is_prepared() const
Definition: mesh_base.h:133
const_iterator end() const
End of the range.
Definition: stored_range.h:266
virtual UniquePtr< 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 ...
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
A class for templated methods that expect output iterator arguments, which adds objects to the Mesh...
Storage for DofConstraint right hand sides for a particular problem.
Definition: dof_map.h:108
std::vector< object_type >::const_iterator const_iterator
Allows an StoredRange to behave like an STL container.
Definition: stored_range.h:58
void heterogenously_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.
void libmesh_ignore(const T &)
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
UniquePtr< FunctionBase< Number > > f
virtual element_iterator active_not_local_elements_end()=0
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.
Definition: dense_matrix.C:85
const Point & point(const unsigned int i) const
Definition: elem.h:1809
void process_constraints(MeshBase &)
Postprocesses any constrained degrees of freedom to be constrained only in terms of unconstrained dof...
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:1793
bool verify(const T &r, const Communicator &comm=Communicator_World)
StoredRange< iterator_type, object_type > & reset(const iterator_type &first, const iterator_type &last)
Resets the StoredRange to contain [first,last).
Definition: stored_range.h:222
virtual bool is_vertex(const unsigned int i) const =0
bool libmesh_isnan(float a)
void add_constraints_to_send_list()
Adds entries to the _send_list vector corresponding to DoFs which are dependencies for constraint equ...
UniquePtr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:30
const Parallel::Communicator & comm() const
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
Access multiple components at once.
std::vector< boundary_id_type > edge_boundary_ids(const Elem *const elem, const unsigned short int edge) const
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
virtual Output component(const FEMContext &, unsigned int i, const Point &p, Real time=0.)
ParallelType type() const
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
std::map< dof_id_type, Real, std::less< dof_id_type >, Threads::scalable_allocator< std::pair< const dof_id_type, Real > > > DofConstraintRow
A row of the Dof constraint matrix.
Definition: dof_map.h:88
virtual unsigned int dim() const =0
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
void heterogenously_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.
UniquePtr< FEMFunctionBase< Gradient > > g_fem
virtual numeric_index_type first_local_index() const =0
void add_dirichlet_boundary(const DirichletBoundary &dirichlet_boundary)
Adds a copy of the specified Dirichlet boundary to the system.
virtual element_iterator active_not_local_elements_begin()=0
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.
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:136
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
Real size() const
Definition: type_vector.h:898
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1831
unsigned int n_comp(const unsigned int s, const unsigned int var) const
Definition: dof_object.h:780
dof_id_type id() const
Definition: dof_object.h:632
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this) * arg.
Definition: dense_matrix.C:382
virtual dof_id_type n_elem() const =0
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:1811
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...
virtual element_iterator active_local_elements_end()=0
UniquePtr< FEMFunctionBase< Number > > f_fem
bool empty() const
Definition: stored_range.h:301
processor_id_type processor_id()
Definition: libmesh_base.h:96
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
For symmetric positive definite (SPD) matrices.
Definition: dense_matrix.C:911
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
virtual void localize(std::vector< T > &v_local) const =0
Creates a copy of the global vector in the local vector v_local.
The constraint matrix storage format.
Definition: dof_map.h:96
UniquePtr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:534
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:939
void enforce_adjoint_constraints_exactly(NumericVector< Number > &v, unsigned int q) const
Heterogenously constrains the numeric vector v, which represents an adjoint solution defined on the m...
Definition: dof_map.h:1820
This class provides single index access to FieldType (i.e.
Definition: raw_accessor.h:93
void constrain_p_dofs(unsigned int var, const Elem *elem, unsigned int s, unsigned int p)
Constrains degrees of freedom on side s of element elem which correspond to variable number var and t...
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
Definition: dof_object.h:694
processor_id_type n_processors()
Definition: libmesh_base.h:88
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.