dof_map_constraints.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // Local Includes
19 #include "libmesh/boundary_info.h" // needed for dirichlet constraints
20 #include "libmesh/dense_matrix.h"
21 #include "libmesh/dense_vector.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"
31 #include "libmesh/mesh_base.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"
37 #include "libmesh/parallel_elem.h"
38 #include "libmesh/parallel_node.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  std::unique_ptr<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
99  FEInterface::compute_constraints (_constraints,
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)
109  FEInterface::compute_periodic_constraints (_constraints,
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  std::unique_ptr<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)
173  FEBase::compute_periodic_node_constraints (_node_constraints,
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
348  std::unique_ptr<FEGenericBase<OutputType>> fe = FEGenericBase<OutputType>::build(dim, fe_type);
349 
350  // Prepare variables for projection
351  std::unique_ptr<QBase> qedgerule (fe_type.default_quadrature_rule(1));
352  std::unique_ptr<QBase> qsiderule (fe_type.default_quadrature_rule(dim-1));
353  std::unique_ptr<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  std::unique_ptr<FEMContext> context;
395  if (f_fem)
396  {
397  libmesh_assert(f_system);
398  if (f_system->current_local_solution->initialized())
399  {
400  context = libmesh_make_unique<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]]);
821  libmesh_assert(std::abs(ui) < TOLERANCE ||
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]]);
941  libmesh_assert(std::abs(ui) < TOLERANCE ||
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]]);
1064  libmesh_assert(std::abs(ui) < TOLERANCE ||
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 
1163 dof_id_type DofMap::n_constrained_dofs() const
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 
1173 dof_id_type DofMap::n_local_constrained_dofs() const
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 
1185 void DofMap::create_dof_constraints(const MeshBase & mesh, Real time)
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 
1342 void DofMap::add_constraint_row (const dof_id_type dof_number,
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 
1555 void DofMap::constrain_element_matrix (DenseMatrix<Number> & matrix,
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 
1625 void DofMap::constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
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 
1704 void DofMap::heterogenously_constrain_element_matrix_and_vector (DenseMatrix<Number> & matrix,
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 
1810 void DofMap::heterogenously_constrain_element_vector (const DenseMatrix<Number> & matrix,
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 
1886 void DofMap::constrain_element_matrix (DenseMatrix<Number> & matrix,
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 
1976 void DofMap::constrain_element_vector (DenseVector<Number> & rhs,
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 
2016 void DofMap::constrain_element_dyad_matrix (DenseVector<Number> & v,
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
2095  std::unique_ptr<NumericVector<Number>> v_built;
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 
2178 void DofMap::enforce_adjoint_constraints_exactly (NumericVector<Number> & v,
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
2190  std::unique_ptr<NumericVector<Number>> v_built;
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>
2282 DofMap::max_constraint_error (const System & system,
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 
2360 void DofMap::build_constraint_matrix (DenseMatrix<Number> & C,
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 
2472 void DofMap::build_constraint_matrix_and_vector (DenseMatrix<Number> & C,
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 
2613 void DofMap::allgather_recursive_constraints(MeshBase & mesh)
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 
3134 void DofMap::process_constraints (MeshBase & mesh)
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  if(_error_on_cyclic_constraint)
3141  {
3142  // Optionally check for cyclic constraints and throw an error
3143  // if they're detected. We always do this check below in dbg/devel
3144  // mode but here we optionally do it in opt mode as well.
3145  check_for_cyclic_constraints();
3146  }
3147 
3148  // Create a set containing the DOFs we already depend on
3149  typedef std::set<dof_id_type> RCSet;
3150  RCSet unexpanded_set;
3151 
3152  for (DofConstraints::iterator i = _dof_constraints.begin();
3153  i != _dof_constraints.end(); ++i)
3154  unexpanded_set.insert(i->first);
3155 
3156  while (!unexpanded_set.empty())
3157  for (RCSet::iterator i = unexpanded_set.begin();
3158  i != unexpanded_set.end(); /* nothing */)
3159  {
3160  // If the DOF is constrained
3161  DofConstraints::iterator
3162  pos = _dof_constraints.find(*i);
3163 
3164  libmesh_assert (pos != _dof_constraints.end());
3165 
3166  DofConstraintRow & constraint_row = pos->second;
3167 
3168  DofConstraintValueMap::iterator rhsit =
3169  _primal_constraint_values.find(*i);
3170  Number constraint_rhs = (rhsit == _primal_constraint_values.end()) ?
3171  0 : rhsit->second;
3172 
3173  std::vector<dof_id_type> constraints_to_expand;
3174 
3175  for (DofConstraintRow::const_iterator
3176  it=constraint_row.begin(); it != constraint_row.end();
3177  ++it)
3178  if (it->first != *i && this->is_constrained_dof(it->first))
3179  {
3180  unexpanded_set.insert(it->first);
3181  constraints_to_expand.push_back(it->first);
3182  }
3183 
3184  for (std::size_t j = 0; j != constraints_to_expand.size();
3185  ++j)
3186  {
3187  dof_id_type expandable = constraints_to_expand[j];
3188 
3189  const Real this_coef = constraint_row[expandable];
3190 
3191  DofConstraints::const_iterator
3192  subpos = _dof_constraints.find(expandable);
3193 
3194  libmesh_assert (subpos != _dof_constraints.end());
3195 
3196  const DofConstraintRow & subconstraint_row = subpos->second;
3197 
3198  for (DofConstraintRow::const_iterator
3199  it=subconstraint_row.begin();
3200  it != subconstraint_row.end(); ++it)
3201  {
3202  // Assert that the constraint is not cyclic.
3203  libmesh_assert(it->first != expandable);
3204 
3205  constraint_row[it->first] += it->second * this_coef;
3206  }
3207  DofConstraintValueMap::const_iterator subrhsit =
3208  _primal_constraint_values.find(expandable);
3209  if (subrhsit != _primal_constraint_values.end())
3210  constraint_rhs += subrhsit->second * this_coef;
3211 
3212  constraint_row.erase(expandable);
3213  }
3214 
3215  if (rhsit == _primal_constraint_values.end())
3216  {
3217  if (constraint_rhs != Number(0))
3218  _primal_constraint_values[*i] = constraint_rhs;
3219  else
3220  _primal_constraint_values.erase(*i);
3221  }
3222  else
3223  {
3224  if (constraint_rhs != Number(0))
3225  rhsit->second = constraint_rhs;
3226  else
3227  _primal_constraint_values.erase(rhsit);
3228  }
3229 
3230  if (constraints_to_expand.empty())
3231  unexpanded_set.erase(i++);
3232  else
3233  ++i;
3234  }
3235 
3236  // In parallel we can't guarantee that nodes/dofs which constrain
3237  // others are on processors which are aware of that constraint, yet
3238  // we need such awareness for sparsity pattern generation. So send
3239  // other processors any constraints they might need to know about.
3240  this->scatter_constraints(mesh);
3241 
3242  // Now that we have our root constraint dependencies sorted out, add
3243  // them to the send_list
3244  this->add_constraints_to_send_list();
3245 }
3246 
3247 
3248 #ifdef LIBMESH_ENABLE_CONSTRAINTS
3249 void DofMap::check_for_cyclic_constraints()
3250 {
3251  // Create a set containing the DOFs we already depend on
3252  typedef std::set<dof_id_type> RCSet;
3253  RCSet unexpanded_set;
3254 
3255  // Use dof_constraints_copy in this method so that we don't
3256  // mess with _dof_constraints.
3257  DofConstraints dof_constraints_copy = _dof_constraints;
3258 
3259  for (DofConstraints::iterator i = dof_constraints_copy.begin();
3260  i != dof_constraints_copy.end(); ++i)
3261  unexpanded_set.insert(i->first);
3262 
3263  while (!unexpanded_set.empty())
3264  for (RCSet::iterator i = unexpanded_set.begin();
3265  i != unexpanded_set.end(); /* nothing */)
3266  {
3267  // If the DOF is constrained
3268  DofConstraints::iterator
3269  pos = dof_constraints_copy.find(*i);
3270 
3271  libmesh_assert (pos != dof_constraints_copy.end());
3272 
3273  DofConstraintRow & constraint_row = pos->second;
3274 
3275  // Comment out "rhs" parts of this method copied from process_constraints
3276  // DofConstraintValueMap::iterator rhsit =
3277  // _primal_constraint_values.find(*i);
3278  // Number constraint_rhs = (rhsit == _primal_constraint_values.end()) ?
3279  // 0 : rhsit->second;
3280 
3281  std::vector<dof_id_type> constraints_to_expand;
3282 
3283  for (DofConstraintRow::const_iterator
3284  it=constraint_row.begin(); it != constraint_row.end();
3285  ++it)
3286  if (it->first != *i && this->is_constrained_dof(it->first))
3287  {
3288  unexpanded_set.insert(it->first);
3289  constraints_to_expand.push_back(it->first);
3290  }
3291 
3292  for (std::size_t j = 0; j != constraints_to_expand.size();
3293  ++j)
3294  {
3295  dof_id_type expandable = constraints_to_expand[j];
3296 
3297  const Real this_coef = constraint_row[expandable];
3298 
3299  DofConstraints::const_iterator
3300  subpos = dof_constraints_copy.find(expandable);
3301 
3302  libmesh_assert (subpos != dof_constraints_copy.end());
3303 
3304  const DofConstraintRow & subconstraint_row = subpos->second;
3305 
3306  for (DofConstraintRow::const_iterator
3307  it=subconstraint_row.begin();
3308  it != subconstraint_row.end(); ++it)
3309  {
3310 
3311  if(it->first == expandable)
3312  {
3313  libmesh_error_msg("Cyclic constraint detected");
3314  }
3315 
3316  constraint_row[it->first] += it->second * this_coef;
3317  }
3318 
3319  // Comment out "rhs" parts of this method copied from process_constraints
3320  // DofConstraintValueMap::const_iterator subrhsit =
3321  // _primal_constraint_values.find(expandable);
3322  // if (subrhsit != _primal_constraint_values.end())
3323  // constraint_rhs += subrhsit->second * this_coef;
3324 
3325  constraint_row.erase(expandable);
3326  }
3327 
3328  // Comment out "rhs" parts of this method copied from process_constraints
3329  // if (rhsit == _primal_constraint_values.end())
3330  // {
3331  // if (constraint_rhs != Number(0))
3332  // _primal_constraint_values[*i] = constraint_rhs;
3333  // else
3334  // _primal_constraint_values.erase(*i);
3335  // }
3336  // else
3337  // {
3338  // if (constraint_rhs != Number(0))
3339  // rhsit->second = constraint_rhs;
3340  // else
3341  // _primal_constraint_values.erase(rhsit);
3342  // }
3343 
3344  if (constraints_to_expand.empty())
3345  unexpanded_set.erase(i++);
3346  else
3347  ++i;
3348  }
3349 }
3350 #else
3351 void DofMap::check_for_cyclic_constraints()
3352 {
3353  // Do nothing
3354 }
3355 #endif
3356 
3357 
3358 void DofMap::scatter_constraints(MeshBase & mesh)
3359 {
3360  // At this point each processor with a constrained node knows
3361  // the corresponding constraint row, but we also need each processor
3362  // with a constrainer node to know the corresponding row(s).
3363 
3364  // This function must be run on all processors at once
3365  parallel_object_only();
3366 
3367  // Return immediately if there's nothing to gather
3368  if (this->n_processors() == 1)
3369  return;
3370 
3371  // We might get to return immediately if none of the processors
3372  // found any constraints
3373  unsigned int has_constraints = !_dof_constraints.empty()
3374 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3375  || !_node_constraints.empty()
3376 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3377  ;
3378  this->comm().max(has_constraints);
3379  if (!has_constraints)
3380  return;
3381 
3382 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3383  std::vector<std::set<dof_id_type>> pushed_node_ids(this->n_processors());
3384 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3385 
3386  std::vector<std::set<dof_id_type>> pushed_ids(this->n_processors());
3387 
3388  // Collect the dof constraints I need to push to each processor
3389  dof_id_type constrained_proc_id = 0;
3390  for (DofConstraints::iterator i = _dof_constraints.begin();
3391  i != _dof_constraints.end(); ++i)
3392  {
3393  const dof_id_type constrained = i->first;
3394  while (constrained >= _end_df[constrained_proc_id])
3395  constrained_proc_id++;
3396 
3397  if (constrained_proc_id != this->processor_id())
3398  continue;
3399 
3400  DofConstraintRow & row = i->second;
3401  for (DofConstraintRow::iterator j = row.begin();
3402  j != row.end(); ++j)
3403  {
3404  const dof_id_type constraining = j->first;
3405 
3406  processor_id_type constraining_proc_id = 0;
3407  while (constraining >= _end_df[constraining_proc_id])
3408  constraining_proc_id++;
3409 
3410  if (constraining_proc_id != this->processor_id() &&
3411  constraining_proc_id != constrained_proc_id)
3412  pushed_ids[constraining_proc_id].insert(constrained);
3413  }
3414  }
3415 
3416 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3417  // Collect the node constraints to push to each processor
3418  for (NodeConstraints::iterator i = _node_constraints.begin();
3419  i != _node_constraints.end(); ++i)
3420  {
3421  const Node * constrained = i->first;
3422 
3423  if (constrained->processor_id() != this->processor_id())
3424  continue;
3425 
3426  NodeConstraintRow & row = i->second.first;
3427  for (NodeConstraintRow::iterator j = row.begin();
3428  j != row.end(); ++j)
3429  {
3430  const Node * constraining = j->first;
3431 
3432  if (constraining->processor_id() != this->processor_id() &&
3433  constraining->processor_id() != constrained->processor_id())
3434  pushed_node_ids[constraining->processor_id()].insert(constrained->id());
3435  }
3436  }
3437 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3438 
3439  // Now trade constraint rows
3440  for (processor_id_type p = 0; p != this->n_processors(); ++p)
3441  {
3442  // Push to processor procup while receiving from procdown
3443  processor_id_type procup =
3444  cast_int<processor_id_type>((this->processor_id() + p) %
3445  this->n_processors());
3446  processor_id_type procdown =
3447  cast_int<processor_id_type>((this->n_processors() +
3448  this->processor_id() - p) %
3449  this->n_processors());
3450 
3451  // Pack the dof constraint rows and rhs's to push to procup
3452  const std::size_t pushed_ids_size = pushed_ids[procup].size();
3453  std::vector<std::vector<dof_id_type>> pushed_keys(pushed_ids_size);
3454  std::vector<std::vector<Real>> pushed_vals(pushed_ids_size);
3455  std::vector<Number> pushed_rhss(pushed_ids_size);
3456 
3457  std::set<dof_id_type>::const_iterator it;
3458  std::size_t push_i;
3459  for (push_i = 0, it = pushed_ids[procup].begin();
3460  it != pushed_ids[procup].end(); ++push_i, ++it)
3461  {
3462  const dof_id_type constrained = *it;
3463  DofConstraintRow & row = _dof_constraints[constrained];
3464  std::size_t row_size = row.size();
3465  pushed_keys[push_i].reserve(row_size);
3466  pushed_vals[push_i].reserve(row_size);
3467  for (DofConstraintRow::iterator j = row.begin();
3468  j != row.end(); ++j)
3469  {
3470  pushed_keys[push_i].push_back(j->first);
3471  pushed_vals[push_i].push_back(j->second);
3472  }
3473 
3474  DofConstraintValueMap::const_iterator rhsit =
3475  _primal_constraint_values.find(constrained);
3476  pushed_rhss[push_i] = (rhsit == _primal_constraint_values.end()) ?
3477  0 : rhsit->second;
3478  }
3479 
3480 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3481  // Pack the node constraint rows to push to procup
3482  const std::size_t pushed_node_ids_size = pushed_node_ids[procup].size();
3483  std::vector<std::vector<dof_id_type>> pushed_node_keys(pushed_node_ids_size);
3484  std::vector<std::vector<Real>> pushed_node_vals(pushed_node_ids_size);
3485  std::vector<Point> pushed_node_offsets(pushed_node_ids_size);
3486  std::set<const Node *> pushed_nodes;
3487 
3488  for (push_i = 0, it = pushed_node_ids[procup].begin();
3489  it != pushed_node_ids[procup].end(); ++push_i, ++it)
3490  {
3491  const Node * constrained = mesh.node_ptr(*it);
3492 
3493  if (constrained->processor_id() != procdown)
3494  pushed_nodes.insert(constrained);
3495 
3496  NodeConstraintRow & row = _node_constraints[constrained].first;
3497  std::size_t row_size = row.size();
3498  pushed_node_keys[push_i].reserve(row_size);
3499  pushed_node_vals[push_i].reserve(row_size);
3500  for (NodeConstraintRow::iterator j = row.begin();
3501  j != row.end(); ++j)
3502  {
3503  const Node * constraining = j->first;
3504 
3505  pushed_node_keys[push_i].push_back(constraining->id());
3506  pushed_node_vals[push_i].push_back(j->second);
3507 
3508  if (constraining->processor_id() != procup)
3509  pushed_nodes.insert(constraining);
3510  }
3511  pushed_node_offsets[push_i] = _node_constraints[constrained].second;
3512  }
3513 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3514 
3515  // Trade pushed dof constraint rows
3516  std::vector<dof_id_type> pushed_ids_from_me
3517  (pushed_ids[procup].begin(), pushed_ids[procup].end());
3518  std::vector<dof_id_type> pushed_ids_to_me;
3519  std::vector<std::vector<dof_id_type>> pushed_keys_to_me;
3520  std::vector<std::vector<Real>> pushed_vals_to_me;
3521  std::vector<Number> pushed_rhss_to_me;
3522  this->comm().send_receive(procup, pushed_ids_from_me,
3523  procdown, pushed_ids_to_me);
3524  this->comm().send_receive(procup, pushed_keys,
3525  procdown, pushed_keys_to_me);
3526  this->comm().send_receive(procup, pushed_vals,
3527  procdown, pushed_vals_to_me);
3528  this->comm().send_receive(procup, pushed_rhss,
3529  procdown, pushed_rhss_to_me);
3530  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_keys_to_me.size());
3531  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_vals_to_me.size());
3532  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_rhss_to_me.size());
3533 
3534 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3535  // Trade pushed node constraint rows
3536  std::vector<dof_id_type> pushed_node_ids_from_me
3537  (pushed_node_ids[procup].begin(), pushed_node_ids[procup].end());
3538  std::vector<dof_id_type> pushed_node_ids_to_me;
3539  std::vector<std::vector<dof_id_type>> pushed_node_keys_to_me;
3540  std::vector<std::vector<Real>> pushed_node_vals_to_me;
3541  std::vector<Point> pushed_node_offsets_to_me;
3542  this->comm().send_receive(procup, pushed_node_ids_from_me,
3543  procdown, pushed_node_ids_to_me);
3544  this->comm().send_receive(procup, pushed_node_keys,
3545  procdown, pushed_node_keys_to_me);
3546  this->comm().send_receive(procup, pushed_node_vals,
3547  procdown, pushed_node_vals_to_me);
3548  this->comm().send_receive(procup, pushed_node_offsets,
3549  procdown, pushed_node_offsets_to_me);
3550 
3551  // Constraining nodes might not even exist on our subset of
3552  // a distributed mesh, so let's make them exist.
3553  if (!mesh.is_serial())
3554  this->comm().send_receive_packed_range
3555  (procup, &mesh, pushed_nodes.begin(), pushed_nodes.end(),
3556  procdown, &mesh, mesh_inserter_iterator<Node>(mesh),
3557  (Node**)libmesh_nullptr);
3558 
3559  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_keys_to_me.size());
3560  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_vals_to_me.size());
3561  libmesh_assert_equal_to (pushed_node_ids_to_me.size(), pushed_node_offsets_to_me.size());
3562 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3563 
3564  // Add the dof constraints that I've been sent
3565  for (std::size_t i = 0; i != pushed_ids_to_me.size(); ++i)
3566  {
3567  libmesh_assert_equal_to (pushed_keys_to_me[i].size(), pushed_vals_to_me[i].size());
3568 
3569  dof_id_type constrained = pushed_ids_to_me[i];
3570 
3571  // If we don't already have a constraint for this dof,
3572  // add the one we were sent
3573  if (!this->is_constrained_dof(constrained))
3574  {
3575  DofConstraintRow & row = _dof_constraints[constrained];
3576  for (std::size_t j = 0; j != pushed_keys_to_me[i].size(); ++j)
3577  {
3578  row[pushed_keys_to_me[i][j]] = pushed_vals_to_me[i][j];
3579  }
3580  if (pushed_rhss_to_me[i] != Number(0))
3581  _primal_constraint_values[constrained] = pushed_rhss_to_me[i];
3582  else
3583  _primal_constraint_values.erase(constrained);
3584  }
3585  }
3586 
3587 #ifdef LIBMESH_ENABLE_NODE_CONSTRAINTS
3588  // Add the node constraints that I've been sent
3589  for (std::size_t i = 0; i != pushed_node_ids_to_me.size(); ++i)
3590  {
3591  libmesh_assert_equal_to (pushed_node_keys_to_me[i].size(), pushed_node_vals_to_me[i].size());
3592 
3593  dof_id_type constrained_id = pushed_node_ids_to_me[i];
3594 
3595  // If we don't already have a constraint for this node,
3596  // add the one we were sent
3597  const Node * constrained = mesh.node_ptr(constrained_id);
3598  if (!this->is_constrained_node(constrained))
3599  {
3600  NodeConstraintRow & row = _node_constraints[constrained].first;
3601  for (std::size_t j = 0; j != pushed_node_keys_to_me[i].size(); ++j)
3602  {
3603  const Node * key_node = mesh.node_ptr(pushed_node_keys_to_me[i][j]);
3604  row[key_node] = pushed_node_vals_to_me[i][j];
3605  }
3606  _node_constraints[constrained].second = pushed_node_offsets_to_me[i];
3607  }
3608  }
3609 #endif // LIBMESH_ENABLE_NODE_CONSTRAINTS
3610  }
3611 
3612  // Next we need to push constraints to processors which don't own
3613  // the constrained dof, don't own the constraining dof, but own an
3614  // element supporting the constraining dof.
3615  //
3616  // We need to be able to quickly look up constrained dof ids by what
3617  // constrains them, so that we can handle the case where we see a
3618  // foreign element containing one of our constraining DoF ids and we
3619  // need to push that constraint.
3620  //
3621  // Getting distributed adaptive sparsity patterns right is hard.
3622 
3623  typedef std::map<dof_id_type, std::set<dof_id_type>> DofConstrainsMap;
3624  DofConstrainsMap dof_id_constrains;
3625 
3626  for (DofConstraints::iterator i = _dof_constraints.begin();
3627  i != _dof_constraints.end(); ++i)
3628  {
3629  const dof_id_type constrained = i->first;
3630  DofConstraintRow & row = i->second;
3631  for (DofConstraintRow::iterator j = row.begin();
3632  j != row.end(); ++j)
3633  {
3634  const dof_id_type constraining = j->first;
3635 
3636  dof_id_type constraining_proc_id = 0;
3637  while (constraining >= _end_df[constraining_proc_id])
3638  constraining_proc_id++;
3639 
3640  if (constraining_proc_id == this->processor_id())
3641  dof_id_constrains[constraining].insert(constrained);
3642  }
3643  }
3644 
3645  // Loop over all foreign elements, find any supporting our
3646  // constrained dof indices.
3647  pushed_ids.clear();
3648  pushed_ids.resize(this->n_processors());
3649 
3652  for (; it != end; ++it)
3653  {
3654  const Elem * elem = *it;
3655 
3656  std::vector<dof_id_type> my_dof_indices;
3657  this->dof_indices (elem, my_dof_indices);
3658 
3659  for (std::size_t i=0; i != my_dof_indices.size(); ++i)
3660  {
3661  DofConstrainsMap::const_iterator dcmi =
3662  dof_id_constrains.find(my_dof_indices[i]);
3663  if (dcmi != dof_id_constrains.end())
3664  {
3665  for (DofConstrainsMap::mapped_type::const_iterator mti =
3666  dcmi->second.begin();
3667  mti != dcmi->second.end(); ++mti)
3668  {
3669  const dof_id_type constrained = *mti;
3670 
3671  dof_id_type the_constrained_proc_id = 0;
3672  while (constrained >= _end_df[the_constrained_proc_id])
3673  the_constrained_proc_id++;
3674 
3675  const dof_id_type elemproc = elem->processor_id();
3676  if (elemproc != the_constrained_proc_id)
3677  pushed_ids[elemproc].insert(constrained);
3678  }
3679  }
3680  }
3681  }
3682 
3683  // One last trade of constraint rows
3684  for (processor_id_type p = 0; p != this->n_processors(); ++p)
3685  {
3686  // Push to processor procup while receiving from procdown
3687  processor_id_type procup =
3688  cast_int<processor_id_type>((this->processor_id() + p) %
3689  this->n_processors());
3690  processor_id_type procdown =
3691  cast_int<processor_id_type>((this->n_processors() +
3692  this->processor_id() - p) %
3693  this->n_processors());
3694 
3695  // Pack the dof constraint rows and rhs's to push to procup
3696  const std::size_t pushed_ids_size = pushed_ids[procup].size();
3697  std::vector<std::vector<dof_id_type>> pushed_keys(pushed_ids_size);
3698  std::vector<std::vector<Real>> pushed_vals(pushed_ids_size);
3699  std::vector<Number> pushed_rhss(pushed_ids_size);
3700 
3701  // As long as we're declaring them outside the loop, let's initialize them too!
3702  std::set<dof_id_type>::const_iterator pushed_ids_iter = pushed_ids[procup].begin();
3703  std::size_t push_i = 0;
3704  for ( ; pushed_ids_iter != pushed_ids[procup].end(); ++push_i, ++pushed_ids_iter)
3705  {
3706  const dof_id_type constrained = *pushed_ids_iter;
3707  DofConstraintRow & row = _dof_constraints[constrained];
3708  std::size_t row_size = row.size();
3709  pushed_keys[push_i].reserve(row_size);
3710  pushed_vals[push_i].reserve(row_size);
3711  for (DofConstraintRow::iterator j = row.begin();
3712  j != row.end(); ++j)
3713  {
3714  pushed_keys[push_i].push_back(j->first);
3715  pushed_vals[push_i].push_back(j->second);
3716  }
3717 
3718  DofConstraintValueMap::const_iterator rhsit =
3719  _primal_constraint_values.find(constrained);
3720  pushed_rhss[push_i] = (rhsit == _primal_constraint_values.end()) ?
3721  0 : rhsit->second;
3722  }
3723 
3724  // Trade pushed dof constraint rows
3725  std::vector<dof_id_type> pushed_ids_from_me
3726  (pushed_ids[procup].begin(), pushed_ids[procup].end());
3727  std::vector<dof_id_type> pushed_ids_to_me;
3728  std::vector<std::vector<dof_id_type>> pushed_keys_to_me;
3729  std::vector<std::vector<Real>> pushed_vals_to_me;
3730  std::vector<Number> pushed_rhss_to_me;
3731  this->comm().send_receive(procup, pushed_ids_from_me,
3732  procdown, pushed_ids_to_me);
3733  this->comm().send_receive(procup, pushed_keys,
3734  procdown, pushed_keys_to_me);
3735  this->comm().send_receive(procup, pushed_vals,
3736  procdown, pushed_vals_to_me);
3737  this->comm().send_receive(procup, pushed_rhss,
3738  procdown, pushed_rhss_to_me);
3739  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_keys_to_me.size());
3740  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_vals_to_me.size());
3741  libmesh_assert_equal_to (pushed_ids_to_me.size(), pushed_rhss_to_me.size());
3742 
3743  // Add the dof constraints that I've been sent
3744  for (std::size_t i = 0; i != pushed_ids_to_me.size(); ++i)
3745  {
3746  libmesh_assert_equal_to (pushed_keys_to_me[i].size(), pushed_vals_to_me[i].size());
3747 
3748  dof_id_type constrained = pushed_ids_to_me[i];
3749 
3750  // If we don't already have a constraint for this dof,
3751  // add the one we were sent
3752  if (!this->is_constrained_dof(constrained))
3753  {
3754  DofConstraintRow & row = _dof_constraints[constrained];
3755  for (std::size_t j = 0; j != pushed_keys_to_me[i].size(); ++j)
3756  {
3757  row[pushed_keys_to_me[i][j]] = pushed_vals_to_me[i][j];
3758  }
3759 
3760  if (pushed_rhss_to_me[i] != Number(0))
3761  _primal_constraint_values[constrained] = pushed_rhss_to_me[i];
3762  else
3763  _primal_constraint_values.erase(constrained);
3764  }
3765  }
3766  }
3767 
3768  // Finally, we need to handle the case of remote dof coupling. If a
3769  // processor's element is coupled to a ghost element, then the
3770  // processor needs to know about all constraints which affect the
3771  // dofs on that ghost element, so we'll have to query the ghost
3772  // element's owner.
3773 
3774  GhostingFunctor::map_type elements_to_couple;
3775 
3776  // Man, I wish we had guaranteed unique_ptr availability...
3777  std::set<CouplingMatrix*> temporary_coupling_matrices;
3778 
3779  this->merge_ghost_functor_outputs
3780  (elements_to_couple,
3781  temporary_coupling_matrices,
3782  this->coupling_functors_begin(),
3783  this->coupling_functors_end(),
3786  this->processor_id());
3787 
3788  // Each ghost-coupled element's owner should get a request for its dofs
3789  std::set<dof_id_type> requested_dofs;
3790 
3791  GhostingFunctor::map_type::iterator etg_it = elements_to_couple.begin();
3792  const GhostingFunctor::map_type::iterator etg_end = elements_to_couple.end();
3793  for (; etg_it != etg_end; ++etg_it)
3794  {
3795  const Elem *elem = etg_it->first;
3796 
3797  // FIXME - optimize for the non-fully-coupled case?
3798  std::vector<dof_id_type> element_dofs;
3799  this->dof_indices(elem, element_dofs);
3800 
3801  for (std::size_t i=0; i != element_dofs.size(); ++i)
3802  requested_dofs.insert(element_dofs[i]);
3803  }
3804 
3805  this->gather_constraints(mesh, requested_dofs, false);
3806 }
3807 
3808 
3809 void DofMap::gather_constraints (MeshBase & /*mesh*/,
3810  std::set<dof_id_type> & unexpanded_dofs,
3811  bool /*look_for_constrainees*/)
3812 {
3813  typedef std::set<dof_id_type> DoF_RCSet;
3814 
3815  // If we have heterogenous adjoint constraints we need to
3816  // communicate those too.
3817  const unsigned int max_qoi_num =
3818  _adjoint_constraint_values.empty() ?
3819  0 : _adjoint_constraint_values.rbegin()->first;
3820 
3821  // We have to keep recursing while the unexpanded set is
3822  // nonempty on *any* processor
3823  bool unexpanded_set_nonempty = !unexpanded_dofs.empty();
3824  this->comm().max(unexpanded_set_nonempty);
3825 
3826  while (unexpanded_set_nonempty)
3827  {
3828  // Let's make sure we don't lose sync in this loop.
3829  parallel_object_only();
3830 
3831  // Request sets
3832  DoF_RCSet dof_request_set;
3833 
3834  // Request sets to send to each processor
3835  std::vector<std::vector<dof_id_type>>
3836  requested_dof_ids(this->n_processors());
3837 
3838  // And the sizes of each
3839  std::vector<dof_id_type>
3840  dof_ids_on_proc(this->n_processors(), 0);
3841 
3842  // Fill (and thereby sort and uniq!) the main request sets
3843  for (DoF_RCSet::iterator i = unexpanded_dofs.begin();
3844  i != unexpanded_dofs.end(); ++i)
3845  {
3846  const dof_id_type unexpanded_dof = *i;
3847  DofConstraints::const_iterator
3848  pos = _dof_constraints.find(unexpanded_dof);
3849 
3850  // If we were asked for a DoF and we don't already have a
3851  // constraint for it, then we need to check for one.
3852  if (pos == _dof_constraints.end())
3853  {
3854  if (((unexpanded_dof < this->first_dof()) ||
3855  (unexpanded_dof >= this->end_dof())) &&
3856  !_dof_constraints.count(unexpanded_dof))
3857  dof_request_set.insert(unexpanded_dof);
3858  }
3859  // If we were asked for a DoF and we already have a
3860  // constraint for it, then we need to check if the
3861  // constraint is recursive.
3862  else
3863  {
3864  const DofConstraintRow & row = pos->second;
3865  for (DofConstraintRow::const_iterator j = row.begin();
3866  j != row.end(); ++j)
3867  {
3868  const dof_id_type constraining_dof = j->first;
3869 
3870  // If it's non-local and we haven't already got a
3871  // constraint for it, we might need to ask for one
3872  if (((constraining_dof < this->first_dof()) ||
3873  (constraining_dof >= this->end_dof())) &&
3874  !_dof_constraints.count(constraining_dof))
3875  dof_request_set.insert(constraining_dof);
3876  }
3877  }
3878  }
3879 
3880  // Clear the unexpanded constraint set; we're about to expand it
3881  unexpanded_dofs.clear();
3882 
3883  // Count requests by processor
3884  processor_id_type proc_id = 0;
3885  for (DoF_RCSet::iterator i = dof_request_set.begin();
3886  i != dof_request_set.end(); ++i)
3887  {
3888  while (*i >= _end_df[proc_id])
3889  proc_id++;
3890  dof_ids_on_proc[proc_id]++;
3891  }
3892 
3893  for (processor_id_type p = 0; p != this->n_processors(); ++p)
3894  {
3895  requested_dof_ids[p].reserve(dof_ids_on_proc[p]);
3896  }
3897 
3898  // Prepare each processor's request set
3899  proc_id = 0;
3900  for (DoF_RCSet::iterator i = dof_request_set.begin();
3901  i != dof_request_set.end(); ++i)
3902  {
3903  while (*i >= _end_df[proc_id])
3904  proc_id++;
3905  requested_dof_ids[proc_id].push_back(*i);
3906  }
3907 
3908  // Now request constraint rows from other processors
3909  for (processor_id_type p=1; p != this->n_processors(); ++p)
3910  {
3911  // Trade my requests with processor procup and procdown
3912  processor_id_type procup =
3913  cast_int<processor_id_type>((this->processor_id() + p) %
3914  this->n_processors());
3915  processor_id_type procdown =
3916  cast_int<processor_id_type>((this->n_processors() +
3917  this->processor_id() - p) %
3918  this->n_processors());
3919  std::vector<dof_id_type> dof_request_to_fill;
3920 
3921  this->comm().send_receive(procup, requested_dof_ids[procup],
3922  procdown, dof_request_to_fill);
3923 
3924  // Fill those requests
3925  std::vector<std::vector<dof_id_type>> dof_row_keys(dof_request_to_fill.size());
3926 
3927  std::vector<std::vector<Real>> dof_row_vals(dof_request_to_fill.size());
3928  std::vector<Number> dof_row_rhss(dof_request_to_fill.size());
3929  std::vector<std::vector<Number>>
3930  dof_adj_rhss(max_qoi_num,
3931  std::vector<Number>(dof_request_to_fill.size()));
3932  for (std::size_t i=0; i != dof_request_to_fill.size(); ++i)
3933  {
3934  dof_id_type constrained = dof_request_to_fill[i];
3935  if (_dof_constraints.count(constrained))
3936  {
3937  DofConstraintRow & row = _dof_constraints[constrained];
3938  std::size_t row_size = row.size();
3939  dof_row_keys[i].reserve(row_size);
3940  dof_row_vals[i].reserve(row_size);
3941  for (DofConstraintRow::iterator j = row.begin();
3942  j != row.end(); ++j)
3943  {
3944  dof_row_keys[i].push_back(j->first);
3945  dof_row_vals[i].push_back(j->second);
3946 
3947  // We should never have a 0 constraint
3948  // coefficient; that's implicit via sparse
3949  // constraint storage
3950  libmesh_assert(j->second);
3951  }
3952  DofConstraintValueMap::const_iterator rhsit =
3953  _primal_constraint_values.find(constrained);
3954  dof_row_rhss[i] = (rhsit == _primal_constraint_values.end()) ?
3955  0 : rhsit->second;
3956 
3957  for (unsigned int q = 0; q != max_qoi_num; ++q)
3958  {
3959  AdjointDofConstraintValues::const_iterator adjoint_map_it =
3960  _adjoint_constraint_values.find(q);
3961 
3962  if (adjoint_map_it == _adjoint_constraint_values.end())
3963  continue;
3964 
3965  const DofConstraintValueMap & constraint_map =
3966  adjoint_map_it->second;
3967 
3968  DofConstraintValueMap::const_iterator rhsit =
3969  constraint_map.find(constrained);
3970  dof_adj_rhss[q][i] = (rhsit == constraint_map.end()) ?
3971  0 : rhsit->second;
3972  }
3973  }
3974  else
3975  {
3976  // Get NaN from Real, where it should exist, not
3977  // from Number, which may be std::complex, in which
3978  // case quiet_NaN() silently returns zero, rather
3979  // than sanely returning NaN or throwing an
3980  // exception or sending Stroustrup hate mail.
3981  dof_row_rhss[i] =
3982  std::numeric_limits<Real>::quiet_NaN();
3983 
3984  // Make sure we don't get caught by "!isnan(NaN)"
3985  // bugs again.
3986  libmesh_assert(libmesh_isnan(dof_row_rhss[i]));
3987  }
3988  }
3989 
3990  // Trade back the results
3991  std::vector<std::vector<dof_id_type>> dof_filled_keys;
3992  std::vector<std::vector<Real>> dof_filled_vals;
3993  std::vector<Number> dof_filled_rhss;
3994  std::vector<std::vector<Number>> adj_filled_rhss;
3995  this->comm().send_receive(procdown, dof_row_keys,
3996  procup, dof_filled_keys);
3997  this->comm().send_receive(procdown, dof_row_vals,
3998  procup, dof_filled_vals);
3999  this->comm().send_receive(procdown, dof_row_rhss,
4000  procup, dof_filled_rhss);
4001  this->comm().send_receive(procdown, dof_adj_rhss,
4002  procup, adj_filled_rhss);
4003 
4004  libmesh_assert_equal_to (dof_filled_keys.size(), requested_dof_ids[procup].size());
4005  libmesh_assert_equal_to (dof_filled_vals.size(), requested_dof_ids[procup].size());
4006  libmesh_assert_equal_to (dof_filled_rhss.size(), requested_dof_ids[procup].size());
4007 #ifndef NDEBUG
4008  for (std::size_t q=0; q != adj_filled_rhss.size(); ++q)
4009  libmesh_assert_equal_to (adj_filled_rhss[q].size(), requested_dof_ids[procup].size());
4010 #endif
4011 
4012  // Add any new constraint rows we've found
4013  for (std::size_t i=0; i != requested_dof_ids[procup].size(); ++i)
4014  {
4015  libmesh_assert_equal_to (dof_filled_keys[i].size(), dof_filled_vals[i].size());
4016  if (!libmesh_isnan(dof_filled_rhss[i]))
4017  {
4018  dof_id_type constrained = requested_dof_ids[procup][i];
4019  DofConstraintRow & row = _dof_constraints[constrained];
4020  for (std::size_t j = 0; j != dof_filled_keys[i].size(); ++j)
4021  row[dof_filled_keys[i][j]] = dof_filled_vals[i][j];
4022  if (dof_filled_rhss[i] != Number(0))
4023  _primal_constraint_values[constrained] = dof_filled_rhss[i];
4024  else
4025  _primal_constraint_values.erase(constrained);
4026 
4027  for (unsigned int q = 0; q != max_qoi_num; ++q)
4028  {
4029  AdjointDofConstraintValues::iterator adjoint_map_it =
4030  _adjoint_constraint_values.find(q);
4031 
4032  if ((adjoint_map_it == _adjoint_constraint_values.end()) &&
4033  adj_filled_rhss[q][constrained] == Number(0))
4034  continue;
4035 
4036  if (adjoint_map_it == _adjoint_constraint_values.end())
4037  adjoint_map_it = _adjoint_constraint_values.insert
4038  (std::make_pair(q,DofConstraintValueMap())).first;
4039 
4040  DofConstraintValueMap & constraint_map =
4041  adjoint_map_it->second;
4042 
4043  if (adj_filled_rhss[q][i] != Number(0))
4044  constraint_map[constrained] =
4045  adj_filled_rhss[q][i];
4046  else
4047  constraint_map.erase(constrained);
4048  }
4049 
4050  // And prepare to check for more recursive constraints
4051  if (!dof_filled_keys[i].empty())
4052  unexpanded_dofs.insert(constrained);
4053  }
4054  }
4055  }
4056 
4057  // We have to keep recursing while the unexpanded set is
4058  // nonempty on *any* processor
4059  unexpanded_set_nonempty = !unexpanded_dofs.empty();
4060  this->comm().max(unexpanded_set_nonempty);
4061  }
4062 }
4063 
4064 void DofMap::add_constraints_to_send_list()
4065 {
4066  // This function must be run on all processors at once
4067  parallel_object_only();
4068 
4069  // Return immediately if there's nothing to gather
4070  if (this->n_processors() == 1)
4071  return;
4072 
4073  // We might get to return immediately if none of the processors
4074  // found any constraints
4075  unsigned int has_constraints = !_dof_constraints.empty();
4076  this->comm().max(has_constraints);
4077  if (!has_constraints)
4078  return;
4079 
4080  for (DofConstraints::iterator i = _dof_constraints.begin();
4081  i != _dof_constraints.end(); ++i)
4082  {
4083  dof_id_type constrained_dof = i->first;
4084 
4085  // We only need the dependencies of our own constrained dofs
4086  if (constrained_dof < this->first_dof() ||
4087  constrained_dof >= this->end_dof())
4088  continue;
4089 
4090  DofConstraintRow & constraint_row = i->second;
4091  for (DofConstraintRow::const_iterator
4092  j=constraint_row.begin(); j != constraint_row.end();
4093  ++j)
4094  {
4095  dof_id_type constraint_dependency = j->first;
4096 
4097  // No point in adding one of our own dofs to the send_list
4098  if (constraint_dependency >= this->first_dof() &&
4099  constraint_dependency < this->end_dof())
4100  continue;
4101 
4102  _send_list.push_back(constraint_dependency);
4103  }
4104  }
4105 }
4106 
4107 
4108 
4109 #endif // LIBMESH_ENABLE_CONSTRAINTS
4110 
4111 
4112 #ifdef LIBMESH_ENABLE_AMR
4113 
4114 void DofMap::constrain_p_dofs (unsigned int var,
4115  const Elem * elem,
4116  unsigned int s,
4117  unsigned int p)
4118 {
4119  // We're constraining dofs on elem which correspond to p refinement
4120  // levels above p - this only makes sense if elem's p refinement
4121  // level is above p.
4122  libmesh_assert_greater (elem->p_level(), p);
4123  libmesh_assert_less (s, elem->n_sides());
4124 
4125  const unsigned int sys_num = this->sys_number();
4126  const unsigned int dim = elem->dim();
4127  ElemType type = elem->type();
4128  FEType low_p_fe_type = this->variable_type(var);
4129  FEType high_p_fe_type = this->variable_type(var);
4130  low_p_fe_type.order = static_cast<Order>(low_p_fe_type.order + p);
4131  high_p_fe_type.order = static_cast<Order>(high_p_fe_type.order +
4132  elem->p_level());
4133 
4134  const unsigned int n_nodes = elem->n_nodes();
4135  for (unsigned int n = 0; n != n_nodes; ++n)
4136  if (elem->is_node_on_side(n, s))
4137  {
4138  const Node & node = elem->node_ref(n);
4139  const unsigned int low_nc =
4140  FEInterface::n_dofs_at_node (dim, low_p_fe_type, type, n);
4141  const unsigned int high_nc =
4142  FEInterface::n_dofs_at_node (dim, high_p_fe_type, type, n);
4143 
4144  // since we may be running this method concurrently
4145  // on multiple threads we need to acquire a lock
4146  // before modifying the _dof_constraints object.
4147  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
4148 
4149  if (elem->is_vertex(n))
4150  {
4151  // Add "this is zero" constraint rows for high p vertex
4152  // dofs
4153  for (unsigned int i = low_nc; i != high_nc; ++i)
4154  {
4155  _dof_constraints[node.dof_number(sys_num,var,i)].clear();
4156  _primal_constraint_values.erase(node.dof_number(sys_num,var,i));
4157  }
4158  }
4159  else
4160  {
4161  const unsigned int total_dofs = node.n_comp(sys_num, var);
4162  libmesh_assert_greater_equal (total_dofs, high_nc);
4163  // Add "this is zero" constraint rows for high p
4164  // non-vertex dofs, which are numbered in reverse
4165  for (unsigned int j = low_nc; j != high_nc; ++j)
4166  {
4167  const unsigned int i = total_dofs - j - 1;
4168  _dof_constraints[node.dof_number(sys_num,var,i)].clear();
4169  _primal_constraint_values.erase(node.dof_number(sys_num,var,i));
4170  }
4171  }
4172  }
4173 }
4174 
4175 #endif // LIBMESH_ENABLE_AMR
4176 
4177 
4178 #ifdef LIBMESH_ENABLE_DIRICHLET
4179 void DofMap::add_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary)
4180 {
4181  _dirichlet_boundaries->push_back(new DirichletBoundary(dirichlet_boundary));
4182 }
4183 
4184 
4185 void DofMap::add_adjoint_dirichlet_boundary (const DirichletBoundary & dirichlet_boundary,
4186  unsigned int qoi_index)
4187 {
4188  unsigned int old_size = cast_int<unsigned int>
4189  (_adjoint_dirichlet_boundaries.size());
4190  for (unsigned int i = old_size; i <= qoi_index; ++i)
4191  _adjoint_dirichlet_boundaries.push_back(new DirichletBoundaries());
4192 
4193  _adjoint_dirichlet_boundaries[qoi_index]->push_back
4194  (new DirichletBoundary(dirichlet_boundary));
4195 }
4196 
4197 
4198 bool DofMap::has_adjoint_dirichlet_boundaries(unsigned int q) const
4199 {
4200  if (_adjoint_dirichlet_boundaries.size() > q)
4201  return true;
4202 
4203  return false;
4204 }
4205 
4206 
4207 const DirichletBoundaries *
4208 DofMap::get_adjoint_dirichlet_boundaries(unsigned int q) const
4209 {
4210  libmesh_assert_greater(_adjoint_dirichlet_boundaries.size(),q);
4211  return _adjoint_dirichlet_boundaries[q];
4212 }
4213 
4214 
4216 DofMap::get_adjoint_dirichlet_boundaries(unsigned int q)
4217 {
4218  unsigned int old_size = cast_int<unsigned int>
4219  (_adjoint_dirichlet_boundaries.size());
4220  for (unsigned int i = old_size; i <= q; ++i)
4221  _adjoint_dirichlet_boundaries.push_back(new DirichletBoundaries());
4222 
4223  return _adjoint_dirichlet_boundaries[q];
4224 }
4225 
4226 
4227 void DofMap::remove_dirichlet_boundary (const DirichletBoundary & boundary_to_remove)
4228 {
4229  // Find a boundary condition matching the one to be removed
4230  std::vector<DirichletBoundary *>::iterator it = _dirichlet_boundaries->begin();
4231  std::vector<DirichletBoundary *>::iterator end = _dirichlet_boundaries->end();
4232  for (; it != end; ++it)
4233  {
4234  DirichletBoundary * bdy = *it;
4235 
4236  if ((bdy->b == boundary_to_remove.b) &&
4237  bdy->variables == boundary_to_remove.variables)
4238  break;
4239  }
4240 
4241  // Delete it and remove it
4242  libmesh_assert (it != end);
4243  delete *it;
4244  _dirichlet_boundaries->erase(it);
4245 }
4246 
4247 
4248 void DofMap::remove_adjoint_dirichlet_boundary (const DirichletBoundary & boundary_to_remove,
4249  unsigned int qoi_index)
4250 {
4251  libmesh_assert_greater(_adjoint_dirichlet_boundaries.size(),
4252  qoi_index);
4253 
4254  // Find a boundary condition matching the one to be removed
4255  std::vector<DirichletBoundary *>::iterator it =
4256  _adjoint_dirichlet_boundaries[qoi_index]->begin();
4257  std::vector<DirichletBoundary *>::iterator end =
4258  _adjoint_dirichlet_boundaries[qoi_index]->end();
4259  for (; it != end; ++it)
4260  {
4261  DirichletBoundary * bdy = *it;
4262 
4263  if ((bdy->b == boundary_to_remove.b) &&
4264  bdy->variables == boundary_to_remove.variables)
4265  break;
4266  }
4267 
4268  // Delete it and remove it
4269  libmesh_assert (it != end);
4270  delete *it;
4271  _adjoint_dirichlet_boundaries[qoi_index]->erase(it);
4272 }
4273 
4274 
4275 DirichletBoundaries::~DirichletBoundaries()
4276 {
4277  for (std::vector<DirichletBoundary *>::iterator it = begin(); it != end(); ++it)
4278  delete *it;
4279 }
4280 
4281 void DofMap::check_dirichlet_bcid_consistency (const MeshBase & mesh,
4282  const DirichletBoundary & boundary) const
4283 {
4284  const std::set<boundary_id_type>& mesh_bcids = mesh.get_boundary_info().get_boundary_ids();
4285  const std::set<boundary_id_type>& dbc_bcids = boundary.b;
4286 
4287  // DirichletBoundary id sets should be consistent across all ranks
4288  libmesh_assert(mesh.comm().verify(dbc_bcids.size()));
4289 
4290  for (std::set<boundary_id_type>::const_iterator bid = dbc_bcids.begin();
4291  bid != dbc_bcids.end(); ++bid)
4292  {
4293  // DirichletBoundary id sets should be consistent across all ranks
4294  libmesh_assert(mesh.comm().verify(*bid));
4295 
4296  bool found_bcid = (mesh_bcids.find(*bid) != mesh_bcids.end());
4297 
4298  // On a distributed mesh, boundary id sets may *not* be
4299  // consistent across all ranks, since not all ranks see all
4300  // boundaries
4301  mesh.comm().max(found_bcid);
4302 
4303  if (!found_bcid)
4304  libmesh_error_msg("Could not find Dirichlet boundary id " << *bid << " in mesh!");
4305  }
4306 }
4307 
4308 #endif // LIBMESH_ENABLE_DIRICHLET
4309 
4310 
4311 #ifdef LIBMESH_ENABLE_PERIODIC
4312 
4313 void DofMap::add_periodic_boundary (const PeriodicBoundaryBase & periodic_boundary)
4314 {
4315  // See if we already have a periodic boundary associated myboundary...
4316  PeriodicBoundaryBase * existing_boundary = _periodic_boundaries->boundary(periodic_boundary.myboundary);
4317 
4318  if (existing_boundary == libmesh_nullptr)
4319  {
4320  // ...if not, clone the input (and its inverse) and add them to the PeriodicBoundaries object
4321  PeriodicBoundaryBase * boundary = periodic_boundary.clone().release();
4322  PeriodicBoundaryBase * inverse_boundary = periodic_boundary.clone(PeriodicBoundaryBase::INVERSE).release();
4323 
4324  // _periodic_boundaries takes ownership of the pointers
4325  _periodic_boundaries->insert(std::make_pair(boundary->myboundary, boundary));
4326  _periodic_boundaries->insert(std::make_pair(inverse_boundary->myboundary, inverse_boundary));
4327  }
4328  else
4329  {
4330  // ...otherwise, merge this object's variable IDs with the existing boundary object's.
4331  existing_boundary->merge(periodic_boundary);
4332 
4333  // Do the same merging process for the inverse boundary. Note: the inverse better already exist!
4334  PeriodicBoundaryBase * inverse_boundary = _periodic_boundaries->boundary(periodic_boundary.pairedboundary);
4335  libmesh_assert(inverse_boundary);
4336  inverse_boundary->merge(periodic_boundary);
4337  }
4338 }
4339 
4340 
4341 
4342 
4343 void DofMap::add_periodic_boundary (const PeriodicBoundaryBase & boundary,
4344  const PeriodicBoundaryBase & inverse_boundary)
4345 {
4346  libmesh_assert_equal_to (boundary.myboundary, inverse_boundary.pairedboundary);
4347  libmesh_assert_equal_to (boundary.pairedboundary, inverse_boundary.myboundary);
4348 
4349  // Allocate copies on the heap. The _periodic_boundaries object will manage this memory.
4350  // Note: this also means that the copy constructor for the PeriodicBoundary (or user class
4351  // derived therefrom) must be implemented!
4352  // PeriodicBoundary * p_boundary = new PeriodicBoundary(boundary);
4353  // PeriodicBoundary * p_inverse_boundary = new PeriodicBoundary(inverse_boundary);
4354 
4355  // We can't use normal copy construction since this leads to slicing with derived classes.
4356  // Use clone() "virtual constructor" instead. But, this *requires* user to override the clone()
4357  // method. Note also that clone() allocates memory. In this case, the _periodic_boundaries object
4358  // takes responsibility for cleanup.
4359  PeriodicBoundaryBase * p_boundary = boundary.clone().release();
4360  PeriodicBoundaryBase * p_inverse_boundary = inverse_boundary.clone().release();
4361 
4362  // Add the periodic boundary and its inverse to the PeriodicBoundaries data structure. The
4363  // PeriodicBoundaries data structure takes ownership of the pointers.
4364  _periodic_boundaries->insert(std::make_pair(p_boundary->myboundary, p_boundary));
4365  _periodic_boundaries->insert(std::make_pair(p_inverse_boundary->myboundary, p_inverse_boundary));
4366 }
4367 
4368 
4369 #endif
4370 
4371 
4372 } // namespace libMesh
std::unique_ptr< FEMFunctionBase< Gradient > > g_fem
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const =0
FEFamily family
Definition: fe_type.h:204
virtual bool closed() const
const BoundaryInfo & get_boundary_info() const
Definition: mesh_base.h:118
double abs(double a)
const FEType & type() const
Definition: variable.h:119
virtual Output component(unsigned int i, const Point &p, Real time=0.)
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
virtual bool is_serial() const
Definition: mesh_base.h:141
bool active() const
Definition: elem.h:2258
unsigned int n() const
virtual numeric_index_type size() const =0
std::unique_ptr< FunctionBase< Number > > f
virtual numeric_index_type last_local_index() const =0
virtual ElemType type() const =0
virtual void zero() libmesh_override
Definition: dense_matrix.h:792
std::unique_ptr< PointLocatorBase > sub_point_locator() const
Definition: mesh_base.C:536
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:2423
void parallel_for(const Range &range, const Body &body)
Definition: threads_none.h:73
Maps between boundary ids and PeriodicBoundaryBase objects.
virtual unsigned int n_edges() const =0
void vector_mult_transpose(DenseVector< T > &dest, const DenseVector< T > &arg) const
virtual element_iterator local_elements_begin()=0
void resize(const unsigned int n)
Definition: dense_vector.h:350
virtual void zero() libmesh_override
Definition: dense_vector.h:374
The base class for all geometric element types.
Definition: elem.h:90
MeshBase & mesh
uint8_t processor_id_type
Definition: id_types.h:99
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
Class for specifying Dirichlet boundary conditions as constraints.
virtual const Node * node_ptr(const dof_id_type i) const =0
std::vector< unsigned int > variables
OrderWrapper order
Definition: fe_type.h:198
static const Real TOLERANCE
IterBase * end
Utility class for defining generic ranges for threading.
Definition: stored_range.h:52
const Variable & variable(const unsigned int c) const
Definition: dof_map.h:1689
const unsigned int n_vars
Definition: tecplot_io.C:68
dof_id_type dof_id
Definition: xdr_io.C:48
const_iterator begin() const
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)
unsigned int first_scalar_number() const
Definition: variable.h:113
void iota(ForwardIter first, ForwardIter last, T value)
Definition: utility.h:57
const std::set< boundary_id_type > & get_boundary_ids() const
Base class for Mesh.
Definition: mesh_base.h:69
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
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:168
std::unique_ptr< FEMFunctionBase< Number > > f_fem
const dof_id_type n_nodes
Definition: tecplot_io.C:67
const MeshBase & get_mesh() const
Definition: system.h:2014
spin_mutex spin_mtx
Definition: threads.C:29
A variable which is solved for in a System of equations.
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
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1875
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
virtual void init_context(const FEMContext &)
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1509
Used by the Mesh to keep track of boundary nodes and elements.
Definition: boundary_info.h:56
bool active_on_subdomain(subdomain_id_type sid) const
Definition: variable.h:136
std::unordered_map< const Elem *, const CouplingMatrix * > map_type
virtual element_iterator local_elements_end()=0
const Node & node_ref(const unsigned int i) const
Definition: elem.h:1897
void libmesh_assert_valid_boundary_ids(const MeshBase &mesh)
Definition: mesh_tools.C:1194
subdomain_id_type subdomain_id() const
Definition: elem.h:1952
virtual unsigned int n_sides() const =0
bool is_prepared() const
Definition: mesh_base.h:134
const_iterator end() const
Definition: stored_range.h:266
virtual void close()=0
An output iterator for use with packed_range functions.
std::vector< object_type >::const_iterator const_iterator
Definition: stored_range.h:58
void libmesh_ignore(const T &)
std::unique_ptr< QBase > default_quadrature_rule(const unsigned int dim, const int extraorder=0) const
Definition: fe_type.C:31
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)
const Point & point(const unsigned int i) const
Definition: elem.h:1810
std::unique_ptr< FunctionBase< Gradient > > g
bool verify(const T &r, const Communicator &comm=Communicator_World)
StoredRange< iterator_type, object_type > & reset(const iterator_type &first, const iterator_type &last)
Definition: stored_range.h:222
virtual bool is_vertex(const unsigned int i) const =0
bool libmesh_isnan(float a)
const Parallel::Communicator & comm() const
virtual void get(const std::vector< numeric_index_type > &index, T *values) const
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:150
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)
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
Definition: dof_map.h:89
virtual unsigned int dim() const =0
virtual numeric_index_type first_local_index() const =0
std::unique_ptr< NumericVector< Number > > current_local_solution
Definition: system.h:1521
virtual element_iterator active_not_local_elements_begin()=0
Base class for all PeriodicBoundary implementations.
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
Definition: dof_map.h:137
virtual void set(const numeric_index_type i, const T value)=0
void cholesky_solve(const DenseVector< T2 > &b, DenseVector< T2 > &x)
Real size() const
Definition: type_vector.h:898
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1832
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
virtual dof_id_type n_elem() const =0
A matrix object used for finite element assembly and numerics.
Definition: dense_matrix.h:54
virtual element_iterator active_local_elements_end()=0
bool empty() const
Definition: stored_range.h:301
processor_id_type processor_id()
Definition: libmesh_base.h:96
A geometric point in (x,y,z) space.
Definition: point.h:38
virtual void localize(std::vector< T > &v_local) const =0
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
Definition: dof_object.h:694
virtual std::unique_ptr< PeriodicBoundaryBase > clone(TransformationType t=FORWARD) const =0
processor_id_type n_processors()
Definition: libmesh_base.h:88
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Definition: dof_map.C:1924