From 0bcf2d26a2cf6155577cc4d07629415c6b7f11f8 Mon Sep 17 00:00:00 2001
From: buerg <buerg@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 6 May 2012 10:03:52 +0000
Subject: [PATCH] Add create_point_source_vector for vector-valued FE.

git-svn-id: https://svn.dealii.org/trunk@25497 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/include/deal.II/numerics/vectors.h    | 57 +++++++++++
 .../deal.II/numerics/vectors.templates.h      | 99 +++++++++++++++++++
 deal.II/source/numerics/vectors.inst.in       | 27 +++++
 3 files changed, 183 insertions(+)

diff --git a/deal.II/include/deal.II/numerics/vectors.h b/deal.II/include/deal.II/numerics/vectors.h
index f2823f9300..f1fc9bf688 100644
--- a/deal.II/include/deal.II/numerics/vectors.h
+++ b/deal.II/include/deal.II/numerics/vectors.h
@@ -1714,6 +1714,63 @@ namespace VectorTools
                                   const Point<spacedim>      &p,
                                   Vector<double>        &rhs_vector);
 
+                                   /**
+                                    * Create a right hand side
+                                    * vector for a point source at point @p p
+                                    * for vector-valued finite elements.
+                                    * Prior content of the
+                                    * given @p rhs_vector vector is
+                                    * deleted.
+                                    *
+                                    * See the general documentation of this
+                                    * class for further information.
+                                    */
+  template <int dim, int spacedim>
+  void create_point_source_vector(const Mapping<dim,spacedim>    &mapping,
+                                  const DoFHandler<dim,spacedim> &dof,
+                                  const Point<spacedim>          &p,
+                                  const Point<dim>               &orientation,
+                                  Vector<double>                 &rhs_vector);
+
+                                   /**
+                                    * Calls the create_point_source_vector()
+                                    * function for vector-valued finite elements,
+                                    * see above, with
+                                    * <tt>mapping=MappingQ1@<dim@>()</tt>.
+                                    */
+  template <int dim, int spacedim>
+  void create_point_source_vector(const DoFHandler<dim,spacedim> &dof,
+                                  const Point<spacedim>          &p,
+                                  const Point<dim>               &orientation,
+                                  Vector<double>                 &rhs_vector);
+
+                                   /**
+                                    * Like the previous set of functions,
+                                    * but for hp objects.
+                                    */
+  template <int dim, int spacedim>
+  void create_point_source_vector(const hp::MappingCollection<dim,spacedim> &mapping,
+                                  const hp::DoFHandler<dim,spacedim>        &dof,
+                                  const Point<spacedim>                     &p,
+                                  const Point<dim>                          &orientation,
+                                  Vector<double>                            &rhs_vector);
+
+                                   /**
+                                    * Like the previous set of functions,
+                                    * but for hp objects. The function uses
+                                    * the default Q1 mapping object. Note
+                                    * that if your hp::DoFHandler uses any
+                                    * active fe index other than zero, then
+                                    * you need to call the function above
+                                    * that provides a mapping object for
+                                    * each active fe index.
+                                    */
+  template <int dim, int spacedim>
+  void create_point_source_vector(const hp::DoFHandler<dim,spacedim> &dof,
+                                  const Point<spacedim>              &p,
+                                  const Point<dim>                   &orientation,
+                                  Vector<double>                     &rhs_vector);
+
                                    /**
                                     * Create a right hand side
                                     * vector from boundary
diff --git a/deal.II/include/deal.II/numerics/vectors.templates.h b/deal.II/include/deal.II/numerics/vectors.templates.h
index 2fb9e657db..72b02f2b7b 100644
--- a/deal.II/include/deal.II/numerics/vectors.templates.h
+++ b/deal.II/include/deal.II/numerics/vectors.templates.h
@@ -1008,6 +1008,105 @@ namespace VectorTools
 
 
 
+
+  template <int dim, int spacedim>
+  void create_point_source_vector (const Mapping<dim, spacedim>       &mapping,
+                                   const DoFHandler<dim,spacedim>     &dof_handler,
+                                   const Point<spacedim>              &p,
+                                   const Point<dim>                   &orientation,
+                                   Vector<double>                     &rhs_vector)
+  {
+    Assert (rhs_vector.size() == dof_handler.n_dofs(),
+            ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
+    Assert (dof_handler.get_fe().n_components() == dim,
+            ExcMessage ("This function only works for vector-valued finite elements."));
+
+    rhs_vector = 0;
+
+    std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
+      cell_point =
+      GridTools::find_active_cell_around_point (mapping, dof_handler, p);
+
+    Quadrature<dim> q(GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
+
+    const FEValuesExtractors::Vector vec (0);
+    FEValues<dim,spacedim> fe_values(mapping, dof_handler.get_fe(),
+                                     q, UpdateFlags(update_values));
+    fe_values.reinit(cell_point.first);
+
+    const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    cell_point.first->get_dof_indices (local_dof_indices);
+
+    for(unsigned int i=0; i<dofs_per_cell; i++)
+      rhs_vector(local_dof_indices[i]) =  orientation * fe_values[vec].value(i,0);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_point_source_vector (const DoFHandler<dim,spacedim>    &dof_handler,
+                                   const Point<spacedim>             &p,
+                                   const Point<dim>                  &orientation,
+                                   Vector<double>                    &rhs_vector)
+  {
+    Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+    create_point_source_vector(StaticMappingQ1<dim,spacedim>::mapping, dof_handler,
+                               p, orientation, rhs_vector);
+  }
+
+
+  template <int dim, int spacedim>
+  void create_point_source_vector (const hp::MappingCollection<dim,spacedim> &mapping,
+                                   const hp::DoFHandler<dim,spacedim>        &dof_handler,
+                                   const Point<spacedim>                     &p,
+                                   const Point<dim>                          &orientation,
+                                   Vector<double>                            &rhs_vector)
+  {
+    Assert (rhs_vector.size() == dof_handler.n_dofs(),
+            ExcDimensionMismatch(rhs_vector.size(), dof_handler.n_dofs()));
+    Assert (dof_handler.get_fe().n_components() == dim,
+            ExcMessage ("This function only works for vector-valued finite elements."));
+
+    rhs_vector = 0;
+
+    std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
+      cell_point =
+      GridTools::find_active_cell_around_point (mapping, dof_handler, p);
+
+    Quadrature<dim> q(GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
+
+    const FEValuesExtractors::Vector vec (0);
+    FEValues<dim> fe_values(mapping[cell_point.first->active_fe_index()],
+                            cell_point.first->get_fe(), q, UpdateFlags(update_values));
+    fe_values.reinit(cell_point.first);
+
+    const unsigned int dofs_per_cell = cell_point.first->get_fe().dofs_per_cell;
+
+    std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+    cell_point.first->get_dof_indices (local_dof_indices);
+
+    for(unsigned int i=0; i<dofs_per_cell; i++)
+      rhs_vector(local_dof_indices[i]) =  orientation * fe_values[vec].value(i,0);
+  }
+
+
+
+  template <int dim, int spacedim>
+  void create_point_source_vector (const hp::DoFHandler<dim,spacedim>   &dof_handler,
+                                   const Point<spacedim>                &p,
+                                   const Point<dim>                     &orientation,
+                                   Vector<double>                       &rhs_vector)
+  {
+    Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
+    create_point_source_vector(hp::StaticMappingQ1<dim>::mapping_collection,
+                               dof_handler,
+                               p, orientation, rhs_vector);
+  }
+
+
+
 // separate implementation for 1D because otherwise we get linker errors since
 // FEFaceValues<1> is not compiled
   template <>
diff --git a/deal.II/source/numerics/vectors.inst.in b/deal.II/source/numerics/vectors.inst.in
index 1062bbc524..c7fd3f27a1 100644
--- a/deal.II/source/numerics/vectors.inst.in
+++ b/deal.II/source/numerics/vectors.inst.in
@@ -462,6 +462,33 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension :  SPACE_DIMENSIONS
         (const hp::DoFHandler<deal_II_dimension> &,
          const Point<deal_II_dimension>      &,
          Vector<double>                      &);
+      template
+        void create_point_source_vector<deal_II_dimension>
+        (const Mapping<deal_II_dimension>    &,
+         const DoFHandler<deal_II_dimension> &,
+         const Point<deal_II_dimension>      &,
+         const Point<deal_II_dimension>      &,
+         Vector<double>                      &);
+      template
+        void create_point_source_vector<deal_II_dimension>
+        (const DoFHandler<deal_II_dimension> &,
+         const Point<deal_II_dimension>      &,
+         const Point<deal_II_dimension>      &,
+         Vector<double>                      &);
+
+      template
+        void create_point_source_vector<deal_II_dimension>
+        (const hp::MappingCollection<deal_II_dimension>    &,
+         const hp::DoFHandler<deal_II_dimension> &,
+         const Point<deal_II_dimension>      &,
+         const Point<deal_II_dimension>      &,
+         Vector<double>                      &);
+      template
+        void create_point_source_vector<deal_II_dimension>
+        (const hp::DoFHandler<deal_II_dimension> &,
+         const Point<deal_II_dimension>      &,
+         const Point<deal_II_dimension>      &,
+         Vector<double>                      &);
 
 #if deal_II_dimension > 1
       template
-- 
2.39.5