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