VectorTools::create_boundary_right_hand_side

git-svn-id: https://svn.dealii.org/trunk@5002 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/deal.II/include/numerics/vectors.h b/deal.II/deal.II/include/numerics/vectors.h
index 49de1319746964fd2fbd6d9219b6f9d1c36a4b07..b3c947646220115265a0cdadb69e7d4cbe1edc28 100644 (file)
--- a/deal.II/deal.II/include/numerics/vectors.h
+++ b/deal.II/deal.II/include/numerics/vectors.h
@@ -14,6 +14,7 @@
 #include <dofs/function_map.h>
 #include <map>
 #include <vector>
+#include <set>
 
 template <int dim> class Function;
 template <int dim> class FunctionMap;
@@ -175,6 +176,15 @@ enum NormType {
  *   @ref{MatrixCreator}@p{::create_*} functions which take a right hand side do,
  *   but without assembling a matrix.
  *
+ * @item Creation of boundary right hand side vectors: The
+ *   @p{create_boundary_right_hand_side} function computes the vector
+ *   $f_i = \int_{\partial\Omega} g(x) \phi_i(x) dx$. This is the
+ *   right hand side contribution of boundary forces when having
+ *   inhomogeneous Neumann boundary values in Laplace's equation or
+ *   other second order operators. This function also takes an
+ *   optional argument denoting over which parts of the boundary the
+ *   integration shall extend.
+ *
  * @item Interpolation of boundary values:
  *   The @ref{MatrixTools}@p{::apply_boundary_values} function takes a list
  *   of boundary nodes and their values. You can get such a list by interpolation
@@ -438,7 +448,10 @@ class VectorTools
                         const bool                project_to_boundary_first = false);
 
                                     /**
-                                     * Create a right hand side vector.
+                                     * Create a right hand side
+                                     * vector. Prior content of the
+                                     * given @p{rhs_vector} vector is
+                                     * deleted.
                                      *
                                      * See the general documentation of this
                                      * class for further information.
@@ -461,6 +474,51 @@ class VectorTools
                                        const Function<dim>   &rhs,
                                        Vector<double>        &rhs_vector);
 
+                                    /**
+                                     * Create a right hand side
+                                     * vector from boundary
+                                     * forces. Prior content of the
+                                     * given @p{rhs_vector} vector is
+                                     * deleted.
+                                     *
+                                     * See the general documentation of this
+                                     * class for further information.
+                                     */
+    template <int dim>
+    static void create_boundary_right_hand_side (const Mapping<dim>      &mapping,
+                                                const DoFHandler<dim>   &dof,
+                                                const Quadrature<dim-1> &q,
+                                                const Function<dim>     &rhs,
+                                                Vector<double>          &rhs_vector,
+                                                const std::set<unsigned char> &boundary_indicators = std::set<unsigned char>());
+
+                                    /**
+                                     * Specialization of above
+                                     * function for 1d. Since the
+                                     * computation is not useful in
+                                     * 1d, this function simply
+                                     * throws an exception.
+                                     */
+    static void create_boundary_right_hand_side (const Mapping<1>    &mapping,
+                                                const DoFHandler<1> &dof,
+                                                const Quadrature<0> &q,
+                                                const Function<1>   &rhs,
+                                                Vector<double>      &rhs_vector,
+                                                const std::set<unsigned char> &boundary_indicators = std::set<unsigned char>());
+    
+                                    /**
+                                     * Calls the
+                                     * @p{create_boundary_right_hand_side}
+                                     * function, see above, with
+                                     * @p{mapping=MappingQ1<dim>()}.
+                                     */
+    template <int dim>
+    static void create_boundary_right_hand_side (const DoFHandler<dim>   &dof,
+                                                const Quadrature<dim-1> &q,
+                                                const Function<dim>     &rhs,
+                                                Vector<double>          &rhs_vector,
+                                                const std::set<unsigned char> &boundary_indicators = std::set<unsigned char>());
+    
                                     /**
                                      * Prepare Dirichlet boundary
                                      * conditions.  Make up the list

diff --git a/deal.II/deal.II/source/numerics/vectors.cc b/deal.II/deal.II/source/numerics/vectors.cc
index 6893354929e675f243df420497a4a4706a2afaed..d13e05133785284a1a6c6374f552327222de480c 100644 (file)
--- a/deal.II/deal.II/source/numerics/vectors.cc
+++ b/deal.II/deal.II/source/numerics/vectors.cc
@@ -545,6 +545,139 @@ void VectorTools::create_right_hand_side (const DoFHandler<dim>    &dof_handler,
 }
 
 
+
+#if deal_II_dimension != 1
+
+template <int dim>
+void
+VectorTools::create_boundary_right_hand_side (const Mapping<dim>      &mapping,
+                                             const DoFHandler<dim>   &dof_handler,
+                                             const Quadrature<dim-1> &quadrature,
+                                             const Function<dim>     &rhs_function,
+                                             Vector<double>          &rhs_vector,
+                                             const std::set<unsigned char> &boundary_indicators)
+{
+  const FiniteElement<dim> &fe  = dof_handler.get_fe();
+  Assert (fe.n_components() == rhs_function.n_components,
+         ExcComponentMismatch());
+  Assert (rhs_vector.size() == dof_handler.n_dofs(),
+         ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
+  rhs_vector.clear();
+  
+  UpdateFlags update_flags = UpdateFlags(update_values   |
+                                        update_q_points |
+                                        update_JxW_values);
+  FEFaceValues<dim> fe_values (mapping, fe, quadrature, update_flags);
+
+  const unsigned int dofs_per_cell = fe_values.dofs_per_cell,
+                    n_q_points    = fe_values.n_quadrature_points,
+                    n_components  = fe.n_components();
+  
+  std::vector<unsigned int> dofs (dofs_per_cell);
+  Vector<double> cell_vector (dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                                endc = dof_handler.end();
+
+  if (n_components==1)
+    {
+      std::vector<double> rhs_values(n_q_points);
+      
+      for (; cell!=endc; ++cell)
+       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+         if (cell->face(face)->at_boundary () &&
+             (boundary_indicators.find (cell->face(face)->boundary_indicator())
+              !=
+              boundary_indicators.end()))
+           {
+             fe_values.reinit(cell, face);
+         
+             const FullMatrix<double>  &values    = fe_values.get_shape_values ();
+             const std::vector<double> &weights   = fe_values.get_JxW_values ();
+             rhs_function.value_list (fe_values.get_quadrature_points(), rhs_values);
+             
+             cell_vector.clear();
+             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] *
+                                   values(i,point) *
+                                   weights[point];
+       
+             cell->get_dof_indices (dofs);
+         
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               rhs_vector(dofs[i]) += cell_vector(i);
+           }
+    }
+  else
+    {
+      std::vector<Vector<double> > rhs_values(n_q_points, Vector<double>(n_components));
+      
+      for (; cell!=endc; ++cell) 
+       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+         if (cell->face(face)->at_boundary () &&
+             (boundary_indicators.find (cell->face(face)->boundary_indicator())
+              !=
+              boundary_indicators.end()))
+           {
+             fe_values.reinit(cell, face);
+         
+             const FullMatrix<double>  &values    = fe_values.get_shape_values ();
+             const std::vector<double> &weights   = fe_values.get_JxW_values ();
+             rhs_function.vector_value_list (fe_values.get_quadrature_points(), rhs_values);
+         
+             cell_vector.clear();
+             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_vector(i) += rhs_values[point](component) *
+                                     values(i,point) *
+                                     weights[point];
+                 }
+             
+             cell->get_dof_indices (dofs);
+             
+             for (unsigned int i=0; i<dofs_per_cell; ++i)
+               rhs_vector(dofs[i]) += cell_vector(i);
+           }
+    }
+};
+
+#else
+
+void
+VectorTools::create_boundary_right_hand_side (const Mapping<1>    &mapping,
+                                             const DoFHandler<1> &dof_handler,
+                                             const Quadrature<0> &quadrature,
+                                             const Function<1>   &rhs_function,
+                                             Vector<double>      &rhs_vector,
+                                             const std::set<unsigned char> &boundary_indicators)
+{
+  Assert (false, ExcImpossibleInDim(1));
+};
+
+#endif
+
+template <int dim>
+void
+VectorTools::create_boundary_right_hand_side (const DoFHandler<dim>   &dof_handler,
+                                             const Quadrature<dim-1> &quadrature,
+                                             const Function<dim>     &rhs_function,
+                                             Vector<double>          &rhs_vector,
+                                             const std::set<unsigned char> &boundary_indicators)
+{
+  Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+  static const MappingQ1<dim> mapping;
+  create_boundary_right_hand_side(mapping, dof_handler, quadrature,
+                                 rhs_function, rhs_vector,
+                                 boundary_indicators);
+}
+
+
+
 #if deal_II_dimension == 1
 
 void
@@ -1361,6 +1494,25 @@ void VectorTools::create_right_hand_side (const DoFHandler<deal_II_dimension> &,
                                          const Quadrature<deal_II_dimension> &,
                                          const Function<deal_II_dimension>   &,
                                          Vector<double>                      &);
+
+#if deal_II_dimension != 1
+template
+void
+VectorTools::create_boundary_right_hand_side (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<unsigned char> &);
+#endif
+
+template
+void
+VectorTools::create_boundary_right_hand_side (const DoFHandler<deal_II_dimension> &,
+                                             const Quadrature<deal_II_dimension-1> &,
+                                             const Function<deal_II_dimension>   &,
+                                             Vector<double>                      &,
+                                             const std::set<unsigned char> &);
 template
 void VectorTools::interpolate_boundary_values (const DoFHandler<deal_II_dimension> &,
                                               const unsigned char,

diff --git a/deal.II/doc/news/2001/c-3-1.html b/deal.II/doc/news/2001/c-3-1.html
index 342dd22fc5825587e157020b5b7ad8080d9581f0..c4eea11daf7d50a16e869c339ecbca6a1d7886a8 100644 (file)
--- a/deal.II/doc/news/2001/c-3-1.html
+++ b/deal.II/doc/news/2001/c-3-1.html
@@ -596,6 +596,13 @@ documentation, etc</a>.
 <h3>deal.II</h3>
 
 <ol>
+  <li> <p> 
+       New: <code class="member">VectorTools::create_boundary_right_hand_side</code>
+       integrates boundary forces for inhomogeneous Neumann boundary values.
+       <br>
+       (WB 2001/09/13)
+       </p>
+
   <li> <p> 
        New: <code class="member">DoFTools::make_flux_sparsity_pattern</code>
        now exists also for 1d.