Allow a template VectorType in VectorTools::create_right_hand_side and VectorTools...

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 1a2aee72cb25fd54a6c7ac925ca9efed1b95abe7..c2fc5e7416546d60b8567ade82dcabb2ba190a59 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -38,6 +38,13 @@ inconvenience this causes.
 </p>
 
 <ol>
+ <li> Changed: VectorTools::create_right_hand_side and
+ VectorTools::create_boundary_right_hand_side now take an additional template
+ parameter VectorType.
+ <br>
+ (Daniel Arndt, 2016/10/25)
+ </li>
+
  <li> Deprecated: ParameterHandler::read_input,
  ParameterHandler::read_input_from_xml, and
  ParameterHandler::read_input_from_string are now deprecated in favor of
@@ -495,7 +502,7 @@ inconvenience this causes.
  (Rajat Arora, 2016/10/29)
  </li>
 
-<li> New: Add MatrixFreeOperators::MassOperator representing a mass matrix.
+ <li> New: Add MatrixFreeOperators::MassOperator representing a mass matrix.
  <br>
  (Daniel Arndt, 2016/10/27)
  </li>
@@ -515,6 +522,12 @@ inconvenience this causes.
  (Martin Kronbichler, 2016/10/27)
  </li>
 
+ <li> >New: VectorTools::create_right_hand_side can be used for parallel
+ Triangulations and homogeneous constraints using a ConstraintMatrix.
+ <br>
+ (Daniel Arndt, 2016/10/25)
+ </li>
+
  <li> Fixed: GridIn::read_vtk() had off-by-one errors in reading face
  boundary indicators from VTK files. Consequently, not all boundary indicators
  were correctly set.

diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h
index 7071137f61ee188c84a105243dd14da25709aaa9..8509ed417b7c0fe86259f52b114a731c1b7d6a85 100644 (file)
--- a/include/deal.II/numerics/vector_tools.h
+++ b/include/deal.II/numerics/vector_tools.h
@@ -739,7 +739,6 @@ namespace VectorTools
                                  const ConstraintMatrix              &constraints,
                                  VectorType                          &u2);
 
-
   /**
    * The same function as above, but takes an InterGridMap object directly as
    * a parameter. Useful for interpolating several vectors at the same time.
@@ -1925,41 +1924,45 @@ namespace VectorTools
    *
    * See the general documentation of this namespace for further information.
    */
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_right_hand_side (const Mapping<dim, spacedim>    &mapping,
                                const DoFHandler<dim,spacedim>  &dof,
                                const Quadrature<dim>           &q,
-                               const Function<spacedim,double> &rhs,
-                               Vector<double>                  &rhs_vector);
+                               const Function<spacedim,typename VectorType::value_type> &rhs,
+                               VectorType                      &rhs_vector,
+                               const ConstraintMatrix          &constraints=ConstraintMatrix());
 
   /**
    * Call the create_right_hand_side() function, see above, with
    * <tt>mapping=MappingQGeneric@<dim@>(1)</tt>.
    */
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_right_hand_side (const DoFHandler<dim,spacedim>  &dof,
                                const Quadrature<dim>           &q,
-                               const Function<spacedim,double> &rhs,
-                               Vector<double>                  &rhs_vector);
+                               const Function<spacedim,typename VectorType::value_type> &rhs,
+                               VectorType                      &rhs_vector,
+                               const ConstraintMatrix          &constraints=ConstraintMatrix());
 
   /**
    * Like the previous set of functions, but for hp objects.
    */
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_right_hand_side (const hp::MappingCollection<dim,spacedim> &mapping,
                                const hp::DoFHandler<dim,spacedim>        &dof,
                                const hp::QCollection<dim>                &q,
-                               const Function<spacedim,double>           &rhs,
-                               Vector<double>                            &rhs_vector);
+                               const Function<spacedim,typename VectorType::value_type> &rhs,
+                               VectorType                                &rhs_vector,
+                               const ConstraintMatrix                    &constraints=ConstraintMatrix());
 
   /**
    * Like the previous set of functions, but for hp objects.
    */
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_right_hand_side (const hp::DoFHandler<dim,spacedim> &dof,
                                const hp::QCollection<dim>         &q,
-                               const Function<spacedim,double>    &rhs,
-                               Vector<double>                     &rhs_vector);
+                               const Function<spacedim,typename VectorType::value_type> &rhs,
+                               VectorType                         &rhs_vector,
+                               const ConstraintMatrix             &constraints=ConstraintMatrix());
 
   /**
    * Create a right hand side vector for a point source at point @p p. In
@@ -2069,12 +2072,12 @@ namespace VectorTools
    * @see
    * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    */
-  template <int dim, int spacedim>
-  void create_boundary_right_hand_side (const Mapping<dim,spacedim>      &mapping,
-                                        const DoFHandler<dim,spacedim>   &dof,
-                                        const Quadrature<dim-1> &q,
-                                        const Function<spacedim,double>     &rhs,
-                                        Vector<double>          &rhs_vector,
+  template <int dim, int spacedim, typename VectorType>
+  void create_boundary_right_hand_side (const Mapping<dim,spacedim>        &mapping,
+                                        const DoFHandler<dim,spacedim>     &dof,
+                                        const Quadrature<dim-1>            &q,
+                                        const Function<spacedim,typename VectorType::value_type> &rhs,
+                                        VectorType                         &rhs_vector,
                                         const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());
 
   /**
@@ -2084,12 +2087,12 @@ namespace VectorTools
    * @see
    * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    */
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_boundary_right_hand_side
   (const DoFHandler<dim,spacedim>     &dof,
    const Quadrature<dim-1>            &q,
-   const Function<spacedim,double>    &rhs,
-   Vector<double>                     &rhs_vector,
+   const Function<spacedim,typename VectorType::value_type> &rhs,
+   VectorType                         &rhs_vector,
    const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());
 
   /**
@@ -2098,13 +2101,13 @@ namespace VectorTools
    * @see
    * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    */
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_boundary_right_hand_side
   (const hp::MappingCollection<dim,spacedim> &mapping,
    const hp::DoFHandler<dim,spacedim>        &dof,
    const hp::QCollection<dim-1>              &q,
-   const Function<spacedim,double>           &rhs,
-   Vector<double>                            &rhs_vector,
+   const Function<spacedim,typename VectorType::value_type> &rhs,
+   VectorType                                &rhs_vector,
    const std::set<types::boundary_id>        &boundary_ids = std::set<types::boundary_id>());
 
   /**
@@ -2115,12 +2118,12 @@ namespace VectorTools
    * @see
    * @ref GlossBoundaryIndicator "Glossary entry on boundary indicators"
    */
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_boundary_right_hand_side
   (const hp::DoFHandler<dim,spacedim> &dof,
    const hp::QCollection<dim-1>       &q,
-   const Function<spacedim,double>    &rhs,
-   Vector<double>                     &rhs_vector,
+   const Function<spacedim,typename VectorType::value_type> &rhs,
+   VectorType                         &rhs_vector,
    const std::set<types::boundary_id> &boundary_ids = std::set<types::boundary_id>());
 
   //@}

diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 6f04aa2ace457292f610d8536da7731a646f9b93..b6bfa29e043012d97b9961500efb06872742d924 100644 (file)
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -777,6 +777,8 @@ namespace VectorTools
       return true;
     }
 
+
+
     template <typename number>
     void invert_mass_matrix(const SparseMatrix<number> &mass_matrix,
                             const Vector<number>       &rhs,
@@ -804,6 +806,7 @@ namespace VectorTools
     }
 
 
+
     /**
      * Generic implementation of the project() function
      */
@@ -1178,6 +1181,7 @@ namespace VectorTools
   }
 
 
+
   template <int dim, typename VectorType, int spacedim>
   void project (const DoFHandler<dim,spacedim> &dof,
                 const ConstraintMatrix         &constraints,
@@ -1228,13 +1232,17 @@ namespace VectorTools
   }
 
 
-  template <int dim, int spacedim>
+
+  template <int dim, int spacedim, typename VectorType>
   void create_right_hand_side (const Mapping<dim, spacedim>   &mapping,
                                const DoFHandler<dim,spacedim> &dof_handler,
                                const Quadrature<dim>          &quadrature,
-                               const Function<spacedim>       &rhs_function,
-                               Vector<double>                 &rhs_vector)
+                               const Function<spacedim, typename VectorType::value_type> &rhs_function,
+                               VectorType                     &rhs_vector,
+                               const ConstraintMatrix         &constraints)
   {
+    typedef typename VectorType::value_type Number;
+
     const FiniteElement<dim,spacedim> &fe  = dof_handler.get_fe();
     Assert (fe.n_components() == rhs_function.n_components,
             ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
@@ -1252,7 +1260,7 @@ namespace VectorTools
                        n_components  = fe.n_components();
 
     std::vector<types::global_dof_index> dofs (dofs_per_cell);
-    Vector<double> cell_vector (dofs_per_cell);
+    Vector<Number> cell_vector (dofs_per_cell);
 
     typename DoFHandler<dim,spacedim>::active_cell_iterator
     cell = dof_handler.begin_active(),
@@ -1260,105 +1268,107 @@ namespace VectorTools
 
     if (n_components==1)
       {
-        std::vector<double> rhs_values(n_q_points);
+        std::vector<Number> rhs_values(n_q_points);
 
         for (; cell!=endc; ++cell)
-          {
-            fe_values.reinit(cell);
-
-            const std::vector<double> &weights   = fe_values.get_JxW_values ();
-            rhs_function.value_list (fe_values.get_quadrature_points(),
-                                     rhs_values);
+          if (cell->is_locally_owned())
+            {
+              fe_values.reinit(cell);
 
-            cell_vector = 0;
-            for (unsigned int point=0; point<n_q_points; ++point)
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                cell_vector(i) += rhs_values[point] *
-                                  fe_values.shape_value(i,point) *
-                                  weights[point];
+              const std::vector<double> &weights   = fe_values.get_JxW_values ();
+              rhs_function.value_list (fe_values.get_quadrature_points(),
+                                       rhs_values);
 
-            cell->get_dof_indices (dofs);
+              cell_vector = 0;
+              for (unsigned int point=0; point<n_q_points; ++point)
+                for (unsigned int i=0; i<dofs_per_cell; ++i)
+                  cell_vector(i) += rhs_values[point] *
+                                    fe_values.shape_value(i,point) *
+                                    weights[point];
 
-            for (unsigned int i=0; i<dofs_per_cell; ++i)
-              rhs_vector(dofs[i]) += cell_vector(i);
-          }
+              cell->get_dof_indices(dofs);
 
+              constraints.distribute_local_to_global(cell_vector, dofs, rhs_vector);
+            }
       }
     else
       {
-        std::vector<Vector<double> > rhs_values(n_q_points,
-                                                Vector<double>(n_components));
+        std::vector<Vector<Number> > rhs_values(n_q_points,
+                                                Vector<Number>(n_components));
 
         for (; cell!=endc; ++cell)
-          {
-            fe_values.reinit(cell);
-
-            const std::vector<double> &weights   = fe_values.get_JxW_values ();
-            rhs_function.vector_value_list (fe_values.get_quadrature_points(),
-                                            rhs_values);
+          if (cell->is_locally_owned())
+            {
+              fe_values.reinit(cell);
 
-            cell_vector = 0;
-            // Use the faster code if the
-            // FiniteElement is primitive
-            if (fe.is_primitive ())
-              {
-                for (unsigned int point=0; point<n_q_points; ++point)
-                  for (unsigned int i=0; i<dofs_per_cell; ++i)
-                    {
-                      const unsigned int component
-                        = fe.system_to_component_index(i).first;
+              const std::vector<double> &weights   = fe_values.get_JxW_values ();
+              rhs_function.vector_value_list (fe_values.get_quadrature_points(),
+                                              rhs_values);
 
-                      cell_vector(i) += rhs_values[point](component) *
-                                        fe_values.shape_value(i,point) *
-                                        weights[point];
-                    }
-              }
-            else
-              {
-                // Otherwise do it the way
-                // proposed for vector valued
-                // elements
-                for (unsigned int point=0; point<n_q_points; ++point)
-                  for (unsigned int i=0; i<dofs_per_cell; ++i)
-                    for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                      if (fe.get_nonzero_components(i)[comp_i])
-                        {
-                          cell_vector(i) += rhs_values[point](comp_i) *
-                                            fe_values.shape_value_component(i,point,comp_i) *
-                                            weights[point];
-                        }
-              }
+              cell_vector = 0;
+              // Use the faster code if the
+              // FiniteElement is primitive
+              if (fe.is_primitive ())
+                {
+                  for (unsigned int point=0; point<n_q_points; ++point)
+                    for (unsigned int i=0; i<dofs_per_cell; ++i)
+                      {
+                        const unsigned int component
+                          = fe.system_to_component_index(i).first;
 
-            cell->get_dof_indices (dofs);
+                        cell_vector(i) += rhs_values[point](component) *
+                                          fe_values.shape_value(i,point) *
+                                          weights[point];
+                      }
+                }
+              else
+                {
+                  // Otherwise do it the way
+                  // proposed for vector valued
+                  // elements
+                  for (unsigned int point=0; point<n_q_points; ++point)
+                    for (unsigned int i=0; i<dofs_per_cell; ++i)
+                      for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                        if (fe.get_nonzero_components(i)[comp_i])
+                          {
+                            cell_vector(i) += rhs_values[point](comp_i) *
+                                              fe_values.shape_value_component(i,point,comp_i) *
+                                              weights[point];
+                          }
+                }
+              cell->get_dof_indices (dofs);
 
-            for (unsigned int i=0; i<dofs_per_cell; ++i)
-              rhs_vector(dofs[i]) += cell_vector(i);
-          }
+              constraints.distribute_local_to_global(cell_vector, dofs, rhs_vector);
+            }
       }
   }
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_right_hand_side (const DoFHandler<dim,spacedim> &dof_handler,
                                const Quadrature<dim>          &quadrature,
-                               const Function<spacedim>       &rhs_function,
-                               Vector<double>                 &rhs_vector)
+                               const Function<spacedim, typename VectorType::value_type> &rhs_function,
+                               VectorType                     &rhs_vector,
+                               const ConstraintMatrix         &constraints)
   {
     create_right_hand_side(StaticMappingQ1<dim,spacedim>::mapping, dof_handler, quadrature,
-                           rhs_function, rhs_vector);
+                           rhs_function, rhs_vector, constraints);
   }
 
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_right_hand_side (const hp::MappingCollection<dim,spacedim> &mapping,
                                const hp::DoFHandler<dim,spacedim>        &dof_handler,
                                const hp::QCollection<dim>                &quadrature,
-                               const Function<spacedim>                  &rhs_function,
-                               Vector<double>                            &rhs_vector)
+                               const Function<spacedim, typename VectorType::value_type> &rhs_function,
+                               VectorType                                &rhs_vector,
+                               const ConstraintMatrix                    &constraints)
   {
+    typedef typename VectorType::value_type Number;
+
     const hp::FECollection<dim,spacedim> &fe  = dof_handler.get_fe();
     Assert (fe.n_components() == rhs_function.n_components,
             ExcDimensionMismatch(fe.n_components(), rhs_function.n_components));
@@ -1374,7 +1384,7 @@ namespace VectorTools
     const unsigned int n_components  = fe.n_components();
 
     std::vector<types::global_dof_index> dofs (fe.max_dofs_per_cell());
-    Vector<double> cell_vector (fe.max_dofs_per_cell());
+    Vector<Number> cell_vector (fe.max_dofs_per_cell());
 
     typename hp::DoFHandler<dim,spacedim>::active_cell_iterator
     cell = dof_handler.begin_active(),
@@ -1382,115 +1392,115 @@ namespace VectorTools
 
     if (n_components==1)
       {
-        std::vector<double> rhs_values;
+        std::vector<Number> rhs_values;
 
         for (; cell!=endc; ++cell)
-          {
-            x_fe_values.reinit(cell);
+          if (cell->is_locally_owned())
+            {
+              x_fe_values.reinit(cell);
 
-            const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
+              const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
 
-            const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                               n_q_points    = fe_values.n_quadrature_points;
-            rhs_values.resize (n_q_points);
-            dofs.resize (dofs_per_cell);
-            cell_vector.reinit (dofs_per_cell);
+              const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                                 n_q_points    = fe_values.n_quadrature_points;
+              rhs_values.resize (n_q_points);
+              dofs.resize (dofs_per_cell);
+              cell_vector.reinit (dofs_per_cell);
 
-            const std::vector<double> &weights   = fe_values.get_JxW_values ();
-            rhs_function.value_list (fe_values.get_quadrature_points(),
-                                     rhs_values);
+              const std::vector<Number> &weights   = fe_values.get_JxW_values ();
+              rhs_function.value_list (fe_values.get_quadrature_points(),
+                                       rhs_values);
 
-            cell_vector = 0;
-            for (unsigned int point=0; point<n_q_points; ++point)
-              for (unsigned int i=0; i<dofs_per_cell; ++i)
-                cell_vector(i) += rhs_values[point] *
-                                  fe_values.shape_value(i,point) *
-                                  weights[point];
+              cell_vector = 0;
+              for (unsigned int point=0; point<n_q_points; ++point)
+                for (unsigned int i=0; i<dofs_per_cell; ++i)
+                  cell_vector(i) += rhs_values[point] *
+                                    fe_values.shape_value(i,point) *
+                                    weights[point];
 
-            cell->get_dof_indices (dofs);
+              cell->get_dof_indices (dofs);
 
-            for (unsigned int i=0; i<dofs_per_cell; ++i)
-              rhs_vector(dofs[i]) += cell_vector(i);
-          }
+              constraints.distribute_local_to_global(cell_vector, dofs, rhs_vector);
+            }
 
       }
     else
       {
-        std::vector<Vector<double> > rhs_values;
+        std::vector<Vector<Number> > rhs_values;
 
         for (; cell!=endc; ++cell)
-          {
-            x_fe_values.reinit(cell);
+          if (cell->is_locally_owned())
+            {
+              x_fe_values.reinit(cell);
 
-            const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
+              const FEValues<dim,spacedim> &fe_values = x_fe_values.get_present_fe_values();
 
-            const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
-                               n_q_points    = fe_values.n_quadrature_points;
-            rhs_values.resize (n_q_points,
-                               Vector<double>(n_components));
-            dofs.resize (dofs_per_cell);
-            cell_vector.reinit (dofs_per_cell);
+              const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                                 n_q_points    = fe_values.n_quadrature_points;
+              rhs_values.resize (n_q_points,
+                                 Vector<Number>(n_components));
+              dofs.resize (dofs_per_cell);
+              cell_vector.reinit (dofs_per_cell);
 
-            const std::vector<double> &weights   = fe_values.get_JxW_values ();
-            rhs_function.vector_value_list (fe_values.get_quadrature_points(),
-                                            rhs_values);
+              const std::vector<Number> &weights   = fe_values.get_JxW_values ();
+              rhs_function.vector_value_list (fe_values.get_quadrature_points(),
+                                              rhs_values);
 
-            cell_vector = 0;
+              cell_vector = 0;
 
-            // Use the faster code if the
-            // FiniteElement is primitive
-            if (cell->get_fe().is_primitive ())
-              {
-                for (unsigned int point=0; point<n_q_points; ++point)
-                  for (unsigned int i=0; i<dofs_per_cell; ++i)
-                    {
-                      const unsigned int component
-                        = cell->get_fe().system_to_component_index(i).first;
+              // Use the faster code if the
+              // FiniteElement is primitive
+              if (cell->get_fe().is_primitive ())
+                {
+                  for (unsigned int point=0; point<n_q_points; ++point)
+                    for (unsigned int i=0; i<dofs_per_cell; ++i)
+                      {
+                        const unsigned int component
+                          = cell->get_fe().system_to_component_index(i).first;
 
-                      cell_vector(i) += rhs_values[point](component) *
-                                        fe_values.shape_value(i,point) *
-                                        weights[point];
-                    }
-              }
-            else
-              {
-                // Otherwise do it the way proposed
-                // for vector valued elements
-                for (unsigned int point=0; point<n_q_points; ++point)
-                  for (unsigned int i=0; i<dofs_per_cell; ++i)
-                    for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
-                      if (cell->get_fe().get_nonzero_components(i)[comp_i])
-                        {
-                          cell_vector(i) += rhs_values[point](comp_i) *
-                                            fe_values.shape_value_component(i,point,comp_i) *
-                                            weights[point];
-                        }
-              }
+                        cell_vector(i) += rhs_values[point](component) *
+                                          fe_values.shape_value(i,point) *
+                                          weights[point];
+                      }
+                }
+              else
+                {
+                  // Otherwise do it the way proposed
+                  // for vector valued elements
+                  for (unsigned int point=0; point<n_q_points; ++point)
+                    for (unsigned int i=0; i<dofs_per_cell; ++i)
+                      for (unsigned int comp_i = 0; comp_i < n_components; ++comp_i)
+                        if (cell->get_fe().get_nonzero_components(i)[comp_i])
+                          {
+                            cell_vector(i) += rhs_values[point](comp_i) *
+                                              fe_values.shape_value_component(i,point,comp_i) *
+                                              weights[point];
+                          }
+                }
 
-            cell->get_dof_indices (dofs);
+              cell->get_dof_indices (dofs);
 
-            for (unsigned int i=0; i<dofs_per_cell; ++i)
-              rhs_vector(dofs[i]) += cell_vector(i);
-          }
+              constraints.distribute_local_to_global(cell_vector, dofs, rhs_vector);
+            }
       }
   }
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void create_right_hand_side (const hp::DoFHandler<dim,spacedim> &dof_handler,
                                const hp::QCollection<dim>         &quadrature,
-                               const Function<spacedim>           &rhs_function,
-                               Vector<double>                     &rhs_vector)
+                               const Function<spacedim, typename VectorType::value_type> &rhs_function,
+                               VectorType                         &rhs_vector,
+                               const ConstraintMatrix             &constraints)
   {
     create_right_hand_side(hp::StaticMappingQ1<dim,spacedim>::mapping_collection,
                            dof_handler, quadrature,
-                           rhs_function, rhs_vector);
+                           rhs_function, rhs_vector, constraints);
   }
 
 
 
-
   template <int dim, int spacedim>
   void create_point_source_vector (const Mapping<dim, spacedim>   &mapping,
                                    const DoFHandler<dim,spacedim> &dof_handler,
@@ -1678,13 +1688,13 @@ namespace VectorTools
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void
   create_boundary_right_hand_side (const Mapping<dim, spacedim>       &mapping,
                                    const DoFHandler<dim,spacedim>     &dof_handler,
                                    const Quadrature<dim-1>            &quadrature,
-                                   const Function<spacedim>           &rhs_function,
-                                   Vector<double>                     &rhs_vector,
+                                   const Function<spacedim,typename VectorType::value_type> &rhs_function,
+                                   VectorType                         &rhs_vector,
                                    const std::set<types::boundary_id> &boundary_ids)
   {
     const FiniteElement<dim> &fe  = dof_handler.get_fe();
@@ -1801,12 +1811,12 @@ namespace VectorTools
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void
   create_boundary_right_hand_side (const DoFHandler<dim,spacedim>     &dof_handler,
                                    const Quadrature<dim-1>            &quadrature,
-                                   const Function<spacedim>           &rhs_function,
-                                   Vector<double>                     &rhs_vector,
+                                   const Function<spacedim,typename VectorType::value_type> &rhs_function,
+                                   VectorType                         &rhs_vector,
                                    const std::set<types::boundary_id> &boundary_ids)
   {
     create_boundary_right_hand_side(StaticMappingQ1<dim>::mapping, dof_handler,
@@ -1817,13 +1827,13 @@ namespace VectorTools
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void
   create_boundary_right_hand_side (const hp::MappingCollection<dim,spacedim> &mapping,
                                    const hp::DoFHandler<dim,spacedim>        &dof_handler,
                                    const hp::QCollection<dim-1>              &quadrature,
-                                   const Function<spacedim>                  &rhs_function,
-                                   Vector<double>                            &rhs_vector,
+                                   const Function<spacedim,typename VectorType::value_type> &rhs_function,
+                                   VectorType                                &rhs_vector,
                                    const std::set<types::boundary_id>        &boundary_ids)
   {
     const hp::FECollection<dim> &fe  = dof_handler.get_fe();
@@ -1952,12 +1962,12 @@ namespace VectorTools
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename VectorType>
   void
   create_boundary_right_hand_side (const hp::DoFHandler<dim,spacedim> &dof_handler,
                                    const hp::QCollection<dim-1>       &quadrature,
-                                   const Function<spacedim>           &rhs_function,
-                                   Vector<double>                     &rhs_vector,
+                                   const Function<spacedim,typename VectorType::value_type> &rhs_function,
+                                   VectorType                         &rhs_vector,
                                    const std::set<types::boundary_id> &boundary_ids)
   {
     create_boundary_right_hand_side(hp::StaticMappingQ1<dim>::mapping_collection,
@@ -2495,7 +2505,7 @@ namespace VectorTools
       // 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);
+      ReductionControl        control(5.*rhs.size(), 0., 1.e-12, false, false);
       GrowingVectorMemory<Vector<number> > memory;
       SolverCG<Vector<number> >            cg(control,memory);
 

diff --git a/source/numerics/vector_tools_rhs.inst.in b/source/numerics/vector_tools_rhs.inst.in
index 2c38f5a6d248d38dc74da78cc6d7caef055b1328..2de583e733bc5fceb956bfa93c0d188be40dbcfe 100644 (file)
--- a/source/numerics/vector_tools_rhs.inst.in
+++ b/source/numerics/vector_tools_rhs.inst.in
@@ -14,25 +14,26 @@
 // ---------------------------------------------------------------------
 
 
-for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
+for (VEC : REAL_SERIAL_VECTORS; deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS)
 {
 #if deal_II_dimension <= deal_II_space_dimension
     namespace VectorTools \{
     template
-    void create_right_hand_side<deal_II_dimension,deal_II_space_dimension>
-    (const Mapping<deal_II_dimension,deal_II_space_dimension>    &,
-     const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
-     const Quadrature<deal_II_dimension> &,
-     const Function<deal_II_space_dimension>   &,
-     Vector<double>                      &);
+    void create_right_hand_side<deal_II_dimension,deal_II_space_dimension,VEC >
+    (const Mapping<deal_II_dimension,deal_II_space_dimension>          &,
+     const DoFHandler<deal_II_dimension,deal_II_space_dimension>       &,
+     const Quadrature<deal_II_dimension>                               &,
+     const Function<deal_II_space_dimension, typename VEC::value_type> &,
+     VEC                                                               &,
+     const ConstraintMatrix                                            &);
 
     template
-    void create_right_hand_side<deal_II_dimension,deal_II_space_dimension>
-    (const DoFHandler<deal_II_dimension,deal_II_space_dimension> &,
-     const Quadrature<deal_II_dimension> &,
-     const Function<deal_II_space_dimension>   &,
-     Vector<double>  &);
-
+    void create_right_hand_side<deal_II_dimension,deal_II_space_dimension,VEC >
+    (const DoFHandler<deal_II_dimension,deal_II_space_dimension>       &,
+     const Quadrature<deal_II_dimension>                               &,
+     const Function<deal_II_space_dimension, typename VEC::value_type> &,
+     VEC                                                               &,
+     const ConstraintMatrix                                            &);
     \}
 #endif
 }
@@ -45,61 +46,63 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
 
 #if deal_II_dimension == deal_II_space_dimension
     template
-    void create_right_hand_side<deal_II_dimension>
-    (const hp::MappingCollection<deal_II_dimension>    &,
-     const hp::DoFHandler<deal_II_dimension> &,
-     const hp::QCollection<deal_II_dimension> &,
-     const Function<deal_II_dimension>   &,
-     Vector<double>                      &);
+    void create_right_hand_side<deal_II_dimension, deal_II_space_dimension, Vector<double> >
+    (const hp::MappingCollection<deal_II_dimension> &,
+     const hp::DoFHandler<deal_II_dimension>        &,
+     const hp::QCollection<deal_II_dimension>       &,
+     const Function<deal_II_space_dimension, typename Vector<double>::value_type> &,
+     Vector<double>                                 &,
+     const ConstraintMatrix                         &);
 
     template
-    void create_right_hand_side<deal_II_dimension>
-    (const hp::DoFHandler<deal_II_dimension> &,
+    void create_right_hand_side<deal_II_dimension, deal_II_space_dimension, Vector<double> >
+    (const hp::DoFHandler<deal_II_dimension>  &,
      const hp::QCollection<deal_II_dimension> &,
-     const Function<deal_II_dimension>   &,
-     Vector<double>                      &);
+     const Function<deal_II_space_dimension, typename Vector<double>::value_type> &,
+     Vector<double>                           &,
+     const ConstraintMatrix                   &);
 
 #if deal_II_dimension > 1
     template
     void
-    create_boundary_right_hand_side<deal_II_dimension>
-    (const Mapping<deal_II_dimension>    &,
-     const DoFHandler<deal_II_dimension> &,
+    create_boundary_right_hand_side<deal_II_dimension, deal_II_space_dimension, Vector<double> >
+    (const Mapping<deal_II_dimension>      &,
+     const DoFHandler<deal_II_dimension>   &,
      const Quadrature<deal_II_dimension-1> &,
-     const Function<deal_II_dimension>   &,
-     Vector<double>                      &,
-     const std::set<types::boundary_id> &);
+     const Function<deal_II_space_dimension, typename Vector<double>::value_type> &,
+     Vector<double>                        &,
+     const std::set<types::boundary_id>    &);
 #endif
 
     template
     void
-    create_boundary_right_hand_side<deal_II_dimension>
-    (const DoFHandler<deal_II_dimension> &,
+    create_boundary_right_hand_side<deal_II_dimension, deal_II_space_dimension, Vector<double> >
+    (const DoFHandler<deal_II_dimension>   &,
      const Quadrature<deal_II_dimension-1> &,
-     const Function<deal_II_dimension>   &,
-     Vector<double>                      &,
-     const std::set<types::boundary_id> &);
+     const Function<deal_II_space_dimension, typename Vector<double>::value_type> &,
+     Vector<double>                        &,
+     const std::set<types::boundary_id>    &);
 
 #if deal_II_dimension > 1
     template
     void
-    create_boundary_right_hand_side<deal_II_dimension>
-    (const hp::MappingCollection<deal_II_dimension>    &,
-     const hp::DoFHandler<deal_II_dimension> &,
-     const hp::QCollection<deal_II_dimension-1> &,
-     const Function<deal_II_dimension>   &,
-     Vector<double>                      &,
-     const std::set<types::boundary_id> &);
+    create_boundary_right_hand_side<deal_II_dimension, deal_II_space_dimension, Vector<double> >
+    (const hp::MappingCollection<deal_II_dimension> &,
+     const hp::DoFHandler<deal_II_dimension>        &,
+     const hp::QCollection<deal_II_dimension-1>     &,
+     const Function<deal_II_space_dimension, typename Vector<double>::value_type> &,
+     Vector<double>                                 &,
+     const std::set<types::boundary_id>             &);
 #endif
 
     template
     void
-    create_boundary_right_hand_side<deal_II_dimension>
-    (const hp::DoFHandler<deal_II_dimension> &,
+    create_boundary_right_hand_side<deal_II_dimension, deal_II_space_dimension, Vector<double> >
+    (const hp::DoFHandler<deal_II_dimension>    &,
      const hp::QCollection<deal_II_dimension-1> &,
-     const Function<deal_II_dimension>   &,
-     Vector<double>                      &,
-     const std::set<types::boundary_id> &);
+     const Function<deal_II_space_dimension, typename Vector<double>::value_type> &,
+     Vector<double>                             &,
+     const std::set<types::boundary_id>         &);
 
 
 #endif