From a74b01abbf61f7b102e0ccd65644f0240a403600 Mon Sep 17 00:00:00 2001
From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Thu, 2 Nov 2017 15:55:04 +0100
Subject: [PATCH] Add test for nearly collinear points

---
 tests/manifold/spherical_manifold_07.cc     | 62 +++++++++++++++++++++
 tests/manifold/spherical_manifold_07.output | 20 +++++++
 2 files changed, 82 insertions(+)
 create mode 100644 tests/manifold/spherical_manifold_07.cc
 create mode 100644 tests/manifold/spherical_manifold_07.output

diff --git a/tests/manifold/spherical_manifold_07.cc b/tests/manifold/spherical_manifold_07.cc
new file mode 100644
index 0000000000..3aa2c56401
--- /dev/null
+++ b/tests/manifold/spherical_manifold_07.cc
@@ -0,0 +1,62 @@
+// ---------------------------------------------------------------------
+//
+// 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 that SphericalManifold also works in a stable manner
+// if two points are almost collinear with the center and another one
+// with a vanishing weight is located in the other hemisphere.
+
+#include "../tests.h"
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/full_matrix.h>
+
+using namespace dealii;
+
+int
+main()
+{
+  initlog();
+  Point<3> center(0.0, 0.0, 0.0);
+  const SphericalManifold<3,3> manifold(center);
+
+  for (unsigned int i=1; i<20; ++i)
+    {
+      const double eps = 1/std::pow(10,i*1.0);
+
+      Point<3> P1( -eps,  std::sqrt(1.-eps*eps), 0.);
+      Point<3> P2( -eps, -std::sqrt(1.-eps*eps), 0.);
+      Point<3> P3(   1.,                     0., 0.);
+
+      std::vector<Point<3> > points1(3);
+      std::vector<double > weights(3);
+
+      points1[0] = P1;
+      points1[1] = P2;
+      points1[2] = P3;
+
+      weights[0] = .5-eps;
+      weights[1] = .5-eps;
+      weights[2] = 2*eps;
+
+      Point<3> Q = manifold.get_new_point(make_array_view(points1),
+                                          make_array_view(weights));
+
+      deallog << Q << std::endl;
+    }
+
+  return 0;
+}
diff --git a/tests/manifold/spherical_manifold_07.output b/tests/manifold/spherical_manifold_07.output
new file mode 100644
index 0000000000..bd37c3774a
--- /dev/null
+++ b/tests/manifold/spherical_manifold_07.output
@@ -0,0 +1,20 @@
+
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
+DEAL::1.00000 0.00000 0.00000
-- 
2.39.5