add geometry utilities functions

diff --git a/doc/news/changes.h b/doc/news/changes.h
index c62ea5a10de492e82b4912af45f7057d752e3c30..3c1111ca72eb76ade55374d4ff4b2593aaaff5dc 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -120,6 +120,13 @@ inconvenience this causes.
 <h3>General</h3>
 
 <ol>
+ <li> New: Add functions to transform Cartesian coordinates to spherical and back:
+ GeometricUtilities::Coordinates::to_spherical and
+ GeometricUtilities::Coordinates::from_spherical.
+ <br>
+ (Denis Davydov, 2016/05/10)
+ </li>
+
  <li> 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 (file)
index 0000000..fc8f39a
--- /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 <deal.II/base/config.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/std_cxx11/array.h>
+
+
+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 <int dim>
+    std_cxx11::array<double,dim>
+    to_spherical(const Point<dim> &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 <std::size_t dim>
+    Point<dim>
+    from_spherical(const std_cxx11::array<double,dim> &scoord);
+
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt
index f13a47d856c36ef2fb5ddefbb81d43b8af4f93ac..587b1dada98c496c88ae2f409c281076c9ed4be1 100644 (file)
--- 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 (file)
index 0000000..427cdf2
--- /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 <deal.II/base/geometric_utilities.h>
+#include <deal.II/base/exceptions.h>
+
+
+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 <int dim>
+    std_cxx11::array<double,dim>
+    to_spherical(const Point<dim> &position)
+    {
+      std_cxx11::array<double,dim> 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<double>::min())
+            scoord[2] = std::acos(position(2)/scoord[0]);
+          else
+            scoord[2] = 0.0;
+        }
+      return scoord;
+    }
+
+    template <std::size_t dim>
+    Point<dim>
+    from_spherical(const std_cxx11::array<double,dim> &scoord)
+    {
+      Point<dim> 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<double,2> &);
+    template Point<3> from_spherical<3>(const std_cxx11::array<double,3> &);
+    template std_cxx11::array<double,2> to_spherical<2>(const Point<2> &);
+    template std_cxx11::array<double,3> 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 (file)
index 0000000..abdbd90
--- /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 <deal.II/base/logstream.h>
+#include <deal.II/base/geometric_utilities.h>
+
+#include <fstream>
+#include <cstdlib>
+
+using namespace dealii;
+using namespace dealii::GeometricUtilities::Coordinates;
+
+template <unsigned int dim,typename T>
+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 <int dim>
+void test ()
+{
+  Assert(dim>1, ExcNotImplemented());
+
+  const Point<dim> origin;
+
+  std_cxx1x::array<double,dim> sorigin;
+  for (unsigned int d= 0; d < dim; d++)
+    sorigin[d] = 0.;
+
+  Point<dim> one;
+  for (unsigned int d=0; d< dim; d++)
+    one[d] = 1.;
+
+  std_cxx1x::array<double,dim> sone;
+  sone[0] = std::sqrt(1.*dim);
+  sone[1] = numbers::PI/4;
+  if (dim==3)
+    sone[2] = std::acos(1/std::sqrt(3.));
+
+  print<dim>(to_spherical(origin),sorigin);
+  print<dim>(origin, from_spherical(sorigin));
+
+  print<dim>(to_spherical(one),sone);
+  print<dim>(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<double,3> sx;
+  sx[0] = 1.;
+  sx[1] = 0;
+  sx[2] = numbers::PI/2;
+  std_cxx1x::array<double,3> sy;
+  sy[0] = 1;
+  sy[1] = numbers::PI/2;
+  sy[2] = numbers::PI/2;
+  std_cxx1x::array<double,3> 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<double,3> 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 (file)
index 0000000..6176048
--- /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