--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
projects / dealii.git / commitdiff