From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Sun, 5 Nov 2017 19:23:36 +0000 (+0100)
Subject: Use PolarManifold in SphericalManifold for 2D
X-Git-Tag: v9.0.0-rc1~818^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5404%2Fhead;p=dealii.git

Use PolarManifold in SphericalManifold for 2D
---

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.