]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Fix mixing interpolate_boundary_values with other types of constraints as e.g. hangin...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 Sep 2010 14:43:35 +0000 (14:43 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 Sep 2010 14:43:35 +0000 (14:43 +0000)
git-svn-id: https://svn.dealii.org/trunk@21992 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/include/numerics/vectors.templates.h
deal.II/doc/news/changes.h
deal.II/examples/step-22/step-22.cc
deal.II/lac/include/lac/constraint_matrix.h

index 2cec72f3becfddd2dfa5a9a22ddf818609081ba3..ccfa873e2ba9b9739b4277aacc45eba70c93f13b 100644 (file)
@@ -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
index f6d43aaf55d0673c5af767479a8eabc0c5478d97..c975d44e566c97f28435dc8edd5bee0e9e4ef9da 100644 (file)
@@ -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
index fb7822eafbc20421416ecd9c5a76260ae899dce3..f03a4e1be85b086309af25a44e1d248b5ad523ac 100644 (file)
@@ -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
index d509a89f398cef82647c85eafb3a432924abab8b..9aa1a95cb0a9c8adab4e1d5c2ba87abd1b329f1b 100644 (file)
@@ -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 ();
index 388c1290bc45d3cc57553ef26cbbae76053930d1..6e454c7370a06ec9575ed42de59cc585ff921743 100644 (file)
@@ -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
 {

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.