From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Wed, 15 Feb 2017 08:38:19 +0000 (+0100)
Subject: Add additional test on cube.
X-Git-Tag: v8.5.0-rc1~17^2~8
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7f71113afcb66d6a8834112134708400360046f6;p=dealii.git

Add additional test on cube.
---

diff --git a/tests/mappings/mapping_q_mixed_manifolds_01.cc b/tests/mappings/mapping_q_mixed_manifolds_01.cc
new file mode 100644
index 0000000000..6c135b2ead
--- /dev/null
+++ b/tests/mappings/mapping_q_mixed_manifolds_01.cc
@@ -0,0 +1,95 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Compute the area of a square where one boundary is deformed by a
+// cylindrical manifold.
+
+#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_tools.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>
+
+
+
+template <int dim>
+void test()
+{
+  Point<dim> direction;
+  direction[dim-1] = 1.;
+
+  std_cxx11::shared_ptr<Manifold<dim> > cylinder_manifold
+  (dim == 2 ? static_cast<Manifold<dim>*>(new SphericalManifold<dim>(Point<dim>())) :
+   static_cast<Manifold<dim>*>(new CylindricalManifold<dim>(direction, Point<dim>())));
+
+  // create cube and shift it to position 0.7
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, -0.5, 0.5);
+  Tensor<1,dim> shift;
+  shift[0] = 1.;
+  GridTools::shift(shift, tria);
+  tria.begin()->face(0)->set_all_manifold_ids(1);
+  tria.set_manifold(1, *cylinder_manifold);
+
+  FE_Nothing<dim> fe;
+  for (unsigned int degree = 6; 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;
+        }
+      // reference area/volume is the area of the cube minus the area of the
+      // circular segment that is missing due to the manifold, whose area is
+      // R^2/2(theta-sin(theta)) with theta 90 degrees and R=sqrt(0.5)
+      const double reference = 1. - 0.25 * (0.5 * numbers::PI - 1.);
+      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.cc b/tests/mappings/mapping_q_mixed_manifolds_01.cc~
similarity index 100%
rename from tests/mappings/mapping_q_mixed_manifolds.cc
rename to tests/mappings/mapping_q_mixed_manifolds_01.cc~
diff --git a/tests/mappings/mapping_q_mixed_manifolds_01.output b/tests/mappings/mapping_q_mixed_manifolds_01.output
new file mode 100644
index 0000000000..f81c33011b
--- /dev/null
+++ b/tests/mappings/mapping_q_mixed_manifolds_01.output
@@ -0,0 +1,3 @@
+
+DEAL::Volume 2D mapping degree 6: 0.857305 error: 5.06446e-06
+DEAL::Volume 3D mapping degree 6: 0.857305 error: 5.06446e-06
diff --git a/tests/mappings/mapping_q_mixed_manifolds_02.cc b/tests/mappings/mapping_q_mixed_manifolds_02.cc
new file mode 100644
index 0000000000..27abc90135
--- /dev/null
+++ b/tests/mappings/mapping_q_mixed_manifolds_02.cc
@@ -0,0 +1,202 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Evaluates the area of a complex geometry when mixing different manifolds on
+// the surfaces than the rest of the cell. The mesh uses two layers of
+// cylinders that finally are changed into a square on the outside of the
+// domain.
+
+#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/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_02.output
similarity index 100%
rename from tests/mappings/mapping_q_mixed_manifolds.output
rename to tests/mappings/mapping_q_mixed_manifolds_02.output