From 4620af6de0f037bf292f4d85262ce6841f359ac6 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sat, 5 Mar 2016 19:36:20 +0100 Subject: [PATCH] Added tangent vector to SphericalManifold. --- include/deal.II/grid/manifold_lib.h | 18 ++++++ source/grid/manifold_lib.cc | 41 +++++++++++++ tests/manifold/spherical_manifold_06.cc | 68 +++++++++++++++++++++ tests/manifold/spherical_manifold_06.output | 16 +++++ 4 files changed, 143 insertions(+) create mode 100644 tests/manifold/spherical_manifold_06.cc create mode 100644 tests/manifold/spherical_manifold_06.output diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h index 78694794e8..21de7b0555 100644 --- a/include/deal.II/grid/manifold_lib.h +++ b/include/deal.II/grid/manifold_lib.h @@ -84,6 +84,24 @@ public: virtual Point push_forward(const Point &chart_point) const; + + /** + * Given a point in the spacedim dimensional Euclidean space, this + * method returns the derivatives of the function $F$ that maps from + * the polar coordinate system to the Euclidean coordinate + * system. In other words, it is a matrix of size + * $\text{spacedim}\times\text{spacedim}$. + * + * This function is used in the computations required by the + * get_tangent_vector() function. + * + * Refer to the general documentation of this class for more information. + */ + virtual + DerivativeForm<1,spacedim,spacedim> + push_forward_gradient(const Point &chart_point) const; + + /** * Let the new point be the average sum of surrounding vertices. * diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index cfe05206f9..bb5b91c861 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -131,6 +131,47 @@ SphericalManifold::pull_back(const Point &space_point) c } +template +DerivativeForm<1,spacedim,spacedim> +SphericalManifold::push_forward_gradient(const Point &spherical_point) const +{ + Assert(spherical_point[0] >=0.0, + ExcMessage("Negative radius for given point.")); + const double rho = spherical_point[0]; + const double theta = spherical_point[1]; + + DerivativeForm<1,spacedim,spacedim> DX; + if (rho > 1e-10) + switch (spacedim) + { + case 2: + DX[0][0] = cos(theta); + DX[0][1] = -rho*sin(theta); + DX[1][0] = sin(theta); + DX[1][1] = rho*cos(theta); + break; + case 3: + { + const double &phi= spherical_point[2]; + DX[0][0] = sin(theta)*cos(phi); + DX[0][1] = rho*cos(theta)*cos(phi); + DX[0][2] = -rho*sin(theta)*sin(phi); + + DX[1][0] = sin(theta)*sin(phi); + DX[1][1] = rho*cos(theta)*sin(phi); + DX[1][2] = rho*sin(theta)*cos(phi); + + DX[2][0] = cos(theta); + DX[2][1] = -rho*sin(theta); + DX[2][2] = 0; + } + break; + default: + Assert(false, ExcInternalError()); + } + return DX; +} + // ============================================================ // CylindricalManifold // ============================================================ diff --git a/tests/manifold/spherical_manifold_06.cc b/tests/manifold/spherical_manifold_06.cc new file mode 100644 index 0000000000..c2800aec13 --- /dev/null +++ b/tests/manifold/spherical_manifold_06.cc @@ -0,0 +1,68 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Check get_tangent_vector for spherical manifold, on simple points. + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include + +#include + +template +void test() +{ + deallog << "dim=" << dim << ", spacedim=" << spacedim << std::endl; + + Point center; + static const SphericalManifold manifold(center); + + // Go from 0,1 to 1,0 + Point p0,p1; + p0[1] = 1.0; + p1[0] = 1.0; + + Tensor<1,spacedim> T = manifold.get_tangent_vector(p0, p1); + + deallog << "P0 : " << p0 << std::endl; + deallog << "P1 : " << p1 << std::endl; + deallog << "T(P0-P1): " << T << std::endl; + deallog << "Error : " << T.norm() - numbers::PI/2 << std::endl; +} + + +int +main() +{ + std::ofstream logfile ("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test<2,2>(); + test<2,3>(); + + test<3,3>(); + + return 0; +} + + + diff --git a/tests/manifold/spherical_manifold_06.output b/tests/manifold/spherical_manifold_06.output new file mode 100644 index 0000000000..d906a07923 --- /dev/null +++ b/tests/manifold/spherical_manifold_06.output @@ -0,0 +1,16 @@ + +DEAL::dim=2, spacedim=2 +DEAL::P0 : 0.00000 1.00000 +DEAL::P1 : 1.00000 0.00000 +DEAL::T(P0-P1): 1.57080 -9.61835e-17 +DEAL::Error : 0 +DEAL::dim=2, spacedim=3 +DEAL::P0 : 0.00000 1.00000 0.00000 +DEAL::P1 : 1.00000 0.00000 0.00000 +DEAL::T(P0-P1): 1.57080 -9.61835e-17 0.00000 +DEAL::Error : 0 +DEAL::dim=3, spacedim=3 +DEAL::P0 : 0.00000 1.00000 0.00000 +DEAL::P1 : 1.00000 0.00000 0.00000 +DEAL::T(P0-P1): 1.57080 -9.61835e-17 0.00000 +DEAL::Error : 0 -- 2.39.5