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