From 24217b0b1cbab235780f35747add6da16f9390d0 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Thu, 16 Sep 2010 14:43:35 +0000
Subject: [PATCH] Fix mixing interpolate_boundary_values with other types of
 constraints as e.g. hanging node constraints: We do not set boundary values
 to dofs that are already constrained.

git-svn-id: https://svn.dealii.org/trunk@21992 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/deal.II/include/numerics/vectors.h    | 117 ++--
 .../include/numerics/vectors.templates.h      | 633 +-----------------
 deal.II/doc/news/changes.h                    |  15 +
 deal.II/examples/step-22/step-22.cc           |  67 +-
 deal.II/lac/include/lac/constraint_matrix.h   |  40 +-
 5 files changed, 185 insertions(+), 687 deletions(-)

diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h
index 2cec72f3be..ccfa873e2b 100644
--- a/deal.II/deal.II/include/numerics/vectors.h
+++ b/deal.II/deal.II/include/numerics/vectors.h
@@ -687,31 +687,43 @@ class VectorTools
 
 
 				     /**
-				      * Insert the (algebraic) constraints
-				      * due to Dirichlet boundary conditions
-				      * to the ConstraintMatrix. This
-				      * function makes up the list of
-				      * degrees of freedom subject to
-				      * Dirichlet boundary conditions and
-				      * the values to be assigned to them,
-				      * by interpolation around the
-				      * boundary. If the ConstraintMatrix @p
-				      * constraints contained values or
-				      * other constraints before, the new
-				      * ones are added, or the old ones
-				      * overwritten if a node of the
-				      * boundary part to be used was already
-				      * in the list of constraints. This is
-				      * handled by using inhomogeneous
-				      * constraints. Please note that when
-				      * combining adaptive meshes and this
-				      * kind of constraints, the Dirichlet
-				      * conditions should be set first, and
-				      * then completed by hanging node
-				      * constraints, in order to make sure
-				      * that the discretization remains
-				      * consistent.
-				      *
+				      * Insert the (algebraic) constraints due
+				      * to Dirichlet boundary conditions into
+				      * a ConstraintMatrix @p
+				      * constraints. This function identifies
+				      * the degrees of freedom subject to
+				      * Dirichlet boundary conditions, adds
+				      * them to the list of constrained DoFs
+				      * in @p constraints and sets the
+				      * respective inhomogeneity to the value
+				      * interpolated around the boundary. If
+				      * this routine encounters a DoF that
+				      * already is constrained (for instance
+				      * by a hanging node constraint, see
+				      * below, or any other type of
+				      * constraint, e.g. from periodic
+				      * boundary conditions), the old setting
+				      * of the constraint (dofs the entry is
+				      * constrained to, inhomogeneities) is
+				      * kept and nothing happens.
+				      *
+				      * Please note that when combining
+				      * adaptively refined meshes with hanging
+				      * node constraints and inhomogeneous
+				      * Dirichlet boundary conditions within
+				      * one ConstraintMatrix object, the
+				      * hanging node constraints should always
+				      * be set first, and then Dirichlet
+				      * boundary conditions should be
+				      * interpolated. This makes sure that the
+				      * discretization remains H<sup>1</sup>
+				      * conforming as is needed e.g. for the
+				      * Laplace equation in 3D, as hanging
+				      * nodes on the boundary should be still
+				      * set to the weighted average of
+				      * neighbors, and not the actual
+				      * Dirichlet value.
+ *				      *
 				      * The parameter @p boundary_component
 				      * corresponds to the number @p
 				      * boundary_indicator of the face.  255
@@ -880,25 +892,42 @@ class VectorTools
 				       std::vector<unsigned int> component_mapping = std::vector<unsigned int>());
 
 				     /**
-				      * Project a function to the boundary
-				      * of the domain, using the given
-				      * quadrature formula for the faces. If
-				      * the ConstraintMatrix @p constraints
-				      * contained values or other
-				      * constraints before, the new ones are
-				      * added, or the old ones overwritten
-				      * if a node of the boundary part to be
-				      * used was already in the list of
-				      * constraints. This is handled by
-				      * using inhomogeneous
-				      * constraints. Please note that when
-				      * combining adaptive meshes and this
-				      * kind of constraints, the Dirichlet
-				      * conditions should be set first, and
-				      * then completed by hanging node
-				      * constraints, in order to make sure
-				      * that the discretization remains
-				      * consistent.
+				      * Project a function to the boundary of
+				      * the domain, using the given quadrature
+				      * formula for the faces. This function
+				      * identifies the degrees of freedom
+				      * subject to Dirichlet boundary
+				      * conditions, adds them to the list of
+				      * constrained DoFs in @p constraints and
+				      * sets the respective inhomogeneity to
+				      * the value resulting from the
+				      * projection operation. If this routine
+				      * encounters a DoF that already is
+				      * constrained (for instance by a hanging
+				      * node constraint, see below, or any
+				      * other type of constraint, e.g. from
+				      * periodic boundary conditions), the old
+				      * setting of the constraint (dofs the
+				      * entry is constrained to,
+				      * inhomogeneities) is kept and nothing
+				      * happens.
+				      *
+				      * Please note that when combining
+				      * adaptively refined meshes with hanging
+				      * node constraints and inhomogeneous
+				      * Dirichlet boundary conditions within
+				      * one ConstraintMatrix object, the
+				      * hanging node constraints should always
+				      * be set first, and then Dirichlet
+				      * boundary conditions should be
+				      * interpolated. This makes sure that the
+				      * discretization remains H<sup>1</sup>
+				      * conforming as is needed e.g. for the
+				      * Laplace equation in 3D, as hanging
+				      * nodes on the boundary should be still
+				      * set to the weighted average of
+				      * neighbors, and not the actual
+				      * Dirichlet value.
 				      *
 				      * If @p component_mapping is empty, it
 				      * is assumed that the number of
diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h
index f6d43aaf55..c975d44e56 100644
--- a/deal.II/deal.II/include/numerics/vectors.templates.h
+++ b/deal.II/deal.II/include/numerics/vectors.templates.h
@@ -1236,6 +1236,10 @@ VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space
 					  std::map<unsigned int,double> &boundary_values,
 					  const std::vector<bool>       &component_mask_)
 {
+  //ConstraintMatrix boundary_constraints();
+  //  interpolate_boundary_values (dof, boundary_component, boundary_function, 
+  //			       boundary_constraints, component_mask);
+  
   const unsigned int dim=DH::dimension;
   const unsigned int spacedim=DH::space_dimension;
 
@@ -1732,132 +1736,6 @@ VectorTools::interpolate_boundary_values (const DH                 &dof,
 
 // ----------- interpolate_boundary_values for ConstraintMatrix --------------
 
-// TODO (M.K.): There is a lot of duplicated code with the above
-// function. We should unify all these interpolate_boundary_values functions
-// in one way or the other.
-
-#if deal_II_dimension == 1
-
-//TODO[?] Actually the Mapping object should be a MappingCollection object for the
-// hp::DoFHandler.
-
-//template <int dim, template <int, int> class DH, int spacedim>
-
-template <class DH>
-void
-VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension> &,
-					  const DH                &dof,
-					  const unsigned char      boundary_component,
-					  const Function<DH::space_dimension> &boundary_function,
-					  ConstraintMatrix        &constraints,
-					  const std::vector<bool> &component_mask_)
-{
-  const unsigned int dim=DH::dimension;
-  const unsigned int spacedim=DH::space_dimension;
-
-  Assert (boundary_component != 255,
-	  ExcInvalidBoundaryIndicator());
-  Assert ((component_mask_.size() == 0) ||
-	  (component_mask_.size() == dof.get_fe().n_components()),
-	  ExcMessage ("The number of components in the mask has to be either "
-		      "zero or equal to the number of components in the finite "
-		      "element."));
-
-				   // check whether boundary values at
-				   // the left or right boundary of
-				   // the line are
-				   // requested. direction denotes
-				   // the neighboring direction in
-				   // which we seek the boundary,
-				   // i.e. 0 is left boundary and 1 is
-				   // right.
-  const unsigned int direction = boundary_component;
-  Assert (direction < 2, ExcInvalidBoundaryIndicator());
-
-				   // first find the outermost active
-				   // cell by first traversing the coarse
-				   // grid to its end and then going
-				   // to the children
-  typename DH::cell_iterator outermost_cell = dof.begin(0);
-  while (outermost_cell->neighbor(direction).state() == IteratorState::valid)
-    outermost_cell = outermost_cell->neighbor(direction);
-
-  while (outermost_cell->has_children())
-    outermost_cell = outermost_cell->child(direction);
-
-                                   // get the FE corresponding to this
-                                   // cell
-  const FiniteElement<dim,spacedim> &fe = outermost_cell->get_fe();
-  Assert (fe.n_components() == boundary_function.n_components,
-	  ExcDimensionMismatch(fe.n_components(), boundary_function.n_components));
-
-				   // set the component mask to either
-				   // the original value or a vector
-				   // of trues
-  const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
-					  std::vector<bool> (fe.n_components(), true) :
-					  component_mask_);
-  Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
-	  ExcNoComponentSelected());
-
-				   // now set the value of the
-				   // outermost degree of
-				   // freedom. setting also
-				   // creates the entry in the map
-				   // if it did not exist
-				   // beforehand
-				   //
-				   // save some time by requesting
-				   // values only once for each point,
-				   // irrespective of the number of
-				   // components of the function
-  Vector<double> function_values (fe.n_components());
-  if (fe.n_components() == 1)
-    function_values(0)
-      = boundary_function.value (outermost_cell->vertex(direction));
-  else
-    boundary_function.vector_value (outermost_cell->vertex(direction),
-				    function_values);
-
-  for (unsigned int i=0; i<fe.dofs_per_vertex; ++i)
-    if (component_mask[fe.face_system_to_component_index(i).first])
-      {
-				   // TODO: should we clear the other
-				   // entries in the line here?
-	const unsigned int row = outermost_cell->vertex_dof_index(direction,i,
-								  outermost_cell->active_fe_index());
-	constraints.add_line (row);
-	constraints.set_inhomogeneity (row,
-		    function_values(fe.face_system_to_component_index(i).first));
-      }
-}
-
-
-//TODO[?] Actually the Mapping object should be a MappingCollection object for the
-// hp::DoFHandler.
-
-// Implementation for 1D
-template <class DH>
-void
-VectorTools::interpolate_boundary_values
-  (const Mapping<DH::dimension, DH::space_dimension>     &mapping,
-   const DH                                              &dof,
-   const typename FunctionMap<DH::space_dimension>::type &function_map,
-   ConstraintMatrix                                      &constraints,
-   const std::vector<bool>                               &component_mask)
-{
-  for (typename FunctionMap<DH::space_dimension>::type::const_iterator i=function_map.begin();
-       i!=function_map.end(); ++i)
-    interpolate_boundary_values (mapping, dof, i->first, *i->second,
-				 constraints, component_mask);
-}
-
-
-//TODO[?] Actually the Mapping object should be a MappingCollection object for the
-// hp::DoFHandler.
-
-
-#else
 
 
 template <class DH>
@@ -1869,317 +1747,20 @@ VectorTools::interpolate_boundary_values
   ConstraintMatrix                                      &constraints,
   const std::vector<bool>                               &component_mask_)
 {
-  const unsigned int dim=DH::dimension;
-
-  Assert ((component_mask_.size() == 0) ||
-	  (component_mask_.size() == dof.get_fe().n_components()),
-	  ExcMessage ("The number of components in the mask has to be either "
-		      "zero or equal to the number of components in the finite "
-		      "element."));
-
-
-				   // if for whatever reason we were
-				   // passed an empty map, return
-				   // immediately
-  if (function_map.size() == 0)
-    return;
-
-  Assert (function_map.find(255) == function_map.end(),
-	  ExcInvalidBoundaryIndicator());
-
-  const unsigned int        n_components = DoFTools::n_components(dof);
-  const bool                fe_is_system = (n_components != 1);
-
-  for (typename FunctionMap<DH::space_dimension>::type::const_iterator i=function_map.begin();
-       i!=function_map.end(); ++i)
-    Assert (n_components == i->second->n_components,
-	    ExcDimensionMismatch(n_components, i->second->n_components));
-
-				   // set the component mask to either
-				   // the original value or a vector
-				   // of trues
-  const std::vector<bool> component_mask ((component_mask_.size() == 0) ?
-					  std::vector<bool> (n_components, true) :
-					  component_mask_);
-  Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
-	  ExcNoComponentSelected());
-
-				   // field to store the indices
-  std::vector<unsigned int> face_dofs;
-  face_dofs.reserve (DoFTools::max_dofs_per_face(dof));
-
-  std::vector<Point<DH::space_dimension> >  dof_locations;
-  dof_locations.reserve (DoFTools::max_dofs_per_face(dof));
-
-				   // array to store the values of the
-				   // boundary function at the boundary
-				   // points. have to arrays for scalar and
-				   // vector functions to use the more
-				   // efficient one respectively
-  std::vector<double>          dof_values_scalar;
-  std::vector<Vector<double> > dof_values_system;
-  dof_values_scalar.reserve (DoFTools::max_dofs_per_face (dof));
-  dof_values_system.reserve (DoFTools::max_dofs_per_face (dof));
-
-				   // before we start with the loop over all
-				   // cells create an hp::FEValues object
-				   // that holds the interpolation points of
-				   // all finite elements that may ever be
-				   // in use
-  hp::FECollection<dim> finite_elements (dof.get_fe());
-  hp::QCollection<dim-1>  q_collection;
-  for (unsigned int f=0; f<finite_elements.size(); ++f)
+  std::map<unsigned int,double> boundary_values;
+  interpolate_boundary_values (mapping, dof, function_map,
+			       boundary_values, component_mask_);
+  std::map<unsigned int,double>::const_iterator boundary_value = 
+    boundary_values.begin();
+  for ( ; boundary_value !=boundary_values.end(); ++boundary_value)
     {
-      const FiniteElement<dim> &fe = finite_elements[f];
-
-				       // generate a quadrature rule on the
-				       // face from the unit support
-				       // points. this will be used to
-				       // obtain the quadrature points on
-				       // the real cell's face
-				       //
-				       // to do this, we check whether
-				       // the FE has support points on
-				       // the face at all:
-      if (fe.has_face_support_points())
-	q_collection.push_back (Quadrature<dim-1>(fe.get_unit_face_support_points()));
-      else
+      if (!constraints.is_constrained(boundary_value->first))
 	{
-					   // if not, then we should try a
-					   // more clever way. the idea is
-					   // that a finite element may not
-					   // offer support points for all
-					   // its shape functions, but maybe
-					   // only some. if it offers
-					   // support points for the
-					   // components we are interested
-					   // in in this function, then
-					   // that's fine. if not, the
-					   // function we call in the finite
-					   // element will raise an
-					   // exception. the support points
-					   // for the other shape functions
-					   // are left uninitialized (well,
-					   // initialized by the default
-					   // constructor), since we don't
-					   // need them anyway.
-					   //
-					   // As a detour, we must make sure
-					   // we only query
-					   // face_system_to_component_index
-					   // if the index corresponds to a
-					   // primitive shape
-					   // function. since we know that
-					   // all the components we are
-					   // interested in are primitive
-					   // (by the above check), we can
-					   // safely put such a check in
-					   // front
-	  std::vector<Point<dim-1> > unit_support_points (fe.dofs_per_face);
-
-	  for (unsigned int i=0; i<fe.dofs_per_face; ++i)
-	    if (fe.is_primitive (fe.face_to_equivalent_cell_index(i)))
-	      if (component_mask[fe.face_system_to_component_index(i).first]
-		  == true)
-		unit_support_points[i] = fe.unit_face_support_point(i);
-
-	  q_collection.push_back (Quadrature<dim-1>(unit_support_points));
-        }
+	  constraints.add_line (boundary_value->first);
+	  constraints.set_inhomogeneity (boundary_value->first,
+					 boundary_value->second);
+	}
     }
-				   // now that we have a q_collection object
-				   // with all the right quadrature points,
-				   // create an hp::FEFaceValues object that
-				   // we can use to evaluate the boundary
-				   // values at
-  hp::MappingCollection<dim> mapping_collection (mapping);
-  hp::FEFaceValues<dim> x_fe_values (mapping_collection, finite_elements, q_collection,
-				     update_quadrature_points);
-
-  typename DH::active_cell_iterator cell = dof.begin_active(),
-				    endc = dof.end();
-  for (; cell!=endc; ++cell)
-    for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
-	 ++face_no)
-      {
-        const FiniteElement<dim,DH::space_dimension> &fe = cell->get_fe();
-
-					 // we can presently deal only with
-					 // primitive elements for boundary
-					 // values. this does not preclude
-					 // us using non-primitive elements
-					 // in components that we aren't
-					 // interested in, however. make
-					 // sure that all shape functions
-					 // that are non-zero for the
-					 // components we are interested in,
-					 // are in fact primitive
-	for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
-	  {
-	    const std::vector<bool> &nonzero_component_array
-	      = cell->get_fe().get_nonzero_components (i);
-	    for (unsigned int c=0; c<n_components; ++c)
-	      if ((nonzero_component_array[c] == true)
-		  &&
-		  (component_mask[c] == true))
-		Assert (cell->get_fe().is_primitive (i),
-			ExcMessage ("This function can only deal with requested boundary "
-				    "values that correspond to primitive (scalar) base "
-				    "elements"));
-	  }
-
-	typename DH::face_iterator face = cell->face(face_no);
-	const unsigned char boundary_component = face->boundary_indicator();
-	if (function_map.find(boundary_component) != function_map.end())
-	  {
-					     // face is of the right
-					     // component
-	    x_fe_values.reinit(cell, face_no);
-	    const FEFaceValues<dim> &fe_values = x_fe_values.get_present_fe_values();
-
-					     // get indices, physical
-					     // location and boundary values
-					     // of dofs on this face
-            face_dofs.resize (fe.dofs_per_face);
-	    face->get_dof_indices (face_dofs, cell->active_fe_index());
-	    const std::vector<Point<DH::space_dimension> > &dof_locations
-              = fe_values.get_quadrature_points ();
-
-	    if (fe_is_system)
-	      {
-                                                 // resize array. avoid
-                                                 // construction of a memory
-                                                 // allocating temporary if
-                                                 // possible
-                if (dof_values_system.size() < fe.dofs_per_face)
-                  dof_values_system.resize (fe.dofs_per_face,
-                                            Vector<double>(fe.n_components()));
-                else
-                  dof_values_system.resize (fe.dofs_per_face);
-
-		function_map.find(boundary_component)->second
-                  ->vector_value_list (dof_locations, dof_values_system);
-
-						 // enter those dofs into
-						 // the list that match the
-						 // component
-						 // signature. avoid the
-						 // usual complication that
-						 // we can't just use
-						 // *_system_to_component_index
-						 // for non-primitive FEs
-		for (unsigned int i=0; i<face_dofs.size(); ++i)
-                  {
-                    unsigned int component;
-                    if (fe.is_primitive())
-                      component = fe.face_system_to_component_index(i).first;
-                    else
-                      {
-                                                         // non-primitive
-                                                         // case. make sure
-                                                         // that this
-                                                         // particular shape
-                                                         // function _is_
-                                                         // primitive, and
-                                                         // get at it's
-                                                         // component. use
-                                                         // usual trick to
-                                                         // transfer face
-                                                         // dof index to
-                                                         // cell dof index
-                        const unsigned int cell_i
-                          = (dim == 1 ?
-                             i
-                             :
-                             (dim == 2 ?
-                              (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
-                              :
-                              (dim == 3 ?
-                               (i<4*fe.dofs_per_vertex ?
-                                i
-                                :
-                                (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
-                                 i+4*fe.dofs_per_vertex
-                                 :
-                                 i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
-                               :
-                               numbers::invalid_unsigned_int)));
-                        Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
-
-                                                         // make sure that
-                                                         // if this is not a
-                                                         // primitive shape
-                                                         // function, then
-                                                         // all the
-                                                         // corresponding
-                                                         // components in
-                                                         // the mask are not
-                                                         // set
-                        if (!fe.is_primitive(cell_i))
-                          for (unsigned int c=0; c<n_components; ++c)
-                            if (fe.get_nonzero_components(cell_i)[c])
-                              Assert (component_mask[c] == false,
-                                      FETools::ExcFENotPrimitive());
-
-                                                         // let's pick the
-                                                         // first of
-                                                         // possibly more
-                                                         // than one
-                                                         // non-zero
-                                                         // components. if
-                                                         // shape function
-                                                         // is
-                                                         // non-primitive,
-                                                         // then we will
-                                                         // ignore the
-                                                         // result in the
-                                                         // following
-                                                         // anyway,
-                                                         // otherwise
-                                                         // there's only one
-                                                         // non-zero
-                                                         // component which
-                                                         // we will use
-                        component = (std::find (fe.get_nonzero_components(cell_i).begin(),
-                                                fe.get_nonzero_components(cell_i).end(),
-                                                true)
-                                     -
-                                     fe.get_nonzero_components(cell_i).begin());
-                      }
-
-                    if (component_mask[component] == true)
-		      {
-				   // TODO: check whether we should clear
-				   // the current line first...
-			constraints.add_line (face_dofs[i]);
-			constraints.set_inhomogeneity (face_dofs[i],
-						       dof_values_system[i](component));
-		      }
-                  }
-	      }
-	    else
-					       // fe has only one component,
-					       // so save some computations
-	      {
-						 // get only the one
-						 // component that this
-						 // function has
-                dof_values_scalar.resize (fe.dofs_per_face);
-		function_map.find(boundary_component)->second
-                  ->value_list (dof_locations, dof_values_scalar, 0);
-
-						 // enter into list
-
-		for (unsigned int i=0; i<face_dofs.size(); ++i)
-				   // TODO: check whether we should clear
-				   // the current line first...
-		  {
-		    constraints.add_line (face_dofs[i]);
-		    constraints.set_inhomogeneity (face_dofs[i],
-						   dof_values_scalar[i]);
-		  }
-	      }
-	  }
-      }
 }
 
 
@@ -2200,7 +1781,6 @@ VectorTools::interpolate_boundary_values
 			       component_mask);
 }
 
-#endif
 
 
 template <class DH>
@@ -2437,26 +2017,6 @@ VectorTools::project_boundary_values (const DoFHandler<dim,spacedim>    &dof,
 
 // ----- implementation for project_boundary_values with ConstraintMatrix -----
 
-#if deal_II_dimension == 1
-
-// Implementation for 1D
-template <int dim, int spacedim>
-void
-VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mapping,
-				      const DoFHandler<dim,spacedim>    &dof,
-				      const typename FunctionMap<spacedim>::type &boundary_functions,
-				      const Quadrature<dim-1>  &,
-				      ConstraintMatrix &constraints,
-				      std::vector<unsigned int> component_mapping)
-{
-  Assert (component_mapping.size() == 0, ExcNotImplemented());
-				   // projection in 1d is equivalent
-				   // to interpolation
-  interpolate_boundary_values (mapping, dof, boundary_functions,
-			       constraints, std::vector<bool>());
-}
-
-#else
 
 
 template <int dim, int spacedim>
@@ -2468,160 +2028,23 @@ VectorTools::project_boundary_values (const Mapping<dim, spacedim>       &mappin
 				      ConstraintMatrix &constraints,
 				      std::vector<unsigned int> component_mapping)
 {
-//TODO:[?] In VectorTools::project_boundary_values, no condensation of sparsity
-//    structures, matrices and right hand sides or distribution of
-//    solution vectors is performed. This is ok for dim<3 because then
-//    there are no constrained nodes on the boundary, but is not
-//    acceptable for higher dimensions. Fix this.
-
-  if (component_mapping.size() == 0)
-    {
-      AssertDimension (dof.get_fe().n_components(), boundary_functions.begin()->second->n_components);
-				       // I still do not see why i
-				       // should create another copy
-				       // here
-      component_mapping.resize(dof.get_fe().n_components());
-      for (unsigned int i= 0 ;i < component_mapping.size() ; ++i)
-	component_mapping[i] = i;
-    }
-  else
-    AssertDimension (dof.get_fe().n_components(), component_mapping.size());
-
-  std::vector<unsigned int> dof_to_boundary_mapping;
-  std::set<unsigned char> selected_boundary_components;
-  for (typename FunctionMap<spacedim>::type::const_iterator i=boundary_functions.begin();
-       i!=boundary_functions.end(); ++i)
-    selected_boundary_components.insert (i->first);
-
-  DoFTools::map_dof_to_boundary_indices (dof, selected_boundary_components,
-					 dof_to_boundary_mapping);
-
-				   // Done if no degrees of freedom on
-				   // the boundary
-  if (dof.n_boundary_dofs (boundary_functions) == 0)
-    return;
-				   // set up sparsity structure
-  SparsityPattern sparsity(dof.n_boundary_dofs (boundary_functions),
-			   dof.max_couplings_between_boundary_dofs());
-  DoFTools::make_boundary_sparsity_pattern (dof,
-					    boundary_functions,
-					    dof_to_boundary_mapping,
-					    sparsity);
-
-				   // note: for three or more dimensions, there
-				   // may be constrained nodes on the boundary
-				   // in this case the boundary mass matrix has
-				   // to be condensed and the solution is to
-				   // be distributed afterwards, which is not
-				   // yet implemented. The reason for this is
-				   // that we cannot simply use the condense
-				   // family of functions, since the matrices
-				   // and vectors do not use the global
-				   // numbering but rather the boundary
-				   // numbering, i.e. the condense function
-				   // needs to use another indirection. There
-				   // should be not many technical problems,
-				   // but it needs to be implemented
-  if (dim>=3)
+  std::map<unsigned int,double> boundary_values;
+  project_boundary_values (mapping, dof, boundary_functions, q,
+			   boundary_values, component_mapping);
+  std::map<unsigned int,double>::const_iterator boundary_value = 
+    boundary_values.begin();
+  for ( ; boundary_value !=boundary_values.end(); ++boundary_value)
     {
-#ifdef DEBUG
-// Assert that there are no hanging nodes at the boundary
-      int level = -1;
-      for (typename DoFHandler<dim,spacedim>::active_cell_iterator cell = dof.begin_active();
-	   cell != dof.end(); ++cell)
-	for (unsigned int f=0;f<GeometryInfo<dim>::faces_per_cell;++f)
-	  {
-	    if (cell->at_boundary(f))
-	      {
-		if (level == -1)
-		  level = cell->level();
-		else
-		  {
-		    Assert (level == cell->level(), ExcNotImplemented());
-		  }
-	      }
-	  }
-#endif
+      if (!constraints.is_constrained(boundary_value->first))
+	{
+	  constraints.add_line (boundary_value->first);
+	  constraints.set_inhomogeneity (boundary_value->first,
+					 boundary_value->second);
+	}
     }
-  sparsity.compress();
-
-
-				   // make mass matrix and right hand side
-  SparseMatrix<double> mass_matrix(sparsity);
-  Vector<double>       rhs(sparsity.n_rows());
-
-
-  MatrixCreator::create_boundary_mass_matrix (mapping, dof, q,
-					      mass_matrix, boundary_functions,
-					      rhs, dof_to_boundary_mapping, (const Function<spacedim>*) 0,
-					      component_mapping);
-
-				   // For certain weird elements,
-				   // there might be degrees of
-				   // freedom on the boundary, but
-				   // their shape functions do not
-				   // have support there. Let's
-				   // eliminate them here.
-
-//TODO: Maybe we should figure out if the element really needs this
-
-  FilteredMatrix<Vector<double> > filtered_mass_matrix(mass_matrix, true);
-  FilteredMatrix<Vector<double> > filtered_precondition;
-  std::vector<bool> excluded_dofs(mass_matrix.m(), false);
-
-  double max_element = 0.;
-  for (unsigned int i=0;i<mass_matrix.m();++i)
-    if (mass_matrix.diag_element(i) > max_element)
-      max_element = mass_matrix.diag_element(i);
-
-  for (unsigned int i=0;i<mass_matrix.m();++i)
-    if (mass_matrix.diag_element(i) < 1.e-8 * max_element)
-      {
-	filtered_mass_matrix.add_constraint(i, 0.);
-	filtered_precondition.add_constraint(i, 0.);
-	mass_matrix.diag_element(i) = 1.;
-	excluded_dofs[i] = true;
-      }
-
-  Vector<double> boundary_projection (rhs.size());
-
-				   // Allow for a maximum of 5*n
-				   // steps to reduce the residual by
-				   // 10^-12. n steps may not be
-				   // sufficient, since roundoff
-				   // errors may accumulate for badly
-				   // conditioned matrices
-  ReductionControl        control(5*rhs.size(), 0., 1.e-12, false, false);
-  GrowingVectorMemory<> memory;
-  SolverCG<>              cg(control,memory);
-
-  PreconditionSSOR<> prec;
-  prec.initialize(mass_matrix, 1.2);
-  filtered_precondition.initialize(prec, true);
-				   // solve
-  cg.solve (filtered_mass_matrix, boundary_projection, rhs, filtered_precondition);
-  filtered_precondition.apply_constraints(boundary_projection, true);
-  filtered_precondition.clear();
-				   // fill in boundary values
-  for (unsigned int i=0; i<dof_to_boundary_mapping.size(); ++i)
-    if (dof_to_boundary_mapping[i] != DoFHandler<dim,spacedim>::invalid_dof_index
-    && ! excluded_dofs[dof_to_boundary_mapping[i]])
-				       // this dof is on one of the
-				       // interesting boundary parts
-				       //
-				       // remember: i is the global dof
-				       // number, dof_to_boundary_mapping[i]
-				       // is the number on the boundary and
-				       // thus in the solution vector
-      {
-				   // TODO: check whether we should clear
-				   // the entries in this line first.
-	constraints.add_line (i);
-	constraints.set_inhomogeneity (i, boundary_projection(dof_to_boundary_mapping[i]));
-      }
 }
 
-#endif
+
 
 template <int dim, int spacedim>
 void
diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index fb7822eafb..f03a4e1be8 100644
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -260,6 +260,21 @@ through DoFHandler::get_tria() and DoFHandler::get_fe(), respectively.
 <h3>deal.II</h3>
 
 <ol>
+  <li>
+  <p>
+  Fixed: The methods VectorTools::interpolate_boundary_values and
+  VectorTools::project_boundary_values with ConstraintMatrix argument added
+  inhomogeneities also to DoFs that were already constrained when entering the
+  call without changing the list of DoFs that those were constrained to. This
+  lead to wrong results for inhomogeneous constraints in case e.g.
+  DoFTools::make_hanging_node_constraints was called before. Now, the correct
+  way for combining hanging node constraints and Dirichlet conditions in one
+  ConstraintMatrix is to first make the hanging node constraints, and then
+  interpolate/project the boundary values.
+  <br>
+  (Luca Heltai, Martin Kronbichler 2010/09/16)
+  </p>
+
   <li>
   <p>
   Fixed: The method FEValuesViews::Vector::curl aborted the program in 2d under certain
diff --git a/deal.II/examples/step-22/step-22.cc b/deal.II/examples/step-22/step-22.cc
index d509a89f39..9aa1a95cb0 100644
--- a/deal.II/examples/step-22/step-22.cc
+++ b/deal.II/examples/step-22/step-22.cc
@@ -573,49 +573,52 @@ void StokesProblem<dim>::setup_dofs ()
 				   // Now comes the implementation of
 				   // Dirichlet boundary conditions, which
 				   // should be evident after the discussion
-				   // in the introduction. All that changed
-				   // is that the function already appears
-				   // in the setup functions, whereas we
-				   // were used to see it in some assembly
-				   // routine. Further down below where we
-				   // set up the mesh, we will associate the
-				   // top boundary where we impose Dirichlet
-				   // boundary conditions with boundary
-				   // indicator 1.  We will have to pass
-				   // this boundary indicator as second
-				   // argument to the function below
-				   // interpolating boundary values.  There
-				   // is one more thing, though.  The
-				   // function describing the Dirichlet
-				   // conditions was defined for all
-				   // components, both velocity and
-				   // pressure. However, the Dirichlet
-				   // conditions are to be set for the
-				   // velocity only.  To this end, we use a
-				   // <code>component_mask</code> that
+				   // in the introduction. All that changed is
+				   // that the function already appears in the
+				   // setup functions, whereas we were used to
+				   // see it in some assembly routine. Further
+				   // down below where we set up the mesh, we
+				   // will associate the top boundary where we
+				   // impose Dirichlet boundary conditions
+				   // with boundary indicator 1.  We will have
+				   // to pass this boundary indicator as
+				   // second argument to the function below
+				   // interpolating boundary values.  There is
+				   // one more thing, though.  The function
+				   // describing the Dirichlet conditions was
+				   // defined for all components, both
+				   // velocity and pressure. However, the
+				   // Dirichlet conditions are to be set for
+				   // the velocity only.  To this end, we use
+				   // a <code>component_mask</code> that
 				   // filters out the pressure component, so
 				   // that the condensation is performed on
-				   // velocity degrees of freedom
-				   // only. Since we use adaptively refined
-				   // grids the constraint matrix needs then
-				   // also be filled with hanging node
-				   // constraints needs to generated from
-				   // the DoF handler. Note the order of the
-				   // two functions &mdash; we first insert
-				   // the boundary values into the
-				   // constraint matrix, and then introduce
-				   // the hanging node constraints.
+				   // velocity degrees of freedom only. Since
+				   // we use adaptively refined grids the
+				   // constraint matrix needs to be first
+				   // filled with hanging node constraints
+				   // generated from the DoF handler. Note the
+				   // order of the two functions &mdash; we
+				   // first compute the hanging node
+				   // constraints, and then insert the
+				   // boundary values into the constraint
+				   // matrix. This makes sure that we respect
+				   // H<sup>1</sup> conformity on boundaries
+				   // with hanging nodes (in three space
+				   // dimensions), where the hanging node
+				   // needs to dominate the Dirichlet boundary
+				   // values.
   {
     constraints.clear ();
     std::vector<bool> component_mask (dim+1, true);
     component_mask[dim] = false;
+    DoFTools::make_hanging_node_constraints (dof_handler,
+					     constraints);
     VectorTools::interpolate_boundary_values (dof_handler,
 					      1,
 					      BoundaryValues<dim>(),
 					      constraints,
 					      component_mask);
-    DoFTools::make_hanging_node_constraints (dof_handler,
-					     constraints);
   }
 
   constraints.close ();
diff --git a/deal.II/lac/include/lac/constraint_matrix.h b/deal.II/lac/include/lac/constraint_matrix.h
index 388c1290bc..6e454c7370 100644
--- a/deal.II/lac/include/lac/constraint_matrix.h
+++ b/deal.II/lac/include/lac/constraint_matrix.h
@@ -96,12 +96,20 @@ namespace internals
  *
  * <h3>Using the %ConstraintMatrix for Dirichlet boundary conditions</h3>
  *
- * The ConstraintMatrix provides an alternative for implementinging
- * Dirichlet boundary conditions (the standard way that is extensively
- * discussed in the tutorial programs is to use
- * MatrixTools::apply_boundary_values). The general principle of Dirichlet
- * conditions are algebraic constraints of the form $x_{i} = b_i$, which
- * fits into the form described above.
+ * The ConstraintMatrix provides an alternative for implementinging Dirichlet
+ * boundary conditions (the standard way that is extensively discussed in the
+ * tutorial programs is to use MatrixTools::apply_boundary_values). The
+ * general principle of Dirichlet conditions are algebraic constraints of the
+ * form $x_{i} = b_i$, which fits into the form described above. 
+ *
+ * Please note that when combining adaptively refined meshes with hanging node
+ * constraints and inhomogeneous Dirichlet boundary conditions within one
+ * ConstraintMatrix object, the hanging node constraints should always be set
+ * first, and then Dirichlet boundary conditions should be interpolated. This
+ * makes sure that the discretization remains H<sup>1</sup> conforming as is
+ * needed e.g. for the Laplace equation in 3D, as hanging nodes on the
+ * boundary should be still set to the weighted average of neighbors, and not
+ * the actual Dirichlet value.
  *
  *
  * <h3>Description of constraints</h3>
@@ -779,6 +787,14 @@ class ConstraintMatrix : public Subscriptor
     const std::vector<std::pair<unsigned int,double> >*
     get_constraint_entries (unsigned int line) const;
 
+                                     /**
+                                      * Returns the value of the inhomogeneity
+                                      * stored in the constrained dof @p
+                                      * line. Unconstrained dofs also return a
+                                      * zero value.
+                                      */
+    double get_inhomogeneity (unsigned int line) const;
+
 				     /**
 				      * Print the constraint lines. Mainly
 				      * for debugging purposes.
@@ -2100,6 +2116,18 @@ ConstraintMatrix::get_constraint_entries (unsigned int line) const
 
 
 
+inline
+double
+ConstraintMatrix::get_inhomogeneity (unsigned int line) const
+{
+  if (is_constrained(line))
+    return lines[lines_cache[calculate_line_index(line)]].inhomogeneity;
+  else
+    return 0;
+}
+
+
+
 inline unsigned int
 ConstraintMatrix::calculate_line_index (const unsigned int line) const
 {
-- 
2.39.5