From: Denis Davydov Date: Tue, 10 May 2016 19:23:24 +0000 (+0200) Subject: add geometry utilities functions X-Git-Tag: v8.5.0-rc1~1040^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=80d0a7181efb11aa8e7f0a0467eac6842dbadf35;p=dealii.git add geometry utilities functions --- diff --git a/doc/news/changes.h b/doc/news/changes.h index c62ea5a10d..3c1111ca72 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -120,6 +120,13 @@ inconvenience this causes.

General

    +
  1. New: Add functions to transform Cartesian coordinates to spherical and back: + GeometricUtilities::Coordinates::to_spherical and + GeometricUtilities::Coordinates::from_spherical. +
    + (Denis Davydov, 2016/05/10) +
    2. +
  2. Improved: The method Triangulation::create_triangulation will now throw an exception if any cells have negative measure. This check is not run if the triangulation keeps track of distorted cells or if the codimension is not zero. diff --git a/include/deal.II/base/geometric_utilities.h b/include/deal.II/base/geometric_utilities.h new file mode 100644 index 0000000000..fc8f39a89a --- /dev/null +++ b/include/deal.II/base/geometric_utilities.h @@ -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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__geometric_utilities_h +#define dealii__geometric_utilities_h + +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +/** + * A namespace for geometric utility functions that are not particularly specific + * to finite element computing or numerical programs, but nevertheless are needed + * in various contexts when writing applications. + * + * @ingroup utilities + * @author Denis Davydov, 2016 + */ +namespace GeometricUtilities +{ + + /** + * A namespace for coordinate transformations. + */ + namespace Coordinates + { + + /** + * Returns spherical coordinates of a Cartesian point @p point. + * The returned array is filled with radius, azimuth angle \f$ \in [0,2 \pi) \f$ + * and polar/inclination angle \f$ \in [0,\pi] \f$ (ommited in 2D). + */ + template + std_cxx11::array + to_spherical(const Point &point); + + /** + * Return the Cartesian coordinates of a spherical point defined by @p scoord + * which is filled with radius \f$ \in [0,\infty) \f$, azimuth angle + * \f$ \in [0,2 \pi) \f$ and polar/inclination angle \f$ \in [0,\pi] \f$ + * (ommited in 2D). + */ + template + Point + from_spherical(const std_cxx11::array &scoord); + + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt index f13a47d856..587b1dada9 100644 --- a/source/base/CMakeLists.txt +++ b/source/base/CMakeLists.txt @@ -32,6 +32,7 @@ SET(_src function_parser.cc function_time.cc geometry_info.cc + geometric_utilities.cc index_set.cc job_identifier.cc logstream.cc diff --git a/source/base/geometric_utilities.cc b/source/base/geometric_utilities.cc new file mode 100644 index 0000000000..427cdf2ca5 --- /dev/null +++ b/source/base/geometric_utilities.cc @@ -0,0 +1,115 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace GeometricUtilities +{ + + namespace Coordinates + { + DeclException1(NegativeRadius, + double, + << "The radius <" << arg1 << "> can not be negative."); + + DeclException1(SphericalAzimuth, + double, + << "The azimuth angle <" << arg1 << "> is not in [0,2Pi)."); + + DeclException1(SphericalPolar, + double, + << "The polar angle <" << arg1 << "> is not in [0,Pi]."); + + + template + std_cxx11::array + to_spherical(const Point &position) + { + std_cxx11::array scoord; + + // radius + scoord[0] = position.norm(); + // azimuth angle + scoord[1] = std::atan2(position(1),position(0)); + // correct to [0,2*pi) + if (scoord[1] < 0.0) + scoord[1] += 2.0*numbers::PI; + + // polar angle + if (dim==3) + { + // acos returns the angle in the range [0,\pi] + if (scoord[0] > std::numeric_limits::min()) + scoord[2] = std::acos(position(2)/scoord[0]); + else + scoord[2] = 0.0; + } + return scoord; + } + + template + Point + from_spherical(const std_cxx11::array &scoord) + { + Point ccoord; + + Assert (scoord[0] >= 0., + NegativeRadius(scoord[0])); + + Assert (scoord[1] >= 0. && scoord[1] < 2.*numbers::PI, + SphericalAzimuth(scoord[1])); + + switch (dim) + { + case 2: + { + ccoord[0] = scoord[0] * std::cos(scoord[1]); + ccoord[1] = scoord[0] * std::sin(scoord[1]); + break; + } + case 3: + { + Assert (scoord[2] >= 0. && scoord[2] <= numbers::PI, + SphericalPolar(scoord[2])); + + ccoord[0] = scoord[0] * std::sin(scoord[2]) * std::cos(scoord[1]); + ccoord[1] = scoord[0] * std::sin(scoord[2]) * std::sin(scoord[1]); + ccoord[2] = scoord[0] * std::cos(scoord[2]); + break; + } + default: + Assert (false, ExcNotImplemented()); + break; + } + + return ccoord; + } + + + template Point<2> from_spherical<2>(const std_cxx11::array &); + template Point<3> from_spherical<3>(const std_cxx11::array &); + template std_cxx11::array to_spherical<2>(const Point<2> &); + template std_cxx11::array to_spherical<3>(const Point<3> &); + + } +} + +DEAL_II_NAMESPACE_CLOSE diff --git a/tests/base/geometric_utilities_01.cc b/tests/base/geometric_utilities_01.cc new file mode 100644 index 0000000000..abdbd90375 --- /dev/null +++ b/tests/base/geometric_utilities_01.cc @@ -0,0 +1,122 @@ +// --------------------------------------------------------------------- +// +// 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 GeometricUtilities::Coordinates::to_spherical and +// GeometricUtilities::Coordinates::from_spherical. + +#include "../tests.h" +#include +#include + +#include +#include + +using namespace dealii; +using namespace dealii::GeometricUtilities::Coordinates; + +template +void print(T point1, T point2) +{ + deallog << std::endl << "Point 1: "; + for (unsigned int i = 0; i < dim; ++i) + deallog << point1[i] << " "; + + deallog << std::endl << "Point 2: "; + for (unsigned int i = 0; i < dim; ++i) + deallog << point2[i] << " "; + + deallog << std::endl; +} + +template +void test () +{ + Assert(dim>1, ExcNotImplemented()); + + const Point origin; + + std_cxx1x::array sorigin; + for (unsigned int d= 0; d < dim; d++) + sorigin[d] = 0.; + + Point one; + for (unsigned int d=0; d< dim; d++) + one[d] = 1.; + + std_cxx1x::array sone; + sone[0] = std::sqrt(1.*dim); + sone[1] = numbers::PI/4; + if (dim==3) + sone[2] = std::acos(1/std::sqrt(3.)); + + print(to_spherical(origin),sorigin); + print(origin, from_spherical(sorigin)); + + print(to_spherical(one),sone); + print(one, from_spherical(sone)); +} + +void test3d() +{ + const dealii::Point<3> x(1,0,0); + const dealii::Point<3> y(0,1,0); + const dealii::Point<3> z(0,0,1); + + std_cxx1x::array sx; + sx[0] = 1.; + sx[1] = 0; + sx[2] = numbers::PI/2; + std_cxx1x::array sy; + sy[0] = 1; + sy[1] = numbers::PI/2; + sy[2] = numbers::PI/2; + std_cxx1x::array sz; + sz[0] = 1.; + sz[1] = 0.; + sz[2] = 0.; + + print<3>(x, from_spherical(sx)); + print<3>(y, from_spherical(sy)); + print<3>(z, from_spherical(sz)); + + print<3>(to_spherical(x),sx); + print<3>(to_spherical(y),sy); + print<3>(to_spherical(z),sz); + + const Point<3> dateline(0,-1,0); + std_cxx1x::array sdateline; + sdateline[0] = 1.; + sdateline[1] = 3*numbers::PI/2; + sdateline[2] = numbers::PI/2; + + print<3>(dateline, from_spherical(sdateline)); + print<3>(to_spherical(dateline),sdateline); +} + + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test<2> (); + test<3> (); + test3d (); + + return 0; +} diff --git a/tests/base/geometric_utilities_01.output b/tests/base/geometric_utilities_01.output new file mode 100644 index 0000000000..61760481af --- /dev/null +++ b/tests/base/geometric_utilities_01.output @@ -0,0 +1,49 @@ + +DEAL:: +DEAL::Point 1: 0 0 +DEAL::Point 2: 0 0 +DEAL:: +DEAL::Point 1: 0 0 +DEAL::Point 2: 0 0 +DEAL:: +DEAL::Point 1: 1.41421 0.785398 +DEAL::Point 2: 1.41421 0.785398 +DEAL:: +DEAL::Point 1: 1.00000 1.00000 +DEAL::Point 2: 1.00000 1.00000 +DEAL:: +DEAL::Point 1: 0 0 0 +DEAL::Point 2: 0 0 0 +DEAL:: +DEAL::Point 1: 0 0 0 +DEAL::Point 2: 0 0 0 +DEAL:: +DEAL::Point 1: 1.73205 0.785398 0.955317 +DEAL::Point 2: 1.73205 0.785398 0.955317 +DEAL:: +DEAL::Point 1: 1.00000 1.00000 1.00000 +DEAL::Point 2: 1.00000 1.00000 1.00000 +DEAL:: +DEAL::Point 1: 1.00000 0 0 +DEAL::Point 2: 1.00000 0 0 +DEAL:: +DEAL::Point 1: 0 1.00000 0 +DEAL::Point 2: 0 1.00000 0 +DEAL:: +DEAL::Point 1: 0 0 1.00000 +DEAL::Point 2: 0 0 1.00000 +DEAL:: +DEAL::Point 1: 1.00000 0 1.57080 +DEAL::Point 2: 1.00000 0 1.57080 +DEAL:: +DEAL::Point 1: 1.00000 1.57080 1.57080 +DEAL::Point 2: 1.00000 1.57080 1.57080 +DEAL:: +DEAL::Point 1: 1.00000 0 0 +DEAL::Point 2: 1.00000 0 0 +DEAL:: +DEAL::Point 1: 0 -1.00000 0 +DEAL::Point 2: 0 -1.00000 0 +DEAL:: +DEAL::Point 1: 1.00000 4.71239 1.57080 +DEAL::Point 2: 1.00000 4.71239 1.57080