From 55471e19722e5c73288a146484bb2193239a5e33 Mon Sep 17 00:00:00 2001
From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Sun, 5 Nov 2017 20:23:36 +0100
Subject: [PATCH] Use PolarManifold in SphericalManifold for 2D

---
 include/deal.II/grid/manifold_lib.h | 7 +++++++
 source/grid/manifold_lib.cc         | 9 ++++++---
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
index bb808be3de..6a0294afe5 100644
--- a/include/deal.II/grid/manifold_lib.h
+++ b/include/deal.II/grid/manifold_lib.h
@@ -248,6 +248,13 @@ public:
    * The center of the spherical coordinate system.
    */
   const Point<spacedim> center;
+
+private:
+
+  /**
+   * A manifold description to be used for get_new_point in 2D.
+   **/
+  const PolarManifold<spacedim> polar_manifold;
 };
 
 
diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
index 5256a5ead9..94d25fc02e 100644
--- a/source/grid/manifold_lib.cc
+++ b/source/grid/manifold_lib.cc
@@ -245,7 +245,8 @@ PolarManifold<dim,spacedim>::push_forward_gradient(const Point<spacedim> &spheri
 
 template <int dim, int spacedim>
 SphericalManifold<dim,spacedim>::SphericalManifold(const Point<spacedim> center):
-  center(center)
+  center(center),
+  polar_manifold(center)
 {}
 
 
@@ -403,8 +404,10 @@ get_new_point (const ArrayView<const Point<spacedim>> &vertices,
     rho /= total_weights;
   }
 
-  if (spacedim<2)
-    return center + rho*candidate;
+  // If not in 3D, just use the implementation from PolarManifold
+  // after we verified that the candidate is not the center.
+  if (spacedim<3)
+    return polar_manifold.get_new_point(vertices, weights);
 
   // Step 2:
   // Do Newton-style iterations to improve the estimate.
-- 
2.39.5