From: Daniel Arndt <d.arndt@math.uni-goettingen.de>
Date: Sun, 3 May 2015 22:42:30 +0000 (+0200)
Subject: Add an entry to changes.h and reduce the number of actual implementations for point_v... 
X-Git-Tag: v8.3.0-rc1~200^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F896%2Fhead;p=dealii.git

Add an entry to changes.h and reduce the number of actual implementations for point_value and point_gradient
---

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 934a491e02..822adabc74 100644
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -472,6 +472,12 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: The function VectorTools::point_gradient has been added to compute
+  the gradient of a given FE function.
+  <br>
+  (Daniel Arndt, 2015/05/03)
+  </li>
+
   <li> New: dealii:Vector, dealii::BlockVector,
   TrilinosWrappers::MPI::Vector, TrilinosWrappers::MPI::BlockVector,
   PETScWrappers::MPI::Vector and PETScWrappers::MPI::BlockVector now have
diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h
index 7ce3c469e6..2da0dec20c 100644
--- a/include/deal.II/numerics/vector_tools.h
+++ b/include/deal.II/numerics/vector_tools.h
@@ -2112,7 +2112,7 @@ namespace VectorTools
     * (vector) gradient of this function through the last argument.
     *
     * This is a wrapper function using a Q1-mapping for cell boundaries to call
-    * the other point_difference() function.
+    * the other point_gradient() function.
     *
     * @note If the cell in which the point is found is not locally owned, an
     * exception of type VectorTools::ExcPointNotAvailableHere is thrown.
diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 08ecab5eff..3ec9ff28b3 100644
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -6729,28 +6729,10 @@ namespace VectorTools
     Assert(fe.n_components() == 1,
            ExcMessage ("Finite element is not scalar as is necessary for this function"));
 
-    // first find the cell in which this point
-    // is, initialize a quadrature rule with
-    // it, and then a FEValues object
-    const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
-    cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
-
-    AssertThrow(cell_point.first->is_locally_owned(),
-                ExcPointNotAvailableHere());
-    Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
-           ExcInternalError());
+    Vector<double> value(1);
+    point_value(mapping, dof, fe_function, point, value);
 
-    const Quadrature<dim>
-    quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
-    FEValues<dim> fe_values(mapping, fe, quadrature, update_values);
-    fe_values.reinit(cell_point.first);
-
-    // then use this to get at the values of
-    // the given fe_function at this point
-    std::vector<double> u_value(1);
-    fe_values.get_function_values(fe_function, u_value);
-
-    return u_value[0];
+    return value(0);
   }
 
 
@@ -6766,29 +6748,10 @@ namespace VectorTools
     Assert(fe.n_components() == 1,
            ExcMessage ("Finite element is not scalar as is necessary for this function"));
 
-    // first find the cell in which this point
-    // is, initialize a quadrature rule with
-    // it, and then a FEValues object
-    const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
-    cell_point =  GridTools::find_active_cell_around_point (mapping, dof, point);
-
-    AssertThrow(cell_point.first->is_locally_owned(),
-                ExcPointNotAvailableHere());
-    Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
-           ExcInternalError());
-
-    const Quadrature<dim>
-    quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
-    hp::FEValues<dim, spacedim> hp_fe_values(mapping, fe, hp::QCollection<dim>(quadrature), update_values);
-    hp_fe_values.reinit(cell_point.first);
-    const FEValues<dim, spacedim> &fe_values = hp_fe_values.get_present_fe_values();
-
-    // then use this to get at the values of
-    // the given fe_function at this point
-    std::vector<double> u_value(1);
-    fe_values.get_function_values(fe_function, u_value);
+    Vector<double> value(1);
+    point_value(mapping, dof, fe_function, point, value);
 
-    return u_value[0];
+    return value(0);
   }
 
 
@@ -6941,28 +6904,10 @@ namespace VectorTools
     Assert(fe.n_components() == 1,
            ExcMessage ("Finite element is not scalar as is necessary for this function"));
 
-    // first find the cell in which this point
-    // is, initialize a quadrature rule with
-    // it, and then a FEValues object
-    const std::pair<typename DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
-    cell_point = GridTools::find_active_cell_around_point (mapping, dof, point);
-
-    AssertThrow(cell_point.first->is_locally_owned(),
-                ExcPointNotAvailableHere());
-    Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
-           ExcInternalError());
+    std::vector<Tensor<1, dim> > gradient(1);
+    point_gradient (mapping, dof, fe_function, point, gradient);
 
-    const Quadrature<dim>
-    quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
-    FEValues<dim> fe_values(mapping, fe, quadrature, update_gradients);
-    fe_values.reinit(cell_point.first);
-
-    // then use this to get the gradients of
-    // the given fe_function at this point
-    std::vector<Tensor<1, dim> > u_gradient(1);
-    fe_values.get_function_gradients(fe_function, u_gradient);
-
-    return u_gradient[0];
+    return gradient[0];
   }
 
 
@@ -6978,29 +6923,10 @@ namespace VectorTools
     Assert(fe.n_components() == 1,
            ExcMessage ("Finite element is not scalar as is necessary for this function"));
 
-    // first find the cell in which this point
-    // is, initialize a quadrature rule with
-    // it, and then a FEValues object
-    const std::pair<typename hp::DoFHandler<dim,spacedim>::active_cell_iterator, Point<spacedim> >
-    cell_point =  GridTools::find_active_cell_around_point (mapping, dof, point);
-
-    AssertThrow(cell_point.first->is_locally_owned(),
-                ExcPointNotAvailableHere());
-    Assert(GeometryInfo<dim>::distance_to_unit_cell(cell_point.second) < 1e-10,
-           ExcInternalError());
-
-    const Quadrature<dim>
-    quadrature (GeometryInfo<dim>::project_to_unit_cell(cell_point.second));
-    hp::FEValues<dim, spacedim> hp_fe_values(mapping, fe, hp::QCollection<dim>(quadrature), update_gradients);
-    hp_fe_values.reinit(cell_point.first);
-    const FEValues<dim, spacedim> &fe_values = hp_fe_values.get_present_fe_values();
-
-    // then use this to get the gradients of
-    // the given fe_function at this point
-    std::vector<Tensor<1, dim > > u_gradient(1);
-    fe_values.get_function_gradients(fe_function, u_gradient);
+    std::vector<Tensor<1, dim> > gradient(1);
+    point_gradient (mapping, dof, fe_function, point, gradient);
 
-    return u_gradient[0];
+    return gradient[0];
   }