--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 SphericalManifold for periodicity issues: check that the
+// spherical manifold finds the right intermediate points independent
+// on the number of surrounding points
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/manifold_lib.h>
+
+
+int
+main()
+{
+ initlog();
+
+ // Center and radius of the Ball
+ Point<2> center(.5, .5);
+ double radius = center.norm();
+
+ const SphericalManifold<2,2> manifold(center);
+
+ // Some points on the circle, that would cross the periodicity
+ // boundary
+ std::vector<Point<2> > points;
+
+ points.push_back(Point<2>(0.0, 0.0));
+ points.push_back(Point<2>(1.0, 0.0));
+
+ // And the weights
+ std::vector<double> weights(2);
+
+ unsigned int n_intermediates=10;
+
+ // First, use only two surrounding points. The second time around,
+ // instead, use four, with additional zero weights
+ for (unsigned int test_no=0; test_no<2; ++test_no)
+ {
+ if (test_no == 1)
+ {
+ // Test that we get the same result with four surrounding points
+ points.push_back(Point<2>(0.0, 1.0));
+ points.push_back(Point<2>(1.0, 1.0));
+ weights.push_back(0.0);
+ weights.push_back(0.0);
+ }
+
+ deallog << "Test " << test_no << std::endl;
+ for (unsigned int i=0; i<n_intermediates; ++i)
+ {
+ weights[0] = (double)i/((double)n_intermediates-1);
+ weights[1] = 1.0-weights[0];
+
+ deallog << manifold.get_new_point(Quadrature<2>(points, weights)) << std::endl;
+ }
+ }
+
+ return 0;
+}
+
+
+
--- /dev/null
+
+DEAL::Test 0
+DEAL::1.00000 -1.11022e-16
+DEAL::0.905580 -0.0792280
+DEAL::0.798836 -0.140856
+DEAL::0.683013 -0.183013
+DEAL::0.561628 -0.204416
+DEAL::0.438372 -0.204416
+DEAL::0.316987 -0.183013
+DEAL::0.201164 -0.140856
+DEAL::0.0944202 -0.0792280
+DEAL::-1.11022e-16 0.00000
+DEAL::Test 1
+DEAL::1.00000 -1.11022e-16
+DEAL::0.905580 -0.0792280
+DEAL::0.798836 -0.140856
+DEAL::0.683013 -0.183013
+DEAL::0.561628 -0.204416
+DEAL::0.438372 -0.204416
+DEAL::0.316987 -0.183013
+DEAL::0.201164 -0.140856
+DEAL::0.0944202 -0.0792280
+DEAL::-1.11022e-16 0.00000
// Compute the points with the faces
deallog << "Face quadratures ("
- << face_quad.size() * GeometryInfo<dim>::faces_per_cell
- << ")" << std::endl;
-
- for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f) {
- std::vector<Point<spacedim> > vertices;
- std::vector<double> weights(GeometryInfo<dim>::vertices_per_face);
-
- for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
- vertices.push_back(cell->face(f)->vertex(i));
-
- for (unsigned int i=0; i<face_quad.size(); ++i)
+ << face_quad.size() * GeometryInfo<dim>::faces_per_cell
+ << ")" << std::endl;
+
+ for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
{
- for(unsigned int v=0; v<weights.size(); ++v)
- weights[v] = face_fe_q.shape_value(v, face_quad.point(i));
+ std::vector<Point<spacedim> > vertices;
+ std::vector<double> weights(GeometryInfo<dim>::vertices_per_face);
- deallog << manifold.get_new_point(Quadrature<spacedim>(vertices, weights))
- << std::endl;
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
+ vertices.push_back(cell->face(f)->vertex(i));
+
+ for (unsigned int i=0; i<face_quad.size(); ++i)
+ {
+ for (unsigned int v=0; v<weights.size(); ++v)
+ weights[v] = face_fe_q.shape_value(v, face_quad.point(i));
+
+ deallog << manifold.get_new_point(Quadrature<spacedim>(vertices, weights))
+ << std::endl;
+ }
}
- }
// Project the face quads
Quadrature<dim> quad = QProjector<dim>::project_to_all_faces(face_quad);
deallog << "Face quadratures projected on cell "
- << "(" << quad.size() << ")" << std::endl;
+ << "(" << quad.size() << ")" << std::endl;
std::vector<Point<spacedim> > vertices;
- for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
+ for (unsigned int i=0; i<GeometryInfo<dim>::vertices_per_cell; ++i)
vertices.push_back(cell->vertex(i));
// Used to compute weights.
FE_Q<dim> fe_q(1);
std::vector<double> weights(GeometryInfo<dim>::vertices_per_cell);
-
+
for (unsigned int i=0; i<quad.size(); ++i)
{
- for(unsigned int v=0; v<weights.size(); ++v)
- weights[v] = fe_q.shape_value(v, quad.point(i));
+ for (unsigned int v=0; v<weights.size(); ++v)
+ weights[v] = fe_q.shape_value(v, quad.point(i));
deallog << manifold.get_new_point(Quadrature<spacedim>(vertices, weights))
- << std::endl;
+ << std::endl;
}
// for(unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f) {
// std::vector<Point<spacedim> > vertices;
// std::vector<double> weights(GeometryInfo<dim>::vertices_per_face);
-
+
// for(unsigned int i=0; i<GeometryInfo<dim>::vertices_per_face; ++i)
// vertices.push_back(cell->face(f)->vertex(i));
-
+
// for (unsigned int i=0; i<face_quad.size(); ++i)
// {
// for(unsigned int v=0; v<weights.size(); ++v)
- // weights[v] = face_fe_q.shape_value(v, face_quad.point(i));
+ // weights[v] = face_fe_q.shape_value(v, face_quad.point(i));
// deallog << manifold.get_new_point(Quadrature<spacedim>(vertices, weights))
- // << std::endl;
+ // << std::endl;
// }
// }
}