-// check ConeBoundary and GridGenerator::truncated_cone
+// check CylindricalManifold and GridGenerator::truncated_cone
#include "../tests.h"
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
template <int dim>
void check ()
{
+ AssertThrow(false, ExcNotImplemented());
+}
+
+template <>
+void check<2> ()
+{
+ constexpr int dim = 2;
+ Triangulation<dim> triangulation;
+ GridGenerator::truncated_cone (triangulation);
+
+ triangulation.refine_global (2);
+
+ GridOut().write_gnuplot (triangulation,
+ deallog.get_file_stream());
+}
+
+template <>
+void check<3> ()
+{
+ constexpr int dim = 3;
Triangulation<dim> triangulation;
GridGenerator::truncated_cone (triangulation);
- Point<dim> p1, p2;
- p1[0] = -1;
- p2[0] = 1;
- static const ConeBoundary<dim> boundary (1, 0.5, p1, p2);
- triangulation.set_boundary (0, boundary);
+ static const CylindricalManifold<dim> boundary;
+ triangulation.set_manifold (0, boundary);
triangulation.refine_global (2);
#include "../tests.h"
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/fe/mapping_c1.h>
-
-
template <int dim>
void check ()
{
+ AssertThrow(false, ExcNotImplemented());
+}
+
+template <>
+void check<2> ()
+{
+ constexpr int dim = 2;
+ deallog << "dim=" << dim << std::endl;
+
+ Triangulation<dim> triangulation;
+ const double r1 = 0.5, r2 = 1.0, halfl = 0.25;
+ GridGenerator::truncated_cone (triangulation, r1, r2, halfl);
+
+ triangulation.refine_global (2);
+
+ GridOut().write_gnuplot (triangulation,
+ deallog.get_file_stream());
+}
+
+template <>
+void check<3> ()
+{
+ constexpr int dim = 3;
deallog << "dim=" << dim << std::endl;
Triangulation<dim> triangulation;
- double r1 = 0.5, r2 = 1.0, halfl = 0.25;
+ const double r1 = 0.5, r2 = 1.0, halfl = 0.25;
GridGenerator::truncated_cone (triangulation, r1, r2, halfl);
- Point<dim> p1, p2;
- p1[0] = -halfl;
- p2[0] = halfl;
- static const ConeBoundary<dim> boundary (r1, r2, p1, p2);
- triangulation.set_boundary (0, boundary);
+ static const CylindricalManifold<dim> boundary;
+ triangulation.set_manifold (0, boundary);
triangulation.refine_global (2);
-// check ConeBoundary
+// check CylindricalManifold
#include "../tests.h"
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/grid/tria.h>
+#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/fe/mapping_c1.h>
-
-
template <int dim>
void check ()
{
+ AssertThrow (false, ExcNotImplemented());
+}
+
+template <>
+void check<2> ()
+{
+ constexpr int dim = 2;
+ Triangulation<dim> triangulation;
+ GridGenerator::truncated_cone (triangulation);
+ static const FlatManifold<dim> boundary;
+ triangulation.set_manifold (0, boundary);
+
+ triangulation.refine_global (2);
+
+
+ const typename Triangulation<dim>::active_cell_iterator cell=triangulation.begin_active();
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ const typename Triangulation<dim>::face_iterator face=cell->face(face_no);
+ if (face->boundary_id() != numbers::internal_face_boundary_id)
+ {
+ deallog << face->boundary_id() << std::endl;
+ typename Manifold<dim>::FaceVertexNormals normals;
+ boundary.get_normals_at_vertices(face, normals);
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+ deallog << face->vertex(v) << ": " << normals[v] << std::endl;
+ }
+ }
+}
+
+template <>
+void check<3> ()
+{
+ constexpr int dim = 3;
Triangulation<dim> triangulation;
GridGenerator::truncated_cone (triangulation);
- Point<dim> p1, p2;
- p1[0] = -1;
- p2[0] = 1;
- static const ConeBoundary<dim> boundary (1, 0.5, p1, p2);
- triangulation.set_boundary (0, boundary);
+ static const CylindricalManifold<dim> boundary;
+ triangulation.set_manifold (0, boundary);
triangulation.refine_global (2);
const typename Triangulation<dim>::active_cell_iterator cell=triangulation.begin_active();
for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
{
- typename Triangulation<dim>::face_iterator face=cell->face(face_no);
- typename Boundary<dim>::FaceVertexNormals normals;
- boundary.get_normals_at_vertices(face, normals);
- for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+ const typename Triangulation<dim>::face_iterator face=cell->face(face_no);
+ if (face->boundary_id() != numbers::internal_face_boundary_id)
{
- deallog << normals[v] << std::endl;
+ deallog << face->boundary_id() << std::endl;
+ typename Manifold<dim>::FaceVertexNormals normals;
+ boundary.get_normals_at_vertices(face, normals);
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+ deallog << face->vertex(v) << ": " << normals[v] << std::endl;
}
}
}
-DEAL::0.00000 -1.00000
-DEAL::0.00000 -1.00000
-DEAL::0.00000 -1.00000
-DEAL::0.00000 -1.00000
-DEAL::0.00000 -1.00000
-DEAL::0.00000 -1.00000
-DEAL::0.00000 -1.00000
-DEAL::0.00000 -1.00000
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 -0.382683 -0.923880
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 -0.382683 -0.923880
-DEAL::0.00000 0.382683 -0.923880
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 0.382683 -0.923880
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 0.382683 -0.923880
-DEAL::0.00000 0.382683 -0.923880
-DEAL::0.00000 -0.382683 -0.923880
-DEAL::0.00000 -0.382683 -0.923880
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 0.382683 -0.923880
-DEAL::0.00000 -0.382683 -0.923880
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 0.00000 -1.00000
-DEAL::0.00000 0.382683 -0.923880
-DEAL::0.00000 -0.382683 -0.923880
-DEAL::0.00000 0.00000 -1.00000
+DEAL::1
+DEAL::-1.00000 -1.00000: 1.00000 0.00000
+DEAL::-1.00000 -0.500000: 1.00000 0.00000
+DEAL::0
+DEAL::-1.00000 -1.00000: 0.242536 -0.970143
+DEAL::-0.500000 -0.875000: 0.242536 -0.970143
+DEAL::0
+DEAL::-1.00000 0.00000 -1.00000: -0.242536 3.81982e-16 0.970143
+DEAL::-1.00000 -0.382683 -0.923880: -0.242536 0.371257 0.896295
+DEAL::-0.500000 -5.55112e-17 -0.875000: -0.242536 3.81982e-16 0.970143
+DEAL::-0.500000 -0.334848 -0.808395: -0.242536 0.371257 0.896295
+DEAL::0
+DEAL::-1.00000 0.00000 -1.00000: -0.242536 -1.66567e-16 0.970143
+DEAL::-0.500000 -5.55112e-17 -0.875000: -0.242536 -7.15117e-16 0.970143
+DEAL::-1.00000 0.382683 -0.923880: -0.242536 -0.371257 0.896295
+DEAL::-0.500000 0.334848 -0.808395: -0.242536 -0.371257 0.896295
+DEAL::1
+DEAL::-1.00000 0.00000 -1.00000: 1.00000 0.00000 0.00000
+DEAL::-1.00000 0.382683 -0.923880: 1.00000 0.00000 0.00000
+DEAL::-1.00000 -0.382683 -0.923880: 1.00000 0.00000 0.00000
+DEAL::-1.00000 0.00000 -0.673880: 1.00000 0.00000 0.00000
-// check ConeBoundary<3>::normal_vector()
+// check CylindricalManifold<3>::normal_vector()
#include "../tests.h"
-#include <deal.II/grid/tria_boundary_lib.h>
-
-
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
void check ()
{
- const Point<3> p1, p2(0,0,2);
- const ConeBoundary<3> boundary (0, 1, p1, p2);
-
- // test a bunch of points that are randomly chosen. first (manually)
- // project them onto the cone (perpendicular to the z-axis) and then
- // compute the normal vector. we test that they are correct by
- // making sure that they (i) have unit length, (ii) are
- // perpendicular to the tangent vector to the cone that runs from
- // the vertex (at the origin) to the projected point, and (iii) are
- // perpendicular to the tangent vector that is tangent to the circle
- // the projected point sits on
- Point<3> points[] = { Point<3>(1,1,1),
- Point<3>(1,0,1),
- Point<3>(0,1,1),
-
- Point<3>(1,1,1.5),
- Point<3>(1,0,1.5),
- Point<3>(0,1,1.5)
- };
- for (unsigned int i=0; i<sizeof(points)/sizeof(points[0]); ++i)
- {
- const Point<2> radial_component (points[i][0],
- points[i][1]);
- const Point<2> projected_radial_component
- = radial_component/radial_component.norm()
- * (points[i][2]/2); // radius of cone at given z
-
- const Point<3> projected_point (projected_radial_component[0],
- projected_radial_component[1],
- points[i][2]);
-
- const Tensor<1,3> tangent_1 = projected_point / projected_point.norm();
- const Tensor<1,3> tangent_2 = cross_product_3d (tangent_1,
- Point<3>(0,0,1));
-
- // get the normal vector and test it
- const Tensor<1,3> normal_vector
- = boundary.normal_vector( /* dummy= */ Triangulation<3>::face_iterator(),
- projected_point);
-
- Assert (std::fabs(normal_vector.norm()-1) < 1e-12,
- ExcInternalError());
- Assert (std::fabs(normal_vector * tangent_1) < 1e-12,
- ExcInternalError());
- Assert (std::fabs(normal_vector * tangent_2) < 1e-12,
- ExcInternalError());
- }
-
+ constexpr int dim = 3;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::truncated_cone (triangulation);
+ GridTools::transform ([](const Point<dim> &p)
+ {
+ return Point<dim> (-p[1], p[0], p[2]);
+ },
+ triangulation);
+ static const CylindricalManifold<dim> boundary(1);
+ triangulation.set_manifold (0, boundary);
+
+ triangulation.refine_global (2);
+
+ for (const auto cell: triangulation.active_cell_iterators())
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ {
+ const auto face=cell->face(face_no);
+ if (face->boundary_id() == 0)
+ for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_face; ++v)
+ {
+ const Point<dim> vertex = face->vertex(v);
+ const Tensor<1,dim> tangent_1 ({-vertex(2), 0., vertex(0)});
+ const Tensor<1,dim> tangent_2 = vertex-Point<dim>(0,3,0);
+
+ // get the normal vector and test it
+ const Tensor<1,3> normal_vector = boundary.normal_vector(face, vertex);
+
+ Assert (std::fabs(normal_vector.norm()-1) < 1e-12,
+ ExcInternalError());
+ Assert (std::fabs(normal_vector * tangent_1) < 1e-12,
+ ExcInternalError());
+ Assert (std::fabs(normal_vector * tangent_2) < 1e-12,
+ ExcInternalError());
+ }
+ }
deallog << "OK" << std::endl;
}