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

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo