From: Wolfgang Bangerth Date: Sun, 28 Feb 2016 18:59:35 +0000 (-0600) Subject: Provide a default implementation of Manifold::get_tangent_vector(). X-Git-Tag: v8.5.0-rc1~1272^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d3973743b91c132f1f5962e72ad84aa0aaada02e;p=dealii.git Provide a default implementation of Manifold::get_tangent_vector(). --- diff --git a/include/deal.II/grid/manifold.h b/include/deal.II/grid/manifold.h index c8390d49ac..e664042ed5 100644 --- a/include/deal.II/grid/manifold.h +++ b/include/deal.II/grid/manifold.h @@ -199,9 +199,11 @@ namespace Manifolds * \frac{\mathbf s(w)-\mathbf s(0)}{w}$ where all we need to do * is compute the new point $\mathbf s(w)$ with weights $w$ and * $1-w$ along the geodesic connecting $\mathbf x_1$ and $\mathbf x_2$. + * The default implementation of the function does this, by evaluating + * the quotient for a small but finite weight $w$. * In practice, however, it is almost always possible to explicitly * compute the direction vector, i.e., without the need to numerically - * approximate the limit process. + * approximate the limit process, and derived classes should do so. * * * @ingroup manifold @@ -365,12 +367,11 @@ public: * $\mathbf s(t)$ must move "faster" if the two points it connects between * arguments $t=0$ and $t=1$ are farther apart. * - * This function is used, among other cases, in computing normal vectors to - * faces or, more generally, surfaces such as the boundary. Since not all - * programs need this functionality, this function has a default - * implementation that just throws an exception. Consequently, derived - * classes only have to implement this function if the program that uses - * them does in fact call it directly or indirectly. + * The default implementation of this function approximates + * $\mathbf s'(0) \approx \frac{$\mathbf s(\epsilon)-\mathbf x_1}{\epsilon}$ + * for a small value of $\epsilon$, and the evaluation of $\mathbf s(\epsilon)$ + * is done by calling get_new_point(). If possible, derived classes should + * override this function by an implement of the exact derivative. * * @param x1 The first point that describes the geodesic, and the one * at which the "direction" is to be evaluated. diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 0bf9d3dcbc..6dbeaec9fc 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -214,11 +214,21 @@ get_new_point_on_hex (const Triangulation<3, 3>::hex_iterator &hex) const template Tensor<1,spacedim> -Manifold::get_tangent_vector(const Point &, - const Point &) const +Manifold::get_tangent_vector(const Point &x1, + const Point &x2) const { - Assert (false, ExcPureFunctionCalled()); - return Tensor<1,spacedim>(); + const double epsilon = 1e-8; + + std::vector > q; + q.push_back(x1); + q.push_back(x2); + + std::vector w; + w.push_back(epsilon); + w.push_back(1.0-epsilon); + + const Tensor<1,spacedim> neighbor_point = get_new_point (Quadrature(q, w)); + return (neighbor_point-x1)/epsilon; } /* -------------------------- FlatManifold --------------------- */