From 310eccbfe3bb91efae8c78e5e40101088aeaaacf Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 7 Dec 2001 14:44:22 +0000 Subject: [PATCH] Cylinder domains git-svn-id: https://svn.dealii.org/trunk@5326 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_generator.h | 41 +++++ .../deal.II/include/grid/tria_boundary_lib.h | 122 ++++++++++++++ deal.II/deal.II/source/grid/grid_generator.cc | 127 ++++++++++++++ deal.II/deal.II/source/grid/grid_out.cc | 80 +-------- .../deal.II/source/grid/tria_boundary_lib.cc | 156 ++++++++++++++++++ 5 files changed, 447 insertions(+), 79 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index 1b010905c5..0e45017fc1 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -251,6 +251,47 @@ class GridGenerator const Point<3> ¢er = Point<3>(), const double radius = 1.); + /** + * Create a cylinder around the + * x-axis. The cylinder extends + * from @p{x=-half_length} to + * @p{x=+half_length} and its + * projection into the + * @p{yz}-plane is a circle of + * radius @p{radius}. + * + * The boundaries are colored + * according to the following + * scheme: 0 for the hull of the + * cylinder, 1 for the left hand + * face and 2 for the right hand + * face. + */ + static void cylinder (Triangulation<3> &tria, + const double radius = 1., + const double half_length = 1.); + + /** + * Projection of the + * three-dimensional cylinder + * into the @p{xy}-plane. + * Therefore, this is simply a square. + * + * Coloring is like in 3D. + */ + static void cylinder (Triangulation<2> &tria, + const double radius = 1., + const double half_length = 1.); + + /** + * Non-implemented dummy for compilation purposes. + */ + static void cylinder (Triangulation<1> &tria, + const double radius, + const double half_length); + + + /** * Initialize the given * triangulation with a hyper-L diff --git a/deal.II/deal.II/include/grid/tria_boundary_lib.h b/deal.II/deal.II/include/grid/tria_boundary_lib.h index 68ea8c8100..11144ef60d 100644 --- a/deal.II/deal.II/include/grid/tria_boundary_lib.h +++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h @@ -18,6 +18,128 @@ #include +/** + * Boundary object for the hull of a cylinder. In three dimensions, + * points are projected on a circular tube along the @p{x}-axis. The + * radius of the tube can be set. Similar to @ref{HyperBallBoundary}, + * new points are projected by dividing the straight line between the + * old two points and adjusting the radius in the @p{yz}-plane. + * + * This class was developed to be used in conjunction with the + * @p{cylinder} function of @ref{GridGenerator}. It should be used for + * the hull of the cylinder only (boundary indicator 0). + * + * This class is derived from @ref{StraightBoundary} rather than from + * @ref{Boundary}, which would seem natural, since this way we can use the + * @ref{StraightBoundary}@p{::in_between(neighbors)} function. + * + * @author Guido Kanschat, 2001 + */ +template +class CylinderBoundary : public StraightBoundary +{ + public: + /** + * Constructor + */ + CylinderBoundary (const double radius = 1.0); + + /** + * Refer to the general documentation of + * this class and the documentation of the + * base class. + */ + virtual Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; + + /** + * Refer to the general documentation of + * this class and the documentation of the + * base class. + */ + virtual Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; + + /** + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + * + * Calls + * @p{get_intermediate_points_between_points}. + */ + virtual void + get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, + typename std::vector > &points) const; + + /** + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + * + * Only implemented for @p{dim=3} + * and for @p{points.size()==1}. + */ + virtual void + get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, + typename std::vector > &points) const; + + /** + * Compute the normals to the + * boundary at the vertices of + * the given face. + * + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + */ + virtual void + get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const; + + /** + * Return the radius of the cylinder. + */ + double get_radius () const; + + /** + * Exception. Thrown by the + * @p{get_radius} if the + * @p{compute_radius_automatically}, + * see below, flag is set true. + */ + DeclException0 (ExcRadiusNotSet); + + + protected: + /** + * Radius of the cylinder. + */ + const double radius; + + private: + + /** + * Called by + * @p{get_intermediate_points_on_line} + * and by + * @p{get_intermediate_points_on_quad}. + * + * Refer to the general + * documentation of + * @p{get_intermediate_points_on_line} + * in the documentation of the + * base class. + */ + void get_intermediate_points_between_points (const Point &p0, const Point &p1, + typename std::vector > &points) const; +}; + + + /** * Specialisation of @ref{Boundary}, which places the new point on * the boundary of a ball in arbitrary dimension. It works by projecting diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index a0f3e382a0..b9130cf121 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -133,6 +133,7 @@ GridGenerator::colorize_hyper_rectangle (Triangulation &tria) cell->face(3)->set_boundary_indicator (1); cell->face(4)->set_boundary_indicator (5); cell->face(5)->set_boundary_indicator (0); + break; default: Assert(false, ExcNotImplemented()); }; @@ -186,6 +187,25 @@ void GridGenerator::hyper_ball (Triangulation<1> &, +void +GridGenerator::cylinder (Triangulation<2> &, + const double, + const double) +{ + Assert (false, ExcInternalError()); +} + + + +void GridGenerator::cylinder (Triangulation<1> &, + const double, + const double) +{ + Assert (false, ExcInternalError()); +}; + + + void GridGenerator::hyper_shell (Triangulation<1> &, const Point<1> &, const double, @@ -413,6 +433,38 @@ void GridGenerator::hyper_shell (Triangulation<2> &tria, +void +GridGenerator::cylinder (Triangulation<2> &tria, + const double radius, + const double half_length) +{ + Point<2> p1 (-half_length, -radius); + Point<2> p2 (half_length, radius); + + hyper_rectangle(tria, p1, p2, true); + + Triangulation<2>::face_iterator f = tria.begin_face(); + Triangulation<2>::face_iterator end = tria.end_face(); + while (f != end) + { + switch (f->boundary_indicator()) + { + case 0: + f->set_boundary_indicator(1); + break; + case 1: + f->set_boundary_indicator(2); + break; + default: + f->set_boundary_indicator(0); + break; + } + ++f; + } +} + + + void GridGenerator::half_hyper_shell (Triangulation<2> &tria, const Point<2> ¢er, @@ -697,6 +749,81 @@ GridGenerator::hyper_ball (Triangulation<3> &tria, +void +GridGenerator::cylinder (Triangulation<3> &tria, + const double radius, + const double half_length) +{ + // Copy the base from hyper_ball<2> + // and transform it to yz + const double d = radius/std::sqrt(2.0); + const double a = d/(1+std::sqrt(2.0)); + const Point<3> vertices[24] = { + Point<3>(-half_length, -d,-d), + Point<3>(-half_length, d,-d), + Point<3>(-half_length, -a,-a), + Point<3>(-half_length, a,-a), + Point<3>(-half_length, -a, a), + Point<3>(-half_length, a, a), + Point<3>(-half_length, -d, d), + Point<3>(-half_length, d, d), + Point<3>(0, -d,-d), + Point<3>(0, d,-d), + Point<3>(0, -a,-a), + Point<3>(0, a,-a), + Point<3>(0, -a, a), + Point<3>(0, a, a), + Point<3>(0, -d, d), + Point<3>(0, d, d), + Point<3>(half_length, -d,-d), + Point<3>(half_length, d,-d), + Point<3>(half_length, -a,-a), + Point<3>(half_length, a,-a), + Point<3>(half_length, -a, a), + Point<3>(half_length, a, a), + Point<3>(half_length, -d, d), + Point<3>(half_length, d, d), + }; + + int cell_vertices[10][8] = { + {0,1,3,2,8,9,11,10}, + {0,2,4,6,8,10,12,14}, + {2,3,5,4,10,11,13,12}, + {1,7,5,3,9,15,13,11}, + {6,4,5,7,14,12,13,15} + }; + for (unsigned int i=0;i<5;++i) + for (unsigned int j=0;j<8;++j) + cell_vertices[i+5][j] = cell_vertices[i][j]+8; + + std::vector > cells (10, CellData<3>()); + + for (unsigned int i=0; i<10; ++i) + { + for (unsigned int j=0; j<8; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + }; + + tria.create_triangulation (std::vector >(&vertices[0], &vertices[24]), + cells, + SubCellData()); // no boundary information + + Triangulation<3>::cell_iterator cell = tria.begin(); + Triangulation<3>::cell_iterator end = tria.end(); + + while (cell != end) + { + if (cell->face(0)->boundary_indicator() != 255) + cell->face(0)->set_boundary_indicator(1); + if (cell->face(1)->boundary_indicator() != 255) + cell->face(1)->set_boundary_indicator(2); + ++cell; + } +}; + + + void GridGenerator::hyper_shell (Triangulation<3> &, const Point<3> &, const double , diff --git a/deal.II/deal.II/source/grid/grid_out.cc b/deal.II/deal.II/source/grid/grid_out.cc index 968a2ffdd4..09706c5908 100644 --- a/deal.II/deal.II/source/grid/grid_out.cc +++ b/deal.II/deal.II/source/grid/grid_out.cc @@ -31,85 +31,7 @@ template void GridOut::write_dx (const Triangulation &tria, std::ostream &out) { - AssertThrow (out, ExcIO()); - - // get the positions of the - // vertices and whether they are - // used. - const std::vector > &vertices = tria.get_vertices(); - const std::vector &vertex_used = tria.get_used_vertices(); - - const unsigned int n_vertices = tria.n_used_vertices(); - - typename Triangulation::active_cell_iterator cell=tria.begin_active(); - const typename Triangulation::active_cell_iterator endc=tria.end(); - - // start with ucd data - out << n_vertices << ' ' - << tria.n_active_cells() + (ucd_flags.write_faces ? - n_boundary_faces(tria) : - 0) - << " 0 0 0" // no data - << std::endl; - - // actually write the vertices. - // note that we shall number them - // with first index 1 instead of 0 - for (unsigned int i=0; i(cell->material_id()) - << " "; - switch (dim) - { - case 1: out << "line "; break; - case 2: out << "quad "; break; - case 3: out << "hex "; break; - default: - Assert (false, ExcNotImplemented()); - }; - - // it follows a list of the - // vertices of each cell. in 1d - // this is simply a list of the - // two vertices, in 2d its counter - // clockwise, as usual in this - // library. in 3d, the same applies - // (special thanks to AVS for - // numbering their vertices in a - // way compatible to deal.II!) - // - // technical reference: - // AVS Developer's Guide, Release 4, - // May, 1992, p. E6 - // - // note: vertex numbers are 1-base - for (unsigned int vertex=0; vertex::vertices_per_cell; - ++vertex) - out << cell->vertex_index(vertex)+1 << ' '; - out << std::endl; - }; - - // write faces with non-zero boundary - // indicator - if (ucd_flags.write_faces) - write_ucd_faces (tria, cell_index, out); - + Assert(false, ExcNotImplemented()); AssertThrow (out, ExcIO()); }; diff --git a/deal.II/deal.II/source/grid/tria_boundary_lib.cc b/deal.II/deal.II/source/grid/tria_boundary_lib.cc index 27b39a1f43..678821e3f4 100644 --- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -20,6 +20,161 @@ #include +template +CylinderBoundary::CylinderBoundary (const double radius) : + radius(radius) +{}; + + + +template +Point +CylinderBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + Point middle = StraightBoundary::get_new_point_on_line (line); + + // project to boundary + if (dim>=3) + middle *= radius / std::sqrt(middle.square()-middle(0)*middle(0)); + return middle; +}; + + +template +Point +CylinderBoundary:: +get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + Point middle = StraightBoundary::get_new_point_on_quad (quad); + + // project to boundary + if (dim>=3) + middle *= radius / std::sqrt(middle.square()-middle(0)*middle(0)); + return middle; +}; + + +template +void +CylinderBoundary::get_intermediate_points_on_line ( + const typename Triangulation::line_iterator &line, + typename std::vector > &points) const +{ + if (points.size()==1) + points[0]=get_new_point_on_line(line); + else + get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points); +} + + +template +void +CylinderBoundary::get_intermediate_points_between_points ( + const Point &v0, const Point &v1, + typename std::vector > &points) const +{ + const unsigned int n=points.size(); + Assert(n>0, ExcInternalError()); + // Do a simple linear interpolation + // followed by projection. If this + // is not sufficient, try to + // understand the sophisticated + // code in HyperBall later. + Point ds = v1-v0; + ds /= n+1; + + for (unsigned int i=0; i +void +CylinderBoundary<3>::get_intermediate_points_on_quad ( + const Triangulation<3>::quad_iterator &quad, + std::vector > &points) const +{ + unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); + Assert(points.size()==m*m, ExcInternalError()); + + std::vector > lp3(m); + std::vector > lp1(m); + + get_intermediate_points_on_line(quad->line(3), lp3); + get_intermediate_points_on_line(quad->line(1), lp1); + + std::vector > lps(m); + for (unsigned int i=0; i +void +CylinderBoundary::get_intermediate_points_on_quad ( + const typename Triangulation::quad_iterator &, + typename std::vector > &) const +{ + Assert(false, Boundary::ExcFunctionNotUseful(dim)); +} + + + +#if deal_II_dimension == 1 + +template <> +void +CylinderBoundary<1>:: +get_normals_at_vertices (const Triangulation<1>::face_iterator &, + Boundary<1>::FaceVertexNormals &) const +{ + Assert (false, Boundary<1>::ExcFunctionNotUseful(1)); +}; + +#endif + + +template +void +CylinderBoundary:: +get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const +{ + for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) + { + face_vertex_normals[vertex] = face->vertex(vertex); + face_vertex_normals[vertex][0] = 0.; + } +}; + + + +template +double +CylinderBoundary::get_radius () const +{ + return radius; +}; + + +//======================================================================// + template HyperBallBoundary::HyperBallBoundary (const Point p, const double radius) : @@ -597,6 +752,7 @@ get_normals_at_vertices (const typename Triangulation::face_iterator &face, // explicit instantiations template class HyperBallBoundary; +template class CylinderBoundary; template class HalfHyperBallBoundary; template class HyperShellBoundary; template class HalfHyperShellBoundary; -- 2.39.5