From 7f9d702dc113f7da371595d7cc83a51e19ec8355 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 17 Aug 2009 17:26:47 +0000 Subject: [PATCH] ConeBoundary and GridGenerator::truncated_cone. Implementation by Markus Buerg. git-svn-id: https://svn.dealii.org/trunk@19288 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/grid/grid_generator.h | 114 +++-- .../deal.II/include/grid/tria_boundary_lib.h | 177 ++++++- deal.II/deal.II/source/grid/grid_generator.cc | 439 +++++++++++------- .../deal.II/source/grid/tria_boundary_lib.cc | 328 +++++++++++-- deal.II/doc/authors.html | 7 +- deal.II/doc/news/changes.h | 9 + 6 files changed, 831 insertions(+), 243 deletions(-) diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index 33b297fe75..39ef246420 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -37,7 +37,7 @@ template class SparseMatrix; * set, parts of the boundary receive different boundary numbers, * allowing them to be distinguished by application programs. See the * documentation of the functions for details. - * + * * Additionally this class provides a function * (@p laplace_transformation) that smoothly transforms a grid * according to given new boundary points. This can be used to @@ -147,7 +147,7 @@ class GridGenerator * adds 2i. For * instance, the center point * (1,-1,1) yields a material id 5. - * + * * @note The triangulation needs to be * void upon calling this * function. @@ -166,7 +166,7 @@ class GridGenerator * dimension @p i, * repetitions[i] cells are * generated. - * + * * To get cells with an aspect * ratio different from that of * the domain, use different @@ -179,7 +179,7 @@ class GridGenerator * of integers denoting the * number of subdivisions in each * coordinate direction. - * + * * If the @p colorize flag is * set, the * @p boundary_indicators of the @@ -279,7 +279,7 @@ class GridGenerator const Point &p, const Table &material_id, const bool colorize=false); - + /** * A parallelogram. The first * corner point is the @@ -304,8 +304,8 @@ class GridGenerator parallelogram(Triangulation& tria, const Tensor<2,dim>& corners, const bool colorize=false); - - + + /** * Hypercube with a layer of * hypercubes around it. The @@ -339,7 +339,7 @@ class GridGenerator const double right= 1., const double thickness = 1., const bool colorize = false); - + /** * Initialize the given * triangulation with a @@ -366,7 +366,7 @@ class GridGenerator * @note The triangulation needs to be * void upon calling this * function. - */ + */ template static void hyper_ball (Triangulation &tria, const Point ¢er = Point(), @@ -381,14 +381,14 @@ class GridGenerator * perpendicular to the * x-axis. * - * The boundary indicators for the final + * The boundary indicators for the final * triangulation are 0 for the curved boundary and - * 1 for the cut plane. + * 1 for the cut plane. * * The appropriate - * boundary class is + * boundary class is * HalfHyperBallBoundary, or HyperBallBoundary. - * + * * @note The triangulation needs to be * void upon calling this * function. @@ -429,7 +429,55 @@ class GridGenerator static void cylinder (Triangulation &tria, const double radius = 1., const double half_length = 1.); - + + /** + * Create a cutted cone around + * the x-axis. The cone extends + * from x=-half_length + * to x=half_length and + * its projection into the @p + * yz-plane is a circle of radius + * @p radius_0 at + * x=-half_length and a + * circle of radius @p radius_1 + * at x=+half_length. + * In between the radius is + * linearly decreasing. + * + * In two dimensions, the cone is + * a trapezoid from + * x=-half_length to + * x=+half_length and + * from y=-radius to + * y=radius_0 at + * x=-half_length and + * from y=-radius_1 to + * y=radius_1 at + * x=+half_length. In + * between the range of + * y is linearly + * decreasing. + * + * The boundaries are colored + * according to the following + * scheme: 0 for the hull of the + * cone, 1 for the left hand + * face and 2 for the right hand + * face. + * + * @note The triangulation needs to be + * void upon calling this + * function. + * + * @author Markus Bürg, 2009 + */ + template + static void + truncated_cone (Triangulation &tria, + const double radius_0 = 1.0, + const double radius_1 = 0.5, + const double half_length = 1.0); + /** * Initialize the given * triangulation with a hyper-L @@ -460,7 +508,7 @@ class GridGenerator static void hyper_L (Triangulation &tria, const double left = -1., const double right= 1.); - + /** * Initialize the given * Triangulation with a hypercube @@ -481,7 +529,7 @@ class GridGenerator * z-direction, such that * a plane cuts the lower half of * a rectangle in two. - + * This function is declared to * exist for triangulations of * all space dimensions, but @@ -497,7 +545,7 @@ class GridGenerator const double left = 0., const double right= 1., const bool colorize = false); - + /** * Produce a hyper-shell, * the region between two @@ -550,7 +598,7 @@ class GridGenerator const double outer_radius, const unsigned int n_cells = 0, bool colorize = false); - + /** * Produce a half hyper-shell, * i.e. the space between two @@ -595,7 +643,7 @@ class GridGenerator const double inner_radius, const double outer_radius, const unsigned int n_cells = 0); - + /** * Produce a domain that is the space * between two cylinders in 3d, with @@ -626,7 +674,7 @@ class GridGenerator const unsigned int n_radial_cells = 0, const unsigned int n_axial_cells = 0); - /** + /** * This class produces a square * on the xy-plane with a * circular hole in the middle, @@ -634,29 +682,29 @@ class GridGenerator * (only in 3d). * * @image html cubes_hole.png - * + * * It is implemented in 2d and * 3d, and takes the following * arguments: - * + * * @arg @p inner_radius: size of the - * internal hole + * internal hole * @arg @p outer_radius: size of the * biggest enclosed cylinder * @arg @p L: extension on the @p z-direction * @arg @p repetitions: number of subdivisions * along the @p z-direction - * @arg @p colorize: wether to assign different + * @arg @p colorize: wether to assign different * boundary indicators to different faces. - * The colors are given in lexicographic - * ordering for the flat faces (0 to 3 in 2d, - * 0 to 5 in 3d) plus the curved hole + * The colors are given in lexicographic + * ordering for the flat faces (0 to 3 in 2d, + * 0 to 5 in 3d) plus the curved hole * (4 in 2d, and 6 in 3d). - * If @p colorize is set to false, then flat faces + * If @p colorize is set to false, then flat faces * get the number 0 and the hole gets number 1. */ template - static void hyper_cube_with_cylindrical_hole (Triangulation &triangulation, + static void hyper_cube_with_cylindrical_hole (Triangulation &triangulation, const double inner_radius = .25, const double outer_radius = .5, const double L = .5, @@ -680,7 +728,7 @@ class GridGenerator const unsigned int n_rotations, const double R, const double r); - + /** * This function transformes the * @p Triangulation @p tria @@ -723,7 +771,7 @@ class GridGenerator */ DeclException1 (ExcInvalidRepetitionsDimension, int, - << "The vector of repetitions must have " + << "The vector of repetitions must have " << arg1 <<" elements."); private: @@ -766,7 +814,7 @@ class GridGenerator const Point& center, const double inner_radius, const double outer_radius); - + /** * Solve the Laplace equation for * @p laplace_transformation 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 bf4ad320c8..62c9c8bae2 100644 --- a/deal.II/deal.II/include/grid/tria_boundary_lib.h +++ b/deal.II/deal.II/include/grid/tria_boundary_lib.h @@ -142,8 +142,8 @@ class CylinderBoundary : public StraightBoundary * see below, flag is set true. */ DeclException0 (ExcRadiusNotSet); - - + + protected: /** * Radius of the cylinder. @@ -188,10 +188,163 @@ class CylinderBoundary : public StraightBoundary /** - * Specialisation of Boundary, which places the new point on + * Boundary object for the hull of a (truncated) cone with two + * different radii at the two ends. If one radius is chosen to be 0 + * the object describes the boundary of a cone. In three dimensions, + * points are projected on an arbitrarily oriented (truncated) cone + * described by the two endpoints and the corresponding radii. Similar + * to HyperBallBoundary, new points are projected by dividing the + * straight line between the old two points and adjusting the radius + * from the axis. + * + * This class is derived from StraightBoundary rather than from + * Boundary, which would seem natural, since this way we can use the + * StraightBoundary::in_between() function. + * + * @author Markus Bürg, 2009 + */ +template +class ConeBoundary : public StraightBoundary +{ + public: + /** + * Constructor. Here the boundary + * object is constructed. The + * points x_0 and + * x_1 describe the + * starting and ending points of + * the axis of the (cutted) + * cone. radius_0 + * denotes the radius + * corresponding to x_0 + * and radius_1 the one + * corresponding to x_1. + */ + ConeBoundary (const double radius_0, + const double radius_1, + const Point x_0, + const Point x_1); + + /** + * Return the radius of the + * (cutted) cone at given point + * x on the axis. + */ + double get_radius (const Point x) const; + + /** + * 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, + std::vector > &points) const; + + /** + * Refer to the general + * documentation of this class + * and the documentation of the + * base class. + * + * Only implemented for + * dim=3 and for + * points.size()==1. + */ + virtual + void + get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &quad, + 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; + + protected: + /** + * First radius of the (cutted) + * cone. + */ + const double radius_0; + + /** + * Second radius of the (cutted) + * cone. + */ + const double radius_1; + + /** + * Starting point of the axis. + */ + const Point x_0; + + /** + * Ending point of the axis. + */ + const Point x_1; + + 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, + std::vector > &points) const; +}; + + + +/** + * Specialization of Boundary, which places the new point on * the boundary of a ball in arbitrary dimension. It works by projecting * the point in the middle of the old points onto the ball. The middle is - * defined as the arithmetic mean of the points. + * defined as the arithmetic mean of the points. * * The center of the ball and its radius may be given upon construction of * an object of this type. They default to the origin and a radius of 1.0. @@ -287,10 +440,10 @@ class HyperBallBoundary : public StraightBoundary * see below, flag is set true. */ DeclException0 (ExcRadiusNotSet); - - + + protected: - + /** * Center point of the hyperball. */ @@ -335,7 +488,7 @@ class HyperBallBoundary : public StraightBoundary * base class. */ void get_intermediate_points_between_points (const Point &p0, const Point &p1, - std::vector > &points) const; + std::vector > &points) const; }; @@ -464,7 +617,7 @@ class HyperShellBoundary : public HyperBallBoundary * @author Wolfgang Bangerth, 2000, 2009 */ template -class HalfHyperShellBoundary : public HyperShellBoundary +class HalfHyperShellBoundary : public HyperShellBoundary { public: /** @@ -482,13 +635,13 @@ class HalfHyperShellBoundary : public HyperShellBoundary HalfHyperShellBoundary (const Point ¢er = Point(), const double inner_radius = -1, const double outer_radius = -1); - + /** * Construct a new point on a line. */ virtual Point - get_new_point_on_line (const typename Triangulation::line_iterator &line) const; - + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; + /** * Construct a new point on a quad. */ diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index 3d1135478a..586f3dcd26 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -70,7 +70,7 @@ namespace Point<3>( 0, 0,-1), Point<3>( 0, 0, 1) }; - + #endif } @@ -86,11 +86,11 @@ GridGenerator::hyper_rectangle (Triangulation &tria, // p1 is lower in all coordinate directions. Point p1(p_1); Point p2(p_2); - + for (unsigned int i=0;i p2(i)) std::swap (p1(i), p2(i)); - + std::vector > vertices (GeometryInfo::vertices_per_cell); switch (dim) { @@ -101,24 +101,24 @@ GridGenerator::hyper_rectangle (Triangulation &tria, case 2: vertices[0] = vertices[1] = p1; vertices[2] = vertices[3] = p2; - + vertices[1](0) = p2(0); vertices[2](0) = p1(0); break; case 3: vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1; vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2; - + vertices[1](0) = p2(0); vertices[2](1) = p2(1); vertices[3](0) = p2(0); vertices[3](1) = p2(1); - + vertices[4](0) = p1(0); vertices[4](1) = p1(1); vertices[5](1) = p1(1); vertices[6](0) = p1(0); - + break; default: Assert (false, ExcNotImplemented ()); @@ -174,7 +174,7 @@ void GridGenerator::hyper_cube (Triangulation &tria, { Assert (left < right, ExcMessage ("Invalid left and right bounds of hypercube")); - + Point p1; Point p2; @@ -204,12 +204,12 @@ GridGenerator::moebius ( Assert (n_cells>4, ExcMessage("More than 4 cells are needed to create a moebius grid.")); Assert (r>0 && R>0, ExcMessage("Outer and inner radius must be positive.")); Assert (R>r, ExcMessage("Outer radius must be greater than inner radius.")); - - + + std::vector > vertices (4*n_cells); double beta_step=n_rotations*numbers::PI/2.0/n_cells; double alpha_step=2.0*numbers::PI/n_cells; - + for (unsigned int i=0; i > cells (n_cells); for (unsigned int i=0; i > vertices (GeometryInfo::vertices_per_cell); vertices[1] = corners[0]; @@ -288,7 +288,7 @@ GridGenerator::subdivided_hyper_cube (Triangulation &tria, Assert (repetitions >= 1, ExcInvalidRepetitions(repetitions)); Assert (left < right, ExcMessage ("Invalid left and right bounds of hypercube")); - + // first generate the necessary // points const double delta = (right-left)/repetitions; @@ -340,7 +340,7 @@ GridGenerator::subdivided_hyper_cube (Triangulation &tria, cells[x].material_id = 0; } break; - + case 2: cells.resize (repetitions*repetitions); for (unsigned int y=0; y &tria, Assert (false, ExcNotImplemented()); } - tria.create_triangulation (points, cells, SubCellData()); + tria.create_triangulation (points, cells, SubCellData()); } @@ -402,14 +402,14 @@ GridGenerator::subdivided_hyper_rectangle ( { // contributed by Joerg R. Weimar // (j.weimar@jweimar.de) 2003 - Assert(repetitions.size() == dim, + Assert(repetitions.size() == dim, ExcInvalidRepetitionsDimension(dim)); // First, normalize input such that // p1 is lower in all coordinate // directions. Point p1(p_1); Point p2(p_2); - + for (unsigned int i=0;i p2(i)) std::swap (p1(i), p2(i)); @@ -419,14 +419,14 @@ GridGenerator::subdivided_hyper_rectangle ( // convert repetitions from double // to int by taking the ceiling. Point delta; - + for (unsigned int i=0; i= 1, ExcInvalidRepetitions(repetitions[i])); - + delta[i] = (p2[i]-p1[i])/repetitions[i]; } - + // then generate the necessary // points std::vector > points; @@ -473,7 +473,7 @@ GridGenerator::subdivided_hyper_rectangle ( } break; } - + case 2: { cells.resize (repetitions[1]*repetitions[0]); @@ -489,12 +489,12 @@ GridGenerator::subdivided_hyper_rectangle ( } break; } - + case 3: { const unsigned int n_x = (repetitions[0]+1); const unsigned int n_xy = (repetitions[0]+1)*(repetitions[1]+1); - + cells.resize (repetitions[2]*repetitions[1]*repetitions[0]); for (unsigned int z=0; z 0, ExcMessage ("The distance between corner points must be positive.")) - + // actual code is external since // 1-D is different from 2/3D. colorize_subdivided_hyper_rectangle (tria, p1, p2, epsilon); @@ -557,7 +557,7 @@ subdivided_hyper_rectangle(Triangulation &tria, // contributed by Joerg R. Weimar // (j.weimar@jweimar.de) 2003 // modified by Yaqi Wang 2006 - Assert(step_sz.size() == dim, + Assert(step_sz.size() == dim, ExcInvalidRepetitionsDimension(dim)); @@ -582,13 +582,13 @@ subdivided_hyper_rectangle(Triangulation &tria, } double x = 0; - for (unsigned int j=0; j > points; @@ -598,7 +598,7 @@ subdivided_hyper_rectangle(Triangulation &tria, { double x=0; for (unsigned int i=0; ; ++i) - { + { points.push_back (Point (p1[0]+x)); // form partial sums. in @@ -608,12 +608,12 @@ subdivided_hyper_rectangle(Triangulation &tria, // and exit early instead if (i == step_sizes[0].size()) break; - + x += step_sizes[0][i]; } break; } - + case 2: { double y=0; @@ -626,13 +626,13 @@ subdivided_hyper_rectangle(Triangulation &tria, p1[1]+y)); if (i == step_sizes[0].size()) break; - + x += step_sizes[0][i]; } - + if (j == step_sizes[1].size()) break; - + y += step_sizes[1][j]; } break; @@ -654,24 +654,24 @@ subdivided_hyper_rectangle(Triangulation &tria, p1[2]+z)); if (i == step_sizes[0].size()) break; - + x += step_sizes[0][i]; } if (j == step_sizes[1].size()) break; - + y += step_sizes[1][j]; } - + if (k == step_sizes[2].size()) break; - + z += step_sizes[2][k]; } break; } - + default: Assert (false, ExcNotImplemented()); } @@ -692,7 +692,7 @@ subdivided_hyper_rectangle(Triangulation &tria, } break; } - + case 2: { cells.resize (step_sizes[1].size()*step_sizes[0].size()); @@ -708,12 +708,12 @@ subdivided_hyper_rectangle(Triangulation &tria, } break; } - + case 3: { const unsigned int n_x = (step_sizes[0].size()+1); const unsigned int n_xy = (step_sizes[0].size()+1)*(step_sizes[1].size()+1); - + cells.resize (step_sizes[2].size()*step_sizes[1].size()*step_sizes[0].size()); for (unsigned int z=0; z &tria, cells[c].material_id = 0; } break; - + } default: Assert (false, ExcNotImplemented()); } - tria.create_triangulation (points, cells, SubCellData()); + tria.create_triangulation (points, cells, SubCellData()); if (colorize) { @@ -747,7 +747,7 @@ subdivided_hyper_rectangle(Triangulation &tria, // faces of all cells and set // boundary indicator to the // correct value if it was 0. - + // use a large epsilon to // compare numbers to avoid // roundoff problems. @@ -758,7 +758,7 @@ subdivided_hyper_rectangle(Triangulation &tria, *std::min_element (step_sizes[i].begin(), step_sizes[i].end())); const double epsilon = 0.01 * min_size; - + // actual code is external since // 1-D is different from 2/3D. colorize_subdivided_hyper_rectangle (tria, p1, p2, epsilon); @@ -779,12 +779,12 @@ subdivided_hyper_rectangle (Triangulation<1>& tria, const bool colorize) { // contributed by Yaqi Wang 2006 - Assert(spacing.size() == 1, + Assert(spacing.size() == 1, ExcInvalidRepetitionsDimension(1)); const unsigned int n_cells = material_id.size(0); - Assert(spacing[0].size() == n_cells, + Assert(spacing[0].size() == n_cells, ExcInvalidRepetitionsDimension(1)); double delta = std::numeric_limits::max(); @@ -792,7 +792,7 @@ subdivided_hyper_rectangle (Triangulation<1>& tria, Assert (spacing[0][i] >= 0, ExcInvalidRepetitions(-1)); delta = std::min (delta, spacing[0][i]); } - + // generate the necessary points std::vector > points; double ax = p[0]; @@ -820,7 +820,7 @@ subdivided_hyper_rectangle (Triangulation<1>& tria, SubCellData t; GridTools::delete_unused_vertices (points, cells, t); - tria.create_triangulation (points, cells, t); + tria.create_triangulation (points, cells, t); // set boundary indicator if (colorize) @@ -841,7 +841,7 @@ subdivided_hyper_rectangle (Triangulation<2>& tria, const bool colorize) { // contributed by Yaqi Wang 2006 - Assert(spacing.size() == 2, + Assert(spacing.size() == 2, ExcInvalidRepetitionsDimension(2)); std::vector repetitions(2); @@ -856,10 +856,10 @@ subdivided_hyper_rectangle (Triangulation<2>& tria, Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1)); delta = std::min (delta, spacing[i][j]); } - Assert(material_id.size(i) == repetitions[i], + Assert(material_id.size(i) == repetitions[i], ExcInvalidRepetitionsDimension(i)); } - + // generate the necessary points std::vector > points; double ay = p[1]; @@ -901,7 +901,7 @@ subdivided_hyper_rectangle (Triangulation<2>& tria, SubCellData t; GridTools::delete_unused_vertices (points, cells, t); - tria.create_triangulation (points, cells, t); + tria.create_triangulation (points, cells, t); // set boundary indicator if (colorize) @@ -943,8 +943,8 @@ subdivided_hyper_rectangle (Triangulation<3>& tria, { // contributed by Yaqi Wang 2006 const unsigned int dim = 3; - - Assert(spacing.size() == dim, + + Assert(spacing.size() == dim, ExcInvalidRepetitionsDimension(dim)); std::vector repetitions(dim); @@ -959,7 +959,7 @@ subdivided_hyper_rectangle (Triangulation<3>& tria, Assert (spacing[i][j] >= 0, ExcInvalidRepetitions(-1)); delta = std::min (delta, spacing[i][j]); } - Assert(material_id.size(i) == repetitions[i], + Assert(material_id.size(i) == repetitions[i], ExcInvalidRepetitionsDimension(i)); } @@ -1018,7 +1018,7 @@ subdivided_hyper_rectangle (Triangulation<3>& tria, SubCellData t; GridTools::delete_unused_vertices (points, cells, t); - tria.create_triangulation (points, cells, t); + tria.create_triangulation (points, cells, t); // set boundary indicator if (colorize) @@ -1077,14 +1077,14 @@ GridGenerator::colorize_subdivided_hyper_rectangle (Triangulation &tria, const Point &p2, const double epsilon) { - + // run through all faces and check // if one of their center coordinates matches // one of the corner points. Comparisons // are made using an epsilon which // should be smaller than the smallest cell // diameter. - + typename Triangulation::raw_face_iterator face = tria.begin_raw_face(), endface = tria.end_face(); for (; face!=endface; ++face) @@ -1110,7 +1110,7 @@ GridGenerator::colorize_subdivided_hyper_rectangle (Triangulation &tria, // but we could not find // on which boundary. Assert (false, ExcInternalError()); - + } } for (typename Triangulation::cell_iterator cell = tria.begin(); @@ -1188,6 +1188,18 @@ void GridGenerator::cylinder (Triangulation &, +// Implementation for 1D only +template +void GridGenerator::truncated_cone (Triangulation &, + const double, + const double, + const double) +{ + Assert (false, ExcNotImplemented()); +} + + + // Implementation for 1D only template void GridGenerator::hyper_shell (Triangulation &, @@ -1276,7 +1288,7 @@ void GridGenerator::enclosed_hyper_cube (Triangulation &tria, 1, 0, 2, 9, 8,10 }; - + std::vector > cells(9); k = 0; for (unsigned int i0=0;i0<3;++i0) @@ -1342,6 +1354,45 @@ GridGenerator::hyper_cube_slit (Triangulation &tria, } + +template +void GridGenerator::truncated_cone (Triangulation &triangulation, + const double radius_0, + const double radius_1, + const double half_length) +{ + Point vertices_tmp[4]; + + vertices_tmp[0] = Point (-half_length, -radius_0); + vertices_tmp[1] = Point (half_length, -radius_1); + vertices_tmp[2] = Point (-half_length, radius_0); + vertices_tmp[3] = Point (half_length, radius_1); + + const std::vector > vertices (&vertices_tmp[0], &vertices_tmp[4]); + unsigned int cell_vertices[1][GeometryInfo::vertices_per_cell]; + + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) + cell_vertices[0][i] = i; + + std::vector > cells (1, CellData ()); + + for (unsigned int i = 0; i < GeometryInfo::vertices_per_cell; ++i) + cells[0].vertices[i] = cell_vertices[0][i]; + + cells[0].material_id = 0; + triangulation.create_triangulation (vertices, cells, SubCellData ()); + + typename Triangulation::cell_iterator cell = triangulation.begin (); + + cell->face (0)->set_boundary_indicator (1); + cell->face (1)->set_boundary_indicator (2); + + for (unsigned int i = 2; i < 4; ++i) + cell->face (i)->set_boundary_indicator (0); +} + + + //TODO: Colorize edges as circumference, left and right radius // Implementation for 2D only template @@ -1363,14 +1414,14 @@ GridGenerator::hyper_L (Triangulation &tria, {3, 4, 6, 7}}; std::vector > cells (3, CellData()); - - for (unsigned int i=0; i<3; ++i) + + for (unsigned int i=0; i<3; ++i) { for (unsigned int j=0; j<4; ++j) cells[i].vertices[j] = cell_vertices[i][j]; cells[i].material_id = 0; }; - + tria.create_triangulation ( std::vector >(&vertices[0], &vertices[8]), cells, @@ -1398,7 +1449,7 @@ GridGenerator::hyper_ball (Triangulation &tria, p+Point(+1,+1)*(radius/std::sqrt(2.0)*a), p+Point(-1,+1)*(radius/std::sqrt(2.0)), p+Point(+1,+1)*(radius/std::sqrt(2.0)) }; - + const int cell_vertices[5][4] = {{0, 1, 2, 3}, {0, 2, 6, 4}, {2, 3, 4, 5}, @@ -1406,14 +1457,14 @@ GridGenerator::hyper_ball (Triangulation &tria, {6, 4, 7, 5}}; std::vector > cells (5, CellData()); - - for (unsigned int i=0; i<5; ++i) + + for (unsigned int i=0; i<5; ++i) { for (unsigned int j=0; j<4; ++j) cells[i].vertices[j] = cell_vertices[i][j]; cells[i].material_id = 0; }; - + tria.create_triangulation ( std::vector >(&vertices[0], &vertices[8]), cells, @@ -1433,9 +1484,9 @@ void GridGenerator::hyper_shell (Triangulation &tria, { Assert ((inner_radius > 0) && (inner_radius < outer_radius), ExcInvalidRadii ()); - + const double pi = numbers::PI; - + // determine the number of cells // for the grid. if not provided by // the user determine it such that @@ -1469,17 +1520,17 @@ void GridGenerator::hyper_shell (Triangulation &tria, }; std::vector > cells (N, CellData()); - - for (unsigned int i=0; i &tria, break; default: f->set_boundary_indicator(0); - break; + break; } ++f; } @@ -1574,21 +1625,21 @@ GridGenerator::half_hyper_ball (Triangulation &tria, p+Point(+1,+1)*(radius/std::sqrt(2.0)*a), p+Point(0,+1)*radius, p+Point(+1,+1)*(radius/std::sqrt(2.0)) }; - + const int cell_vertices[5][4] = {{0, 1, 2, 3}, {2, 3, 4, 5}, {1, 7, 3, 5}, {6, 4, 7, 5}}; std::vector > cells (4, CellData()); - - for (unsigned int i=0; i<4; ++i) + + for (unsigned int i=0; i<4; ++i) { for (unsigned int j=0; j<4; ++j) cells[i].vertices[j] = cell_vertices[i][j]; cells[i].material_id = 0; }; - + tria.create_triangulation ( std::vector >(&vertices[0], &vertices[8]), cells, @@ -1626,7 +1677,7 @@ GridGenerator::half_hyper_shell (Triangulation &tria, { Assert ((inner_radius > 0) && (inner_radius < outer_radius), ExcInvalidRadii ()); - + const double pi = numbers::PI; // determine the number of cells // for the grid. if not provided by @@ -1670,17 +1721,17 @@ GridGenerator::half_hyper_shell (Triangulation &tria, std::vector > cells (N, CellData()); - - for (unsigned int i=0; i& tria, { const double rl2=(right+left)/2; const double len = (right-left)/2.; - + const Point vertices[20] = { Point(left, left , -len/2.), Point(rl2, left , -len/2.), @@ -1784,7 +1835,7 @@ void GridGenerator::enclosed_hyper_cube (Triangulation &tria, 33,32,34, 41,40,42 }; - + std::vector > cells(27); k = 0; for (unsigned int z=0;z<3;++z) @@ -1811,6 +1862,80 @@ void GridGenerator::enclosed_hyper_cube (Triangulation &tria, +template +void GridGenerator::truncated_cone (Triangulation &triangulation, + const double radius_0, + const double radius_1, + const double half_length) +{ + Assert (dim <= 3, ExcNotImplemented()); + // Determine number of cells and vertices + const unsigned int n_cells = std::floor (half_length / std::max (radius_0, radius_1) + 0.5); + const unsigned int n_vertices = 4 * (n_cells + 1); + std::vector > vertices_tmp(n_vertices); + + vertices_tmp[0] = Point (-half_length, 0, -radius_0); + vertices_tmp[1] = Point (-half_length, radius_0, 0); + vertices_tmp[2] = Point (-half_length, -radius_0, 0); + vertices_tmp[3] = Point (-half_length, 0, radius_0); + + const double dx = 2 * half_length / n_cells; + + for (unsigned int i = 0; i < n_cells; ++i) { + vertices_tmp[4 * (i + 1)] + = vertices_tmp[4 * i] + + Point (dx, 0, 0.5 * (radius_0 - radius_1) * dx / half_length); + vertices_tmp[4 * i + 5] + = vertices_tmp[4 * i + 1] + + Point (dx, 0.5 * (radius_1 - radius_0) * dx / half_length, 0); + vertices_tmp[4 * i + 6] + = vertices_tmp[4 * i + 2] + + Point (dx, 0.5 * (radius_0 - radius_1) * dx / half_length, 0); + vertices_tmp[4 * i + 7] + = vertices_tmp[4 * i + 3] + + Point (dx, 0, 0.5 * (radius_1 - radius_0) * dx / half_length); + } + + const std::vector > vertices (&vertices_tmp[0], &vertices_tmp[n_vertices]); + Table<2,unsigned int> cell_vertices(n_cells,GeometryInfo::vertices_per_cell); + + for (unsigned int i = 0; i < n_cells; ++i) + for (unsigned int j = 0; j < GeometryInfo::vertices_per_cell; ++j) + cell_vertices[i][j] = 4 * i + j; + + std::vector > cells (n_cells, CellData ()); + + for (unsigned int i = 0; i < n_cells; ++i) { + for (unsigned int j = 0; j < GeometryInfo::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + + cells[i].material_id = 0; + } + + triangulation.create_triangulation (vertices, cells, SubCellData ()); + + for (typename Triangulation::cell_iterator cell = triangulation.begin (); + cell != triangulation.end (); ++cell) { + if (cell->vertex (0) (0) == -half_length) { + cell->face (4)->set_boundary_indicator (1); + + for (unsigned int i = 0; i < 4; ++i) + cell->line (i)->set_boundary_indicator (0); + } + + if (cell->vertex (4) (0) == half_length) { + cell->face (5)->set_boundary_indicator (2); + + for (unsigned int i = 4; i < 8; ++i) + cell->line (i)->set_boundary_indicator (0); + } + + for (unsigned int i = 0; i < 4; ++i) + cell->face (i)->set_boundary_indicator (0); + } +} + + // Implementation for 3D only template void @@ -1862,8 +1987,8 @@ GridGenerator::hyper_L (Triangulation &tria, {12, 13, 21, 22, 15, 16, 24, 25}}; std::vector > cells (7, CellData()); - - for (unsigned int i=0; i<7; ++i) + + for (unsigned int i=0; i<7; ++i) { for (unsigned int j=0; j<8; ++j) cells[i].vertices[j] = cell_vertices[i][j]; @@ -1923,10 +2048,10 @@ GridGenerator::hyper_ball (Triangulation &tria, {8, 0, 12, 4, 11, 3, 15, 7}, // left {8, 9, 0, 1, 11, 10, 3, 2}, // front {12, 4, 13, 5, 15, 7, 14, 6}}; // back - + std::vector > cells (n_cells, CellData()); - - for (unsigned int i=0; i::vertices_per_cell; ++j) cells[i].vertices[j] = cell_vertices[i][j]; @@ -1949,7 +2074,7 @@ GridGenerator::cylinder (Triangulation &tria, const double half_length) { Assert (dim <= 3, ExcNotImplemented()); - + // Copy the base from hyper_ball // and transform it to yz const double d = radius/std::sqrt(2.0); @@ -1987,7 +2112,7 @@ GridGenerator::cylinder (Triangulation &tria, vertices[i](1) = -vertices[i](0); vertices[i](0) = h; } - + int cell_vertices[10][8] = { {0, 1, 8, 9, 2, 3, 10, 11}, {0, 2, 8, 10, 6, 4, 14, 12}, @@ -1998,10 +2123,10 @@ GridGenerator::cylinder (Triangulation &tria, 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()); - - for (unsigned int i=0; i<10; ++i) + + for (unsigned int i=0; i<10; ++i) { for (unsigned int j=0; j<8; ++j) cells[i].vertices[j] = cell_vertices[i][j]; @@ -2025,7 +2150,7 @@ GridGenerator::cylinder (Triangulation &tria, // above typename Triangulation::cell_iterator cell = tria.begin(); typename Triangulation::cell_iterator end = tria.end(); - + for (; cell != end; ++cell) for (unsigned int i=0; i::faces_per_cell; ++i) if (cell->at_boundary(i)) @@ -2105,7 +2230,7 @@ GridGenerator::half_hyper_ball (Triangulation& tria, std::vector > cells (6, CellData()); - for (unsigned int i=0; i<6; ++i) + for (unsigned int i=0; i<6; ++i) { for (unsigned int j=0; j<8; ++j) cells[i].vertices[j] = cell_vertices[i][j]; @@ -2131,7 +2256,7 @@ GridGenerator::half_hyper_ball (Triangulation& tria, // element if (cell->face(i)->center()(0) < center(0)+1.e-5) { cell->face(i)->set_boundary_indicator(1); - for (unsigned int j=0;j::lines_per_face;++j) + for (unsigned int j=0;j::lines_per_face;++j) cell->face(i)->line(j)->set_boundary_indicator(1); } } @@ -2144,7 +2269,7 @@ GridGenerator::half_hyper_ball (Triangulation& tria, // If the center is not on the plane x=0, this is a curvilinear // element if (cell->face(i)->center()(0) > center(0)+1.e-5) { - for (unsigned int j=0;j::lines_per_face;++j) + for (unsigned int j=0;j::lines_per_face;++j) cell->face(i)->line(j)->set_boundary_indicator(0); } } @@ -2168,7 +2293,7 @@ void GridGenerator::hyper_shell (Triangulation& tria, const double orad = outer_radius/std::sqrt(3.0); std::vector > vertices; std::vector > cells; - + // Start with the shell bounded by // two nested cubes if (n <= 6) @@ -2177,7 +2302,7 @@ void GridGenerator::hyper_shell (Triangulation& tria, vertices.push_back(p+hexahedron[i]*irad); for (unsigned int i=0;i<8;++i) vertices.push_back(p+hexahedron[i]*orad); - + const unsigned int n_cells = 6; const int cell_vertices[n_cells][8] = {{8, 9, 10, 11, 0, 1, 2, 3}, // bottom @@ -2186,10 +2311,10 @@ void GridGenerator::hyper_shell (Triangulation& tria, {8, 0, 10, 2, 12, 4, 14, 6}, // left {8, 9, 0, 1, 12, 13, 4, 5}, // front {10, 2, 11, 3, 14, 6, 15, 7}}; // back - + cells.resize(n_cells, CellData()); - - for (unsigned int i=0; i::vertices_per_cell; ++j) cells[i].vertices[j] = cell_vertices[i][j]; @@ -2224,10 +2349,10 @@ void GridGenerator::hyper_shell (Triangulation& tria, { 2, 8, 11, 6}, { 12, 0, 2, 8}, { 1, 10, 12, 0}}; - + cells.resize(n_cells, CellData()); - - for (unsigned int i=0; i& tria, { Assert(false, ExcIndexRange(n, 1, 7)); } - + tria.create_triangulation (vertices, cells, SubCellData()); // no boundary // information - + if (colorize) colorize_hyper_shell(tria, p, inner_radius, outer_radius); } @@ -2271,7 +2396,7 @@ GridGenerator::colorize_hyper_shell ( (++cell)->face(2)->set_boundary_indicator(1); (++cell)->face(0)->set_boundary_indicator(1); (++cell)->face(2)->set_boundary_indicator(1); - (++cell)->face(0)->set_boundary_indicator(1); + (++cell)->face(0)->set_boundary_indicator(1); } else // For higher polyhedra, this is regular. @@ -2337,7 +2462,7 @@ GridGenerator::half_hyper_shell (Triangulation& tria, std::vector > cells (5, CellData()); - for (unsigned int i=0; i<5; ++i) + for (unsigned int i=0; i<5; ++i) { for (unsigned int j=0; j<8; ++j) cells[i].vertices[j] = cell_vertices[i][j]; @@ -2369,9 +2494,9 @@ void GridGenerator::cylinder_shell (Triangulation &tria, { Assert ((inner_radius > 0) && (inner_radius < outer_radius), ExcInvalidRadii ()); - + const double pi = numbers::PI; - + // determine the number of cells // for the grid. if not provided by // the user determine it such that @@ -2416,11 +2541,11 @@ void GridGenerator::cylinder_shell (Triangulation &tria, j*length/N_z); vertices_3d.push_back (v); } - + std::vector > cells (N_r*N_z, CellData()); - + for (unsigned int j=0; j &tria, cells[i+j*N_r].vertices[5] = N_r+((i+1)%N_r) + (j+1)*2*N_r; cells[i+j*N_r].vertices[6] = N_r+i + j*2*N_r; cells[i+j*N_r].vertices[7] = N_r+((i+1)%N_r) + j*2*N_r; - + cells[i+j*N_r].material_id = 0; } - + tria.create_triangulation ( vertices_3d, cells, SubCellData()); } @@ -2459,13 +2584,13 @@ void GridGenerator::laplace_solve (const SparseMatrix &S, PreconditionJacobi > prec; prec.initialize(S, 1.2); FilteredMatrix > PF (prec); - + SolverControl control (n_dofs, 1.e-10, false, false); GrowingVectorMemory > mem; SolverCG > solver (control, mem); - + Vector f(n_dofs); - + SF.add_constraints(m); SF.apply_constraints (f, true); solver.solve(SF, u, f, PF); @@ -2491,21 +2616,21 @@ void GridGenerator::laplace_transformation (Triangulation &tria, { // first provide everything that is // needed for solving a Laplace - // equation. - MappingQ1 mapping_q1; + // equation. + MappingQ1 mapping_q1; FE_Q q1(1); - + DoFHandler dof_handler(tria); dof_handler.distribute_dofs(q1); SparsityPattern sparsity_pattern (dof_handler.n_dofs (), dof_handler.n_dofs (), dof_handler.max_couplings_between_dofs()); DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); sparsity_pattern.compress (); - + SparseMatrix S(sparsity_pattern); - + QGauss4 quadrature; - + MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S); // set up the boundary values for @@ -2539,20 +2664,20 @@ void GridGenerator::laplace_transformation (Triangulation &tria, } } } - + // solve the dim problems with // different right hand sides. Vector us[dim]; for (unsigned int i=0; i tasks; for (unsigned int i=0; i &, #if deal_II_dimension == 2 template -void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation &triangulation, +void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation &triangulation, const double inner_radius, const double outer_radius, const double, // width, @@ -2632,7 +2757,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation &triang } } } - } + } } double eps = 1e-3 * outer_radius; cell = triangulation.begin_active(); @@ -2650,7 +2775,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation &triang cell->face(f)->set_boundary_indicator(2); else if(std::abs(dy - outer_radius) < eps) cell->face(f)->set_boundary_indicator(3); - else + else cell->face(f)->set_boundary_indicator(4); } else { double d = (cell->face(f)->center() - center).norm(); @@ -2660,7 +2785,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation &triang cell->face(f)->set_boundary_indicator(0); } } - } + } } #endif @@ -2668,7 +2793,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation &triang #if deal_II_dimension == 3 template -void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation &triangulation, +void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation &triangulation, const double inner_radius, const double outer_radius, const double L, @@ -2677,15 +2802,15 @@ void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation &triangu { Assert(inner_radius < outer_radius, ExcMessage("outer_radius has to be bigger than inner_radius.")); - Assert(L > 0, + Assert(L > 0, ExcMessage("Must give positive extension L")); Assert(Nz >= 1, ExcLowerRange(1, Nz)); - + GridGenerator::cylinder_shell (triangulation, L, inner_radius, outer_radius, 8, Nz); - + typename Triangulation::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); @@ -2718,7 +2843,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation &triangu } } } - } + } } double eps = 1e-3 * outer_radius; cell = triangulation.begin_active(); @@ -2728,26 +2853,26 @@ void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation &triangu double dx = cell->face(f)->center()(0); double dy = cell->face(f)->center()(1); double dz = cell->face(f)->center()(2); - + if(colorize) { if(std::abs(dx + outer_radius) < eps) cell->face(f)->set_boundary_indicator(0); - + else if(std::abs(dx - outer_radius) < eps) cell->face(f)->set_boundary_indicator(1); - + else if(std::abs(dy + outer_radius) < eps) cell->face(f)->set_boundary_indicator(2); - + else if(std::abs(dy - outer_radius) < eps) cell->face(f)->set_boundary_indicator(3); - + else if(std::abs(dz) < eps) cell->face(f)->set_boundary_indicator(4); else if(std::abs(dz - L) < eps) cell->face(f)->set_boundary_indicator(5); - + else { cell->face(f)->set_boundary_indicator(6); for(unsigned int l=0; l::lines_per_face; ++l) @@ -2766,7 +2891,7 @@ void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation &triangu cell->face(f)->set_boundary_indicator(0); } } - } + } } #endif @@ -2789,14 +2914,14 @@ GridGenerator::hyper_rectangle ( const bool); template void -GridGenerator::subdivided_hyper_rectangle +GridGenerator::subdivided_hyper_rectangle (Triangulation &, const std::vector&, const Point&, const Point&, bool); template void -GridGenerator::subdivided_hyper_rectangle +GridGenerator::subdivided_hyper_rectangle (Triangulation &, const std::vector >&, const Point&, @@ -2822,6 +2947,10 @@ GridGenerator::cylinder ( Triangulation&, double, double); +template void +GridGenerator::truncated_cone ( + Triangulation&, double, double, double); + template void GridGenerator::hyper_L ( Triangulation&, double, double); @@ -2851,7 +2980,7 @@ GridGenerator::half_hyper_shell ( const Point&, double, double, unsigned int); -template void +template void GridGenerator::hyper_cube_with_cylindrical_hole ( Triangulation &, const double, const double, const double, const unsigned int, bool); 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 803508b3b7..871bf7d7b1 100644 --- a/deal.II/deal.II/source/grid/tria_boundary_lib.cc +++ b/deal.II/deal.II/source/grid/tria_boundary_lib.cc @@ -96,7 +96,7 @@ CylinderBoundary<3>:: get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const { const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad); - + // same algorithm as above const unsigned int dim = 3; const Point vector_from_axis = (middle-point_on_axis) - @@ -135,7 +135,7 @@ CylinderBoundary::get_intermediate_points_on_line ( get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points); } - + template void CylinderBoundary::get_intermediate_points_between_points ( @@ -181,18 +181,18 @@ CylinderBoundary<3>::get_intermediate_points_on_quad ( { unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); Assert(points.size()==m*m, ExcInternalError()); - + std::vector > lp0(m); std::vector > lp1(m); get_intermediate_points_on_line(quad->line(0), lp0); get_intermediate_points_on_line(quad->line(1), lp1); - + std::vector > lps(m); for (unsigned int i=0; i::face_iterator &face, const Point vector_from_axis = (vertex-point_on_axis) - ((vertex-point_on_axis) * direction) * direction; - + face_vertex_normals[v] = (vector_from_axis / vector_from_axis.norm()); } } @@ -248,12 +248,257 @@ get_normals_at_vertices (const typename Triangulation::face_iterator &face, template double -CylinderBoundary::get_radius () const +CylinderBoundary::get_radius () const { return radius; } +//======================================================================// + +template +ConeBoundary::ConeBoundary (const double radius_0, + const double radius_1, + const Point x_0, + const Point x_1) + : + radius_0 (radius_0), + radius_1 (radius_1), + x_0 (x_0), + x_1 (x_1) +{} + + + +template +double ConeBoundary::get_radius (Point x) const +{ + for (unsigned int i = 0; i < dim; ++i) + if ((x_1 (i) - x_0 (i)) != 0) + return (radius_1 - radius_0) * (x (i) - x_0 (i)) / (x_1 (i) - x_0 (i)) + radius_0; + + return 0; +} + + + +template +void +ConeBoundary:: +get_intermediate_points_between_points (const Point &p0, + const Point &p1, + std::vector > &points) const +{ + const unsigned int n = points.size (); + const Point axis = x_1 - x_0; + const Point dx = (p1 - p0) / (n + 1); + + Assert (n > 0, ExcInternalError ()); + + for (unsigned int i = 0; i < n; ++i) + { + // Compute the current point. + const Point x_i = p0 + (i + 1) * dx; + // To project this point on the + // boundary of the cone we first + // compute the orthogonal + // projection of this point onto + // the axis of the cone. + const double c = (x_i - x_0) * axis / axis.square (); + const Point x_ip = x_0 + c * axis; + // If the middle point is located + // on the axis hand over the + // current point itself, else + // hand over the projection of it + // on the boundary. + if (x_ip.norm () <= 1e-10 * x_i.norm ()) + points[i] = x_i; + else { + // Compute the projection of + // the middle point on the + // boundary of the cone. + points[i] = x_ip + get_radius (x_ip) * (x_i - x_ip) / (x_i - x_ip).norm (); + } + } +} + +template +Point +ConeBoundary:: +get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + const Point axis = x_1 - x_0; + // Compute the middle point of the line. + const Point middle = StraightBoundary::get_new_point_on_line (line); + // To project it on the boundary of + // the cone we first compute the + // orthogonal projection of the + // middle point onto the axis of + // the cone. + const double c = (middle - x_0) * axis / axis.square (); + const Point middle_p = x_0 + c * axis; + // If the middle point is located + // on the axis return the middle + // point itself, else return the + // projection of it on the + // boundary. + if (middle_p.norm () <= 1e-10 * middle.norm ()) + return middle; + else { + // Compute the projection of the + // middle point on the boundary + // of the cone. + return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm (); + } +} + + + +template<> +Point<3> +ConeBoundary<3>::get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const +{ + const Point<3> axis = x_1 - x_0; + // Compute the middle point of the + // quad. + const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad); + // Same algorithm as above: To + // project it on the boundary of + // the cone we first compute the + // orthogonal projection of the + // middle point onto the axis of + // the cone. + const double c = (middle - x_0) * axis / axis.square (); + const Point<3> middle_p = x_0 + c * axis; + // If the middle point is located + // on the axis return the middle + // point itself, else return the + // projection of it on the + // boundary. + if (middle_p.norm () <= 1e-10 * middle.norm ()) + return middle; + else { + // Compute the projection of the + // middle point on the boundary + // of the cone. + return middle_p + get_radius (middle_p) * (middle - middle_p) / (middle - middle_p).norm (); + } +} + + + +template +Point +ConeBoundary:: +get_new_point_on_quad (const typename Triangulation::quad_iterator &) const +{ + Assert (false, ExcImpossibleInDim (dim)); + + return Point(); +} + + + +template +void +ConeBoundary:: +get_intermediate_points_on_line (const typename Triangulation::line_iterator &line, + 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); +} + + +#if deal_II_dimension == 3 + +template<> +void +ConeBoundary<3>:: +get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &quad, + std::vector > &points) const +{ + if (points.size () == 1) + points[0] = get_new_point_on_quad (quad); + else + { + unsigned int n = static_cast (std::sqrt (static_cast (points.size ()))); + + Assert (points.size () == n * n, ExcInternalError ()); + + std::vector > points_line_0 (n); + std::vector > points_line_1 (n); + + get_intermediate_points_on_line (quad->line (0), points_line_0); + get_intermediate_points_on_line (quad->line (1), points_line_1); + + std::vector > points_line_segment (n); + + for (unsigned int i = 0; i < n; ++i) + { + get_intermediate_points_between_points (points_line_0[i], + points_line_1[i], + points_line_segment); + + for (unsigned int j = 0; j < n; ++j) + points[i * n + j] = points_line_segment[j]; + } + } +} + +#endif + +template +void +ConeBoundary:: +get_intermediate_points_on_quad (const typename Triangulation::quad_iterator &, + std::vector > &) const +{ + Assert (false, ExcImpossibleInDim (dim)); +} + + +#if deal_II_dimension == 1 + +template<> +void +ConeBoundary<1>:: +get_normals_at_vertices (const Triangulation<1>::face_iterator &, + Boundary<1>::FaceVertexNormals &) const +{ + Assert (false, ExcImpossibleInDim (1)); +} + +#endif + +template +void +ConeBoundary:: +get_normals_at_vertices (const typename Triangulation::face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const +{ + const Point axis = x_1 - x_0; + + for (unsigned int vertex = 0; vertex < GeometryInfo::vertices_per_cell; ++vertex) + { + // Compute the orthogonal + // projection of the vertex onto + // the axis of the cone. + const double c = (face->vertex (vertex) - x_0) * axis / axis.square (); + const Point vertex_p = x_0 + c * axis; + // Then compute the vector + // pointing from the point + // vertex_p on the axis + // to the vertex. + const Point axis_to_vertex = face->vertex (vertex) - vertex_p; + + face_vertex_normals[vertex] = axis_to_vertex / axis_to_vertex.norm (); + } +} + + //======================================================================// template @@ -269,7 +514,7 @@ Point HyperBallBoundary::get_new_point_on_line (const typename Triangulation::line_iterator &line) const { Point middle = StraightBoundary::get_new_point_on_line (line); - + middle -= center; double r=0; @@ -282,7 +527,7 @@ HyperBallBoundary::get_new_point_on_line (const typename Triangula r=radius; // project to boundary middle *= r / std::sqrt(middle.square()); - + middle += center; return middle; } @@ -319,7 +564,7 @@ HyperBallBoundary:: get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { Point middle = StraightBoundary::get_new_point_on_quad (quad); - + middle -= center; double r=0; @@ -332,7 +577,7 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator r=radius; // project to boundary middle *= r / std::sqrt(middle.square()); - + middle += center; return middle; } @@ -364,7 +609,7 @@ HyperBallBoundary::get_intermediate_points_on_line ( get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points); } - + template void HyperBallBoundary::get_intermediate_points_between_points ( @@ -377,7 +622,7 @@ HyperBallBoundary::get_intermediate_points_between_points ( const Point v0=p0-center, v1=p1-center; const double length=std::sqrt((v1-v0).square()); - + double eps=1e-12; double r=0; if (compute_radius_automatically) @@ -388,17 +633,17 @@ HyperBallBoundary::get_intermediate_points_between_points ( else r=radius; - + const double r2=r*r; Assert(std::fabs(v0.square()-r2) pm=0.5*(v0+v1); - + const double h=std::sqrt(pm.square()); - + double beta=0; unsigned int left_index=0, right_index=0; if ((n+1)%2==0) @@ -419,7 +664,7 @@ HyperBallBoundary::get_intermediate_points_between_points ( right_index=n/2; left_index=n/2-1; } - + // n even: m=n/2, // n odd: m=(n-1)/2 const unsigned int m=n/2; @@ -427,12 +672,12 @@ HyperBallBoundary::get_intermediate_points_between_points ( { Assert(right_index=1, ExcInternalError()); - + double d=h*std::tan(beta); points[right_index]=pm+d/length*(v1-v0); points[left_index]=pm-d/length*(v1-v0); } - + // project the points from the // straight line to the @@ -458,18 +703,18 @@ HyperBallBoundary<3>::get_intermediate_points_on_quad ( { unsigned int m=static_cast (std::sqrt(static_cast(points.size()))); Assert(points.size()==m*m, ExcInternalError()); - + std::vector > lp0(m); std::vector > lp1(m); - + get_intermediate_points_on_line(quad->line(0), lp0); get_intermediate_points_on_line(quad->line(1), lp1); - + std::vector > lps(m); for (unsigned int i=0; i::face_iterat template Point -HyperBallBoundary::get_center () const +HyperBallBoundary::get_center () const { return center; } @@ -536,7 +781,7 @@ HyperBallBoundary::get_center () const template double -HyperBallBoundary::get_radius () const +HyperBallBoundary::get_radius () const { Assert(!compute_radius_automatically, ExcRadiusNotSet()); return radius; @@ -557,7 +802,7 @@ HalfHyperBallBoundary::HalfHyperBallBoundary (const Point center, template Point HalfHyperBallBoundary:: -get_new_point_on_line (const typename Triangulation::line_iterator &line) const +get_new_point_on_line (const typename Triangulation::line_iterator &line) const { // check whether center of object is // at x==0, since then it belongs @@ -577,7 +822,7 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c template <> Point<1> HalfHyperBallBoundary<1>:: -get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const +get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const { Assert (false, ExcInternalError()); return Point<1>(); @@ -590,7 +835,7 @@ get_new_point_on_quad (const Triangulation<1>::quad_iterator &) const template Point HalfHyperBallBoundary:: -get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const { const Point quad_center = quad->center(); if (quad_center(0) == this->center(0)) @@ -728,7 +973,7 @@ HalfHyperShellBoundary::HalfHyperShellBoundary (const Point ¢er, template Point HalfHyperShellBoundary:: -get_new_point_on_line (const typename Triangulation::line_iterator &line) const +get_new_point_on_line (const typename Triangulation::line_iterator &line) const { switch (dim) { @@ -748,7 +993,7 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c // as in the base class return HyperShellBoundary::get_new_point_on_line (line); } - + // in 3d, a line is a straight // line if it is on the symmetry // plane and if not both of its @@ -756,7 +1001,7 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c // inner or outer sphere case 3: { - + if (((line->vertex(0)(0) == this->center(0)) && (line->vertex(1)(0) == this->center(0))) @@ -783,7 +1028,7 @@ get_new_point_on_line (const typename Triangulation::line_iterator &line) c default: Assert (false, ExcNotImplemented()); } - + return Point(); } @@ -822,7 +1067,7 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) c quad->vertex(2) + quad->vertex(3) )/4; const Point quad_center_offset = quad_center - this->center; - + if (std::fabs (quad->line(0)->center().distance(this->center) - quad->line(1)->center().distance(this->center)) < 1e-12 * outer_radius) @@ -830,7 +1075,7 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) c // lines 0 and 1 are radial const double needed_radius = quad->line(0)->center().distance(this->center); - + return (this->center + quad_center_offset/quad_center_offset.norm() * needed_radius); } @@ -841,14 +1086,14 @@ get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) c // lines 2 and 3 are radial const double needed_radius = quad->line(2)->center().distance(this->center); - + return (this->center + quad_center_offset/quad_center_offset.norm() * needed_radius); } else Assert (false, ExcInternalError()); } - + // otherwise we are on the outer or // inner part of the shell. proceed // as in the base class @@ -881,7 +1126,7 @@ get_intermediate_points_on_line (const typename Triangulation::line_iterato // as in the base class HyperShellBoundary::get_intermediate_points_on_line (line, points); } - + // in 3d, a line is a straight // line if it is on the symmetry // plane and if not both of its @@ -889,7 +1134,7 @@ get_intermediate_points_on_line (const typename Triangulation::line_iterato // inner or outer sphere case 3: { - + if (((line->vertex(0)(0) == this->center(0)) && (line->vertex(1)(0) == this->center(0))) @@ -927,7 +1172,7 @@ get_intermediate_points_on_quad (const typename Triangulation::quad_iterato std::vector > &points) const { Assert (dim < 3, ExcNotImplemented()); - + // check whether center of object is // at x==0, since then it belongs // to the plane part of the @@ -986,8 +1231,9 @@ get_normals_at_vertices (const typename Triangulation::face_iterator &face, // explicit instantiations -template class HyperBallBoundary; template class CylinderBoundary; +template class ConeBoundary; +template class HyperBallBoundary; template class HalfHyperBallBoundary; template class HyperShellBoundary; template class HalfHyperShellBoundary; diff --git a/deal.II/doc/authors.html b/deal.II/doc/authors.html index 8a817b40c2..ffc62c805f 100644 --- a/deal.II/doc/authors.html +++ b/deal.II/doc/authors.html @@ -60,6 +60,9 @@ contact with them, send email to authors at dealii.org.
  • Michael Anderson: Linear complexity grid reordering algorithm in 2d and 3d. +
  • Markus Bürg: + Conical meshes and geometries. +
  • John Burnell: Configuration on Microsoft Windows systems. @@ -81,7 +84,7 @@ contact with them, send email to authors at dealii.org. cubit journal file to export to ucd mesh format;   FEFieldFunction and ParsedFunction classes;   work on the codimension-one meshes, DoFHandler, and finite - elements;   + elements;   singular integration;   Step-34 tutorial program;   random bug fixes and enhancements. @@ -170,7 +173,7 @@ contact with them, send email to authors at dealii.org.
  • Yaqi Wang: The step-28 tutorial program; -   +   some grid generation functions.
  • Joshua White: diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index b5944929ba..e9e97e08a0 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -359,6 +359,15 @@ inconvenience this causes.

    deal.II

      +
    1. +

      + New: The function GridGenerator::truncated_cone() and the class ConeBoundary + can now be used to describe conical objects. +
      + (Markus Bürg 2009/08/17) +

      +
    2. +
    3. New: Instead of asking face or edge iterators for their boundary indicator -- 2.39.5