]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add functions interpolate_boundary_values and project_boundary_values for a Constrain...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 10 Mar 2009 16:55:43 +0000 (16:55 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 10 Mar 2009 16:55:43 +0000 (16:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@18474 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/vectors.h
deal.II/deal.II/include/numerics/vectors.templates.h
deal.II/deal.II/source/numerics/vectors.cc

index 92c14cad838923b6f7acc645441ad7bc82a324db..e63d0b2e5f41feb06595ac3aea723e15fa2f8754 100644 (file)
@@ -682,53 +682,179 @@ class VectorTools
                               const typename FunctionMap<DH::space_dimension>::type &function_map,
                               std::map<unsigned int,double>         &boundary_values,
                               const std::vector<bool>               &component_mask = std::vector<bool>());
+
+
+                                    /**
+                                     * 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.
+                                     *
+                                     * The parameter @p boundary_component
+                                     * corresponds to the number @p
+                                     * boundary_indicator of the face.  255
+                                     * is an illegal value, since it is
+                                     * reserved for interior faces.
+                                     *
+                                     * The flags in the last parameter, @p
+                                     * component_mask denote which
+                                     * components of the finite element
+                                     * space shall be interpolated. If it
+                                     * is left as specified by the default
+                                     * value (i.e. an empty array), all
+                                     * components are interpolated. If it
+                                     * is different from the default value,
+                                     * it is assumed that the number of
+                                     * entries equals the number of
+                                     * components in the boundary functions
+                                     * and the finite element, and those
+                                     * components in the given boundary
+                                     * function will be used for which the
+                                     * respective flag was set in the
+                                     * component mask.
+                                     *
+                                     * It is assumed that the number of
+                                     * components of the function in @p
+                                     * boundary_function matches that of
+                                     * the finite element used by @p dof.
+                                     *
+                                     * If the finite element used has shape
+                                     * functions that are non-zero in more
+                                     * than one component (in deal.II
+                                     * speak: they are non-primitive), then
+                                     * these components can presently not
+                                     * be used for interpolating boundary
+                                     * values. Thus, the elements in the
+                                     * component mask corresponding to the
+                                     * components of these non-primitive
+                                     * shape functions must be @p false.
+                                     *
+                                     * See the general doc for more
+                                     * information.
+                                     */
+  template <class DH>
+  static
+  void
+  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 = std::vector<bool>());
   
+                                    /**
+                                     * @deprecated This function is there
+                                     * mainly for backward compatibility.
+                                     *
+                                     * Same function as above, but taking
+                                     * only one pair of boundary indicator
+                                     * and corresponding boundary
+                                     * function. Calls the other function
+                                     * with remapped arguments.
+                                     *
+                                     */
+  template <class DH>
+  static
+  void
+  interpolate_boundary_values (const Mapping<DH::dimension,DH::space_dimension> &mapping,
+                              const DH                            &dof,
+                              const unsigned char                  boundary_component,
+                              const Function<DH::space_dimension> &boundary_function,
+                              ConstraintMatrix                    &constraints,
+                              const std::vector<bool>             &component_mask = std::vector<bool>());
+
+                                    /**
+                                     * Calls the other
+                                     * interpolate_boundary_values()
+                                     * function, see above, with
+                                     * <tt>mapping=MappingQ1@<dim@>()</tt>.
+                                     */
+  template <class DH>
+  static
+  void
+  interpolate_boundary_values (const DH                            &dof,
+                              const unsigned char                  boundary_component,
+                              const Function<DH::space_dimension> &boundary_function,
+                              ConstraintMatrix                    &constraints,
+                              const std::vector<bool>             &component_mask = std::vector<bool>());
+
     
+                                    /**
+                                     * Calls the other
+                                     * interpolate_boundary_values()
+                                     * function, see above, with
+                                     * <tt>mapping=MappingQ1@<dim@>()</tt>.
+                                     */
+  template <class DH>
+  static
+  void
+  interpolate_boundary_values (const DH                &dof,
+                              const typename FunctionMap<DH::space_dimension>::type &function_map,
+                              ConstraintMatrix        &constraints,
+                              const std::vector<bool> &component_mask = std::vector<bool>());
+
+
                                     /**
                                      * Project a function to the boundary
-                                     * of the domain, using the given quadrature
-                                     * formula for the faces. If the
-                                     * @p boundary_values contained values
-                                     * before, the new ones are added, or
-                                     * the old one overwritten if a node
-                                     * of the boundary part to be projected
-                                     * on already was in the variable.
-                                     *
-                                     * If @p component_mapping is
-                                     * empty, it is assumed that the
-                                     * number of components of @p
-                                     * boundary_function matches that
-                                     * of the finite element used by
-                                     * @p dof.
+                                     * of the domain, using the given
+                                     * quadrature formula for the faces. If
+                                     * the @p boundary_values contained
+                                     * values before, the new ones are
+                                     * added, or the old one overwritten if
+                                     * a node of the boundary part to be
+                                     * projected on already was in the
+                                     * variable.
+                                     *
+                                     * If @p component_mapping is empty, it
+                                     * is assumed that the number of
+                                     * components of @p boundary_function
+                                     * matches that of the finite element
+                                     * used by @p dof.
                                      *
                                      * In 1d, projection equals
                                      * interpolation. Therefore,
                                      * interpolate_boundary_values is
                                      * called.
                                      *
-                                     * @arg @p boundary_values: the
-                                     * result of this function, a map
-                                     * containing all indices of
-                                     * degrees of freedom at the
-                                     * boundary (as covered by the
+                                     * @arg @p boundary_values: the result
+                                     * of this function, a map containing
+                                     * all indices of degrees of freedom at
+                                     * the boundary (as covered by the
                                      * boundary parts in @p
-                                     * boundary_functions) and the
-                                     * computed dof value for this
-                                     * degree of freedom.
-                                     *
-                                     * @arg @p component_mapping: if
-                                     * the components in @p
-                                     * boundary_functions and @p dof
-                                     * do not coincide, this vector
-                                     * allows them to be remapped. If
-                                     * the vector is not empty, it
-                                     * has to have one entry for each
-                                     * component in @p dof. This
-                                     * entry is the component number
-                                     * in @p boundary_functions that
-                                     * should be used for this
-                                     * component in @p dof. By
-                                     * default, no remapping is
+                                     * boundary_functions) and the computed
+                                     * dof value for this degree of
+                                     * freedom.
+                                     *
+                                     * @arg @p component_mapping: if the
+                                     * components in @p boundary_functions
+                                     * and @p dof do not coincide, this
+                                     * vector allows them to be
+                                     * remapped. If the vector is not
+                                     * empty, it has to have one entry for
+                                     * each component in @p dof. This entry
+                                     * is the component number in @p
+                                     * boundary_functions that should be
+                                     * used for this component in @p
+                                     * dof. By default, no remapping is
                                      * applied.
                                      */
   template <int dim, int spacedim>
@@ -751,6 +877,72 @@ class VectorTools
                                       std::map<unsigned int,double> &boundary_values,
                                       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.
+                                     *
+                                     * If @p component_mapping is empty, it
+                                     * is assumed that the number of
+                                     * components of @p boundary_function
+                                     * matches that of the finite element
+                                     * used by @p dof.
+                                     *
+                                     * In 1d, projection equals
+                                     * interpolation. Therefore,
+                                     * interpolate_boundary_values is
+                                     * called.
+                                     *
+                                     * @arg @p component_mapping: if the
+                                     * components in @p boundary_functions
+                                     * and @p dof do not coincide, this
+                                     * vector allows them to be
+                                     * remapped. If the vector is not
+                                     * empty, it has to have one entry for
+                                     * each component in @p dof. This entry
+                                     * is the component number in @p
+                                     * boundary_functions that should be
+                                     * used for this component in @p
+                                     * dof. By default, no remapping is
+                                     * applied.
+                                     */
+  template <int dim, int spacedim>
+  static void project_boundary_values (const Mapping<dim, spacedim>   &mapping,
+                                      const DoFHandler<dim,spacedim> &dof,
+                                      const typename FunctionMap<spacedim>::type &boundary_functions,
+                                      const Quadrature<dim-1>        &q,
+                                      ConstraintMatrix               &constraints,
+                                      std::vector<unsigned int>       component_mapping = std::vector<unsigned int>());
+
+                                    /**
+                                     * Calls the project_boundary_values()
+                                     * function, see above, with
+                                     * <tt>mapping=MappingQ1@<dim@>()</tt>.
+                                     */
+  template <int dim, int spacedim>
+  static void project_boundary_values (const DoFHandler<dim,spacedim> &dof,
+                                      const typename FunctionMap<spacedim>::type &boundary_function,
+                                      const Quadrature<dim-1>        &q,
+                                      ConstraintMatrix               &constraints,
+                                      std::vector<unsigned int>       component_mapping = std::vector<unsigned int>());
+
+
                                     /**
                                      * Compute the constraints that
                                      * correspond to boundary conditions of
index 1677cbd99f0badf038bdfc77f3ab2364b200f2c9..350ad1504fdc991aa2bcf559056ad4fd919d4a8f 100644 (file)
@@ -1213,6 +1213,7 @@ VectorTools::create_boundary_right_hand_side (const hp::DoFHandler<dim,spacedim>
 
 
 
+// ----------- interpolate_boundary_values for std::map --------------------
 
 #if deal_II_dimension == 1
 
@@ -1223,10 +1224,10 @@ VectorTools::create_boundary_right_hand_side (const hp::DoFHandler<dim,spacedim>
 
 template <class DH>
 void
-VectorTools::interpolate_boundary_values (const Mapping<DH::dimension, DH::space_dimension>         &,
+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,
+                                         const unsigned char    boundary_component,
+                                         const Function<DH::space_dimension> &boundary_function,
                                          std::map<unsigned int,double> &boundary_values,
                                          const std::vector<bool>       &component_mask_)
 {
@@ -1720,6 +1721,516 @@ 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);
+       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>
+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_)
+{
+  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)
+    {
+      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 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));
+        }    
+    }
+                                  // 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<dim> > &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]);
+                 }
+             }
+         }
+      }
+}
+
+
+
+template <class DH>
+void
+VectorTools::interpolate_boundary_values 
+  (const Mapping<DH::dimension, DH::space_dimension> &mapping,
+   const DH                                          &dof,
+   const unsigned char                                boundary_component,
+   const Function<DH::space_dimension>               &boundary_function,
+   ConstraintMatrix                                  &constraints,
+   const std::vector<bool>                           &component_mask)
+{
+  typename FunctionMap<DH::space_dimension>::type function_map;
+  function_map[boundary_component] = &boundary_function;
+  interpolate_boundary_values (mapping, dof, function_map, constraints,
+                              component_mask);
+}
+
+#endif
+  
+
+template <class DH>
+void
+VectorTools::interpolate_boundary_values 
+  (const DH                            &dof,
+   const unsigned char                  boundary_component,
+   const Function<DH::space_dimension> &boundary_function,
+   ConstraintMatrix                    &constraints,
+   const std::vector<bool>             &component_mask)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping, 
+                             dof, boundary_component,
+                             boundary_function, constraints, component_mask);
+}
+
+
+
+template <class DH>
+void
+VectorTools::interpolate_boundary_values 
+  (const DH                                              &dof,
+   const typename FunctionMap<DH::space_dimension>::type &function_map,
+   ConstraintMatrix                                      &constraints,
+   const std::vector<bool>                               &component_mask)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping, 
+                             dof, function_map,
+                             constraints, component_mask);
+}
+
+
+
+
+// -------- implementation for project_boundary_values with std::map --------
+
 #if deal_II_dimension == 1
 
 // Implementation for 1D
@@ -1916,6 +2427,210 @@ 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>
+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>  &q,
+                                     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)
+    {
+#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
+    }
+  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
+VectorTools::project_boundary_values (const DoFHandler<dim,spacedim>    &dof,
+                                     const typename FunctionMap<spacedim>::type &boundary_functions,
+                                     const Quadrature<dim-1>  &q,
+                                     ConstraintMatrix &constraints,
+                                     std::vector<unsigned int> component_mapping)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  project_boundary_values(StaticMappingQ1<dim,spacedim>::mapping, dof, boundary_functions, q,
+                         constraints, component_mapping);
+}
+
+
+
+
 namespace internal
 {
   namespace VectorTools
index 542ffbc45d1b109a5e43116dc3449ad7f0cdf9a8..0f1ffa8e85129de6ecb3d9404d246231a7339a98 100644 (file)
@@ -137,6 +137,30 @@ void VectorTools::interpolate_boundary_values (
   std::map<unsigned int,double>       &,
   const std::vector<bool>    &);
 
+template
+void VectorTools::interpolate_boundary_values ( 
+  const DoFHandler<deal_II_dimension> &,
+  const unsigned char,
+  const Function<deal_II_dimension>   &,
+  ConstraintMatrix                    &,
+  const std::vector<bool>    &);
+
+template
+void VectorTools::interpolate_boundary_values (
+  const hp::DoFHandler<deal_II_dimension> &,
+  const unsigned char,
+  const Function<deal_II_dimension>   &,
+  ConstraintMatrix                    &,
+  const std::vector<bool>    &);
+
+template
+void VectorTools::interpolate_boundary_values (
+  const MGDoFHandler<deal_II_dimension> &,
+  const unsigned char,
+  const Function<deal_II_dimension>   &,
+  ConstraintMatrix                    &,
+  const std::vector<bool>    &);
+
 template
 void VectorTools::project_boundary_values<deal_II_dimension>
 (const Mapping<deal_II_dimension>     &,
@@ -153,6 +177,22 @@ void VectorTools::project_boundary_values<deal_II_dimension>
  std::map<unsigned int,double>&, std::vector<unsigned int>);
 
 
+template
+void VectorTools::project_boundary_values<deal_II_dimension>
+(const Mapping<deal_II_dimension>     &,
+ const DoFHandler<deal_II_dimension>  &,
+ const FunctionMap<deal_II_dimension>::type &,
+ const Quadrature<deal_II_dimension-1>&,
+ ConstraintMatrix&, std::vector<unsigned int>);
+
+template
+void VectorTools::project_boundary_values<deal_II_dimension>
+(const DoFHandler<deal_II_dimension>  &,
+ const FunctionMap<deal_II_dimension>::type &,
+ const Quadrature<deal_II_dimension-1>&,
+ ConstraintMatrix&, std::vector<unsigned int>);
+
+
 #if deal_II_dimension != 1
 template
 void

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.