Fix support points on quad for MappingQGeneric and StraightBoundary.

diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc
index a7a48963d66102426f1726c6b5bfc24b853e226e..b9c6c65dc5ef86b9b671a52dab1145755bde643e 100644 (file)
--- a/source/fe/mapping_q_generic.cc
+++ b/source/fe/mapping_q_generic.cc
@@ -3775,9 +3775,9 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell,
 {
   const unsigned int faces_per_cell    = GeometryInfo<3>::faces_per_cell;
 
-  static const StraightBoundary<3> straight_boundary;
   // used if face quad at boundary or entirely in the interior of the domain
   std::vector<Point<3> > quad_points ((polynomial_degree-1)*(polynomial_degree-1));
+  std::vector<Point<3> > tmp_points;
 
   // loop over all faces and collect points on them
   for (unsigned int face_no=0; face_no<faces_per_cell; ++face_no)
@@ -3810,19 +3810,45 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell,
                ExcInternalError());
 #endif
 
-      // ask the boundary/manifold object to return intermediate points on it
-      get_intermediate_points_on_object(face->get_manifold(), line_support_points,
-                                        face, quad_points);
-
-      // in 3D, the orientation, flip and rotation of the face might not
-      // match what we expect here, namely the standard orientation. thus
-      // reorder points accordingly. since a Mapping uses the same shape
-      // function as an FE_Q, we can ask a FE_Q to do the reordering for us.
-      for (unsigned int i=0; i<quad_points.size(); ++i)
-        a.push_back(quad_points[fe_q->adjust_quad_dof_index_for_face_orientation(i,
-                                face_orientation,
-                                face_flip,
-                                face_rotation)]);
+      // On a quad, we have to check whether the manifold should determine the
+      // point distribution or rather a weighted sum should be created. This
+      // is the same logic as in the compute_mapping_support_points
+      // function below.
+      if (dynamic_cast<const Boundary<3,3> *>(&face->get_manifold()) == NULL)
+        {
+          // ask the boundary/manifold object to return intermediate points on it
+          get_intermediate_points_on_object(face->get_manifold(), line_support_points,
+                                            face, quad_points);
+
+          // in 3D, the orientation, flip and rotation of the face might not
+          // match what we expect here, namely the standard orientation. thus
+          // reorder points accordingly. since a Mapping uses the same shape
+          // function as an FE_Q, we can ask a FE_Q to do the reordering for us.
+          for (unsigned int i=0; i<quad_points.size(); ++i)
+            a.push_back(quad_points[fe_q->adjust_quad_dof_index_for_face_orientation(i,
+                                    face_orientation,
+                                    face_flip,
+                                    face_rotation)]);
+        }
+      else
+        {
+          // need to extract the points surrounding a quad from the points
+          // already computed. First get the 4 vertices and then the points on
+          // the four lines
+          tmp_points.resize(4 + 4*(polynomial_degree-1));
+          for (unsigned int v=0; v<GeometryInfo<2>::vertices_per_cell; ++v)
+            tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no,v)];
+          if (polynomial_degree > 1)
+            for (unsigned int line=0; line<GeometryInfo<2>::lines_per_cell; ++line)
+              for (unsigned int i=0; i<polynomial_degree-1; ++i)
+                tmp_points[4+line*(polynomial_degree-1)+i] =
+                  a[GeometryInfo<3>::vertices_per_cell +
+                    (polynomial_degree-1)*
+                    GeometryInfo<3>::face_to_cell_lines(face_no,line) + i];
+          add_weighted_interior_points (support_point_weights_on_quad, tmp_points);
+          a.insert(a.end(), tmp_points.begin()+4+4*(polynomial_degree-1),
+                   tmp_points.end());
+        }
     }
 }
 

diff --git a/tests/mappings/mapping_q_mixed_manifolds.cc b/tests/mappings/mapping_q_mixed_manifolds.cc
new file mode 100644 (file)
index 0000000..fa0840c
--- /dev/null
+++ b/tests/mappings/mapping_q_mixed_manifolds.cc
@@ -0,0 +1,201 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Shows the shape functions implemented and computes the area of cells.
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/fe/mapping_q_generic.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <vector>
+#include <fstream>
+#include <iomanip>
+#include <string>
+#include <sstream>
+
+const double D = 0.1;
+const double R = D/2.0;
+const double R_1 = 1.2*R;
+const double R_2 = 1.7*R;
+const double H = 0.41;
+const double X_0 = 0.0;
+const double X_1 = 0.3;
+const double X_C = 0.5; // center
+const double X_2 = 0.7;
+
+const double Y_0 = 0.0;
+const double Y_C = 0.2; // center
+
+const unsigned int MANIFOLD_ID = 1;
+
+
+void create_triangulation(Triangulation<2> &tria)
+{
+  AssertThrow(std::abs((X_2-X_1) - 2.0*(X_C-X_1))<1.0e-12, ExcMessage("Geometry parameters X_1,X_2,X_C invalid!"));
+  SphericalManifold<2> spherical_manifold(Point<2>(X_C,Y_C));
+
+  Triangulation<2> circle_1, circle_2, circle_tmp, middle, middle_tmp, middle_tmp2, tmp_3D;
+  std::vector<unsigned int> ref_1(2, 2);
+  ref_1[1] = 2;
+
+  // create middle part first as a hyper shell
+  const double outer_radius = (X_2-X_1)/2.0;
+  const unsigned int n_cells = 4;
+  GridGenerator::hyper_shell(middle, Point<2>(X_C, Y_C), R_2, outer_radius, n_cells, true);
+  middle.set_all_manifold_ids(MANIFOLD_ID);
+  middle.set_manifold(MANIFOLD_ID, spherical_manifold);
+  middle.refine_global(1);
+
+  // two inner circles in order to refine towards the cylinder surface
+  const unsigned int n_cells_circle = 8;
+  GridGenerator::hyper_shell(circle_1, Point<2>(X_C, Y_C), R, R_1, n_cells_circle, true);
+  circle_1.set_all_manifold_ids(MANIFOLD_ID);
+  circle_1.set_manifold(MANIFOLD_ID,spherical_manifold);
+
+  GridGenerator::hyper_shell(circle_2, Point<2>(X_C, Y_C), R_1, R_2, n_cells_circle, true);
+  circle_2.set_all_manifold_ids(MANIFOLD_ID);
+  circle_2.set_manifold(MANIFOLD_ID,spherical_manifold);
+
+  // then move the vertices to the points where we want them to be to create a slightly asymmetric cube with a hole
+  for (Triangulation<2>::cell_iterator cell = middle.begin(); cell != middle.end(); ++cell)
+    {
+      for (unsigned int v=0; v < GeometryInfo<2>::vertices_per_cell; ++v)
+        {
+          Point<2> &vertex = cell->vertex(v);
+          if (std::abs(vertex[0] - X_2) < 1e-10 && std::abs(vertex[1] - Y_C) < 1e-10)
+            {
+              vertex = Point<2>(X_2, H/2.0);
+            }
+          else if (std::abs(vertex[0] - (X_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10)
+            {
+              vertex = Point<2>(X_2, H);
+            }
+          else if (std::abs(vertex[0] - (X_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10)
+            {
+              vertex = Point<2>(X_2, Y_0);
+            }
+          else if (std::abs(vertex[0] - X_C) < 1e-10 && std::abs(vertex[1] - (Y_C +(X_2-X_1)/2.0)) < 1e-10)
+            {
+              vertex = Point<2>(X_C, H);
+            }
+          else if (std::abs(vertex[0] - X_C) < 1e-10 && std::abs(vertex[1] - (Y_C-(X_2-X_1)/2.0)) < 1e-10)
+            {
+              vertex = Point<2>(X_C, Y_0);
+            }
+          else if (std::abs(vertex[0] - (X_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10)
+            {
+              vertex = Point<2>(X_1, H);
+            }
+          else if (std::abs(vertex[0] - (X_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10)
+            {
+              vertex = Point<2>(X_1, Y_0);
+            }
+          else if (std::abs(vertex[0] - X_1) < 1e-10 && std::abs(vertex[1] - Y_C) < 1e-10)
+            {
+              vertex = Point<2>(X_1, H/2.0);
+            }
+        }
+    }
+
+  // must copy the triangulation because we cannot merge triangulations with refinement...
+  GridGenerator::flatten_triangulation(middle, middle_tmp);
+  GridGenerator::merge_triangulations(circle_1,circle_2,circle_tmp);
+  GridGenerator::merge_triangulations(middle_tmp,circle_tmp,tria);
+
+  // Set the cylinder boundary  to 2, outflow to 1, the rest to 0.
+  //tria.set_all_manifold_ids(0);
+  for (Triangulation<2>::active_cell_iterator cell=tria.begin(); cell != tria.end(); ++cell)
+    {
+      if (Point<2>(X_C,Y_C).distance(cell->center())<= R_2)
+        cell->set_all_manifold_ids(MANIFOLD_ID);
+    }
+}
+
+void create_triangulation(Triangulation<3> &tria)
+{
+  Triangulation<2> tria_2d;
+  create_triangulation(tria_2d);
+  GridGenerator::extrude_triangulation(tria_2d, 3, H, tria);
+
+  // Set the cylinder boundary  to 2, outflow to 1, the rest to 0.
+  tria.set_all_manifold_ids(0);
+  for (Triangulation<3>::active_cell_iterator cell=tria.begin(); cell != tria.end(); ++cell)
+    {
+      if (Point<3>(X_C,Y_C,cell->center()[2]).distance(cell->center())<= R_2)
+        cell->set_all_manifold_ids(MANIFOLD_ID);
+    }
+}
+
+
+
+template <int dim>
+void test()
+{
+  Point<dim> center;
+  center[0] = X_C;
+  center[1] = Y_C;
+  Point<dim> direction;
+  direction[dim-1] = 1.;
+
+  std_cxx11::shared_ptr<Manifold<dim> > cylinder_manifold
+  (dim == 2 ? static_cast<Manifold<dim>*>(new SphericalManifold<dim>(center)) :
+   static_cast<Manifold<dim>*>(new CylindricalManifold<dim>(direction, center)));
+  Triangulation<dim> tria;
+  create_triangulation(tria);
+  tria.set_manifold(MANIFOLD_ID, *cylinder_manifold);
+
+  FE_Nothing<dim> fe;
+  for (unsigned int degree = 1; degree < 7; ++degree)
+    {
+      MappingQGeneric<dim> mapping(degree);
+      QGauss<dim> quad(degree+1);
+      FEValues<dim> fe_values(mapping, fe, quad, update_JxW_values);
+      double sum = 0.;
+      for (typename Triangulation<dim>::active_cell_iterator cell=tria.begin_active();
+           cell != tria.end(); ++cell)
+        {
+          fe_values.reinit(cell);
+          double local_sum = 0;
+          for (unsigned int q=0; q<quad.size(); ++q)
+            local_sum += fe_values.JxW(q);
+          sum += local_sum;
+        }
+      const double reference = (dim==2 ? 1. : H) * (H*(X_2-X_1)-D*D*numbers::PI*0.25);
+      deallog << "Volume " << dim << "D mapping degree " << degree << ": "
+              << sum << " error: " << (sum-reference)/reference << std::endl;
+    }
+}
+
+
+int main()
+{
+  std::ofstream logfile ("output");
+  deallog.attach(logfile);
+  deallog.threshold_double(1.e-14);
+
+  test<2>();
+  test<3>();
+}

diff --git a/tests/mappings/mapping_q_mixed_manifolds.output b/tests/mappings/mapping_q_mixed_manifolds.output
new file mode 100644 (file)
index 0000000..19c6438
--- /dev/null
+++ b/tests/mappings/mapping_q_mixed_manifolds.output
@@ -0,0 +1,13 @@
+
+DEAL::Volume 2D mapping degree 1: 0.156929 error: 0.00501399
+DEAL::Volume 2D mapping degree 2: 0.156152 error: 3.91474e-05
+DEAL::Volume 2D mapping degree 3: 0.156147 error: 6.04868e-06
+DEAL::Volume 2D mapping degree 4: 0.156146 error: 5.88206e-08
+DEAL::Volume 2D mapping degree 5: 0.156146 error: 8.80543e-09
+DEAL::Volume 2D mapping degree 6: 0.156146 error: 8.99033e-11
+DEAL::Volume 3D mapping degree 1: 0.0643409 error: 0.00501399
+DEAL::Volume 3D mapping degree 2: 0.0640224 error: 3.91474e-05
+DEAL::Volume 3D mapping degree 3: 0.0640203 error: 6.04868e-06
+DEAL::Volume 3D mapping degree 4: 0.0640199 error: 5.88206e-08
+DEAL::Volume 3D mapping degree 5: 0.0640199 error: 8.80543e-09
+DEAL::Volume 3D mapping degree 6: 0.0640199 error: 8.99045e-11