From 9c2abb2432a2d08f232f931ab1a2025370be527b Mon Sep 17 00:00:00 2001
From: Luca Heltai <luca.heltai@sissa.it>
Date: Sat, 5 Mar 2016 16:53:37 +0100
Subject: [PATCH] Fixed quadrature points on faces for MappingManifold.

---
 source/fe/mapping_manifold.cc       |  6 +-
 tests/fe/mapping_manifold_07.cc     | 91 +++++++++++++++++++++++++++++
 tests/fe/mapping_manifold_07.output |  7 +++
 3 files changed, 99 insertions(+), 5 deletions(-)
 create mode 100644 tests/fe/mapping_manifold_07.cc
 create mode 100644 tests/fe/mapping_manifold_07.output

diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc
index b71f22269b..c8ce69650e 100644
--- a/source/fe/mapping_manifold.cc
+++ b/source/fe/mapping_manifold.cc
@@ -396,14 +396,11 @@ namespace internal
 
       if (update_flags & update_quadrature_points)
         {
-          AssertDimension(quadrature_points.size(),
-                          data.cell_manifold_quadrature_weights.size());
-
           for (unsigned int point=0; point<quadrature_points.size(); ++point)
             {
               quadrature_points[point] = data.manifold->
                                          get_new_point(Quadrature<spacedim>(data.vertices,
-                                                       data.cell_manifold_quadrature_weights[point]));
+                                                       data.cell_manifold_quadrature_weights[point+data_set]));
             }
         }
     }
@@ -1329,7 +1326,6 @@ fill_fe_face_values (const typename Triangulation<dim,spacedim>::cell_iterator &
                                     quadrature,
                                     data,
                                     output_data);
-  Assert(false, ExcNotImplemented());
 }
 
 
diff --git a/tests/fe/mapping_manifold_07.cc b/tests/fe/mapping_manifold_07.cc
new file mode 100644
index 0000000000..ecaef6b656
--- /dev/null
+++ b/tests/fe/mapping_manifold_07.cc
@@ -0,0 +1,91 @@
+// ---------------------------------------------------------------------
+//
+// 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 MappingManifold on faces. Create a sphere around a one cell
+// Triangulation, and make sure all points are on the circle.
+
+#include "../tests.h"
+
+#include <deal.II/base/utilities.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/fe/mapping_q_generic.h>
+#include <deal.II/fe/mapping_manifold.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/grid/manifold_lib.h>
+
+template<int dim, int spacedim>
+void test()
+{
+  deallog << "dim=" << dim << ", spacedim=" << spacedim << std::endl;
+
+  Triangulation<dim, spacedim>   triangulation;
+
+  GridGenerator::hyper_cube (triangulation, 2.0, 3.0);
+
+  typename Triangulation<dim,spacedim>::active_cell_iterator
+  cell = triangulation.begin_active();
+
+  // Center and radius of the Ball
+  Point<spacedim> center = cell->center();
+  double radius = center.distance(cell->vertex(0));
+
+  static  const SphericalManifold<dim,spacedim> manifold(cell->center());
+
+  triangulation.set_all_manifold_ids_on_boundary(0);
+  triangulation.set_manifold (0, manifold);
+
+  MappingManifold<dim,spacedim> map_manifold;
+  FE_Q<dim,spacedim> fe(1);
+  const QGauss<dim-1> quad(5);
+
+  FEFaceValues<dim,spacedim> fe_v(map_manifold, fe, quad,
+                                  update_quadrature_points);
+
+
+
+  for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+    {
+      fe_v.reinit(cell, f);
+      const std::vector<Point<spacedim> > &qps = fe_v.get_quadrature_points();
+      for (unsigned int i=0; i<qps.size(); ++i)
+        if ( std::abs(qps[i].distance(center) - radius) > 1e-10)
+          {
+            deallog << "Expected radius: " << radius << ", got: "
+                    <<  qps[i].distance(center)
+                    << ", on point " << qps[i] << std::endl;
+          }
+    }
+  deallog << "OK" << std::endl;
+}
+
+
+int
+main()
+{
+  std::ofstream logfile ("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-10);
+
+  test<2,2>();
+  test<2,3>();
+
+  test<3,3>();
+
+  return 0;
+}
+
+
+
diff --git a/tests/fe/mapping_manifold_07.output b/tests/fe/mapping_manifold_07.output
new file mode 100644
index 0000000000..72b2a69be5
--- /dev/null
+++ b/tests/fe/mapping_manifold_07.output
@@ -0,0 +1,7 @@
+
+DEAL::dim=2, spacedim=2
+DEAL::OK
+DEAL::dim=2, spacedim=3
+DEAL::OK
+DEAL::dim=3, spacedim=3
+DEAL::OK
-- 
2.39.5