From 4003814b3c5aa11405f8b375d80050e62ced97dc Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 24 Jul 2013 12:58:44 +0000 Subject: [PATCH] Change GridGenerator from a class to a namespace. Undo introduction of local typedef size_type. git-svn-id: https://svn.dealii.org/trunk@30149 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 17 +- deal.II/include/deal.II/grid/grid_generator.h | 256 +- deal.II/source/grid/grid_generator.cc | 6659 +++++++++-------- deal.II/source/grid/grid_generator.inst.in | 32 +- 4 files changed, 3455 insertions(+), 3509 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 59758c9601..3af5c44c42 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -44,7 +44,22 @@ inconvenience this causes.

Specific improvements

    - +
  1. Changed: GridGenerator used to be a class with only static members + but is now a namespace, like all other similar constructs in deal.II. +
    + (Wolfgang Bangerth, 2013/07/24) +
    2. + +
  2. Changed: In GridGenerator, several functions had erroneously been changed + to take an argument of type size_type rather than unsigned + int. GridGenerator::size_type was a typedef to + types::global_dof_index, which for most users was unsigned int + anyway, but could also be set to be a 64-bit integer type. In any case, the + change has been reverted and these functions take just a regular + unsigned int again. +
    + (Wolfgang Bangerth, 2013/07/24) +
diff --git a/deal.II/include/deal.II/grid/grid_generator.h b/deal.II/include/deal.II/grid/grid_generator.h index 3967669860..58d4d3e2cf 100644 --- a/deal.II/include/deal.II/grid/grid_generator.h +++ b/deal.II/include/deal.II/grid/grid_generator.h @@ -33,34 +33,25 @@ template class SparseMatrix; /** - * This class provides a collection of functions for generating basic - * triangulations. Below, we try to provide some pictures in order to - * illustrate at least the more complex ones. + * This namespace provides a collection of functions for generating + * triangulations for some basic geometries. * * Some of these functions receive a flag @p colorize. If this is - * 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 - * transform (simple-shaped) grids to a more complicated ones, like a - * shell onto a grid of an airfoil, for example. - * - * No meshes for the codimension one case are provided at the moment. + * set, parts of the boundary receive different boundary indicators + * (@ref GlossBoundaryIndicator), + * allowing them to be distinguished for the purpose of attaching geometry + * objects and evaluating different boundary conditions. * + * This namespace also provides a function + * GridGenerator::laplace_transformation that smoothly transforms a domain + * into another one. This can be used to + * transform basic geometries to more complicated ones, like a + * shell to a grid of an airfoil, for example. * * @ingroup grid */ -class GridGenerator +namespace GridGenerator { -public: - /** - * Declare type for number of cell. - */ - typedef types::global_dof_index size_type; - /** * Initialize the given triangulation with a hypercube (line in 1D, square * in 2D, etc) consisting of exactly one cell. The hypercube volume is the @@ -79,9 +70,9 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void hyper_cube (Triangulation &tria, - const double left = 0., - const double right= 1.); + void hyper_cube (Triangulation &tria, + const double left = 0., + const double right= 1.); /** * Same as hyper_cube(), but with the difference that not only one cell is @@ -97,10 +88,10 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void subdivided_hyper_cube (Triangulation &tria, - const unsigned int repetitions, - const double left = 0., - const double right= 1.); + void subdivided_hyper_cube (Triangulation &tria, + const unsigned int repetitions, + const double left = 0., + const double right= 1.); /** * Create a coordinate-parallel brick from the two diagonally opposite @@ -118,10 +109,10 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void hyper_rectangle (Triangulation &tria, - const Point &p1, - const Point &p2, - const bool colorize = false); + void hyper_rectangle (Triangulation &tria, + const Point &p1, + const Point &p2, + const bool colorize = false); /** * Create a coordinate-parallel parallelepiped from the two diagonally @@ -155,7 +146,6 @@ public: * program. */ template - static void subdivided_hyper_rectangle (Triangulation &tria, const std::vector &repetitions, @@ -179,7 +169,6 @@ public: * specified by the points @p p1 and @p p2. */ template - static void subdivided_hyper_rectangle(Triangulation &tria, const std::vector > &step_sizes, @@ -196,7 +185,6 @@ public: * i.e. the domain will have a void there. */ template - static void subdivided_hyper_rectangle (Triangulation &tria, const std::vector< std::vector > &spacing, @@ -215,7 +203,6 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void parallelogram(Triangulation &tria, const Point (&corners)[dim], @@ -225,7 +212,6 @@ public: * @deprecated Use the other function of same name. */ template - static void parallelogram(Triangulation &tria, const Tensor<2,dim> &corners, @@ -249,7 +235,6 @@ public: * function. */ template - static void parallelepiped (Triangulation &tria, const Point (&corners) [dim], @@ -268,10 +253,9 @@ public: * function. */ template - static void subdivided_parallelepiped (Triangulation &tria, - const size_type n_subdivisions, + const unsigned int n_subdivisions, const Point (&corners) [dim], const bool colorize = false); @@ -284,10 +268,9 @@ public: * function. */ template - static void subdivided_parallelepiped (Triangulation &tria, - const size_type ( n_subdivisions) [dim], + const unsigned int ( n_subdivisions) [dim], const Point (&corners) [dim], const bool colorize = false); @@ -307,11 +290,11 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void enclosed_hyper_cube (Triangulation &tria, - const double left = 0., - const double right= 1., - const double thickness = 1., - const bool colorize = false); + void enclosed_hyper_cube (Triangulation &tria, + const double left = 0., + const double right= 1., + const double thickness = 1., + const bool colorize = false); /** * Initialize the given triangulation with a hyperball, i.e. a circle or a @@ -328,9 +311,9 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void hyper_ball (Triangulation &tria, - const Point ¢er = Point(), - const double radius = 1.); + void hyper_ball (Triangulation &tria, + const Point ¢er = Point(), + const double radius = 1.); /** * This class produces a half hyper-ball around center, which @@ -346,9 +329,9 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void half_hyper_ball (Triangulation &tria, - const Point ¢er = Point(), - const double radius = 1.); + void half_hyper_ball (Triangulation &tria, + const Point ¢er = Point(), + const double radius = 1.); /** * Create a cylinder around the x-axis. The cylinder extends from @@ -366,9 +349,9 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void cylinder (Triangulation &tria, - const double radius = 1., - const double half_length = 1.); + 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 @@ -395,7 +378,7 @@ public: * @author Markus Bürg, 2009 */ template - static void + void truncated_cone (Triangulation &tria, const double radius_0 = 1.0, const double radius_1 = 0.5, @@ -417,9 +400,9 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void hyper_L (Triangulation &tria, - const double left = -1., - const double right= 1.); + void hyper_L (Triangulation &tria, + const double left = -1., + const double right= 1.); /** * Initialize the given Triangulation with a hypercube with a slit. In each @@ -437,10 +420,10 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void hyper_cube_slit (Triangulation &tria, - const double left = 0., - const double right= 1., - const bool colorize = false); + void hyper_cube_slit (Triangulation &tria, + const double left = 0., + const double right= 1., + const bool colorize = false); /** * Produce a hyper-shell, the region between two spheres around @@ -491,12 +474,12 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void hyper_shell (Triangulation &tria, - const Point ¢er, - const double inner_radius, - const double outer_radius, - const size_type n_cells = 0, - bool colorize = false); + void hyper_shell (Triangulation &tria, + const Point ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells = 0, + bool colorize = false); /** * Produce a half hyper-shell, i.e. the space between two circles in two @@ -520,12 +503,12 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void half_hyper_shell (Triangulation &tria, - const Point ¢er, - const double inner_radius, - const double outer_radius, - const size_type n_cells = 0, - const bool colorize = false); + void half_hyper_shell (Triangulation &tria, + const Point ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells = 0, + const bool colorize = false); /** @@ -549,12 +532,12 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void quarter_hyper_shell (Triangulation &tria, - const Point ¢er, - const double inner_radius, - const double outer_radius, - const size_type n_cells = 0, - const bool colorize = false); + void quarter_hyper_shell (Triangulation &tria, + const Point ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells = 0, + const bool colorize = false); /** * Produce a domain that is the space between two cylinders in 3d, with @@ -569,12 +552,12 @@ public: * @note The triangulation needs to be void upon calling this function. */ template - static void cylinder_shell (Triangulation &tria, - const double length, - const double inner_radius, - const double outer_radius, - const size_type n_radial_cells = 0, - const size_type n_axial_cells = 0); + void cylinder_shell (Triangulation &tria, + const double length, + const double inner_radius, + const double outer_radius, + const unsigned int n_radial_cells = 0, + const unsigned int n_axial_cells = 0); @@ -592,9 +575,9 @@ public: * torus. */ - static void torus (Triangulation<2,3> &tria, - const double R, - const double r); + void torus (Triangulation<2,3> &tria, + const double R, + const double r); /** @@ -622,12 +605,12 @@ public: * get the number 0 and the hole gets number 1. */ template - static void hyper_cube_with_cylindrical_hole (Triangulation &triangulation, - const double inner_radius = .25, - const double outer_radius = .5, - const double L = .5, - const size_type repetition = 1, - const bool colorize = false); + void hyper_cube_with_cylindrical_hole (Triangulation &triangulation, + const double inner_radius = .25, + const double outer_radius = .5, + const double L = .5, + const unsigned int repetition = 1, + const bool colorize = false); /** * Produce a ring of cells in 3D that is cut open, twisted and glued @@ -639,11 +622,11 @@ public: * @param R The radius of the circle, which forms the middle line of the torus containing the loop of cells. Must be greater than @p r. * @param r The radius of the cylinder bend together as loop. */ - static void moebius (Triangulation<3,3> &tria, - const size_type n_cells, - const unsigned int n_rotations, - const double R, - const double r); + void moebius (Triangulation<3,3> &tria, + const unsigned int n_cells, + const unsigned int n_rotations, + const double R, + const double r); /** * Given the two triangulations specified as the first two arguments, create @@ -676,7 +659,6 @@ public: * GridTools::create_union_triangulation . */ template - static void merge_triangulations (const Triangulation &triangulation_1, const Triangulation &triangulation_2, @@ -690,10 +672,9 @@ public: * to the corresponding side walls in z direction. The bottom and top * get the next two free boundary indicators. */ - static void extrude_triangulation(const Triangulation<2, 2> &input, - const size_type n_slices, + const unsigned int n_slices, const double height, Triangulation<3,3> &result); @@ -712,9 +693,8 @@ public: * @deprecated This function has been moved to GridTools::laplace_transform */ template - static void laplace_transformation (Triangulation &tria, - const std::map > &new_points) DEAL_II_DEPRECATED; + const std::map > &new_points) DEAL_II_DEPRECATED; /** * Exception @@ -734,67 +714,7 @@ public: int, << "The vector of repetitions must have " << arg1 <<" elements."); - -private: - /** - * Perform the action specified by the @p colorize flag of the - * hyper_rectangle() function of this class. - */ - template - static - void - colorize_hyper_rectangle (Triangulation &tria); - - /** - * Perform the action specified by the @p colorize flag of the - * subdivided_hyper_rectangle() function of this class. This function is - * singled out because it is dimension specific. - */ - template - static - void - colorize_subdivided_hyper_rectangle (Triangulation &tria, - const Point &p1, - const Point &p2, - const double epsilon); - - /** - * Assign boundary number zero to the inner shell boundary and 1 to the - * outer. - */ - template - static - void - colorize_hyper_shell (Triangulation &tria, - const Point ¢er, - const double inner_radius, - const double outer_radius); - - - /** - * Assign boundary number zero the inner shell boundary, one to the outer - * shell boundary, two to the face with x=0, three to the face with y=0, - * four to the face with z=0. - */ - template - static - void - colorize_quarter_hyper_shell(Triangulation &tria, - const Point ¢er, - const double inner_radius, - const double outer_radius); - - /** - * Solve the Laplace equation for @p laplace_transformation function for one - * of the @p dim space dimensions. Externalized into a function of its own - * in order to allow parallel execution. - */ - static - void - laplace_solve (const SparseMatrix &S, - const std::map &m, - Vector &u); -}; +} diff --git a/deal.II/source/grid/grid_generator.cc b/deal.II/source/grid/grid_generator.cc index e834ca6334..443cc9d479 100644 --- a/deal.II/source/grid/grid_generator.cc +++ b/deal.II/source/grid/grid_generator.cc @@ -46,3784 +46,3793 @@ DEAL_II_NAMESPACE_OPEN -namespace +namespace GridGenerator { - // Corner points of the cube [-1,1]^3 - const Point<3> hexahedron[8] = + namespace { - Point<3>(-1,-1,-1), - Point<3>(+1,-1,-1), - Point<3>(-1,+1,-1), - Point<3>(+1,+1,-1), - Point<3>(-1,-1,+1), - Point<3>(+1,-1,+1), - Point<3>(-1,+1,+1), - Point<3>(+1,+1,+1) - }; - - // Octahedron inscribed in the cube - // [-1,1]^3 - const Point<3> octahedron[6] = - { - Point<3>(-1, 0, 0), - Point<3>( 1, 0, 0), - Point<3>( 0,-1, 0), - Point<3>( 0, 1, 0), - Point<3>( 0, 0,-1), - Point<3>( 0, 0, 1) - }; -} - + // Corner points of the cube [-1,1]^3 + const Point<3> hexahedron[8] = + { + Point<3>(-1,-1,-1), + Point<3>(+1,-1,-1), + Point<3>(-1,+1,-1), + Point<3>(+1,+1,-1), + Point<3>(-1,-1,+1), + Point<3>(+1,-1,+1), + Point<3>(-1,+1,+1), + Point<3>(+1,+1,+1) + }; -template -void -GridGenerator::hyper_rectangle (Triangulation &tria, - const Point &p_1, - const Point &p_2, - const bool colorize) -{ - // First, normalize input such that - // p1 is lower in all coordinate directions. - Point p1(p_1); - Point p2(p_2); + // Octahedron inscribed in the cube + // [-1,1]^3 + const Point<3> octahedron[6] = + { + Point<3>(-1, 0, 0), + Point<3>( 1, 0, 0), + Point<3>( 0,-1, 0), + Point<3>( 0, 1, 0), + Point<3>( 0, 0,-1), + Point<3>( 0, 0, 1) + }; - for (unsigned int i=0; i p2(i)) - std::swap (p1(i), p2(i)); - std::vector > vertices (GeometryInfo::vertices_per_cell); - switch (dim) + /** + * Perform the action specified by the @p colorize flag of the + * hyper_rectangle() function of this class. + */ + template + void + colorize_hyper_rectangle (Triangulation &tria) { - case 1: - vertices[0] = p1; - vertices[1] = p2; - break; - 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; + // there is nothing to do in 1d + if (dim > 1) + { + // there is only one cell, so + // simple task + const typename Triangulation::cell_iterator + cell = tria.begin(); + for (unsigned int f=0; f::faces_per_cell; ++f) + cell->face(f)->set_boundary_indicator (f); + } + } - 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 ()); + void + colorize_subdivided_hyper_rectangle (Triangulation<1> &tria, + const Point<1> &, + const Point<1> &, + const double) + { + for (Triangulation<1>::cell_iterator cell = tria.begin(); + cell != tria.end(); ++cell) + if (cell->center()(0) > 0) + cell->set_material_id(1); + // boundary indicators are set to + // 0 (left) and 1 (right) by default. } - // Prepare cell data - std::vector > cells (1); - for (unsigned int i=0; i::vertices_per_cell; ++i) - cells[0].vertices[i] = i; - cells[0].material_id = 0; - tria.create_triangulation (vertices, cells, SubCellData()); - // Assign boundary indicators - if (colorize) - colorize_hyper_rectangle (tria); -} + template + void + colorize_subdivided_hyper_rectangle (Triangulation &tria, + const Point &p1, + 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::face_iterator face = tria.begin_face(), + endface = tria.end_face(); + for (; face!=endface; ++face) + { + if (face->boundary_indicator() == 0) + { + const Point center (face->center()); + if (std::abs(center(0)-p1[0]) < epsilon) + face->set_boundary_indicator(0); + else if (std::abs(center(0) - p2[0]) < epsilon) + face->set_boundary_indicator(1); + else if (dim > 1 && std::abs(center(1) - p1[1]) < epsilon) + face->set_boundary_indicator(2); + else if (dim > 1 && std::abs(center(1) - p2[1]) < epsilon) + face->set_boundary_indicator(3); + else if (dim > 2 && std::abs(center(2) - p1[2]) < epsilon) + face->set_boundary_indicator(4); + else if (dim > 2 && std::abs(center(2) - p2[2]) < epsilon) + face->set_boundary_indicator(5); + else + // triangulation says it + // is on the boundary, + // but we could not find + // on which boundary. + Assert (false, ExcInternalError()); -template -void -GridGenerator::colorize_hyper_rectangle (Triangulation &tria) -{ - // there is nothing to do in 1d - if (dim > 1) - { - // there is only one cell, so - // simple task - const typename Triangulation::cell_iterator - cell = tria.begin(); - for (unsigned int f=0; f::faces_per_cell; ++f) - cell->face(f)->set_boundary_indicator (f); + } + } + for (typename Triangulation::cell_iterator cell = tria.begin(); + cell != tria.end(); ++cell) + { + char id = 0; + for (unsigned int d=0; dcenter()(d) > 0) id += 1 << d; + cell->set_material_id(id); + } } -} - -template -void GridGenerator::hyper_cube (Triangulation &tria, - const double left, - const double right) -{ - Assert (left < right, - ExcMessage ("Invalid left-to-right bounds of hypercube")); - - Point p1; - Point p2; + /** + * Assign boundary number zero to the inner shell boundary and 1 to the + * outer. + */ + void colorize_hyper_shell (Triangulation<1> &, + const Point<1> &, + const double, + const double) + { + Assert (false, ExcNotImplemented()); + } - p1(spacedim-1) = 0; - p2(spacedim-1) = 0; - for (unsigned int i=0; i &tria, + const Point<2> &, const double, const double) { - p1(i) = left; - p2(i) = right; + // In spite of receiving geometrical + // data, we do this only based on + // topology. + + // For the mesh based on cube, + // this is highly irregular + for (Triangulation<2>::cell_iterator cell = tria.begin(); + cell != tria.end(); ++cell) + { + Assert (cell->face(2)->at_boundary(), ExcInternalError()); + cell->face(2)->set_boundary_indicator(1); + } } - hyper_rectangle (tria, p1, p2); -} + /** + * Assign boundary number zero to the inner shell boundary and 1 to the + * outer. + */ + void + colorize_hyper_shell (Triangulation<3> &tria, + const Point<3> &, + const double, + const double) + { + // the following uses a good amount + // of knowledge about the + // orientation of cells. this is + // probably not good style... + if (tria.n_cells() == 6) + { + Triangulation<3>::cell_iterator cell = tria.begin(); -void -GridGenerator::moebius ( - Triangulation<3> &tria, - const size_type n_cells, - const unsigned int n_rotations, - const double R, - const double r) -{ - const unsigned int dim=3; - 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.")); + cell->face(4)->set_boundary_indicator(1); + Assert (cell->face(4)->at_boundary(), ExcInternalError()); + (++cell)->face(2)->set_boundary_indicator(1); + Assert (cell->face(2)->at_boundary(), ExcInternalError()); - 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; + (++cell)->face(2)->set_boundary_indicator(1); + Assert (cell->face(2)->at_boundary(), ExcInternalError()); - for (size_type i=0; iface(0)->set_boundary_indicator(1); + Assert (cell->face(0)->at_boundary(), ExcInternalError()); - size_type offset=0; + (++cell)->face(2)->set_boundary_indicator(1); + Assert (cell->face(2)->at_boundary(), ExcInternalError()); - std::vector > cells (n_cells); - for (size_type i=0; iface(0)->set_boundary_indicator(1); + Assert (cell->face(0)->at_boundary(), ExcInternalError()); + } + else if (tria.n_cells() == 12) + { + // again use some internal + // knowledge + for (Triangulation<3>::cell_iterator cell = tria.begin(); + cell != tria.end(); ++cell) + { + Assert (cell->face(5)->at_boundary(), ExcInternalError()); + cell->face(5)->set_boundary_indicator(1); + } + } + else if (tria.n_cells() == 96) { - cells[i].vertices[0+4*j]=offset+0+4*j; - cells[i].vertices[1+4*j]=offset+3+4*j; - cells[i].vertices[2+4*j]=offset+2+4*j; - cells[i].vertices[3+4*j]=offset+1+4*j; + // the 96-cell hypershell is + // based on a once refined + // 12-cell mesh. consequently, + // since the outer faces all + // are face_no==5 above, so + // they are here (unless they + // are in the interior). Use + // this to assign boundary + // indicators, but also make + // sure that we encounter + // exactly 48 such faces + unsigned int count = 0; + for (Triangulation<3>::cell_iterator cell = tria.begin(); + cell != tria.end(); ++cell) + if (cell->face(5)->at_boundary()) + { + cell->face(5)->set_boundary_indicator(1); + ++count; + } + Assert (count == 48, ExcInternalError()); } - offset+=4; - cells[i].material_id=0; + else + Assert (false, ExcNotImplemented()); } - // now correct the last four vertices - cells[n_cells-1].vertices[4]=(0+n_rotations)%4; - cells[n_cells-1].vertices[5]=(3+n_rotations)%4; - cells[n_cells-1].vertices[6]=(2+n_rotations)%4; - cells[n_cells-1].vertices[7]=(1+n_rotations)%4; - GridReordering::invert_all_cells_of_negative_grid(vertices,cells); - tria.create_triangulation_compatibility (vertices, cells, SubCellData()); -} + /** + * Assign boundary number zero the inner shell boundary, one to the outer + * shell boundary, two to the face with x=0, three to the face with y=0, + * four to the face with z=0. + */ + void + colorize_quarter_hyper_shell(Triangulation<3> &tria, + const Point<3> ¢er, + const double inner_radius, + const double outer_radius) + { + if (tria.n_cells() != 3) + AssertThrow (false, ExcNotImplemented()); + double middle = (outer_radius-inner_radius)/2e0 + inner_radius; + double eps = 1e-3*middle; + Triangulation<3>::cell_iterator cell = tria.begin(); -void -GridGenerator::torus (Triangulation<2,3> &tria, - const double R, - const double r) -{ - Assert (R>r, ExcMessage("Outer radius must be greater than inner radius.")); - - const unsigned int dim=2; - const unsigned int spacedim=3; - std::vector > vertices (16); - - vertices[0]=Point(R-r,0,0); - vertices[1]=Point(R,-r,0); - vertices[2]=Point(R+r,0,0); - vertices[3]=Point(R, r,0); - vertices[4]=Point(0,0,R-r); - vertices[5]=Point(0,-r,R); - vertices[6]=Point(0,0,R+r); - vertices[7]=Point(0,r,R); - vertices[8]=Point(-(R-r),0,0); - vertices[9]=Point(-R,-r,0); - vertices[10]=Point(-(R+r),0,0); - vertices[11]=Point(-R, r,0); - vertices[12]=Point(0,0,-(R-r)); - vertices[13]=Point(0,-r,-R); - vertices[14]=Point(0,0,-(R+r)); - vertices[15]=Point(0,r,-R); - - std::vector > cells (16); - //Right Hand Orientation - cells[0].vertices[0] = 0; - cells[0].vertices[1] = 4; - cells[0].vertices[2] = 7; - cells[0].vertices[3] = 3; - cells[0].material_id = 0; - - cells[1].vertices[0] = 1; - cells[1].vertices[1] = 5; - cells[1].vertices[2] = 4; - cells[1].vertices[3] = 0; - cells[1].material_id = 0; - - cells[2].vertices[0] = 2; - cells[2].vertices[1] = 6; - cells[2].vertices[2] = 5; - cells[2].vertices[3] = 1; - cells[2].material_id = 0; - - cells[3].vertices[0] = 3; - cells[3].vertices[1] = 7; - cells[3].vertices[2] = 6; - cells[3].vertices[3] = 2; - cells[3].material_id = 0; - - cells[4].vertices[0] = 4; - cells[4].vertices[1] = 8; - cells[4].vertices[2] = 11; - cells[4].vertices[3] = 7; - cells[4].material_id = 0; - - cells[5].vertices[0] = 5; - cells[5].vertices[1] = 9; - cells[5].vertices[2] = 8; - cells[5].vertices[3] = 4; - cells[5].material_id = 0; - - cells[6].vertices[0] = 6; - cells[6].vertices[1] = 10; - cells[6].vertices[2] = 9; - cells[6].vertices[3] = 5; - cells[6].material_id = 0; - - cells[7].vertices[0] = 7; - cells[7].vertices[1] = 11; - cells[7].vertices[2] = 10; - cells[7].vertices[3] = 6; - cells[7].material_id = 0; - - cells[8].vertices[0] = 8; - cells[8].vertices[1] = 12; - cells[8].vertices[2] = 15; - cells[8].vertices[3] = 11; - cells[8].material_id = 0; - - cells[9].vertices[0] = 9; - cells[9].vertices[1] = 13; - cells[9].vertices[2] = 12; - cells[9].vertices[3] = 8; - cells[9].material_id = 0; - - cells[10].vertices[0] = 10; - cells[10].vertices[1] = 14; - cells[10].vertices[2] = 13; - cells[10].vertices[3] = 9; - cells[10].material_id = 0; - - cells[11].vertices[0] = 11; - cells[11].vertices[1] = 15; - cells[11].vertices[2] = 14; - cells[11].vertices[3] = 10; - cells[11].material_id = 0; - - cells[12].vertices[0] = 12; - cells[12].vertices[1] = 0; - cells[12].vertices[2] = 3; - cells[12].vertices[3] = 15; - cells[12].material_id = 0; - - cells[13].vertices[0] = 13; - cells[13].vertices[1] = 1; - cells[13].vertices[2] = 0; - cells[13].vertices[3] = 12; - cells[13].material_id = 0; - - cells[14].vertices[0] = 14; - cells[14].vertices[1] = 2; - cells[14].vertices[2] = 1; - cells[14].vertices[3] = 13; - cells[14].material_id = 0; - - cells[15].vertices[0] = 15; - cells[15].vertices[1] = 3; - cells[15].vertices[2] = 2; - cells[15].vertices[3] = 14; - cells[15].material_id = 0; - - // Must call this to be able to create a - // correct triangulation in dealii, read - // GridReordering<> doc - GridReordering::reorder_cells (cells); - tria.create_triangulation_compatibility (vertices, cells, SubCellData()); -} + for (; cell!=tria.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + { + if (!cell->face(f)->at_boundary()) + continue; + double radius = cell->face(f)->center().norm() - center.norm(); + if (std::fabs(cell->face(f)->center()(0)) < eps ) // x = 0 set boundary 2 + { + cell->face(f)->set_boundary_indicator(2); + for (unsigned int j=0; j::lines_per_face; ++j) + if (cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(2); + } + else if (std::fabs(cell->face(f)->center()(1)) < eps) // y = 0 set boundary 3 + { + cell->face(f)->set_boundary_indicator(3); + for (unsigned int j=0; j::lines_per_face; ++j) + if (cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(3); + } + else if (std::fabs(cell->face(f)->center()(2)) < eps ) // z = 0 set boundary 4 + { + cell->face(f)->set_boundary_indicator(4); + for (unsigned int j=0; j::lines_per_face; ++j) + if (cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) + cell->face(f)->line(j)->set_boundary_indicator(4); + } + else if (radius < middle) // inner radius set boundary 0 + { + cell->face(f)->set_boundary_indicator(0); + for (unsigned int j=0; j::lines_per_face; ++j) + if (cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps) + cell->face(f)->line(j)->set_boundary_indicator(0); + } + else if (radius > middle) // outer radius set boundary 1 + { + cell->face(f)->set_boundary_indicator(1); + for (unsigned int j=0; j::lines_per_face; ++j) + if (cell->face(f)->line(j)->at_boundary()) + if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps) + cell->face(f)->line(j)->set_boundary_indicator(1); + } + else + AssertThrow (false, ExcInternalError()); + } + } -// Implementation for 2D only -template<> -void -GridGenerator::parallelogram (Triangulation<2> &tria, - const Point<2> (&corners)[2], - const bool colorize) -{ - std::vector > vertices (GeometryInfo<2>::vertices_per_cell); - - vertices[1] = corners[0]; - vertices[2] = corners[1]; - vertices[3] = vertices[1] + vertices[2]; - // Prepare cell data - std::vector > cells (1); - for (unsigned int i=0; i::vertices_per_cell; ++i) - cells[0].vertices[i] = i; - cells[0].material_id = 0; - - tria.create_triangulation (vertices, cells, SubCellData()); - - // Assign boundary indicators - if (colorize) - colorize_hyper_rectangle (tria); -} + } -template<> -void -GridGenerator::parallelogram (Triangulation<2> &tria, - const Tensor<2,2> &corners, - const bool colorize) -{ - // simply pass everything to the other function of same name - const Point<2> x[2] = { corners[0], corners[1] }; - parallelogram (tria, x, colorize); -} + template + void + hyper_rectangle (Triangulation &tria, + const Point &p_1, + const Point &p_2, + const bool colorize) + { + // 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)); + std::vector > vertices (GeometryInfo::vertices_per_cell); + switch (dim) + { + case 1: + vertices[0] = p1; + vertices[1] = p2; + break; + 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 ()); + } -// Parallelepiped implementation in 1d, 2d, and 3d. @note The -// implementation in 1d is similar to hyper_rectangle(), and in 2d is -// similar to parallelogram(). -// -// The GridReordering::reorder_grid is made use of towards the end of -// this function. Thus the triangulation is explicitly constructed for -// all dim here since it is slightly different in that respect -// (cf. hyper_rectangle(), parallelogram()). -template -void -GridGenerator::parallelepiped (Triangulation &tria, - const Point (&corners) [dim], - const bool colorize) -{ - // Check that none of the user defined vertices overlap - for (unsigned int i=0; i > cells (1); + for (unsigned int i=0; i::vertices_per_cell; ++i) + cells[0].vertices[i] = i; + cells[0].material_id = 0; - // Note: vertex[0] is the origin and is initialised as so here: - std::vector > vertices (GeometryInfo::vertices_per_cell); + tria.create_triangulation (vertices, cells, SubCellData()); - switch (dim) - { - // A line (1d parallelepiped) - case 1: - vertices[1] = corners[0]; - break; + // Assign boundary indicators + if (colorize) + colorize_hyper_rectangle (tria); + } - // A parallelogram (2d parallelepiped) - case 2: - // assign corners to vertices: - vertices[1] = corners[0]; - vertices[2] = corners[1]; - // compose the remaining vertex: - vertices[3] = vertices[1] + vertices[2]; - break; + template + void hyper_cube (Triangulation &tria, + const double left, + const double right) + { + Assert (left < right, + ExcMessage ("Invalid left-to-right bounds of hypercube")); - // A parallelepiped (3d parallelepiped) - case 3: - // assign corners to vertices: - vertices[1] = corners[0]; - vertices[2] = corners[1]; - vertices[4] = corners[2]; - - // compose the remaining vertices: - vertices[3] = vertices[1] + vertices[2]; - vertices[5] = vertices[1] + vertices[4]; - vertices[6] = vertices[2] + vertices[4]; - vertices[7] = vertices[1] + vertices[2] + vertices[4]; - break; + Point p1; + Point p2; - default: - Assert (false, ExcNotImplemented()); - } + p1(spacedim-1) = 0; + p2(spacedim-1) = 0; - // Prepare cell data and wipe material identity - std::vector > cells (1); - for (unsigned int i=0; i::vertices_per_cell; ++i) - cells[0].vertices[i] = i; - cells[0].material_id = 0; + for (unsigned int i=0; i::reorder_cells (cells); - tria.create_triangulation (vertices, cells, SubCellData()); - // Finally assign boundary indicators according to hyper_rectangle - if (colorize) - colorize_hyper_rectangle (tria); -} -template -void -GridGenerator::subdivided_parallelepiped (Triangulation &tria, - const size_type n_subdivisions, - const Point (&corners) [dim], - const bool colorize) -{ - // Equalise number of subdivisions in each dim-direction, heir - // validity will be checked later - size_type (n_subdivisions_) [dim]; - for (unsigned int i=0; i &tria, + const unsigned int n_cells, + const unsigned int n_rotations, + const double R, + const double r) + { + const unsigned int dim=3; + 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.")); -template -void -GridGenerator::subdivided_parallelepiped (Triangulation &tria, - const size_type ( n_subdivisions) [dim], - const Point (&corners) [dim], - const bool colorize) -{ - // Zero n_subdivisions is the origin only, which makes no sense - for (unsigned int i=0; i0, ExcInvalidRepetitions(n_subdivisions[i])); - // Check corners do not overlap (unique) - for (unsigned int i=0; i > vertices (4*n_cells); + double beta_step=n_rotations*numbers::PI/2.0/n_cells; + double alpha_step=2.0*numbers::PI/n_cells; - // Create a list of points - Point (delta) [dim]; + for (unsigned int i=0; i > points; + unsigned int offset=0; - switch (dim) - { - case 1: - for (size_type x=0; x<=n_subdivisions[0]; ++x) - points.push_back (Point (x*delta[0])); - break; + std::vector > cells (n_cells); + for (unsigned int i=0; i (x*delta[0] + y*delta[1])); - break; + // now correct the last four vertices + cells[n_cells-1].vertices[4]=(0+n_rotations)%4; + cells[n_cells-1].vertices[5]=(3+n_rotations)%4; + cells[n_cells-1].vertices[6]=(2+n_rotations)%4; + cells[n_cells-1].vertices[7]=(1+n_rotations)%4; - case 3: - for (size_type z=0; z<=n_subdivisions[2]; ++z) - for (size_type y=0; y<=n_subdivisions[1]; ++y) - for (size_type x=0; x<=n_subdivisions[0]; ++x) - points.push_back (Point (x*delta[0] + y*delta[1] + z*delta[2])); - break; + GridReordering::invert_all_cells_of_negative_grid(vertices,cells); + tria.create_triangulation_compatibility (vertices, cells, SubCellData()); + } - default: - Assert (false, ExcNotImplemented()); - } - // Prepare cell data - size_type n_cells = 1; - for (unsigned int i=0; i > cells (n_cells); - // Create fixed ordering of - switch (dim) - { - case 1: - for (size_type x=0; x &tria, + const double R, + const double r) + { + Assert (R>r, ExcMessage("Outer radius must be greater than inner radius.")); + + const unsigned int dim=2; + const unsigned int spacedim=3; + std::vector > vertices (16); + + vertices[0]=Point(R-r,0,0); + vertices[1]=Point(R,-r,0); + vertices[2]=Point(R+r,0,0); + vertices[3]=Point(R, r,0); + vertices[4]=Point(0,0,R-r); + vertices[5]=Point(0,-r,R); + vertices[6]=Point(0,0,R+r); + vertices[7]=Point(0,r,R); + vertices[8]=Point(-(R-r),0,0); + vertices[9]=Point(-R,-r,0); + vertices[10]=Point(-(R+r),0,0); + vertices[11]=Point(-R, r,0); + vertices[12]=Point(0,0,-(R-r)); + vertices[13]=Point(0,-r,-R); + vertices[14]=Point(0,0,-(R+r)); + vertices[15]=Point(0,r,-R); + + std::vector > cells (16); + //Right Hand Orientation + cells[0].vertices[0] = 0; + cells[0].vertices[1] = 4; + cells[0].vertices[2] = 7; + cells[0].vertices[3] = 3; + cells[0].material_id = 0; + + cells[1].vertices[0] = 1; + cells[1].vertices[1] = 5; + cells[1].vertices[2] = 4; + cells[1].vertices[3] = 0; + cells[1].material_id = 0; + + cells[2].vertices[0] = 2; + cells[2].vertices[1] = 6; + cells[2].vertices[2] = 5; + cells[2].vertices[3] = 1; + cells[2].material_id = 0; + + cells[3].vertices[0] = 3; + cells[3].vertices[1] = 7; + cells[3].vertices[2] = 6; + cells[3].vertices[3] = 2; + cells[3].material_id = 0; + + cells[4].vertices[0] = 4; + cells[4].vertices[1] = 8; + cells[4].vertices[2] = 11; + cells[4].vertices[3] = 7; + cells[4].material_id = 0; + + cells[5].vertices[0] = 5; + cells[5].vertices[1] = 9; + cells[5].vertices[2] = 8; + cells[5].vertices[3] = 4; + cells[5].material_id = 0; + + cells[6].vertices[0] = 6; + cells[6].vertices[1] = 10; + cells[6].vertices[2] = 9; + cells[6].vertices[3] = 5; + cells[6].material_id = 0; + + cells[7].vertices[0] = 7; + cells[7].vertices[1] = 11; + cells[7].vertices[2] = 10; + cells[7].vertices[3] = 6; + cells[7].material_id = 0; + + cells[8].vertices[0] = 8; + cells[8].vertices[1] = 12; + cells[8].vertices[2] = 15; + cells[8].vertices[3] = 11; + cells[8].material_id = 0; + + cells[9].vertices[0] = 9; + cells[9].vertices[1] = 13; + cells[9].vertices[2] = 12; + cells[9].vertices[3] = 8; + cells[9].material_id = 0; + + cells[10].vertices[0] = 10; + cells[10].vertices[1] = 14; + cells[10].vertices[2] = 13; + cells[10].vertices[3] = 9; + cells[10].material_id = 0; + + cells[11].vertices[0] = 11; + cells[11].vertices[1] = 15; + cells[11].vertices[2] = 14; + cells[11].vertices[3] = 10; + cells[11].material_id = 0; + + cells[12].vertices[0] = 12; + cells[12].vertices[1] = 0; + cells[12].vertices[2] = 3; + cells[12].vertices[3] = 15; + cells[12].material_id = 0; + + cells[13].vertices[0] = 13; + cells[13].vertices[1] = 1; + cells[13].vertices[2] = 0; + cells[13].vertices[3] = 12; + cells[13].material_id = 0; + + cells[14].vertices[0] = 14; + cells[14].vertices[1] = 2; + cells[14].vertices[2] = 1; + cells[14].vertices[3] = 13; + cells[14].material_id = 0; + + cells[15].vertices[0] = 15; + cells[15].vertices[1] = 3; + cells[15].vertices[2] = 2; + cells[15].vertices[3] = 14; + cells[15].material_id = 0; + + // Must call this to be able to create a + // correct triangulation in dealii, read + // GridReordering<> doc + GridReordering::reorder_cells (cells); + tria.create_triangulation_compatibility (vertices, cells, SubCellData()); + } - // wipe material id - cells[x].material_id = 0; - } - break; - case 2: - { - // Shorthand - const size_type n_dy = n_subdivisions[1]; - const size_type n_dx = n_subdivisions[0]; +// Implementation for 2D only + template<> + void + parallelogram (Triangulation<2> &tria, + const Point<2> (&corners)[2], + const bool colorize) + { + std::vector > vertices (GeometryInfo<2>::vertices_per_cell); + + vertices[1] = corners[0]; + vertices[2] = corners[1]; + vertices[3] = vertices[1] + vertices[2]; + // Prepare cell data + std::vector > cells (1); + for (unsigned int i=0; i::vertices_per_cell; ++i) + cells[0].vertices[i] = i; + cells[0].material_id = 0; + + tria.create_triangulation (vertices, cells, SubCellData()); + + // Assign boundary indicators + if (colorize) + colorize_hyper_rectangle (tria); + } + + + template<> + void + parallelogram (Triangulation<2> &tria, + const Tensor<2,2> &corners, + const bool colorize) + { + // simply pass everything to the other function of same name + const Point<2> x[2] = { corners[0], corners[1] }; + parallelogram (tria, x, colorize); + } + + + +// Parallelepiped implementation in 1d, 2d, and 3d. @note The +// implementation in 1d is similar to hyper_rectangle(), and in 2d is +// similar to parallelogram(). +// +// The GridReordering::reorder_grid is made use of towards the end of +// this function. Thus the triangulation is explicitly constructed for +// all dim here since it is slightly different in that respect +// (cf. hyper_rectangle(), parallelogram()). + template + void + parallelepiped (Triangulation &tria, + const Point (&corners) [dim], + const bool colorize) + { + // Check that none of the user defined vertices overlap + for (unsigned int i=0; i > vertices (GeometryInfo::vertices_per_cell); + + switch (dim) + { + // A line (1d parallelepiped) + case 1: + vertices[1] = corners[0]; + break; + + // A parallelogram (2d parallelepiped) + case 2: + // assign corners to vertices: + vertices[1] = corners[0]; + vertices[2] = corners[1]; + + // compose the remaining vertex: + vertices[3] = vertices[1] + vertices[2]; + break; + + // A parallelepiped (3d parallelepiped) + case 3: + // assign corners to vertices: + vertices[1] = corners[0]; + vertices[2] = corners[1]; + vertices[4] = corners[2]; + + // compose the remaining vertices: + vertices[3] = vertices[1] + vertices[2]; + vertices[5] = vertices[1] + vertices[4]; + vertices[6] = vertices[2] + vertices[4]; + vertices[7] = vertices[1] + vertices[2] + vertices[4]; + break; + + default: + Assert (false, ExcNotImplemented()); + } + + // Prepare cell data and wipe material identity + std::vector > cells (1); + for (unsigned int i=0; i::vertices_per_cell; ++i) + cells[0].vertices[i] = i; + cells[0].material_id = 0; + + // Check ordering of vertices and create triangulation + GridReordering::reorder_cells (cells); + tria.create_triangulation (vertices, cells, SubCellData()); + + // Finally assign boundary indicators according to hyper_rectangle + if (colorize) + colorize_hyper_rectangle (tria); + } + + template + void + subdivided_parallelepiped (Triangulation &tria, + const unsigned int n_subdivisions, + const Point (&corners) [dim], + const bool colorize) + { + // Equalise number of subdivisions in each dim-direction, heir + // validity will be checked later + unsigned int (n_subdivisions_) [dim]; + for (unsigned int i=0; i + void + subdivided_parallelepiped (Triangulation &tria, + const unsigned int ( n_subdivisions) [dim], + const Point (&corners) [dim], + const bool colorize) + { + // Zero n_subdivisions is the origin only, which makes no sense + for (unsigned int i=0; i0, ExcInvalidRepetitions(n_subdivisions[i])); - for (size_type y=0; y (delta) [dim]; + + for (unsigned int i=0; i > points; + + switch (dim) + { + case 1: + for (unsigned int x=0; x<=n_subdivisions[0]; ++x) + points.push_back (Point (x*delta[0])); + break; + + case 2: + for (unsigned int y=0; y<=n_subdivisions[1]; ++y) + for (unsigned int x=0; x<=n_subdivisions[0]; ++x) + points.push_back (Point (x*delta[0] + y*delta[1])); + break; + + case 3: + for (unsigned int z=0; z<=n_subdivisions[2]; ++z) + for (unsigned int y=0; y<=n_subdivisions[1]; ++y) + for (unsigned int x=0; x<=n_subdivisions[0]; ++x) + points.push_back (Point (x*delta[0] + y*delta[1] + z*delta[2])); + break; + + default: + Assert (false, ExcNotImplemented()); + } + + // Prepare cell data + unsigned int n_cells = 1; + for (unsigned int i=0; i > cells (n_cells); + + // Create fixed ordering of + switch (dim) + { + case 1: + for (unsigned int x=0; x -void -GridGenerator::subdivided_hyper_cube (Triangulation &tria, - const unsigned int repetitions, - const double left, - const double right) -{ - Assert (repetitions >= 1, ExcInvalidRepetitions(repetitions)); - Assert (left < right, - ExcMessage ("Invalid left-to-right bounds of hypercube")); - - // first generate the necessary - // points - const double delta = (right-left)/repetitions; - std::vector > points; - switch (dim) - { - case 1: - for (unsigned int x=0; x<=repetitions; ++x) - points.push_back (Point (left+x*delta)); - break; - case 2: - for (unsigned int y=0; y<=repetitions; ++y) + template + void + subdivided_hyper_cube (Triangulation &tria, + const unsigned int repetitions, + const double left, + const double right) + { + Assert (repetitions >= 1, ExcInvalidRepetitions(repetitions)); + Assert (left < right, + ExcMessage ("Invalid left-to-right bounds of hypercube")); + + // first generate the necessary + // points + const double delta = (right-left)/repetitions; + std::vector > points; + switch (dim) + { + case 1: for (unsigned int x=0; x<=repetitions; ++x) - points.push_back (Point (left+x*delta, - left+y*delta)); - break; + points.push_back (Point (left+x*delta)); + break; - case 3: - for (unsigned int z=0; z<=repetitions; ++z) + case 2: for (unsigned int y=0; y<=repetitions; ++y) for (unsigned int x=0; x<=repetitions; ++x) points.push_back (Point (left+x*delta, - left+y*delta, - left+z*delta)); - break; - - default: - Assert (false, ExcNotImplemented()); - } - - // next create the cells - // Prepare cell data - std::vector > cells; - // Define these as abbreviations - // for the step sizes below. The - // number of points in a single - // direction is repetitions+1 - const unsigned int dy = repetitions+1; - const unsigned int dz = dy*dy; - switch (dim) - { - case 1: - cells.resize (repetitions); - for (unsigned int x=0; x (left+x*delta, + left+y*delta, + left+z*delta)); + break; + + default: + Assert (false, ExcNotImplemented()); + } - case 2: - cells.resize (repetitions*repetitions); - for (unsigned int y=0; y > cells; + // Define these as abbreviations + // for the step sizes below. The + // number of points in a single + // direction is repetitions+1 + const unsigned int dy = repetitions+1; + const unsigned int dz = dy*dy; + switch (dim) + { + case 1: + cells.resize (repetitions); for (unsigned int x=0; x -void -GridGenerator::subdivided_hyper_rectangle ( - Triangulation &tria, - const std::vector &repetitions, - const Point &p_1, - const Point &p_2, - const bool colorize) -{ - // contributed by Joerg R. Weimar - // (j.weimar@jweimar.de) 2003 - 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)); - - // then check that all repetitions - // are >= 1, and calculate deltas - // convert repetitions from double - // to int by taking the ceiling. - Point delta; - - for (unsigned int i=0; i= 1, ExcInvalidRepetitions(repetitions[i])); + template + void + subdivided_hyper_rectangle ( + Triangulation &tria, + const std::vector &repetitions, + const Point &p_1, + const Point &p_2, + const bool colorize) + { + // contributed by Joerg R. Weimar + // (j.weimar@jweimar.de) 2003 + 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)); - delta[i] = (p2[i]-p1[i])/repetitions[i]; - } + // then check that all repetitions + // are >= 1, and calculate deltas + // convert repetitions from double + // to int by taking the ceiling. + Point delta; - // then generate the necessary - // points - std::vector > points; - switch (dim) - { - case 1: - for (unsigned int x=0; x<=repetitions[0]; ++x) - points.push_back (Point (p1[0]+x*delta[0])); - break; + for (unsigned int i=0; i= 1, ExcInvalidRepetitions(repetitions[i])); - case 2: - for (unsigned int y=0; y<=repetitions[1]; ++y) + delta[i] = (p2[i]-p1[i])/repetitions[i]; + } + + // then generate the necessary + // points + std::vector > points; + switch (dim) + { + case 1: for (unsigned int x=0; x<=repetitions[0]; ++x) - points.push_back (Point (p1[0]+x*delta[0], - p1[1]+y*delta[1])); - break; + points.push_back (Point (p1[0]+x*delta[0])); + break; - case 3: - for (unsigned int z=0; z<=repetitions[2]; ++z) + case 2: for (unsigned int y=0; y<=repetitions[1]; ++y) for (unsigned int x=0; x<=repetitions[0]; ++x) points.push_back (Point (p1[0]+x*delta[0], - p1[1]+y*delta[1], - p1[2]+z*delta[2])); - break; - - default: - Assert (false, ExcNotImplemented()); - } - - // next create the cells - // Prepare cell data - std::vector > cells; - switch (dim) - { - case 1: - { - cells.resize (repetitions[0]); - for (unsigned int x=0; x (p1[0]+x*delta[0], + p1[1]+y*delta[1], + p1[2]+z*delta[2])); + break; + + default: + Assert (false, ExcNotImplemented()); + } - case 2: - { - cells.resize (repetitions[1]*repetitions[0]); - for (unsigned int y=0; y > cells; + switch (dim) + { + case 1: + { + cells.resize (repetitions[0]); for (unsigned int x=0; x 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); - } -} + default: + Assert (false, ExcNotImplemented()); + } + tria.create_triangulation (points, cells, SubCellData()); + if (colorize) + { + // to colorize, run through all + // 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. + const double epsilon + = 0.01 * *std::min_element (&delta[0], &delta[0]+dim); + Assert (epsilon > 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); + } + } -template -void -GridGenerator::subdivided_hyper_rectangle( - Triangulation &tria, - const std::vector > &step_sz, - const Point &p_1, - const Point &p_2, - const bool colorize) -{ - // contributed by Joerg R. Weimar - // (j.weimar@jweimar.de) 2003 - // modified by Yaqi Wang 2006 - Assert(step_sz.size() == dim, - ExcInvalidRepetitionsDimension(dim)); - - - // First, normalize input such that - // p1 is lower in all coordinate - // directions. - - // and check the consistency of - // step sizes, i.e. that they all - // add up to the sizes specified by - // p_1 and p_2 - Point p1(p_1); - Point p2(p_2); - std::vector< std::vector > step_sizes(step_sz); - - for (unsigned int i=0; i p2(i)) - { - std::swap (p1(i), p2(i)); - std::reverse (step_sizes[i].begin(), step_sizes[i].end()); - } - double x = 0; - for (size_type j=0; j + void + subdivided_hyper_rectangle( + Triangulation &tria, + const std::vector > &step_sz, + const Point &p_1, + const Point &p_2, + const bool colorize) + { + // contributed by Joerg R. Weimar + // (j.weimar@jweimar.de) 2003 + // modified by Yaqi Wang 2006 + Assert(step_sz.size() == dim, + ExcInvalidRepetitionsDimension(dim)); + + + // First, normalize input such that + // p1 is lower in all coordinate + // directions. + + // and check the consistency of + // step sizes, i.e. that they all + // add up to the sizes specified by + // p_1 and p_2 + Point p1(p_1); + Point p2(p_2); + std::vector< std::vector > step_sizes(step_sz); + + for (unsigned int i=0; i p2(i)) + { + std::swap (p1(i), p2(i)); + std::reverse (step_sizes[i].begin(), step_sizes[i].end()); + } - // then generate the necessary - // points - std::vector > points; - switch (dim) - { - case 1: - { - double x=0; - for (size_type i=0; ; ++i) - { - points.push_back (Point (p1[0]+x)); - - // form partial sums. in - // the last run through - // avoid accessing - // non-existent values - // and exit early instead - if (i == step_sizes[0].size()) - break; + double x = 0; + for (unsigned int j=0; j (p1[0]+x, - p1[1]+y)); - if (i == step_sizes[0].size()) - break; + // then generate the necessary + // points + std::vector > points; + switch (dim) + { + case 1: + { + double x=0; + for (unsigned int i=0; ; ++i) + { + points.push_back (Point (p1[0]+x)); - x += step_sizes[0][i]; - } + // form partial sums. in + // the last run through + // avoid accessing + // non-existent values + // and exit early instead + if (i == step_sizes[0].size()) + break; - if (j == step_sizes[1].size()) - break; + x += step_sizes[0][i]; + } + break; + } - y += step_sizes[1][j]; - } - break; + case 2: + { + double y=0; + for (unsigned int j=0; ; ++j) + { + double x=0; + for (unsigned int i=0; ; ++i) + { + points.push_back (Point (p1[0]+x, + p1[1]+y)); + if (i == step_sizes[0].size()) + break; - } - case 3: - { - double z=0; - for (size_type k=0; ; ++k) - { - double y=0; - for (size_type j=0; ; ++j) - { - double x=0; - for (size_type i=0; ; ++i) - { - points.push_back (Point (p1[0]+x, - p1[1]+y, - p1[2]+z)); - if (i == step_sizes[0].size()) - break; + x += step_sizes[0][i]; + } - x += step_sizes[0][i]; - } + if (j == step_sizes[1].size()) + break; - if (j == step_sizes[1].size()) - break; + y += step_sizes[1][j]; + } + break; - y += step_sizes[1][j]; - } + } + case 3: + { + double z=0; + for (unsigned int k=0; ; ++k) + { + double y=0; + for (unsigned int j=0; ; ++j) + { + double x=0; + for (unsigned int i=0; ; ++i) + { + points.push_back (Point (p1[0]+x, + p1[1]+y, + p1[2]+z)); + if (i == step_sizes[0].size()) + break; - if (k == step_sizes[2].size()) - break; + x += step_sizes[0][i]; + } - z += step_sizes[2][k]; - } - break; - } + if (j == step_sizes[1].size()) + break; - default: - Assert (false, ExcNotImplemented()); - } + y += step_sizes[1][j]; + } - // next create the cells - // Prepare cell data - std::vector > cells; - switch (dim) - { - case 1: - { - cells.resize (step_sizes[0].size()); - for (size_type x=0; x > cells; + switch (dim) + { + case 1: + { + cells.resize (step_sizes[0].size()); + for (unsigned int x=0; x -void -GridGenerator::subdivided_hyper_rectangle ( - Triangulation<1> &tria, - const std::vector< std::vector > &spacing, - const Point<1> &p, - const Table<1,types::material_id> &material_id, - const bool colorize) -{ - // contributed by Yaqi Wang 2006 - Assert(spacing.size() == 1, - ExcInvalidRepetitionsDimension(1)); - const size_type n_cells = material_id.size(0); - Assert(spacing[0].size() == n_cells, - ExcInvalidRepetitionsDimension(1)); + template <> + void + subdivided_hyper_rectangle ( + Triangulation<1> &tria, + const std::vector< std::vector > &spacing, + const Point<1> &p, + const Table<1,types::material_id> &material_id, + const bool colorize) + { + // contributed by Yaqi Wang 2006 + Assert(spacing.size() == 1, + ExcInvalidRepetitionsDimension(1)); - double delta = std::numeric_limits::max(); - for (size_type i=0; i= 0, ExcInvalidRepetitions(-1)); - delta = std::min (delta, spacing[0][i]); - } + const unsigned int n_cells = material_id.size(0); - // generate the necessary points - std::vector > points; - double ax = p[0]; - for (size_type x=0; x<=n_cells; ++x) - { - points.push_back (Point<1> (ax)); - if (x > cells(n_val_cells); - size_type id = 0; - for (size_type x=0; x::max(); + for (unsigned int i=0; i= 0, ExcInvalidRepetitions(-1)); + delta = std::min (delta, spacing[0][i]); } - // create triangulation - SubCellData t; - GridTools::delete_unused_vertices (points, cells, t); - tria.create_triangulation (points, cells, t); + // generate the necessary points + std::vector > points; + double ax = p[0]; + for (unsigned int x=0; x<=n_cells; ++x) + { + points.push_back (Point<1> (ax)); + if (x > cells(n_val_cells); + unsigned int id = 0; + for (unsigned int x=0; x -void -GridGenerator::subdivided_hyper_rectangle ( - Triangulation<2> &tria, - const std::vector< std::vector > &spacing, - const Point<2> &p, - const Table<2,types::material_id> &material_id, - const bool colorize) -{ - // contributed by Yaqi Wang 2006 - Assert(spacing.size() == 2, - ExcInvalidRepetitionsDimension(2)); - - std::vector repetitions(2); - size_type n_cells = 1; - double delta = std::numeric_limits::max(); - for (unsigned int i=0; i<2; i++) - { - repetitions[i] = spacing[i].size(); - n_cells *= repetitions[i]; - for (size_type j=0; j= 0, ExcInvalidRepetitions(-1)); - delta = std::min (delta, spacing[i][j]); - } - Assert(material_id.size(i) == repetitions[i], - ExcInvalidRepetitionsDimension(i)); - } - // generate the necessary points - std::vector > points; - double ay = p[1]; - for (size_type y=0; y<=repetitions[1]; ++y) - { - double ax = p[0]; - for (size_type x=0; x<=repetitions[0]; ++x) - { - points.push_back (Point<2> (ax,ay)); - if (x + void + subdivided_hyper_rectangle ( + Triangulation<2> &tria, + const std::vector< std::vector > &spacing, + const Point<2> &p, + const Table<2,types::material_id> &material_id, + const bool colorize) + { + // contributed by Yaqi Wang 2006 + Assert(spacing.size() == 2, + ExcInvalidRepetitionsDimension(2)); + + std::vector repetitions(2); + unsigned int n_cells = 1; + double delta = std::numeric_limits::max(); + for (unsigned int i=0; i<2; i++) + { + repetitions[i] = spacing[i].size(); + n_cells *= repetitions[i]; + for (unsigned int j=0; j= 0, ExcInvalidRepetitions(-1)); + delta = std::min (delta, spacing[i][j]); + } + Assert(material_id.size(i) == repetitions[i], + ExcInvalidRepetitionsDimension(i)); + } - // create the cells - size_type n_val_cells = 0; - for (size_type i=0; i > cells(n_val_cells); - size_type id = 0; - for (size_type y=0; y > points; + double ay = p[1]; + for (unsigned int y=0; y<=repetitions[1]; ++y) + { + double ax = p[0]; + for (unsigned int x=0; x<=repetitions[0]; ++x) + { + points.push_back (Point<2> (ax,ay)); + if (x > cells(n_val_cells); + unsigned int id = 0; + for (unsigned int y=0; y::cell_iterator cell = tria.begin(), - endc = tria.end(); - for (; cell !=endc; ++cell) - { - Point<2> cell_center = cell->center(); - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->boundary_indicator() == 0) - { - Point<2> face_center = cell->face(f)->center(); - for (unsigned int i=0; i<2; ++i) - { - if (face_center[i]face(f)->set_boundary_indicator(i*2); - if (face_center[i]>cell_center[i]+eps) - cell->face(f)->set_boundary_indicator(i*2+1); - } - } - } - } -} + // create triangulation + SubCellData t; + GridTools::delete_unused_vertices (points, cells, t); + tria.create_triangulation (points, cells, t); -template <> -void -GridGenerator::subdivided_hyper_rectangle ( - Triangulation<3> &tria, - const std::vector< std::vector > &spacing, - const Point<3> &p, - const Table<3,types::material_id> &material_id, - const bool colorize) -{ - // contributed by Yaqi Wang 2006 - const unsigned int dim = 3; + // set boundary indicator + if (colorize) + { + double eps = 0.01 * delta; + Triangulation<2>::cell_iterator cell = tria.begin(), + endc = tria.end(); + for (; cell !=endc; ++cell) + { + Point<2> cell_center = cell->center(); + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->boundary_indicator() == 0) + { + Point<2> face_center = cell->face(f)->center(); + for (unsigned int i=0; i<2; ++i) + { + if (face_center[i]face(f)->set_boundary_indicator(i*2); + if (face_center[i]>cell_center[i]+eps) + cell->face(f)->set_boundary_indicator(i*2+1); + } + } + } + } + } - Assert(spacing.size() == dim, - ExcInvalidRepetitionsDimension(dim)); - std::vector repetitions(dim); - size_type n_cells = 1; - double delta = std::numeric_limits::max(); - for (unsigned int i=0; i= 0, ExcInvalidRepetitions(-1)); - delta = std::min (delta, spacing[i][j]); - } - Assert(material_id.size(i) == repetitions[i], - ExcInvalidRepetitionsDimension(i)); - } + template <> + void + subdivided_hyper_rectangle ( + Triangulation<3> &tria, + const std::vector< std::vector > &spacing, + const Point<3> &p, + const Table<3,types::material_id> &material_id, + const bool colorize) + { + // contributed by Yaqi Wang 2006 + const unsigned int dim = 3; - // generate the necessary points - std::vector > points; - double az = p[2]; - for (size_type z=0; z<=repetitions[2]; ++z) - { - double ay = p[1]; - for (size_type y=0; y<=repetitions[1]; ++y) - { - double ax = p[0]; - for (size_type x=0; x<=repetitions[0]; ++x) + Assert(spacing.size() == dim, + ExcInvalidRepetitionsDimension(dim)); + + std::vector repetitions(dim); + unsigned int n_cells = 1; + double delta = std::numeric_limits::max(); + for (unsigned int i=0; i= 0, ExcInvalidRepetitions(-1)); + delta = std::min (delta, spacing[i][j]); + } + Assert(material_id.size(i) == repetitions[i], + ExcInvalidRepetitionsDimension(i)); + } + + // generate the necessary points + std::vector > points; + double az = p[2]; + for (unsigned int z=0; z<=repetitions[2]; ++z) + { + double ay = p[1]; + for (unsigned int y=0; y<=repetitions[1]; ++y) + { + double ax = p[0]; + for (unsigned int x=0; x<=repetitions[0]; ++x) + { + points.push_back (Point (ax,ay,az)); + if (x > cells(n_val_cells); + unsigned int id = 0; + const unsigned int n_x = (repetitions[0]+1); + const unsigned int n_xy = (repetitions[0]+1)*(repetitions[1]+1); + for (unsigned int z=0; z (ax,ay,az)); - if (x > cells(n_val_cells); - size_type id = 0; - const size_type n_x = (repetitions[0]+1); - const size_type n_xy = (repetitions[0]+1)*(repetitions[1]+1); - for (size_type z=0; z1) + { + double eps = 0.01 * delta; + Triangulation::cell_iterator cell = tria.begin(), + endc = tria.end(); + for (; cell !=endc; ++cell) { - cells[id].vertices[0] = z*n_xy + y*n_x + x; - cells[id].vertices[1] = z*n_xy + y*n_x + x+1; - cells[id].vertices[2] = z*n_xy + (y+1)*n_x + x; - cells[id].vertices[3] = z*n_xy + (y+1)*n_x + x+1; - cells[id].vertices[4] = (z+1)*n_xy + y*n_x + x; - cells[id].vertices[5] = (z+1)*n_xy + y*n_x + x+1; - cells[id].vertices[6] = (z+1)*n_xy + (y+1)*n_x + x; - cells[id].vertices[7] = (z+1)*n_xy + (y+1)*n_x + x+1; - cells[id].material_id = material_id[x][y][z]; - id++; + Point cell_center = cell->center(); + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->boundary_indicator() == 0) + { + Point face_center = cell->face(f)->center(); + for (unsigned int i=0; iface(f)->set_boundary_indicator(i*2); + if (face_center[i]>cell_center[i]+eps) + cell->face(f)->set_boundary_indicator(i*2+1); + } + } } + } + } - // create triangulation - SubCellData t; - GridTools::delete_unused_vertices (points, cells, t); - tria.create_triangulation (points, cells, t); - // set boundary indicator - if (colorize && dim>1) - { - double eps = 0.01 * delta; - Triangulation::cell_iterator cell = tria.begin(), - endc = tria.end(); - for (; cell !=endc; ++cell) - { - Point cell_center = cell->center(); - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->boundary_indicator() == 0) - { - Point face_center = cell->face(f)->center(); - for (unsigned int i=0; iface(f)->set_boundary_indicator(i*2); - if (face_center[i]>cell_center[i]+eps) - cell->face(f)->set_boundary_indicator(i*2+1); - } - } - } - } -} + template <> + void hyper_cube_slit (Triangulation<1> &, + const double, + const double, + const bool) + { + Assert (false, ExcNotImplemented()); + } -template <> -void -GridGenerator::colorize_subdivided_hyper_rectangle ( - Triangulation<1> &tria, - const Point<1> &, - const Point<1> &, - const double) -{ - for (Triangulation<1>::cell_iterator cell = tria.begin(); - cell != tria.end(); ++cell) - if (cell->center()(0) > 0) - cell->set_material_id(1); - // boundary indicators are set to - // 0 (left) and 1 (right) by default. -} + template <> + void enclosed_hyper_cube (Triangulation<1> &, + const double, + const double, + const double, + const bool) + { + Assert (false, ExcNotImplemented()); + } -template -void -GridGenerator::colorize_subdivided_hyper_rectangle (Triangulation &tria, - const Point &p1, - const Point &p2, - const double epsilon) -{ + template <> + void hyper_L (Triangulation<1> &, + const double, + const double) + { + Assert (false, ExcNotImplemented()); + } - // 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::face_iterator face = tria.begin_face(), - endface = tria.end_face(); - for (; face!=endface; ++face) - { - if (face->boundary_indicator() == 0) - { - const Point center (face->center()); - if (std::abs(center(0)-p1[0]) < epsilon) - face->set_boundary_indicator(0); - else if (std::abs(center(0) - p2[0]) < epsilon) - face->set_boundary_indicator(1); - else if (dim > 1 && std::abs(center(1) - p1[1]) < epsilon) - face->set_boundary_indicator(2); - else if (dim > 1 && std::abs(center(1) - p2[1]) < epsilon) - face->set_boundary_indicator(3); - else if (dim > 2 && std::abs(center(2) - p1[2]) < epsilon) - face->set_boundary_indicator(4); - else if (dim > 2 && std::abs(center(2) - p2[2]) < epsilon) - face->set_boundary_indicator(5); - else - // triangulation says it - // is on the boundary, - // but we could not find - // on which boundary. - Assert (false, ExcInternalError()); - } - } - for (typename Triangulation::cell_iterator cell = tria.begin(); - cell != tria.end(); ++cell) - { - char id = 0; - for (unsigned int d=0; dcenter()(d) > 0) id += 1 << d; - cell->set_material_id(id); - } -} + template <> + void hyper_ball (Triangulation<1> &, + const Point<1> &, + const double) + { + Assert (false, ExcNotImplemented()); + } -template <> -void GridGenerator::hyper_cube_slit (Triangulation<1> &, - const double, - const double, - const bool) -{ - Assert (false, ExcNotImplemented()); -} + template <> + void cylinder (Triangulation<1> &, + const double, + const double) + { + Assert (false, ExcNotImplemented()); + } -template <> -void GridGenerator::enclosed_hyper_cube (Triangulation<1> &, - const double, - const double, - const double, - const bool) -{ - Assert (false, ExcNotImplemented()); -} + template <> + void truncated_cone (Triangulation<1> &, + const double, + const double, + const double) + { + Assert (false, ExcNotImplemented()); + } -template <> -void GridGenerator::hyper_L (Triangulation<1> &, - const double, - const double) -{ - Assert (false, ExcNotImplemented()); -} + template <> + void hyper_shell (Triangulation<1> &, + const Point<1> &, + const double, + const double, + const unsigned int , + const bool) + { + Assert (false, ExcNotImplemented()); + } -template <> -void GridGenerator::hyper_ball (Triangulation<1> &, - const Point<1> &, - const double) -{ - Assert (false, ExcNotImplemented()); -} + template <> + void cylinder_shell (Triangulation<1> &, + const double, + const double, + const double, + const unsigned int , + const unsigned int ) + { + Assert (false, ExcNotImplemented()); + } + template <> + void + half_hyper_ball (Triangulation<1> &, + const Point<1> &, + const double) + { + Assert (false, ExcNotImplemented()); + } -template <> -void GridGenerator::cylinder (Triangulation<1> &, - const double, - const double) -{ - Assert (false, ExcNotImplemented()); -} + template <> + void + half_hyper_shell (Triangulation<1> &, + const Point<1> &, + const double, + const double, + const unsigned int , + const bool) + { + Assert (false, ExcNotImplemented()); + } + + template <> + void quarter_hyper_shell (Triangulation<1> &, + const Point<1> &, + const double, + const double, + const unsigned int , + const bool) + { + Assert (false, ExcNotImplemented()); + } + + template <> + void enclosed_hyper_cube (Triangulation<2> &tria, + const double left, + const double right, + const double thickness, + const bool colorize) + { + Assert(left > vertices(16); + double coords[4]; + coords[0] = left-thickness; + coords[1] = left; + coords[2] = right; + coords[3] = right+thickness; + + unsigned int k=0; + for (unsigned int i0=0; i0<4; ++i0) + for (unsigned int i1=0; i1<4; ++i1) + vertices[k++] = Point<2>(coords[i1], coords[i0]); + + const types::material_id materials[9] = { 5, 4, 6, + 1, 0, 2, + 9, 8,10 + }; + + std::vector > cells(9); + k = 0; + for (unsigned int i0=0; i0<3; ++i0) + for (unsigned int i1=0; i1<3; ++i1) + { + cells[k].vertices[0] = i1+4*i0; + cells[k].vertices[1] = i1+4*i0+1; + cells[k].vertices[2] = i1+4*i0+4; + cells[k].vertices[3] = i1+4*i0+5; + if (colorize) + cells[k].material_id = materials[k]; + ++k; + } + tria.create_triangulation (vertices, + cells, + SubCellData()); // no boundary information + } -template <> -void GridGenerator::truncated_cone (Triangulation<1> &, - const double, - const double, - const double) -{ - Assert (false, ExcNotImplemented()); -} +// Implementation for 2D only + template <> + void + hyper_cube_slit (Triangulation<2> &tria, + const double left, + const double right, + const bool colorize) + { + const double rl2=(right+left)/2; + const Point<2> vertices[10] = { Point<2>(left, left ), + Point<2>(rl2, left ), + Point<2>(rl2, rl2 ), + Point<2>(left, rl2 ), + Point<2>(right,left ), + Point<2>(right,rl2 ), + Point<2>(rl2, right), + Point<2>(left, right), + Point<2>(right,right), + Point<2>(rl2, left ) + }; + const int cell_vertices[4][4] = { { 0,1,3,2 }, + { 9,4,2,5 }, + { 3,2,7,6 }, + { 2,5,6,8 } + }; + std::vector > cells (4, CellData<2>()); + 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[10]), + cells, + SubCellData()); // no boundary information + if (colorize) + { + Triangulation<2>::cell_iterator cell = tria.begin(); + cell->face(1)->set_boundary_indicator(1); + ++cell; + cell->face(3)->set_boundary_indicator(2); + } + } -template <> -void GridGenerator::hyper_shell (Triangulation<1> &, - const Point<1> &, - const double, - const double, - const size_type , - const bool) -{ - Assert (false, ExcNotImplemented()); -} -template <> -void GridGenerator::colorize_hyper_shell (Triangulation<1> &, - const Point<1> &, - const double, - const double) -{ - Assert (false, ExcNotImplemented()); -} + template <> + void truncated_cone (Triangulation<2> &triangulation, + const double radius_0, + const double radius_1, + const double half_length) + { + Point<2> vertices_tmp[4]; + vertices_tmp[0] = Point<2> (-half_length, -radius_0); + vertices_tmp[1] = Point<2> (half_length, -radius_1); + vertices_tmp[2] = Point<2> (-half_length, radius_0); + vertices_tmp[3] = Point<2> (half_length, radius_1); -template <> -void GridGenerator::cylinder_shell (Triangulation<1> &, - const double, - const double, - const double, - const size_type , - const size_type ) -{ - Assert (false, ExcNotImplemented()); -} + const std::vector > vertices (&vertices_tmp[0], &vertices_tmp[4]); + unsigned int cell_vertices[1][GeometryInfo<2>::vertices_per_cell]; + for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i) + cell_vertices[0][i] = i; -template <> -void -GridGenerator::half_hyper_ball (Triangulation<1> &, - const Point<1> &, - const double) -{ - Assert (false, ExcNotImplemented()); -} + std::vector > cells (1, CellData<2> ()); + for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i) + cells[0].vertices[i] = cell_vertices[0][i]; -template <> -void -GridGenerator::half_hyper_shell (Triangulation<1> &, - const Point<1> &, - const double, - const double, - const size_type , - const bool) -{ - Assert (false, ExcNotImplemented()); -} + cells[0].material_id = 0; + triangulation.create_triangulation (vertices, cells, SubCellData ()); -template <> -void GridGenerator::quarter_hyper_shell (Triangulation<1> &, - const Point<1> &, - const double, - const double, - const size_type , - const bool) -{ - Assert (false, ExcNotImplemented()); -} + Triangulation<2>::cell_iterator cell = triangulation.begin (); -template <> -void GridGenerator::enclosed_hyper_cube (Triangulation<2> &tria, - const double left, - const double right, - const double thickness, - const bool colorize) -{ - Assert(left > vertices(16); - double coords[4]; - coords[0] = left-thickness; - coords[1] = left; - coords[2] = right; - coords[3] = right+thickness; - - unsigned int k=0; - for (unsigned int i0=0; i0<4; ++i0) - for (unsigned int i1=0; i1<4; ++i1) - vertices[k++] = Point<2>(coords[i1], coords[i0]); - - const types::material_id materials[9] = { 5, 4, 6, - 1, 0, 2, - 9, 8,10 - }; - - std::vector > cells(9); - k = 0; - for (unsigned int i0=0; i0<3; ++i0) - for (unsigned int i1=0; i1<3; ++i1) - { - cells[k].vertices[0] = i1+4*i0; - cells[k].vertices[1] = i1+4*i0+1; - cells[k].vertices[2] = i1+4*i0+4; - cells[k].vertices[3] = i1+4*i0+5; - if (colorize) - cells[k].material_id = materials[k]; - ++k; - } - tria.create_triangulation (vertices, - cells, - SubCellData()); // no boundary information -} + 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 <> -void -GridGenerator::hyper_cube_slit (Triangulation<2> &tria, - const double left, - const double right, - const bool colorize) -{ - const double rl2=(right+left)/2; - const Point<2> vertices[10] = { Point<2>(left, left ), - Point<2>(rl2, left ), - Point<2>(rl2, rl2 ), - Point<2>(left, rl2 ), - Point<2>(right,left ), - Point<2>(right,rl2 ), - Point<2>(rl2, right), - Point<2>(left, right), - Point<2>(right,right), - Point<2>(rl2, left ) - }; - const int cell_vertices[4][4] = { { 0,1,3,2 }, - { 9,4,2,5 }, - { 3,2,7,6 }, - { 2,5,6,8 } - }; - std::vector > cells (4, CellData<2>()); - 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; + template <> + void + hyper_L (Triangulation<2> &tria, + const double a, + const double b) + { + const Point<2> vertices[8] = { Point<2> (a,a), + Point<2> ((a+b)/2,a), + Point<2> (b,a), + Point<2> (a,(a+b)/2), + Point<2> ((a+b)/2,(a+b)/2), + Point<2> (b,(a+b)/2), + Point<2> (a,b), + Point<2> ((a+b)/2,b) + }; + const int cell_vertices[3][4] = {{0, 1, 3, 4}, + {1, 2, 4, 5}, + {3, 4, 6, 7} }; - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[10]), - cells, - SubCellData()); // no boundary information - if (colorize) - { - Triangulation<2>::cell_iterator cell = tria.begin(); - cell->face(1)->set_boundary_indicator(1); - ++cell; - cell->face(3)->set_boundary_indicator(2); - } -} + std::vector > cells (3, CellData<2>()); + 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, + SubCellData()); // no boundary information + } -template <> -void GridGenerator::truncated_cone (Triangulation<2> &triangulation, - const double radius_0, - const double radius_1, - const double half_length) -{ - Point<2> vertices_tmp[4]; - vertices_tmp[0] = Point<2> (-half_length, -radius_0); - vertices_tmp[1] = Point<2> (half_length, -radius_1); - vertices_tmp[2] = Point<2> (-half_length, radius_0); - vertices_tmp[3] = Point<2> (half_length, radius_1); - const std::vector > vertices (&vertices_tmp[0], &vertices_tmp[4]); - unsigned int cell_vertices[1][GeometryInfo<2>::vertices_per_cell]; +// Implementation for 2D only + template <> + void + hyper_ball (Triangulation<2> &tria, + const Point<2> &p, + const double radius) + { + // equilibrate cell sizes at + // transition from the inner part + // to the radial cells + const double a = 1./(1+std::sqrt(2.0)); + const Point<2> vertices[8] = { p+Point<2>(-1,-1) *(radius/std::sqrt(2.0)), + p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)), + p+Point<2>(-1,-1) *(radius/std::sqrt(2.0)*a), + p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)*a), + p+Point<2>(-1,+1) *(radius/std::sqrt(2.0)*a), + p+Point<2>(+1,+1) *(radius/std::sqrt(2.0)*a), + p+Point<2>(-1,+1) *(radius/std::sqrt(2.0)), + p+Point<2>(+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}, + {1, 7, 3, 5}, + {6, 4, 7, 5} + }; - for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i) - cell_vertices[0][i] = i; + std::vector > cells (5, CellData<2>()); - std::vector > cells (1, CellData<2> ()); + 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; + }; - for (unsigned int i = 0; i < GeometryInfo<2>::vertices_per_cell; ++i) - cells[0].vertices[i] = cell_vertices[0][i]; + tria.create_triangulation ( + std::vector >(&vertices[0], &vertices[8]), + cells, + SubCellData()); // no boundary information + } - cells[0].material_id = 0; - triangulation.create_triangulation (vertices, cells, SubCellData ()); - Triangulation<2>::cell_iterator cell = triangulation.begin (); - cell->face (0)->set_boundary_indicator (1); - cell->face (1)->set_boundary_indicator (2); + template <> + void hyper_shell (Triangulation<2> &tria, + const Point<2> ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells, + const bool colorize) + { + 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 + // the length of each cell on the + // median (in the middle between + // the two circles) is equal to its + // radial extent (which is the + // difference between the two + // radii) + const unsigned int N = (n_cells == 0 ? + static_cast + (std::ceil((2*pi* (outer_radius + inner_radius)/2) / + (outer_radius - inner_radius))) : + n_cells); + + // set up N vertices on the + // outer and N vertices on + // the inner circle. the + // first N ones are on the + // outer one, and all are + // numbered counter-clockwise + std::vector > vertices(2*N); + for (unsigned int i=0; i( std::cos(2*pi * i/N), + std::sin(2*pi * i/N)) * outer_radius; + vertices[i+N] = vertices[i] * (inner_radius/outer_radius); - for (unsigned int i = 2; i < 4; ++i) - cell->face (i)->set_boundary_indicator (0); -} + vertices[i] += center; + vertices[i+N] += center; + }; + std::vector > cells (N, CellData<2>()); + for (unsigned int i=0; i -void -GridGenerator::hyper_L (Triangulation<2> &tria, - const double a, - const double b) -{ - const Point<2> vertices[8] = { Point<2> (a,a), - Point<2> ((a+b)/2,a), - Point<2> (b,a), - Point<2> (a,(a+b)/2), - Point<2> ((a+b)/2,(a+b)/2), - Point<2> (b,(a+b)/2), - Point<2> (a,b), - Point<2> ((a+b)/2,b) - }; - const int cell_vertices[3][4] = {{0, 1, 3, 4}, - {1, 2, 4, 5}, - {3, 4, 6, 7} - }; - - std::vector > cells (3, CellData<2>()); - - 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; - }; + cells[i].material_id = 0; + }; - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[8]), - cells, - SubCellData()); // no boundary information -} + tria.create_triangulation ( + vertices, cells, SubCellData()); + if (colorize) + colorize_hyper_shell(tria, center, inner_radius, outer_radius); + } // Implementation for 2D only -template <> -void -GridGenerator::hyper_ball (Triangulation<2> &tria, - const Point<2> &p, - const double radius) -{ - // equilibrate cell sizes at - // transition from the inner part - // to the radial cells - const double a = 1./(1+std::sqrt(2.0)); - const Point<2> vertices[8] = { p+Point<2>(-1,-1) *(radius/std::sqrt(2.0)), - p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)), - p+Point<2>(-1,-1) *(radius/std::sqrt(2.0)*a), - p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)*a), - p+Point<2>(-1,+1) *(radius/std::sqrt(2.0)*a), - p+Point<2>(+1,+1) *(radius/std::sqrt(2.0)*a), - p+Point<2>(-1,+1) *(radius/std::sqrt(2.0)), - p+Point<2>(+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}, - {1, 7, 3, 5}, - {6, 4, 7, 5} - }; - - std::vector > cells (5, CellData<2>()); - - 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; - }; + template <> + void + cylinder (Triangulation<2> &tria, + const double radius, + const double half_length) + { + Point<2> p1 (-half_length, -radius); + Point<2> p2 (half_length, radius); - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[8]), - cells, - SubCellData()); // no boundary information -} + 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; + } + } // Implementation for 2D only -template<> -void -GridGenerator::colorize_hyper_shell ( - Triangulation<2> &tria, - const Point<2> &, const double, const double) -{ - // In spite of receiving geometrical - // data, we do this only based on - // topology. - - // For the mesh based on cube, - // this is highly irregular - for (Triangulation<2>::cell_iterator cell = tria.begin(); - cell != tria.end(); ++cell) - { - Assert (cell->face(2)->at_boundary(), ExcInternalError()); - cell->face(2)->set_boundary_indicator(1); - } -} - + template <> + void cylinder_shell (Triangulation<2> &, + const double, + const double, + const double, + const unsigned int, + const unsigned int) + { + Assert (false, ExcNotImplemented()); + } -template <> -void GridGenerator::hyper_shell (Triangulation<2> &tria, - const Point<2> ¢er, - const double inner_radius, - const double outer_radius, - const size_type n_cells, - const bool colorize) -{ - 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 - // the length of each cell on the - // median (in the middle between - // the two circles) is equal to its - // radial extent (which is the - // difference between the two - // radii) - const size_type N = (n_cells == 0 ? - static_cast - (std::ceil((2*pi* (outer_radius + inner_radius)/2) / - (outer_radius - inner_radius))) : - n_cells); - - // set up N vertices on the - // outer and N vertices on - // the inner circle. the - // first N ones are on the - // outer one, and all are - // numbered counter-clockwise - std::vector > vertices(2*N); - for (size_type i=0; i( std::cos(2*pi * i/N), - std::sin(2*pi * i/N)) * outer_radius; - vertices[i+N] = vertices[i] * (inner_radius/outer_radius); - vertices[i] += center; - vertices[i+N] += center; + template <> + void + half_hyper_ball (Triangulation<2> &tria, + const Point<2> &p, + const double radius) + { + // equilibrate cell sizes at + // transition from the inner part + // to the radial cells + const double a = 1./(1+std::sqrt(2.0)); + const Point<2> vertices[8] = { p+Point<2>(0,-1) *radius, + p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)), + p+Point<2>(0,-1) *(radius/std::sqrt(2.0)*a), + p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)*a), + p+Point<2>(0,+1) *(radius/std::sqrt(2.0)*a), + p+Point<2>(+1,+1) *(radius/std::sqrt(2.0)*a), + p+Point<2>(0,+1) *radius, + p+Point<2>(+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 (N, CellData<2>()); - - for (size_type i=0; i > cells (4, CellData<2>()); - tria.create_triangulation ( - vertices, cells, SubCellData()); + 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; + }; - if (colorize) - colorize_hyper_shell(tria, center, inner_radius, outer_radius); -} + tria.create_triangulation ( + std::vector >(&vertices[0], &vertices[8]), + cells, + SubCellData()); // no boundary information + Triangulation<2>::cell_iterator cell = tria.begin(); + Triangulation<2>::cell_iterator end = tria.end(); -// Implementation for 2D only -template <> -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); + while (cell != end) + { + for (unsigned int i=0; i::faces_per_cell; ++i) + { + if (cell->face(i)->boundary_indicator() == numbers::internal_face_boundary_id) + continue; - 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; - } -} + // If x is zero, then this is part of the plane + if (cell->face(i)->center()(0) < p(0)+1.e-5) + cell->face(i)->set_boundary_indicator(1); + } + ++cell; + } + } // Implementation for 2D only -template <> -void GridGenerator::cylinder_shell (Triangulation<2> &, - const double, - const double, - const double, - const size_type, - const size_type) -{ - Assert (false, ExcNotImplemented()); -} + template <> + void + half_hyper_shell (Triangulation<2> &tria, + const Point<2> ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells, + const bool colorize) + { + 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 + // the length of each cell on the + // median (in the middle between + // the two circles) is equal to its + // radial extent (which is the + // difference between the two + // radii) + const unsigned int N = (n_cells == 0 ? + static_cast + (std::ceil((pi* (outer_radius + inner_radius)/2) / + (outer_radius - inner_radius))) : + n_cells); + + // set up N+1 vertices on the + // outer and N+1 vertices on + // the inner circle. the + // first N+1 ones are on the + // outer one, and all are + // numbered counter-clockwise + std::vector > vertices(2*(N+1)); + for (unsigned int i=0; i<=N; ++i) + { + // enforce that the x-coordinates + // of the first and last point of + // each half-circle are exactly + // zero (contrary to what we may + // compute using the imprecise + // value of pi) + vertices[i] = Point<2>( ( (i==0) || (i==N) ? + 0 : + std::cos(pi * i/N - pi/2) ), + std::sin(pi * i/N - pi/2)) * outer_radius; + vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius); + + vertices[i] += center; + vertices[i+N+1] += center; + }; -template <> -void -GridGenerator::half_hyper_ball (Triangulation<2> &tria, - const Point<2> &p, - const double radius) -{ - // equilibrate cell sizes at - // transition from the inner part - // to the radial cells - const double a = 1./(1+std::sqrt(2.0)); - const Point<2> vertices[8] = { p+Point<2>(0,-1) *radius, - p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)), - p+Point<2>(0,-1) *(radius/std::sqrt(2.0)*a), - p+Point<2>(+1,-1) *(radius/std::sqrt(2.0)*a), - p+Point<2>(0,+1) *(radius/std::sqrt(2.0)*a), - p+Point<2>(+1,+1) *(radius/std::sqrt(2.0)*a), - p+Point<2>(0,+1) *radius, - p+Point<2>(+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<2>()); - - 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; - }; + std::vector > cells (N, CellData<2>()); - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[8]), - cells, - SubCellData()); // no boundary information + for (unsigned int i=0; i::cell_iterator cell = tria.begin(); - Triangulation<2>::cell_iterator end = tria.end(); + cells[i].material_id = 0; + }; + tria.create_triangulation (vertices, cells, SubCellData()); - while (cell != end) - { - for (unsigned int i=0; i::faces_per_cell; ++i) - { - if (cell->face(i)->boundary_indicator() == numbers::internal_face_boundary_id) - continue; + if (colorize) + { + Triangulation<2>::cell_iterator cell = tria.begin(); + for (; cell!=tria.end(); ++cell) + { + cell->face(2)->set_boundary_indicator(1); + } + tria.begin()->face(0)->set_boundary_indicator(3); - // If x is zero, then this is part of the plane - if (cell->face(i)->center()(0) < p(0)+1.e-5) - cell->face(i)->set_boundary_indicator(1); - } - ++cell; - } -} + tria.last()->face(1)->set_boundary_indicator(2); + } + } + template <> + void quarter_hyper_shell (Triangulation<2> &tria, + const Point<2> ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n_cells, + const bool colorize) + { + 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 + // the length of each cell on the + // median (in the middle between + // the two circles) is equal to its + // radial extent (which is the + // difference between the two + // radii) + const unsigned int N = (n_cells == 0 ? + static_cast + (std::ceil((pi* (outer_radius + inner_radius)/4) / + (outer_radius - inner_radius))) : + n_cells); + + // set up N+1 vertices on the + // outer and N+1 vertices on + // the inner circle. the + // first N+1 ones are on the + // outer one, and all are + // numbered counter-clockwise + std::vector > vertices(2*(N+1)); + for (unsigned int i=0; i<=N; ++i) + { + // enforce that the x-coordinates + // of the last point is exactly + // zero (contrary to what we may + // compute using the imprecise + // value of pi) + vertices[i] = Point<2>( ( (i==N) ? + 0 : + std::cos(pi * i/N/2) ), + std::sin(pi * i/N/2)) * outer_radius; + vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius); + + vertices[i] += center; + vertices[i+N+1] += center; + }; -// Implementation for 2D only -template <> -void -GridGenerator::half_hyper_shell (Triangulation<2> &tria, - const Point<2> ¢er, - const double inner_radius, - const double outer_radius, - const size_type n_cells, - const bool colorize) -{ - 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 - // the length of each cell on the - // median (in the middle between - // the two circles) is equal to its - // radial extent (which is the - // difference between the two - // radii) - const size_type N = (n_cells == 0 ? - static_cast - (std::ceil((pi* (outer_radius + inner_radius)/2) / - (outer_radius - inner_radius))) : - n_cells); - - // set up N+1 vertices on the - // outer and N+1 vertices on - // the inner circle. the - // first N+1 ones are on the - // outer one, and all are - // numbered counter-clockwise - std::vector > vertices(2*(N+1)); - for (size_type i=0; i<=N; ++i) - { - // enforce that the x-coordinates - // of the first and last point of - // each half-circle are exactly - // zero (contrary to what we may - // compute using the imprecise - // value of pi) - vertices[i] = Point<2>( ( (i==0) || (i==N) ? - 0 : - std::cos(pi * i/N - pi/2) ), - std::sin(pi * i/N - pi/2)) * outer_radius; - vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius); - - vertices[i] += center; - vertices[i+N+1] += center; - }; + std::vector > cells (N, CellData<2>()); - std::vector > cells (N, CellData<2>()); + for (unsigned int i=0; i::cell_iterator cell = tria.begin(); + for (; cell!=tria.end(); ++cell) + { + cell->face(2)->set_boundary_indicator(1); + } + tria.begin()->face(0)->set_boundary_indicator(3); - if (colorize) - { - Triangulation<2>::cell_iterator cell = tria.begin(); - for (; cell!=tria.end(); ++cell) - { - cell->face(2)->set_boundary_indicator(1); - } - tria.begin()->face(0)->set_boundary_indicator(3); + tria.last()->face(1)->set_boundary_indicator(2); + } + } - tria.last()->face(1)->set_boundary_indicator(2); - } -} -template <> -void GridGenerator::quarter_hyper_shell (Triangulation<2> &tria, - const Point<2> ¢er, - const double inner_radius, - const double outer_radius, - const size_type n_cells, - const bool colorize) -{ - 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 - // the length of each cell on the - // median (in the middle between - // the two circles) is equal to its - // radial extent (which is the - // difference between the two - // radii) - const size_type N = (n_cells == 0 ? - static_cast - (std::ceil((pi* (outer_radius + inner_radius)/4) / - (outer_radius - inner_radius))) : - n_cells); - - // set up N+1 vertices on the - // outer and N+1 vertices on - // the inner circle. the - // first N+1 ones are on the - // outer one, and all are - // numbered counter-clockwise - std::vector > vertices(2*(N+1)); - for (size_type i=0; i<=N; ++i) - { - // enforce that the x-coordinates - // of the last point is exactly - // zero (contrary to what we may - // compute using the imprecise - // value of pi) - vertices[i] = Point<2>( ( (i==N) ? - 0 : - std::cos(pi * i/N/2) ), - std::sin(pi * i/N/2)) * outer_radius; - vertices[i+N+1] = vertices[i] * (inner_radius/outer_radius); - - vertices[i] += center; - vertices[i+N+1] += center; +// Implementation for 3D only + template <> + void hyper_cube_slit (Triangulation<3> &tria, + const double left, + const double right, + const bool colorize) + { + const double rl2=(right+left)/2; + const double len = (right-left)/2.; + + const Point<3> vertices[20] = + { + Point<3>(left, left , -len/2.), + Point<3>(rl2, left , -len/2.), + Point<3>(rl2, rl2 , -len/2.), + Point<3>(left, rl2 , -len/2.), + Point<3>(right,left , -len/2.), + Point<3>(right,rl2 , -len/2.), + Point<3>(rl2, right, -len/2.), + Point<3>(left, right, -len/2.), + Point<3>(right,right, -len/2.), + Point<3>(rl2, left , -len/2.), + Point<3>(left, left , len/2.), + Point<3>(rl2, left , len/2.), + Point<3>(rl2, rl2 , len/2.), + Point<3>(left, rl2 , len/2.), + Point<3>(right,left , len/2.), + Point<3>(right,rl2 , len/2.), + Point<3>(rl2, right, len/2.), + Point<3>(left, right, len/2.), + Point<3>(right,right, len/2.), + Point<3>(rl2, left , len/2.) }; - - - std::vector > cells (N, CellData<2>()); - - for (size_type i=0; i > cells (4, CellData<3>()); + for (unsigned int i=0; i<4; ++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[20]), + cells, + SubCellData()); // no boundary information - tria.create_triangulation (vertices, cells, SubCellData()); - - if (colorize) - { - Triangulation<2>::cell_iterator cell = tria.begin(); - for (; cell!=tria.end(); ++cell) - { - cell->face(2)->set_boundary_indicator(1); - } - tria.begin()->face(0)->set_boundary_indicator(3); - - tria.last()->face(1)->set_boundary_indicator(2); - } -} + if (colorize) + { + Assert(false, ExcNotImplemented()); + Triangulation<3>::cell_iterator cell = tria.begin(); + cell->face(1)->set_boundary_indicator(1); + ++cell; + cell->face(3)->set_boundary_indicator(2); + } + } // Implementation for 3D only -template <> -void GridGenerator::hyper_cube_slit (Triangulation<3> &tria, - const double left, - const double right, - const bool colorize) -{ - const double rl2=(right+left)/2; - const double len = (right-left)/2.; - - const Point<3> vertices[20] = + template <> + void enclosed_hyper_cube (Triangulation<3> &tria, + const double left, + const double right, + const double thickness, + const bool colorize) { - Point<3>(left, left , -len/2.), - Point<3>(rl2, left , -len/2.), - Point<3>(rl2, rl2 , -len/2.), - Point<3>(left, rl2 , -len/2.), - Point<3>(right,left , -len/2.), - Point<3>(right,rl2 , -len/2.), - Point<3>(rl2, right, -len/2.), - Point<3>(left, right, -len/2.), - Point<3>(right,right, -len/2.), - Point<3>(rl2, left , -len/2.), - Point<3>(left, left , len/2.), - Point<3>(rl2, left , len/2.), - Point<3>(rl2, rl2 , len/2.), - Point<3>(left, rl2 , len/2.), - Point<3>(right,left , len/2.), - Point<3>(right,rl2 , len/2.), - Point<3>(rl2, right, len/2.), - Point<3>(left, right, len/2.), - Point<3>(right,right, len/2.), - Point<3>(rl2, left , len/2.) - }; - const int cell_vertices[4][8] = { { 0,1,3,2, 10, 11, 13, 12 }, - { 9,4,2,5, 19,14, 12, 15 }, - { 3,2,7,6,13,12,17,16 }, - { 2,5,6,8,12,15,16,18 } - }; - std::vector > cells (4, CellData<3>()); - for (unsigned int i=0; i<4; ++i) - { - for (unsigned int j=0; j<8; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; + Assert(left > vertices(64); + double coords[4]; + coords[0] = left-thickness; + coords[1] = left; + coords[2] = right; + coords[3] = right+thickness; + + unsigned int k=0; + for (unsigned int z=0; z<4; ++z) + for (unsigned int y=0; y<4; ++y) + for (unsigned int x=0; x<4; ++x) + vertices[k++] = Point<3>(coords[x], coords[y], coords[z]); + + const types::material_id materials[27] = + { + 21,20,22, + 17,16,18, + 25,24,26, + 5 , 4, 6, + 1 , 0, 2, + 9 , 8,10, + 37,36,38, + 33,32,34, + 41,40,42 }; - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[20]), - cells, - SubCellData()); // no boundary information - if (colorize) - { - Assert(false, ExcNotImplemented()); - Triangulation<3>::cell_iterator cell = tria.begin(); - cell->face(1)->set_boundary_indicator(1); - ++cell; - cell->face(3)->set_boundary_indicator(2); - } -} + std::vector > cells(27); + k = 0; + for (unsigned int z=0; z<3; ++z) + for (unsigned int y=0; y<3; ++y) + for (unsigned int x=0; x<3; ++x) + { + cells[k].vertices[0] = x+4*y+16*z; + cells[k].vertices[1] = x+4*y+16*z+1; + cells[k].vertices[2] = x+4*y+16*z+4; + cells[k].vertices[3] = x+4*y+16*z+5; + cells[k].vertices[4] = x+4*y+16*z+16; + cells[k].vertices[5] = x+4*y+16*z+17; + cells[k].vertices[6] = x+4*y+16*z+20; + cells[k].vertices[7] = x+4*y+16*z+21; + if (colorize) + cells[k].material_id = materials[k]; + ++k; + } + tria.create_triangulation ( + vertices, + cells, + SubCellData()); // no boundary information + } -// Implementation for 3D only -template <> -void GridGenerator::enclosed_hyper_cube (Triangulation<3> &tria, - const double left, - const double right, - const double thickness, - const bool colorize) -{ - Assert(left > vertices(64); - double coords[4]; - coords[0] = left-thickness; - coords[1] = left; - coords[2] = right; - coords[3] = right+thickness; - - unsigned int k=0; - for (unsigned int z=0; z<4; ++z) - for (unsigned int y=0; y<4; ++y) - for (unsigned int x=0; x<4; ++x) - vertices[k++] = Point<3>(coords[x], coords[y], coords[z]); - - const types::material_id materials[27] = + template <> + void truncated_cone (Triangulation<3> &triangulation, + const double radius_0, + const double radius_1, + const double half_length) { - 21,20,22, - 17,16,18, - 25,24,26, - 5 , 4, 6, - 1 , 0, 2, - 9 , 8,10, - 37,36,38, - 33,32,34, - 41,40,42 - }; - - std::vector > cells(27); - k = 0; - for (unsigned int z=0; z<3; ++z) - for (unsigned int y=0; y<3; ++y) - for (unsigned int x=0; x<3; ++x) - { - cells[k].vertices[0] = x+4*y+16*z; - cells[k].vertices[1] = x+4*y+16*z+1; - cells[k].vertices[2] = x+4*y+16*z+4; - cells[k].vertices[3] = x+4*y+16*z+5; - cells[k].vertices[4] = x+4*y+16*z+16; - cells[k].vertices[5] = x+4*y+16*z+17; - cells[k].vertices[6] = x+4*y+16*z+20; - cells[k].vertices[7] = x+4*y+16*z+21; - if (colorize) - cells[k].material_id = materials[k]; - ++k; - } - tria.create_triangulation ( - vertices, - cells, - SubCellData()); // no boundary information -} - - - -template <> -void GridGenerator::truncated_cone (Triangulation<3> &triangulation, - const double radius_0, - const double radius_1, - const double half_length) -{ - // Determine number of cells and vertices - const size_type - n_cells = static_cast(std::floor (half_length / - std::max (radius_0, - radius_1) + - 0.5)); - const size_type n_vertices = 4 * (n_cells + 1); - std::vector > vertices_tmp(n_vertices); - - vertices_tmp[0] = Point<3> (-half_length, 0, -radius_0); - vertices_tmp[1] = Point<3> (-half_length, radius_0, 0); - vertices_tmp[2] = Point<3> (-half_length, -radius_0, 0); - vertices_tmp[3] = Point<3> (-half_length, 0, radius_0); - - const double dx = 2 * half_length / n_cells; - - for (size_type i = 0; i < n_cells; ++i) - { - vertices_tmp[4 * (i + 1)] - = vertices_tmp[4 * i] + - Point<3> (dx, 0, 0.5 * (radius_0 - radius_1) * dx / half_length); - vertices_tmp[4 * i + 5] - = vertices_tmp[4 * i + 1] + - Point<3> (dx, 0.5 * (radius_1 - radius_0) * dx / half_length, 0); - vertices_tmp[4 * i + 6] - = vertices_tmp[4 * i + 2] + - Point<3> (dx, 0.5 * (radius_0 - radius_1) * dx / half_length, 0); - vertices_tmp[4 * i + 7] - = vertices_tmp[4 * i + 3] + - Point<3> (dx, 0, 0.5 * (radius_1 - radius_0) * dx / half_length); - } + // Determine number of cells and vertices + const unsigned int + n_cells = static_cast(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<3> (-half_length, 0, -radius_0); + vertices_tmp[1] = Point<3> (-half_length, radius_0, 0); + vertices_tmp[2] = Point<3> (-half_length, -radius_0, 0); + vertices_tmp[3] = Point<3> (-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<3> (dx, 0, 0.5 * (radius_0 - radius_1) * dx / half_length); + vertices_tmp[4 * i + 5] + = vertices_tmp[4 * i + 1] + + Point<3> (dx, 0.5 * (radius_1 - radius_0) * dx / half_length, 0); + vertices_tmp[4 * i + 6] + = vertices_tmp[4 * i + 2] + + Point<3> (dx, 0.5 * (radius_0 - radius_1) * dx / half_length, 0); + vertices_tmp[4 * i + 7] + = vertices_tmp[4 * i + 3] + + Point<3> (dx, 0, 0.5 * (radius_1 - radius_0) * dx / half_length); + } - const std::vector > vertices (vertices_tmp.begin(), - vertices_tmp.end()); - Table<2,unsigned int> cell_vertices(n_cells,GeometryInfo<3>::vertices_per_cell); + const std::vector > vertices (vertices_tmp.begin(), + vertices_tmp.end()); + Table<2,unsigned int> cell_vertices(n_cells,GeometryInfo<3>::vertices_per_cell); - for (size_type i = 0; i < n_cells; ++i) - for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j) - cell_vertices[i][j] = 4 * i + j; + for (unsigned int i = 0; i < n_cells; ++i) + for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j) + cell_vertices[i][j] = 4 * i + j; - std::vector > cells (n_cells, CellData<3> ()); + std::vector > cells (n_cells, CellData<3> ()); - for (size_type i = 0; i < n_cells; ++i) - { - for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; + for (unsigned int i = 0; i < n_cells; ++i) + { + for (unsigned int j = 0; j < GeometryInfo<3>::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; - } + cells[i].material_id = 0; + } - triangulation.create_triangulation (vertices, cells, SubCellData ()); + triangulation.create_triangulation (vertices, cells, SubCellData ()); - for (Triangulation<3>::cell_iterator cell = triangulation.begin (); - cell != triangulation.end (); ++cell) - { - if (cell->vertex (0) (0) == -half_length) - { - cell->face (4)->set_boundary_indicator (1); + for (Triangulation<3>::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); - } + 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); + 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 = 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); - } -} + for (unsigned int i = 0; i < 4; ++i) + cell->face (i)->set_boundary_indicator (0); + } + } // Implementation for 3D only -template <> -void -GridGenerator::hyper_L (Triangulation<3> &tria, - const double a, - const double b) -{ - // we slice out the top back right - // part of the cube - const Point<3> vertices[26] - = + template <> + void + hyper_L (Triangulation<3> &tria, + const double a, + const double b) { - // front face of the big cube - Point<3> (a, a,a), - Point<3> ((a+b)/2,a,a), - Point<3> (b, a,a), - Point<3> (a, a,(a+b)/2), - Point<3> ((a+b)/2,a,(a+b)/2), - Point<3> (b, a,(a+b)/2), - Point<3> (a, a,b), - Point<3> ((a+b)/2,a,b), - Point<3> (b, a,b), - // middle face of the big cube - Point<3> (a, (a+b)/2,a), - Point<3> ((a+b)/2,(a+b)/2,a), - Point<3> (b, (a+b)/2,a), - Point<3> (a, (a+b)/2,(a+b)/2), - Point<3> ((a+b)/2,(a+b)/2,(a+b)/2), - Point<3> (b, (a+b)/2,(a+b)/2), - Point<3> (a, (a+b)/2,b), - Point<3> ((a+b)/2,(a+b)/2,b), - Point<3> (b, (a+b)/2,b), - // back face of the big cube - // last (top right) point is missing - Point<3> (a, b,a), - Point<3> ((a+b)/2,b,a), - Point<3> (b, b,a), - Point<3> (a, b,(a+b)/2), - Point<3> ((a+b)/2,b,(a+b)/2), - Point<3> (b, b,(a+b)/2), - Point<3> (a, b,b), - Point<3> ((a+b)/2,b,b) - }; - const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13}, - {1, 2, 10, 11, 4, 5, 13, 14}, - {3, 4, 12, 13, 6, 7, 15, 16}, - {4, 5, 13, 14, 7, 8, 16, 17}, - {9, 10, 18, 19, 12, 13, 21, 22}, - {10, 11, 19, 20, 13, 14, 22, 23}, - {12, 13, 21, 22, 15, 16, 24, 25} - }; - - std::vector > cells (7, CellData<3>()); - - for (unsigned int i=0; i<7; ++i) - { - for (unsigned int j=0; j<8; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; + // we slice out the top back right + // part of the cube + const Point<3> vertices[26] + = + { + // front face of the big cube + Point<3> (a, a,a), + Point<3> ((a+b)/2,a,a), + Point<3> (b, a,a), + Point<3> (a, a,(a+b)/2), + Point<3> ((a+b)/2,a,(a+b)/2), + Point<3> (b, a,(a+b)/2), + Point<3> (a, a,b), + Point<3> ((a+b)/2,a,b), + Point<3> (b, a,b), + // middle face of the big cube + Point<3> (a, (a+b)/2,a), + Point<3> ((a+b)/2,(a+b)/2,a), + Point<3> (b, (a+b)/2,a), + Point<3> (a, (a+b)/2,(a+b)/2), + Point<3> ((a+b)/2,(a+b)/2,(a+b)/2), + Point<3> (b, (a+b)/2,(a+b)/2), + Point<3> (a, (a+b)/2,b), + Point<3> ((a+b)/2,(a+b)/2,b), + Point<3> (b, (a+b)/2,b), + // back face of the big cube + // last (top right) point is missing + Point<3> (a, b,a), + Point<3> ((a+b)/2,b,a), + Point<3> (b, b,a), + Point<3> (a, b,(a+b)/2), + Point<3> ((a+b)/2,b,(a+b)/2), + Point<3> (b, b,(a+b)/2), + Point<3> (a, b,b), + Point<3> ((a+b)/2,b,b) + }; + const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13}, + {1, 2, 10, 11, 4, 5, 13, 14}, + {3, 4, 12, 13, 6, 7, 15, 16}, + {4, 5, 13, 14, 7, 8, 16, 17}, + {9, 10, 18, 19, 12, 13, 21, 22}, + {10, 11, 19, 20, 13, 14, 22, 23}, + {12, 13, 21, 22, 15, 16, 24, 25} }; - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[26]), - cells, - SubCellData()); // no boundary information -} + std::vector > cells (7, CellData<3>()); + + for (unsigned int i=0; i<7; ++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[26]), + cells, + SubCellData()); // no boundary information + } // Implementation for 3D only -template <> -void -GridGenerator::hyper_ball (Triangulation<3> &tria, - const Point<3> &p, - const double radius) -{ - const double a = 1./(1+std::sqrt(3.0)); // equilibrate cell sizes at transition - // from the inner part to the radial - // cells - const unsigned int n_vertices = 16; - const Point<3> vertices[n_vertices] - = + template <> + void + hyper_ball (Triangulation<3> &tria, + const Point<3> &p, + const double radius) { - // first the vertices of the inner - // cell - p+Point<3>(-1,-1,-1) *(radius/std::sqrt(3.0)*a), - p+Point<3>(+1,-1,-1) *(radius/std::sqrt(3.0)*a), - p+Point<3>(+1,-1,+1) *(radius/std::sqrt(3.0)*a), - p+Point<3>(-1,-1,+1) *(radius/std::sqrt(3.0)*a), - p+Point<3>(-1,+1,-1) *(radius/std::sqrt(3.0)*a), - p+Point<3>(+1,+1,-1) *(radius/std::sqrt(3.0)*a), - p+Point<3>(+1,+1,+1) *(radius/std::sqrt(3.0)*a), - p+Point<3>(-1,+1,+1) *(radius/std::sqrt(3.0)*a), - // now the eight vertices at - // the outer sphere - p+Point<3>(-1,-1,-1) *(radius/std::sqrt(3.0)), - p+Point<3>(+1,-1,-1) *(radius/std::sqrt(3.0)), - p+Point<3>(+1,-1,+1) *(radius/std::sqrt(3.0)), - p+Point<3>(-1,-1,+1) *(radius/std::sqrt(3.0)), - p+Point<3>(-1,+1,-1) *(radius/std::sqrt(3.0)), - p+Point<3>(+1,+1,-1) *(radius/std::sqrt(3.0)), - p+Point<3>(+1,+1,+1) *(radius/std::sqrt(3.0)), - p+Point<3>(-1,+1,+1) *(radius/std::sqrt(3.0)), - }; - - // one needs to draw the seven cubes to - // understand what's going on here - const unsigned int n_cells = 7; - const int cell_vertices[n_cells][8] = {{0, 1, 4, 5, 3, 2, 7, 6}, // center - {8, 9, 12, 13, 0, 1, 4, 5}, // bottom - {9, 13, 1, 5, 10, 14, 2, 6}, // right - {11, 10, 3, 2, 15, 14, 7, 6}, // top - {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<3>()); - - for (unsigned int i=0; i::vertices_per_cell; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; + const double a = 1./(1+std::sqrt(3.0)); // equilibrate cell sizes at transition + // from the inner part to the radial + // cells + const unsigned int n_vertices = 16; + const Point<3> vertices[n_vertices] + = + { + // first the vertices of the inner + // cell + p+Point<3>(-1,-1,-1) *(radius/std::sqrt(3.0)*a), + p+Point<3>(+1,-1,-1) *(radius/std::sqrt(3.0)*a), + p+Point<3>(+1,-1,+1) *(radius/std::sqrt(3.0)*a), + p+Point<3>(-1,-1,+1) *(radius/std::sqrt(3.0)*a), + p+Point<3>(-1,+1,-1) *(radius/std::sqrt(3.0)*a), + p+Point<3>(+1,+1,-1) *(radius/std::sqrt(3.0)*a), + p+Point<3>(+1,+1,+1) *(radius/std::sqrt(3.0)*a), + p+Point<3>(-1,+1,+1) *(radius/std::sqrt(3.0)*a), + // now the eight vertices at + // the outer sphere + p+Point<3>(-1,-1,-1) *(radius/std::sqrt(3.0)), + p+Point<3>(+1,-1,-1) *(radius/std::sqrt(3.0)), + p+Point<3>(+1,-1,+1) *(radius/std::sqrt(3.0)), + p+Point<3>(-1,-1,+1) *(radius/std::sqrt(3.0)), + p+Point<3>(-1,+1,-1) *(radius/std::sqrt(3.0)), + p+Point<3>(+1,+1,-1) *(radius/std::sqrt(3.0)), + p+Point<3>(+1,+1,+1) *(radius/std::sqrt(3.0)), + p+Point<3>(-1,+1,+1) *(radius/std::sqrt(3.0)), }; - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[n_vertices]), - cells, - SubCellData()); // no boundary information -} + // one needs to draw the seven cubes to + // understand what's going on here + const unsigned int n_cells = 7; + const int cell_vertices[n_cells][8] = {{0, 1, 4, 5, 3, 2, 7, 6}, // center + {8, 9, 12, 13, 0, 1, 4, 5}, // bottom + {9, 13, 1, 5, 10, 14, 2, 6}, // right + {11, 10, 3, 2, 15, 14, 7, 6}, // top + {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<3>()); + + for (unsigned int i=0; i::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + }; + + tria.create_triangulation ( + std::vector >(&vertices[0], &vertices[n_vertices]), + cells, + SubCellData()); // no boundary information + } // Implementation for 3D only -template <> -void -GridGenerator::cylinder (Triangulation<3> &tria, - const double radius, - const double half_length) -{ - // Copy the base from hyper_ball<3> - // and transform it to yz - const double d = radius/std::sqrt(2.0); - const double a = d/(1+std::sqrt(2.0)); - Point<3> vertices[24] = + template <> + void + cylinder (Triangulation<3> &tria, + const double radius, + const double half_length) { - Point<3>(-d, -half_length,-d), - Point<3>( d, -half_length,-d), - Point<3>(-a, -half_length,-a), - Point<3>( a, -half_length,-a), - Point<3>(-a, -half_length, a), - Point<3>( a, -half_length, a), - Point<3>(-d, -half_length, d), - Point<3>( d, -half_length, d), - Point<3>(-d, 0,-d), - Point<3>( d, 0,-d), - Point<3>(-a, 0,-a), - Point<3>( a, 0,-a), - Point<3>(-a, 0, a), - Point<3>( a, 0, a), - Point<3>(-d, 0, d), - Point<3>( d, 0, d), - Point<3>(-d, half_length,-d), - Point<3>( d, half_length,-d), - Point<3>(-a, half_length,-a), - Point<3>( a, half_length,-a), - Point<3>(-a, half_length, a), - Point<3>( a, half_length, a), - Point<3>(-d, half_length, d), - Point<3>( d, half_length, d), - }; - // Turn cylinder such that y->x - for (unsigned int i=0; i<24; ++i) - { - const double h = vertices[i](1); - vertices[i](1) = -vertices[i](0); - vertices[i](0) = h; - } + // Copy the base from hyper_ball<3> + // and transform it to yz + const double d = radius/std::sqrt(2.0); + const double a = d/(1+std::sqrt(2.0)); + Point<3> vertices[24] = + { + Point<3>(-d, -half_length,-d), + Point<3>( d, -half_length,-d), + Point<3>(-a, -half_length,-a), + Point<3>( a, -half_length,-a), + Point<3>(-a, -half_length, a), + Point<3>( a, -half_length, a), + Point<3>(-d, -half_length, d), + Point<3>( d, -half_length, d), + Point<3>(-d, 0,-d), + Point<3>( d, 0,-d), + Point<3>(-a, 0,-a), + Point<3>( a, 0,-a), + Point<3>(-a, 0, a), + Point<3>( a, 0, a), + Point<3>(-d, 0, d), + Point<3>( d, 0, d), + Point<3>(-d, half_length,-d), + Point<3>( d, half_length,-d), + Point<3>(-a, half_length,-a), + Point<3>( a, half_length,-a), + Point<3>(-a, half_length, a), + Point<3>( a, half_length, a), + Point<3>(-d, half_length, d), + Point<3>( d, half_length, d), + }; + // Turn cylinder such that y->x + for (unsigned int i=0; i<24; ++i) + { + const double h = vertices[i](1); + 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}, - {2, 3, 10, 11, 4, 5, 12, 13}, - {1, 7, 9, 15, 3, 5, 11, 13}, - {6, 4, 14, 12, 7, 5, 15, 13} - }; - 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) + int cell_vertices[10][8] = { - for (unsigned int j=0; j<8; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; + {0, 1, 8, 9, 2, 3, 10, 11}, + {0, 2, 8, 10, 6, 4, 14, 12}, + {2, 3, 10, 11, 4, 5, 12, 13}, + {1, 7, 9, 15, 3, 5, 11, 13}, + {6, 4, 14, 12, 7, 5, 15, 13} }; + 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; - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[24]), - cells, - SubCellData()); // no boundary information - - // set boundary indicators for the - // faces at the ends to 1 and 2, - // respectively. note that we also - // have to deal with those lines - // that are purely in the interior - // of the ends. we determine whether - // an edge is purely in the - // interior if one of its vertices - // is at coordinates '+-a' as set - // above - Triangulation<3>::cell_iterator cell = tria.begin(); - Triangulation<3>::cell_iterator end = tria.end(); - - for (; cell != end; ++cell) - for (unsigned int i=0; i::faces_per_cell; ++i) - if (cell->at_boundary(i)) - { - if (cell->face(i)->center()(0) > half_length-1.e-5) - { - cell->face(i)->set_boundary_indicator(2); - - for (unsigned int e=0; e::lines_per_face; ++e) - if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a)) - cell->face(i)->line(e)->set_boundary_indicator(2); - } - else if (cell->face(i)->center()(0) < -half_length+1.e-5) - { - cell->face(i)->set_boundary_indicator(1); + std::vector > cells (10, CellData<3>()); - for (unsigned int e=0; e::lines_per_face; ++e) - if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) || - (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a)) - cell->face(i)->line(e)->set_boundary_indicator(1); - } - } -} + 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 + + // set boundary indicators for the + // faces at the ends to 1 and 2, + // respectively. note that we also + // have to deal with those lines + // that are purely in the interior + // of the ends. we determine whether + // an edge is purely in the + // interior if one of its vertices + // is at coordinates '+-a' as set + // above + Triangulation<3>::cell_iterator cell = tria.begin(); + Triangulation<3>::cell_iterator end = tria.end(); + + for (; cell != end; ++cell) + for (unsigned int i=0; i::faces_per_cell; ++i) + if (cell->at_boundary(i)) + { + if (cell->face(i)->center()(0) > half_length-1.e-5) + { + cell->face(i)->set_boundary_indicator(2); + + for (unsigned int e=0; e::lines_per_face; ++e) + if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a)) + cell->face(i)->line(e)->set_boundary_indicator(2); + } + else if (cell->face(i)->center()(0) < -half_length+1.e-5) + { + cell->face(i)->set_boundary_indicator(1); + + for (unsigned int e=0; e::lines_per_face; ++e) + if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) || + (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a)) + cell->face(i)->line(e)->set_boundary_indicator(1); + } + } + } // Implementation for 3D only -template <> -void -GridGenerator::half_hyper_ball (Triangulation<3> &tria, - const Point<3> ¢er, - const double radius) -{ - // These are for the two lower squares - const double d = radius/std::sqrt(2.0); - const double a = d/(1+std::sqrt(2.0)); - // These are for the two upper square - const double b = a/2.0; - const double c = d/2.0; - // And so are these - const double hb = radius*std::sqrt(3.0)/4.0; - const double hc = radius*std::sqrt(3.0)/2.0; - - Point<3> vertices[16] = + template <> + void + half_hyper_ball (Triangulation<3> &tria, + const Point<3> ¢er, + const double radius) { - center+Point<3>( 0, d, -d), - center+Point<3>( 0, -d, -d), - center+Point<3>( 0, a, -a), - center+Point<3>( 0, -a, -a), - center+Point<3>( 0, a, a), - center+Point<3>( 0, -a, a), - center+Point<3>( 0, d, d), - center+Point<3>( 0, -d, d), - - center+Point<3>(hc, c, -c), - center+Point<3>(hc, -c, -c), - center+Point<3>(hb, b, -b), - center+Point<3>(hb, -b, -b), - center+Point<3>(hb, b, b), - center+Point<3>(hb, -b, b), - center+Point<3>(hc, c, c), - center+Point<3>(hc, -c, c), - }; - - int cell_vertices[6][8] = - { - {0, 1, 8, 9, 2, 3, 10, 11}, - {0, 2, 8, 10, 6, 4, 14, 12}, - {2, 3, 10, 11, 4, 5, 12, 13}, - {1, 7, 9, 15, 3, 5, 11, 13}, - {6, 4, 14, 12, 7, 5, 15, 13}, - {8, 10, 9, 11, 14, 12, 15, 13} - }; - - std::vector > cells (6, CellData<3>()); + // These are for the two lower squares + const double d = radius/std::sqrt(2.0); + const double a = d/(1+std::sqrt(2.0)); + // These are for the two upper square + const double b = a/2.0; + const double c = d/2.0; + // And so are these + const double hb = radius*std::sqrt(3.0)/4.0; + const double hc = radius*std::sqrt(3.0)/2.0; + + Point<3> vertices[16] = + { + center+Point<3>( 0, d, -d), + center+Point<3>( 0, -d, -d), + center+Point<3>( 0, a, -a), + center+Point<3>( 0, -a, -a), + center+Point<3>( 0, a, a), + center+Point<3>( 0, -a, a), + center+Point<3>( 0, d, d), + center+Point<3>( 0, -d, d), + + center+Point<3>(hc, c, -c), + center+Point<3>(hc, -c, -c), + center+Point<3>(hb, b, -b), + center+Point<3>(hb, -b, -b), + center+Point<3>(hb, b, b), + center+Point<3>(hb, -b, b), + center+Point<3>(hc, c, c), + center+Point<3>(hc, -c, c), + }; - for (unsigned int i=0; i<6; ++i) + int cell_vertices[6][8] = { - for (unsigned int j=0; j<8; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; + {0, 1, 8, 9, 2, 3, 10, 11}, + {0, 2, 8, 10, 6, 4, 14, 12}, + {2, 3, 10, 11, 4, 5, 12, 13}, + {1, 7, 9, 15, 3, 5, 11, 13}, + {6, 4, 14, 12, 7, 5, 15, 13}, + {8, 10, 9, 11, 14, 12, 15, 13} }; - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[16]), - cells, - SubCellData()); // no boundary information + std::vector > cells (6, CellData<3>()); - Triangulation<3>::cell_iterator cell = tria.begin(); - Triangulation<3>::cell_iterator end = tria.end(); - - // go over all faces. for the ones on the flat face, set boundary - // indicator for face and edges to one; the rest will remain at - // zero but we have to pay attention to those edges that are - // at the perimeter of the flat face since they should not be - // set to one - while (cell != end) - { - for (unsigned int i=0; i::faces_per_cell; ++i) - { - if (!cell->at_boundary(i)) - continue; + for (unsigned int i=0; i<6; ++i) + { + for (unsigned int j=0; j<8; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + }; - // If the center is on the plane x=0, this is a planar element. set - // its boundary indicator. also set the boundary indicators of the - // bounding faces unless both vertices are on the perimeter - if (cell->face(i)->center()(0) < center(0)+1.e-5*radius) - { - cell->face(i)->set_boundary_indicator(1); - for (unsigned int j=0; j::lines_per_face; ++j) - { - const Point<3> vertices[2] - = { cell->face(i)->line(j)->vertex(0), - cell->face(i)->line(j)->vertex(1) - }; - if ((std::fabs(vertices[0].distance(center)-radius) > - 1e-5*radius) - || - (std::fabs(vertices[1].distance(center)-radius) > - 1e-5*radius)) - cell->face(i)->line(j)->set_boundary_indicator(1); - } - } - } - ++cell; - } -} + tria.create_triangulation ( + std::vector >(&vertices[0], &vertices[16]), + cells, + SubCellData()); // no boundary information -// Implementation for 3D only -template<> -void -GridGenerator:: -colorize_hyper_shell (Triangulation<3> &tria, - const Point<3> &, - const double, - const double) -{ - // the following uses a good amount - // of knowledge about the - // orientation of cells. this is - // probably not good style... - if (tria.n_cells() == 6) - { - Triangulation<3>::cell_iterator cell = tria.begin(); + Triangulation<3>::cell_iterator cell = tria.begin(); + Triangulation<3>::cell_iterator end = tria.end(); - cell->face(4)->set_boundary_indicator(1); - Assert (cell->face(4)->at_boundary(), ExcInternalError()); - - (++cell)->face(2)->set_boundary_indicator(1); - Assert (cell->face(2)->at_boundary(), ExcInternalError()); - - (++cell)->face(2)->set_boundary_indicator(1); - Assert (cell->face(2)->at_boundary(), ExcInternalError()); - - (++cell)->face(0)->set_boundary_indicator(1); - Assert (cell->face(0)->at_boundary(), ExcInternalError()); - - (++cell)->face(2)->set_boundary_indicator(1); - Assert (cell->face(2)->at_boundary(), ExcInternalError()); - - (++cell)->face(0)->set_boundary_indicator(1); - Assert (cell->face(0)->at_boundary(), ExcInternalError()); - } - else if (tria.n_cells() == 12) - { - // again use some internal - // knowledge - for (Triangulation<3>::cell_iterator cell = tria.begin(); - cell != tria.end(); ++cell) - { - Assert (cell->face(5)->at_boundary(), ExcInternalError()); - cell->face(5)->set_boundary_indicator(1); - } - } - else if (tria.n_cells() == 96) - { - // the 96-cell hypershell is - // based on a once refined - // 12-cell mesh. consequently, - // since the outer faces all - // are face_no==5 above, so - // they are here (unless they - // are in the interior). Use - // this to assign boundary - // indicators, but also make - // sure that we encounter - // exactly 48 such faces - unsigned int count = 0; - for (Triangulation<3>::cell_iterator cell = tria.begin(); - cell != tria.end(); ++cell) - if (cell->face(5)->at_boundary()) + // go over all faces. for the ones on the flat face, set boundary + // indicator for face and edges to one; the rest will remain at + // zero but we have to pay attention to those edges that are + // at the perimeter of the flat face since they should not be + // set to one + while (cell != end) + { + for (unsigned int i=0; i::faces_per_cell; ++i) { - cell->face(5)->set_boundary_indicator(1); - ++count; - } - Assert (count == 48, ExcInternalError()); - } - else - Assert (false, ExcNotImplemented()); -} - + if (!cell->at_boundary(i)) + continue; + // If the center is on the plane x=0, this is a planar element. set + // its boundary indicator. also set the boundary indicators of the + // bounding faces unless both vertices are on the perimeter + if (cell->face(i)->center()(0) < center(0)+1.e-5*radius) + { + cell->face(i)->set_boundary_indicator(1); + for (unsigned int j=0; j::lines_per_face; ++j) + { + const Point<3> vertices[2] + = { cell->face(i)->line(j)->vertex(0), + cell->face(i)->line(j)->vertex(1) + }; + if ((std::fabs(vertices[0].distance(center)-radius) > + 1e-5*radius) + || + (std::fabs(vertices[1].distance(center)-radius) > + 1e-5*radius)) + cell->face(i)->line(j)->set_boundary_indicator(1); + } + } + } + ++cell; + } + } -template <> -void -GridGenerator::hyper_shell (Triangulation<3> &tria, - const Point<3> &p, - const double inner_radius, - const double outer_radius, - const size_type n, - const bool colorize) -{ - Assert ((inner_radius > 0) && (inner_radius < outer_radius), - ExcInvalidRadii ()); - const double irad = inner_radius/std::sqrt(3.0); - const double orad = outer_radius/std::sqrt(3.0); - std::vector > vertices; - std::vector > cells; + template <> + void + hyper_shell (Triangulation<3> &tria, + const Point<3> &p, + const double inner_radius, + const double outer_radius, + const unsigned int n, + const bool colorize) + { + Assert ((inner_radius > 0) && (inner_radius < outer_radius), + ExcInvalidRadii ()); - // Start with the shell bounded by - // two nested cubes - if (n == 6) - { - for (unsigned int i=0; i<8; ++i) - vertices.push_back(p+hexahedron[i]*irad); - for (unsigned int i=0; i<8; ++i) - vertices.push_back(p+hexahedron[i]*orad); + const double irad = inner_radius/std::sqrt(3.0); + const double orad = outer_radius/std::sqrt(3.0); + std::vector > vertices; + std::vector > cells; - const unsigned int n_cells = 6; - const int cell_vertices[n_cells][8] = + // Start with the shell bounded by + // two nested cubes + if (n == 6) { - {8, 9, 10, 11, 0, 1, 2, 3}, // bottom - {9, 11, 1, 3, 13, 15, 5, 7}, // right - {12, 13, 4, 5, 14, 15, 6, 7}, // top - {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 + for (unsigned int i=0; i<8; ++i) + vertices.push_back(p+hexahedron[i]*irad); + for (unsigned int i=0; i<8; ++i) + vertices.push_back(p+hexahedron[i]*orad); - cells.resize(n_cells, CellData<3>()); - - for (unsigned int i=0; i::vertices_per_cell; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; - } + {8, 9, 10, 11, 0, 1, 2, 3}, // bottom + {9, 11, 1, 3, 13, 15, 5, 7}, // right + {12, 13, 4, 5, 14, 15, 6, 7}, // top + {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 - tria.create_triangulation (vertices, cells, SubCellData()); - } - // A more regular subdivision can - // be obtained by two nested - // rhombic dodecahedra - else if (n == 12) - { - for (unsigned int i=0; i<8; ++i) - vertices.push_back(p+hexahedron[i]*irad); - for (unsigned int i=0; i<6; ++i) - vertices.push_back(p+octahedron[i]*inner_radius); - for (unsigned int i=0; i<8; ++i) - vertices.push_back(p+hexahedron[i]*orad); - for (unsigned int i=0; i<6; ++i) - vertices.push_back(p+octahedron[i]*outer_radius); - - const unsigned int n_cells = 12; - const unsigned int rhombi[n_cells][4] = - { - { 10, 4, 0, 8}, - { 4, 13, 8, 6}, - { 10, 5, 4, 13}, - { 1, 9, 10, 5}, - { 9, 7, 5, 13}, - { 7, 11, 13, 6}, - { 9, 3, 7, 11}, - { 1, 12, 9, 3}, - { 12, 2, 3, 11}, - { 2, 8, 11, 6}, - { 12, 0, 2, 8}, - { 1, 10, 12, 0} - }; + cells.resize(n_cells, CellData<3>()); - cells.resize(n_cells, CellData<3>()); + for (unsigned int i=0; i::vertices_per_cell; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; + } - for (unsigned int i=0; i boundary (p); - Triangulation<3> tmp; - GridGenerator::hyper_shell (tmp, p, inner_radius, outer_radius, 12); - tmp.set_boundary(0, boundary); - tmp.set_boundary(1, boundary); - tmp.refine_global (1); - - // let's determine the distance at - // which the interior nodes should be - // from the center. let's say we - // measure distances in multiples of - // outer_radius and call - // r=inner_radius. - // - // then note - // that we now have 48 faces on the - // inner and 48 on the outer sphere, - // each with an area of approximately - // 4*pi/48*r^2 and 4*pi/48, for - // a face edge length of approximately - // sqrt(pi/12)*r and sqrt(pi/12) - // - // let's say we put the interior nodes - // at a distance rho, then a measure of - // deformation for the inner cells - // would be - // di=max(sqrt(pi/12)*r/(rho-r), - // (rho-r)/sqrt(pi/12)/r) - // and for the outer cells - // do=max(sqrt(pi/12)/(1-rho), - // (1-rho)/sqrt(pi/12)) - // - // we now seek a rho so that the - // deformation of cells on the inside - // and outside is equal. there are in - // principle four possibilities for one - // of the branches of do== one of the - // branches of di, though not all of - // them satisfy do==di, of - // course. however, we are not - // interested in cases where the inner - // cell is long and skinny and the - // outer one tall -- yes, they have the - // same aspect ratio, but in different - // space directions. - // - // so it only boils down to the - // following two possibilities: the - // first branch of each max(.,.) - // functions are equal, or the second - // one are. on the other hand, since - // they two branches are reciprocals of - // each other, if one pair of branches - // is equal, so is the other - // - // this yields the following equation - // for rho: - // sqrt(pi/12)*r/(rho-r) - // == sqrt(pi/12)/(1-rho) - // with solution rho=2r/(1+r) - const double r = inner_radius / outer_radius; - const double rho = 2*r/(1+r); - - // then this is the distance of the - // interior nodes from the center: - const double middle_radius = rho * outer_radius; - - // mark vertices we've already moved or - // that we want to ignore: we don't - // want to move vertices at the inner - // or outer boundaries - std::vector vertex_already_treated (tmp.n_vertices(), false); - for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active(); - cell != tmp.end(); ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->at_boundary(f)) - for (unsigned int v=0; v::vertices_per_face; ++v) - vertex_already_treated[cell->face(f)->vertex_index(v)] = true; - - // now move the remaining vertices - for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active(); - cell != tmp.end(); ++cell) - for (unsigned int v=0; v::vertices_per_cell; ++v) - if (vertex_already_treated[cell->vertex_index(v)] == false) - { - // this is a new interior - // vertex. mesh refinement may - // have placed it at a number - // of places in radial - // direction and oftentimes not - // in a particularly good - // one. move it to halfway - // between inner and outer - // sphere - const Point<3> old_distance = cell->vertex(v) - p; - const double old_radius = cell->vertex(v).distance(p); - cell->vertex(v) = p + old_distance * (middle_radius / old_radius); - - vertex_already_treated[cell->vertex_index(v)] = true; - } + cells.resize(n_cells, CellData<3>()); - // now copy the resulting level 1 cells - // into the new triangulation, - cells.resize(tmp.n_active_cells(), CellData<3>()); + for (unsigned int i=0; i::active_cell_iterator cell = tmp.begin_active(); - cell != tmp.end(); ++cell, ++index) - { + tria.create_triangulation (vertices, cells, SubCellData()); + } + else if (n == 96) + { + // create a triangulation based on the + // 12-cell one where we refine the mesh + // once and then re-arrange all + // interior nodes so that the mesh is + // the least distorted + HyperShellBoundary<3> boundary (p); + Triangulation<3> tmp; + hyper_shell (tmp, p, inner_radius, outer_radius, 12); + tmp.set_boundary(0, boundary); + tmp.set_boundary(1, boundary); + tmp.refine_global (1); + + // let's determine the distance at + // which the interior nodes should be + // from the center. let's say we + // measure distances in multiples of + // outer_radius and call + // r=inner_radius. + // + // then note + // that we now have 48 faces on the + // inner and 48 on the outer sphere, + // each with an area of approximately + // 4*pi/48*r^2 and 4*pi/48, for + // a face edge length of approximately + // sqrt(pi/12)*r and sqrt(pi/12) + // + // let's say we put the interior nodes + // at a distance rho, then a measure of + // deformation for the inner cells + // would be + // di=max(sqrt(pi/12)*r/(rho-r), + // (rho-r)/sqrt(pi/12)/r) + // and for the outer cells + // do=max(sqrt(pi/12)/(1-rho), + // (1-rho)/sqrt(pi/12)) + // + // we now seek a rho so that the + // deformation of cells on the inside + // and outside is equal. there are in + // principle four possibilities for one + // of the branches of do== one of the + // branches of di, though not all of + // them satisfy do==di, of + // course. however, we are not + // interested in cases where the inner + // cell is long and skinny and the + // outer one tall -- yes, they have the + // same aspect ratio, but in different + // space directions. + // + // so it only boils down to the + // following two possibilities: the + // first branch of each max(.,.) + // functions are equal, or the second + // one are. on the other hand, since + // they two branches are reciprocals of + // each other, if one pair of branches + // is equal, so is the other + // + // this yields the following equation + // for rho: + // sqrt(pi/12)*r/(rho-r) + // == sqrt(pi/12)/(1-rho) + // with solution rho=2r/(1+r) + const double r = inner_radius / outer_radius; + const double rho = 2*r/(1+r); + + // then this is the distance of the + // interior nodes from the center: + const double middle_radius = rho * outer_radius; + + // mark vertices we've already moved or + // that we want to ignore: we don't + // want to move vertices at the inner + // or outer boundaries + std::vector vertex_already_treated (tmp.n_vertices(), false); + for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active(); + cell != tmp.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->at_boundary(f)) + for (unsigned int v=0; v::vertices_per_face; ++v) + vertex_already_treated[cell->face(f)->vertex_index(v)] = true; + + // now move the remaining vertices + for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active(); + cell != tmp.end(); ++cell) for (unsigned int v=0; v::vertices_per_cell; ++v) - cells[index].vertices[v] = cell->vertex_index(v); - cells[index].material_id = 0; - } - - tria.create_triangulation (tmp.get_vertices(), cells, SubCellData()); - } - else - { - Assert(false, ExcMessage ("Invalid number of coarse mesh cells.")); - } + if (vertex_already_treated[cell->vertex_index(v)] == false) + { + // this is a new interior + // vertex. mesh refinement may + // have placed it at a number + // of places in radial + // direction and oftentimes not + // in a particularly good + // one. move it to halfway + // between inner and outer + // sphere + const Point<3> old_distance = cell->vertex(v) - p; + const double old_radius = cell->vertex(v).distance(p); + cell->vertex(v) = p + old_distance * (middle_radius / old_radius); + + vertex_already_treated[cell->vertex_index(v)] = true; + } - if (colorize) - colorize_hyper_shell(tria, p, inner_radius, outer_radius); -} + // now copy the resulting level 1 cells + // into the new triangulation, + cells.resize(tmp.n_active_cells(), CellData<3>()); + unsigned int index = 0; + for (Triangulation<3>::active_cell_iterator cell = tmp.begin_active(); + cell != tmp.end(); ++cell, ++index) + { + for (unsigned int v=0; v::vertices_per_cell; ++v) + cells[index].vertices[v] = cell->vertex_index(v); + cells[index].material_id = 0; + } + tria.create_triangulation (tmp.get_vertices(), cells, SubCellData()); + } + else + { + Assert(false, ExcMessage ("Invalid number of coarse mesh cells.")); + } + if (colorize) + colorize_hyper_shell(tria, p, inner_radius, outer_radius); + } -// Implementation for 3D only -template <> -void -GridGenerator::half_hyper_shell (Triangulation<3> &tria, - const Point<3> ¢er, - const double inner_radius, - const double outer_radius, - const size_type n, - const bool colorize) -{ - Assert ((inner_radius > 0) && (inner_radius < outer_radius), - ExcInvalidRadii ()); - Assert(colorize == false, ExcNotImplemented()); - if (n <= 5) - { - // These are for the two lower squares - const double d = outer_radius/std::sqrt(2.0); - const double a = inner_radius/std::sqrt(2.0); - // These are for the two upper square - const double b = a/2.0; - const double c = d/2.0; - // And so are these - const double hb = inner_radius*std::sqrt(3.0)/2.0; - const double hc = outer_radius*std::sqrt(3.0)/2.0; - - Point<3> vertices[16] = - { - center+Point<3>( 0, d, -d), - center+Point<3>( 0, -d, -d), - center+Point<3>( 0, a, -a), - center+Point<3>( 0, -a, -a), - center+Point<3>( 0, a, a), - center+Point<3>( 0, -a, a), - center+Point<3>( 0, d, d), - center+Point<3>( 0, -d, d), - - center+Point<3>(hc, c, -c), - center+Point<3>(hc, -c, -c), - center+Point<3>(hb, b, -b), - center+Point<3>(hb, -b, -b), - center+Point<3>(hb, b, b), - center+Point<3>(hb, -b, b), - center+Point<3>(hc, c, c), - center+Point<3>(hc, -c, c), - }; - int cell_vertices[5][8] = - { - {0, 1, 8, 9, 2, 3, 10, 11}, - {0, 2, 8, 10, 6, 4, 14, 12}, - {1, 7, 9, 15, 3, 5, 11, 13}, - {6, 4, 14, 12, 7, 5, 15, 13}, - {8, 10, 9, 11, 14, 12, 15, 13} - }; - std::vector > cells (5, CellData<3>()); +// Implementation for 3D only + template <> + void + half_hyper_shell (Triangulation<3> &tria, + const Point<3> ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n, + const bool colorize) + { + Assert ((inner_radius > 0) && (inner_radius < outer_radius), + ExcInvalidRadii ()); + Assert(colorize == false, ExcNotImplemented()); - for (unsigned int i=0; i<5; ++i) + if (n <= 5) + { + // These are for the two lower squares + const double d = outer_radius/std::sqrt(2.0); + const double a = inner_radius/std::sqrt(2.0); + // These are for the two upper square + const double b = a/2.0; + const double c = d/2.0; + // And so are these + const double hb = inner_radius*std::sqrt(3.0)/2.0; + const double hc = outer_radius*std::sqrt(3.0)/2.0; + + Point<3> vertices[16] = { - for (unsigned int j=0; j<8; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; + center+Point<3>( 0, d, -d), + center+Point<3>( 0, -d, -d), + center+Point<3>( 0, a, -a), + center+Point<3>( 0, -a, -a), + center+Point<3>( 0, a, a), + center+Point<3>( 0, -a, a), + center+Point<3>( 0, d, d), + center+Point<3>( 0, -d, d), + + center+Point<3>(hc, c, -c), + center+Point<3>(hc, -c, -c), + center+Point<3>(hb, b, -b), + center+Point<3>(hb, -b, -b), + center+Point<3>(hb, b, b), + center+Point<3>(hb, -b, b), + center+Point<3>(hc, c, c), + center+Point<3>(hc, -c, c), }; - tria.create_triangulation ( - std::vector >(&vertices[0], &vertices[16]), - cells, - SubCellData()); // no boundary information - } - else - { - Assert(false, ExcIndexRange(n, 0, 5)); - } + int cell_vertices[5][8] = + { + {0, 1, 8, 9, 2, 3, 10, 11}, + {0, 2, 8, 10, 6, 4, 14, 12}, + {1, 7, 9, 15, 3, 5, 11, 13}, + {6, 4, 14, 12, 7, 5, 15, 13}, + {8, 10, 9, 11, 14, 12, 15, 13} + }; -} + std::vector > cells (5, CellData<3>()); -// Implementation for 3D only -template <> -void -GridGenerator::colorize_quarter_hyper_shell(Triangulation<3> &tria, - const Point<3> ¢er, - const double inner_radius, - const double outer_radius) -{ + for (unsigned int i=0; i<5; ++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[16]), + cells, + SubCellData()); // no boundary information + } + else + { + Assert(false, ExcIndexRange(n, 0, 5)); + } - if (tria.n_cells() != 3) - AssertThrow (false, ExcNotImplemented()); + } - double middle = (outer_radius-inner_radius)/2e0 + inner_radius; - double eps = 1e-3*middle; - Triangulation<3>::cell_iterator cell = tria.begin(); - for (; cell!=tria.end(); ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) +// Implementation for 3D only + template <> + void quarter_hyper_shell (Triangulation<3> &tria, + const Point<3> ¢er, + const double inner_radius, + const double outer_radius, + const unsigned int n, + const bool colorize) + { + Assert ((inner_radius > 0) && (inner_radius < outer_radius), + ExcInvalidRadii ()); + if (n == 0 || n == 3) { - if (!cell->face(f)->at_boundary()) - continue; + const double a = inner_radius*std::sqrt(2.0)/2e0; + const double b = outer_radius*std::sqrt(2.0)/2e0; + const double c = a*std::sqrt(3.0)/2e0; + const double d = b*std::sqrt(3.0)/2e0; + const double e = outer_radius/2e0; + const double h = inner_radius/2e0; + + std::vector > vertices; + + vertices.push_back (center+Point<3>( 0, inner_radius, 0)); //0 + vertices.push_back (center+Point<3>( a, a, 0)); //1 + vertices.push_back (center+Point<3>( b, b, 0)); //2 + vertices.push_back (center+Point<3>( 0, outer_radius, 0)); //3 + vertices.push_back (center+Point<3>( 0, a , a)); //4 + vertices.push_back (center+Point<3>( c, c , h)); //5 + vertices.push_back (center+Point<3>( d, d , e)); //6 + vertices.push_back (center+Point<3>( 0, b , b)); //7 + vertices.push_back (center+Point<3>( inner_radius, 0 , 0)); //8 + vertices.push_back (center+Point<3>( outer_radius, 0 , 0)); //9 + vertices.push_back (center+Point<3>( a, 0 , a)); //10 + vertices.push_back (center+Point<3>( b, 0 , b)); //11 + vertices.push_back (center+Point<3>( 0, 0 , inner_radius)); //12 + vertices.push_back (center+Point<3>( 0, 0 , outer_radius)); //13 + + const int cell_vertices[3][8] = + { + {0, 1, 3, 2, 4, 5, 7, 6}, + {1, 8, 2, 9, 5, 10, 6, 11}, + {4, 5, 7, 6, 12, 10, 13, 11}, + }; + std::vector > cells(3); - double radius = cell->face(f)->center().norm() - center.norm(); - if (std::fabs(cell->face(f)->center()(0)) < eps ) // x = 0 set boundary 2 - { - cell->face(f)->set_boundary_indicator(2); - for (unsigned int j=0; j::lines_per_face; ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) - cell->face(f)->line(j)->set_boundary_indicator(2); - } - else if (std::fabs(cell->face(f)->center()(1)) < eps) // y = 0 set boundary 3 - { - cell->face(f)->set_boundary_indicator(3); - for (unsigned int j=0; j::lines_per_face; ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) - cell->face(f)->line(j)->set_boundary_indicator(3); - } - else if (std::fabs(cell->face(f)->center()(2)) < eps ) // z = 0 set boundary 4 - { - cell->face(f)->set_boundary_indicator(4); - for (unsigned int j=0; j::lines_per_face; ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) > eps) - cell->face(f)->line(j)->set_boundary_indicator(4); - } - else if (radius < middle) // inner radius set boundary 0 + for (unsigned int i=0; i<3; ++i) { - cell->face(f)->set_boundary_indicator(0); - for (unsigned int j=0; j::lines_per_face; ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps) - cell->face(f)->line(j)->set_boundary_indicator(0); + for (unsigned int j=0; j<8; ++j) + cells[i].vertices[j] = cell_vertices[i][j]; + cells[i].material_id = 0; } - else if (radius > middle) // outer radius set boundary 1 - { - cell->face(f)->set_boundary_indicator(1); - for (unsigned int j=0; j::lines_per_face; ++j) - if (cell->face(f)->line(j)->at_boundary()) - if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() - cell->face(f)->line(j)->vertex(1).norm()) < eps) - cell->face(f)->line(j)->set_boundary_indicator(1); - } - else - AssertThrow (false, ExcInternalError()); + + tria.create_triangulation ( vertices, cells, SubCellData()); // no boundary information + } + else + { + AssertThrow(false, ExcNotImplemented()); } -} + if (colorize) + colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius); + } // Implementation for 3D only -template <> -void GridGenerator::quarter_hyper_shell (Triangulation<3> &tria, - const Point<3> ¢er, - const double inner_radius, - const double outer_radius, - const size_type n, - const bool colorize) -{ - Assert ((inner_radius > 0) && (inner_radius < outer_radius), - ExcInvalidRadii ()); - if (n == 0 || n == 3) - { - const double a = inner_radius*std::sqrt(2.0)/2e0; - const double b = outer_radius*std::sqrt(2.0)/2e0; - const double c = a*std::sqrt(3.0)/2e0; - const double d = b*std::sqrt(3.0)/2e0; - const double e = outer_radius/2e0; - const double h = inner_radius/2e0; - - std::vector > vertices; - - vertices.push_back (center+Point<3>( 0, inner_radius, 0)); //0 - vertices.push_back (center+Point<3>( a, a, 0)); //1 - vertices.push_back (center+Point<3>( b, b, 0)); //2 - vertices.push_back (center+Point<3>( 0, outer_radius, 0)); //3 - vertices.push_back (center+Point<3>( 0, a , a)); //4 - vertices.push_back (center+Point<3>( c, c , h)); //5 - vertices.push_back (center+Point<3>( d, d , e)); //6 - vertices.push_back (center+Point<3>( 0, b , b)); //7 - vertices.push_back (center+Point<3>( inner_radius, 0 , 0)); //8 - vertices.push_back (center+Point<3>( outer_radius, 0 , 0)); //9 - vertices.push_back (center+Point<3>( a, 0 , a)); //10 - vertices.push_back (center+Point<3>( b, 0 , b)); //11 - vertices.push_back (center+Point<3>( 0, 0 , inner_radius)); //12 - vertices.push_back (center+Point<3>( 0, 0 , outer_radius)); //13 - - const int cell_vertices[3][8] = - { - {0, 1, 3, 2, 4, 5, 7, 6}, - {1, 8, 2, 9, 5, 10, 6, 11}, - {4, 5, 7, 6, 12, 10, 13, 11}, + template <> + void cylinder_shell (Triangulation<3> &tria, + const double length, + const double inner_radius, + const double outer_radius, + const unsigned int n_radial_cells, + const unsigned int n_axial_cells) + { + 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 + // the length of each cell on the + // median (in the middle between + // the two circles) is equal to its + // radial extent (which is the + // difference between the two + // radii) + const unsigned int N_r = (n_radial_cells == 0 ? + static_cast + (std::ceil((2*pi* (outer_radius + inner_radius)/2) / + (outer_radius - inner_radius))) : + n_radial_cells); + const unsigned int N_z = (n_axial_cells == 0 ? + static_cast + (std::ceil (length / + (2*pi*(outer_radius + inner_radius)/2/N_r))) : + n_axial_cells); + + // set up N vertices on the + // outer and N vertices on + // the inner circle. the + // first N ones are on the + // outer one, and all are + // numbered counter-clockwise + std::vector > vertices_2d(2*N_r); + for (unsigned int i=0; i( std::cos(2*pi * i/N_r), + std::sin(2*pi * i/N_r)) * outer_radius; + vertices_2d[i+N_r] = vertices_2d[i] * (inner_radius/outer_radius); }; - std::vector > cells(3); - for (unsigned int i=0; i<3; ++i) + std::vector > vertices_3d; + vertices_3d.reserve (2*N_r*(N_z+1)); + for (unsigned int j=0; j<=N_z; ++j) + for (unsigned int i=0; i<2*N_r; ++i) { - for (unsigned int j=0; j<8; ++j) - cells[i].vertices[j] = cell_vertices[i][j]; - cells[i].material_id = 0; + const Point<3> v (vertices_2d[i][0], + vertices_2d[i][1], + j*length/N_z); + vertices_3d.push_back (v); } - tria.create_triangulation ( vertices, cells, SubCellData()); // no boundary information - } - else - { - AssertThrow(false, ExcNotImplemented()); - } - - if (colorize) - colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius); -} + std::vector > cells (N_r*N_z, CellData<3>()); + for (unsigned int j=0; j -void GridGenerator::cylinder_shell (Triangulation<3> &tria, - const double length, - const double inner_radius, - const double outer_radius, - const size_type n_radial_cells, - const size_type n_axial_cells) -{ - 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 - // the length of each cell on the - // median (in the middle between - // the two circles) is equal to its - // radial extent (which is the - // difference between the two - // radii) - const size_type N_r = (n_radial_cells == 0 ? - static_cast - (std::ceil((2*pi* (outer_radius + inner_radius)/2) / - (outer_radius - inner_radius))) : - n_radial_cells); - const size_type N_z = (n_axial_cells == 0 ? - static_cast - (std::ceil (length / - (2*pi*(outer_radius + inner_radius)/2/N_r))) : - n_axial_cells); - - // set up N vertices on the - // outer and N vertices on - // the inner circle. the - // first N ones are on the - // outer one, and all are - // numbered counter-clockwise - std::vector > vertices_2d(2*N_r); - for (size_type i=0; i( std::cos(2*pi * i/N_r), - std::sin(2*pi * i/N_r)) * outer_radius; - vertices_2d[i+N_r] = vertices_2d[i] * (inner_radius/outer_radius); - }; + cells[i+j*N_r].vertices[4] = N_r+i + (j+1)*2*N_r; + 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; - std::vector > vertices_3d; - vertices_3d.reserve (2*N_r*(N_z+1)); - for (size_type j=0; j<=N_z; ++j) - for (size_type i=0; i<2*N_r; ++i) - { - const Point<3> v (vertices_2d[i][0], - vertices_2d[i][1], - j*length/N_z); - vertices_3d.push_back (v); - } + cells[i+j*N_r].material_id = 0; + } - std::vector > cells (N_r*N_z, CellData<3>()); + tria.create_triangulation ( + vertices_3d, cells, SubCellData()); + } - for (size_type j=0; j + void + merge_triangulations (const Triangulation &triangulation_1, + const Triangulation &triangulation_2, + Triangulation &result) + { + Assert (triangulation_1.n_levels() == 1, + ExcMessage ("The input triangulations must be coarse meshes.")); + Assert (triangulation_2.n_levels() == 1, + ExcMessage ("The input triangulations must be coarse meshes.")); + + // get the union of the set of vertices + std::vector > vertices = triangulation_1.get_vertices(); + vertices.insert (vertices.end(), + triangulation_2.get_vertices().begin(), + triangulation_2.get_vertices().end()); + + // now form the union of the set of cells + std::vector > cells; + cells.reserve (triangulation_1.n_cells() + triangulation_2.n_cells()); + for (typename Triangulation::cell_iterator + cell = triangulation_1.begin(); cell != triangulation_1.end(); ++cell) + { + CellData this_cell; + for (unsigned int v=0; v::vertices_per_cell; ++v) + this_cell.vertices[v] = cell->vertex_index(v); + this_cell.material_id = cell->material_id(); + cells.push_back (this_cell); } - tria.create_triangulation ( - vertices_3d, cells, SubCellData()); -} - - - -template -void -GridGenerator:: -merge_triangulations (const Triangulation &triangulation_1, - const Triangulation &triangulation_2, - Triangulation &result) -{ - Assert (triangulation_1.n_levels() == 1, - ExcMessage ("The input triangulations must be coarse meshes.")); - Assert (triangulation_2.n_levels() == 1, - ExcMessage ("The input triangulations must be coarse meshes.")); - - // get the union of the set of vertices - std::vector > vertices = triangulation_1.get_vertices(); - vertices.insert (vertices.end(), - triangulation_2.get_vertices().begin(), - triangulation_2.get_vertices().end()); - - // now form the union of the set of cells - std::vector > cells; - cells.reserve (triangulation_1.n_cells() + triangulation_2.n_cells()); - for (typename Triangulation::cell_iterator - cell = triangulation_1.begin(); cell != triangulation_1.end(); ++cell) - { - CellData this_cell; - for (unsigned int v=0; v::vertices_per_cell; ++v) - this_cell.vertices[v] = cell->vertex_index(v); - this_cell.material_id = cell->material_id(); - cells.push_back (this_cell); - } - - // now do the same for the other other mesh. note that we have to - // translate the vertex indices - for (typename Triangulation::cell_iterator - cell = triangulation_2.begin(); cell != triangulation_2.end(); ++cell) - { - CellData this_cell; - for (unsigned int v=0; v::vertices_per_cell; ++v) - this_cell.vertices[v] = cell->vertex_index(v) + triangulation_1.n_vertices(); - this_cell.material_id = cell->material_id(); - cells.push_back (this_cell); - } - - // throw out duplicated vertices from the two meshes - // and create the triangulation - SubCellData subcell_data; - std::vector considered_vertices; - GridTools::delete_duplicated_vertices (vertices, cells, subcell_data, considered_vertices); - result.clear (); - result.create_triangulation (vertices, cells, subcell_data); -} - -void -GridGenerator:: -extrude_triangulation(const Triangulation<2, 2> &input, - const size_type n_slices, - const double height, - Triangulation<3,3> &result) -{ - Assert (input.n_levels() == 1, - ExcMessage ("The input triangulations must be coarse meshes.")); - Assert(result.n_cells()==0, ExcMessage("resultin Triangulation need to be empty upon calling extrude_triangulation.")); - Assert(height>0, ExcMessage("The height in extrude_triangulation needs to be positive.")); - Assert(n_slices>=2, ExcMessage("The number of slices in extrude_triangulation needs to be at least 2.")); - - std::vector > points(n_slices*input.n_vertices()); - std::vector > cells; - cells.reserve((n_slices-1)*input.n_active_cells()); + // now do the same for the other other mesh. note that we have to + // translate the vertex indices + for (typename Triangulation::cell_iterator + cell = triangulation_2.begin(); cell != triangulation_2.end(); ++cell) + { + CellData this_cell; + for (unsigned int v=0; v::vertices_per_cell; ++v) + this_cell.vertices[v] = cell->vertex_index(v) + triangulation_1.n_vertices(); + this_cell.material_id = cell->material_id(); + cells.push_back (this_cell); + } - for (size_type slice=0; slice considered_vertices; + GridTools::delete_duplicated_vertices (vertices, cells, subcell_data, considered_vertices); + result.clear (); + result.create_triangulation (vertices, cells, subcell_data); + } + + + void + extrude_triangulation(const Triangulation<2, 2> &input, + const unsigned int n_slices, + const double height, + Triangulation<3,3> &result) + { + Assert (input.n_levels() == 1, + ExcMessage ("The input triangulations must be coarse meshes.")); + Assert(result.n_cells()==0, ExcMessage("resultin Triangulation need to be empty upon calling extrude_triangulation.")); + Assert(height>0, ExcMessage("The height in extrude_triangulation needs to be positive.")); + Assert(n_slices>=2, ExcMessage("The number of slices in extrude_triangulation needs to be at least 2.")); - { - const Point<2> &v = input.get_vertices()[i]; - points[i+slice*input.n_vertices()](0) = v(0); - points[i+slice*input.n_vertices()](1) = v(1); - points[i+slice*input.n_vertices()](2) = height * slice / (n_slices-1); - } - } + std::vector > points(n_slices*input.n_vertices()); + std::vector > cells; + cells.reserve((n_slices-1)*input.n_active_cells()); - for (Triangulation<2,2>::cell_iterator - cell = input.begin(); cell != input.end(); ++cell) - { - for (size_type slice=0; slice this_cell; - for (unsigned int v=0; v::vertices_per_cell; ++v) - { - this_cell.vertices[v] - = cell->vertex_index(v)+slice*input.n_vertices(); - this_cell.vertices[v+GeometryInfo<2>::vertices_per_cell] - = cell->vertex_index(v)+(slice+1)*input.n_vertices(); - } + for (unsigned int slice=0; slicematerial_id(); - cells.push_back(this_cell); - } - } + { + const Point<2> &v = input.get_vertices()[i]; + points[i+slice*input.n_vertices()](0) = v(0); + points[i+slice*input.n_vertices()](1) = v(1); + points[i+slice*input.n_vertices()](2) = height * slice / (n_slices-1); + } + } - SubCellData s; - types::boundary_id bid=0; - s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2); - for (Triangulation<2,2>::cell_iterator - cell = input.begin(); cell != input.end(); ++cell) - { - CellData<2> quad; - for (unsigned int f=0; f<4; ++f) - if (cell->at_boundary(f)) + for (Triangulation<2,2>::cell_iterator + cell = input.begin(); cell != input.end(); ++cell) + { + for (unsigned int slice=0; sliceface(f)->boundary_indicator(); - bid = std::max(bid, quad.boundary_id); - for (size_type slice=0; slice this_cell; + for (unsigned int v=0; v::vertices_per_cell; ++v) { - quad.vertices[0] = cell->face(f)->vertex_index(0)+slice*input.n_vertices(); - quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_vertices(); - quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_vertices(); - quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_vertices(); - s.boundary_quads.push_back(quad); + this_cell.vertices[v] + = cell->vertex_index(v)+slice*input.n_vertices(); + this_cell.vertices[v+GeometryInfo<2>::vertices_per_cell] + = cell->vertex_index(v)+(slice+1)*input.n_vertices(); } - } - } - for (Triangulation<2,2>::cell_iterator - cell = input.begin(); cell != input.end(); ++cell) - { - CellData<2> quad; - quad.boundary_id = bid + 1; - quad.vertices[0] = cell->vertex_index(0); - quad.vertices[1] = cell->vertex_index(1); - quad.vertices[2] = cell->vertex_index(2); - quad.vertices[3] = cell->vertex_index(3); - s.boundary_quads.push_back(quad); - - quad.boundary_id = bid + 2; - for (int i=0; i<4; ++i) - quad.vertices[i] += (n_slices-1)*input.n_vertices(); - s.boundary_quads.push_back(quad); - } + this_cell.material_id = cell->material_id(); + cells.push_back(this_cell); + } + } + SubCellData s; + types::boundary_id bid=0; + s.boundary_quads.reserve(input.n_active_lines()*(n_slices-1) + input.n_active_cells()*2); + for (Triangulation<2,2>::cell_iterator + cell = input.begin(); cell != input.end(); ++cell) + { + CellData<2> quad; + for (unsigned int f=0; f<4; ++f) + if (cell->at_boundary(f)) + { + quad.boundary_id = cell->face(f)->boundary_indicator(); + bid = std::max(bid, quad.boundary_id); + for (unsigned int slice=0; sliceface(f)->vertex_index(0)+slice*input.n_vertices(); + quad.vertices[1] = cell->face(f)->vertex_index(1)+slice*input.n_vertices(); + quad.vertices[2] = cell->face(f)->vertex_index(0)+(slice+1)*input.n_vertices(); + quad.vertices[3] = cell->face(f)->vertex_index(1)+(slice+1)*input.n_vertices(); + s.boundary_quads.push_back(quad); + } + } + } + for (Triangulation<2,2>::cell_iterator + cell = input.begin(); cell != input.end(); ++cell) + { + CellData<2> quad; + quad.boundary_id = bid + 1; + quad.vertices[0] = cell->vertex_index(0); + quad.vertices[1] = cell->vertex_index(1); + quad.vertices[2] = cell->vertex_index(2); + quad.vertices[3] = cell->vertex_index(3); + s.boundary_quads.push_back(quad); + + quad.boundary_id = bid + 2; + for (int i=0; i<4; ++i) + quad.vertices[i] += (n_slices-1)*input.n_vertices(); + s.boundary_quads.push_back(quad); + } - result.create_triangulation ( - points, - cells, - s); -} + result.create_triangulation (points, + cells, + s); + } -// make the following function inline. this is so that it becomes marked -// internal/weak for the linker and we don't get multiply defined errors -// when linking with more than one dimension at a time. Usually we used -// the trick of putting these functions in a .all_dimensions.cc file, but -// this is not necessary here as this is an internal only function. -inline -void GridGenerator::laplace_solve (const SparseMatrix &S, - const std::map &m, - Vector &u) -{ - const size_type n_dofs=S.n(); - FilteredMatrix > SF (S); - PreconditionJacobi > prec; - prec.initialize(S, 1.2); - FilteredMatrix > PF (prec); + /** + * Solve the Laplace equation for @p laplace_transformation function for one + * of the @p dim space dimensions. Factorized into a function of its own + * in order to allow parallel execution. + */ + void laplace_solve (const SparseMatrix &S, + const std::map &m, + Vector &u) + { + const unsigned int n_dofs=S.n(); + FilteredMatrix > SF (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); + SolverControl control (n_dofs, 1.e-10, false, false); + GrowingVectorMemory > mem; + SolverCG > solver (control, mem); - Vector f(n_dofs); + Vector f(n_dofs); - SF.add_constraints(m); - SF.apply_constraints (f, true); - solver.solve(SF, u, f, PF); -} + SF.add_constraints(m); + SF.apply_constraints (f, true); + solver.solve(SF, u, f, PF); + } // Implementation for 1D only -template <> -void GridGenerator::laplace_transformation (Triangulation<1> &, - const std::map > &) -{ - Assert(false, ExcNotImplemented()); -} + template <> + void laplace_transformation (Triangulation<1> &, + const std::map > &) + { + Assert(false, ExcNotImplemented()); + } // Implementation for dimensions except 1 -template -void GridGenerator::laplace_transformation (Triangulation &tria, - const std::map > &new_points) -{ - // first provide everything that is - // needed for solving a Laplace - // equation. - MappingQ1 mapping_q1; - FE_Q q1(1); - - DoFHandler dof_handler(tria); - dof_handler.distribute_dofs(q1); - - CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs (), - dof_handler.n_dofs ()); - DoFTools::make_sparsity_pattern (dof_handler, c_sparsity_pattern); - c_sparsity_pattern.compress (); - - SparsityPattern sparsity_pattern; - sparsity_pattern.copy_from (c_sparsity_pattern); - sparsity_pattern.compress (); - - SparseMatrix S(sparsity_pattern); - - QGauss quadrature(4); - - MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S); - - // set up the boundary values for - // the laplace problem - std::vector > m(dim); - typename std::map >::const_iterator map_iter; - typename std::map >::const_iterator map_end=new_points.end(); - - // fill these maps using the data - // given by new_points - typename DoFHandler::cell_iterator cell=dof_handler.begin_active(), - endc=dof_handler.end(); - typename DoFHandler::face_iterator face; - for (; cell!=endc; ++cell) - { - if (cell->at_boundary()) - for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) - { - face=cell->face(face_no); - if (face->at_boundary()) - for (unsigned int vertex_no=0; - vertex_no::vertices_per_face; ++vertex_no) - { - const size_type vertex_index=face->vertex_index(vertex_no); - map_iter=new_points.find(vertex_index); - - if (map_iter!=map_end) - for (unsigned int i=0; i ( - face->vertex_dof_index(vertex_no, 0), map_iter->second(i))); - } - } - } + template + void laplace_transformation (Triangulation &tria, + const std::map > &new_points) + { + // first provide everything that is + // needed for solving a Laplace + // equation. + MappingQ1 mapping_q1; + FE_Q q1(1); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(q1); + + CompressedSparsityPattern c_sparsity_pattern (dof_handler.n_dofs (), + dof_handler.n_dofs ()); + DoFTools::make_sparsity_pattern (dof_handler, c_sparsity_pattern); + c_sparsity_pattern.compress (); + + SparsityPattern sparsity_pattern; + sparsity_pattern.copy_from (c_sparsity_pattern); + sparsity_pattern.compress (); + + SparseMatrix S(sparsity_pattern); + + QGauss quadrature(4); + + MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, S); + + // set up the boundary values for + // the laplace problem + std::vector > m(dim); + typename std::map >::const_iterator map_iter; + typename std::map >::const_iterator map_end=new_points.end(); + + // fill these maps using the data + // given by new_points + typename DoFHandler::cell_iterator cell=dof_handler.begin_active(), + endc=dof_handler.end(); + typename DoFHandler::face_iterator face; + for (; cell!=endc; ++cell) + { + if (cell->at_boundary()) + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + face=cell->face(face_no); + if (face->at_boundary()) + for (unsigned int vertex_no=0; + vertex_no::vertices_per_face; ++vertex_no) + { + const unsigned int vertex_index=face->vertex_index(vertex_no); + map_iter=new_points.find(vertex_index); - // 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::vertices_per_cell; ++vertex_no) - { - Point &v=cell->vertex(vertex_no); - const size_type dof_index=cell->vertex_dof_index(vertex_no, 0); - for (unsigned int i=0; i ( + face->vertex_dof_index(vertex_no, 0), map_iter->second(i))); + } + } } -} + // 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::vertices_per_cell; ++vertex_no) + { + Point &v=cell->vertex(vertex_no); + const unsigned int dof_index=cell->vertex_dof_index(vertex_no, 0); + for (unsigned int i=0; i -void GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<1> &, - const double, - const double, - const double, - const size_type, - bool) -{ - Assert(false, ExcNotImplemented()); -} + template <> + void hyper_cube_with_cylindrical_hole (Triangulation<1> &, + const double, + const double, + const double, + const unsigned int, + bool) + { + Assert(false, ExcNotImplemented()); + } -template <> -void -GridGenerator::hyper_cube_with_cylindrical_hole (Triangulation<2> &triangulation, - const double inner_radius, - const double outer_radius, - const double, // width, - const size_type, // width_repetition, - bool colorize) -{ - const int dim = 2; - - Assert(inner_radius < outer_radius, - ExcMessage("outer_radius has to be bigger than inner_radius.")); - - Point center; - // We create an hyper_shell in two dimensions, and then we modify it. - GridGenerator::hyper_shell (triangulation, - center, inner_radius, outer_radius, - 8); - Triangulation::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); - std::vector treated_vertices(triangulation.n_vertices(), false); - for (; cell != endc; ++cell) - { - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->at_boundary()) - { - for (unsigned int v=0; v < GeometryInfo::vertices_per_face; ++v) - { - size_type vv = cell->face(f)->vertex_index(v); - if (treated_vertices[vv] == false) - { - treated_vertices[vv] = true; - switch (vv) - { - case 1: - cell->face(f)->vertex(v) = center+Point(outer_radius,outer_radius); - break; - case 3: - cell->face(f)->vertex(v) = center+Point(-outer_radius,outer_radius); - break; - case 5: - cell->face(f)->vertex(v) = center+Point(-outer_radius,-outer_radius); - break; - case 7: - cell->face(f)->vertex(v) = center+Point(outer_radius,-outer_radius); - default: - break; - } - } - } - } - } - double eps = 1e-3 * outer_radius; - cell = triangulation.begin_active(); - for (; cell != endc; ++cell) - { - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->at_boundary()) - { - double dx = cell->face(f)->center()(0) - center(0); - double dy = cell->face(f)->center()(1) - center(1); - 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 - cell->face(f)->set_boundary_indicator(4); - } - else - { - double d = (cell->face(f)->center() - center).norm(); - if (d-inner_radius < 0) - cell->face(f)->set_boundary_indicator(1); - else - cell->face(f)->set_boundary_indicator(0); - } - } - } -} -template <> -void GridGenerator::hyper_cube_with_cylindrical_hole(Triangulation<3> &triangulation, - const double inner_radius, - const double outer_radius, - const double L, - const size_type Nz, - bool colorize) -{ - const int dim = 3; - - Assert(inner_radius < outer_radius, - ExcMessage("outer_radius has to be bigger than inner_radius.")); - 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); - - Triangulation::active_cell_iterator - cell = triangulation.begin_active(), - endc = triangulation.end(); - std::vector treated_vertices(triangulation.n_vertices(), false); - for (; cell != endc; ++cell) - { - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->at_boundary()) - { - for (unsigned int v=0; v < GeometryInfo::vertices_per_face; ++v) - { - size_type vv = cell->face(f)->vertex_index(v); - if (treated_vertices[vv] == false) - { - treated_vertices[vv] = true; - for (size_type i=0; i<=Nz; ++i) - { - double d = ((double) i)*L/((double) Nz); - switch (vv-i*16) - { - case 1: - cell->face(f)->vertex(v) = Point(outer_radius,outer_radius,d); - break; - case 3: - cell->face(f)->vertex(v) = Point(-outer_radius,outer_radius,d); - break; - case 5: - cell->face(f)->vertex(v) = Point(-outer_radius,-outer_radius,d); - break; - case 7: - cell->face(f)->vertex(v) = Point(outer_radius,-outer_radius,d); - break; - default: - break; - } - } - } - } - } - } - double eps = 1e-3 * outer_radius; - cell = triangulation.begin_active(); - for (; cell != endc; ++cell) - { - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->at_boundary()) - { - double dx = cell->face(f)->center()(0); - double dy = cell->face(f)->center()(1); - double dz = cell->face(f)->center()(2); + template <> + void + hyper_cube_with_cylindrical_hole (Triangulation<2> &triangulation, + const double inner_radius, + const double outer_radius, + const double, // width, + const unsigned int, // width_repetition, + bool colorize) + { + const int dim = 2; + + Assert(inner_radius < outer_radius, + ExcMessage("outer_radius has to be bigger than inner_radius.")); + + Point center; + // We create an hyper_shell in two dimensions, and then we modify it. + hyper_shell (triangulation, + center, inner_radius, outer_radius, + 8); + Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + std::vector treated_vertices(triangulation.n_vertices(), false); + for (; cell != endc; ++cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary()) + { + for (unsigned int v=0; v < GeometryInfo::vertices_per_face; ++v) + { + unsigned int vv = cell->face(f)->vertex_index(v); + if (treated_vertices[vv] == false) + { + treated_vertices[vv] = true; + switch (vv) + { + case 1: + cell->face(f)->vertex(v) = center+Point(outer_radius,outer_radius); + break; + case 3: + cell->face(f)->vertex(v) = center+Point(-outer_radius,outer_radius); + break; + case 5: + cell->face(f)->vertex(v) = center+Point(-outer_radius,-outer_radius); + break; + case 7: + cell->face(f)->vertex(v) = center+Point(outer_radius,-outer_radius); + default: + break; + } + } + } + } + } + double eps = 1e-3 * outer_radius; + cell = triangulation.begin_active(); + for (; cell != endc; ++cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary()) + { + double dx = cell->face(f)->center()(0) - center(0); + double dy = cell->face(f)->center()(1) - center(1); + 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 + cell->face(f)->set_boundary_indicator(4); + } + else + { + double d = (cell->face(f)->center() - center).norm(); + if (d-inner_radius < 0) + cell->face(f)->set_boundary_indicator(1); + else + cell->face(f)->set_boundary_indicator(0); + } + } + } + } - 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); + template <> + void hyper_cube_with_cylindrical_hole(Triangulation<3> &triangulation, + const double inner_radius, + const double outer_radius, + const double L, + const unsigned int Nz, + bool colorize) + { + const int dim = 3; + + Assert(inner_radius < outer_radius, + ExcMessage("outer_radius has to be bigger than inner_radius.")); + Assert(L > 0, + ExcMessage("Must give positive extension L")); + Assert(Nz >= 1, ExcLowerRange(1, Nz)); + + cylinder_shell (triangulation, + L, inner_radius, outer_radius, + 8, + Nz); + + Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + std::vector treated_vertices(triangulation.n_vertices(), false); + for (; cell != endc; ++cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary()) + { + for (unsigned int v=0; v < GeometryInfo::vertices_per_face; ++v) + { + unsigned int vv = cell->face(f)->vertex_index(v); + if (treated_vertices[vv] == false) + { + treated_vertices[vv] = true; + for (unsigned int i=0; i<=Nz; ++i) + { + double d = ((double) i)*L/((double) Nz); + switch (vv-i*16) + { + case 1: + cell->face(f)->vertex(v) = Point(outer_radius,outer_radius,d); + break; + case 3: + cell->face(f)->vertex(v) = Point(-outer_radius,outer_radius,d); + break; + case 5: + cell->face(f)->vertex(v) = Point(-outer_radius,-outer_radius,d); + break; + case 7: + cell->face(f)->vertex(v) = Point(outer_radius,-outer_radius,d); + break; + default: + break; + } + } + } + } + } + } + double eps = 1e-3 * outer_radius; + cell = triangulation.begin_active(); + for (; cell != endc; ++cell) + { + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->at_boundary()) + { + double dx = cell->face(f)->center()(0); + double dy = cell->face(f)->center()(1); + double dz = cell->face(f)->center()(2); - else if (std::abs(dy - outer_radius) < eps) - cell->face(f)->set_boundary_indicator(3); + if (colorize) + { + if (std::abs(dx + outer_radius) < eps) + cell->face(f)->set_boundary_indicator(0); - else if (std::abs(dz) < eps) - cell->face(f)->set_boundary_indicator(4); + else if (std::abs(dx - outer_radius) < eps) + cell->face(f)->set_boundary_indicator(1); - else if (std::abs(dz - L) < eps) - cell->face(f)->set_boundary_indicator(5); + else if (std::abs(dy + outer_radius) < eps) + cell->face(f)->set_boundary_indicator(2); - else - { - cell->face(f)->set_boundary_indicator(6); - for (unsigned int l=0; l::lines_per_face; ++l) - cell->face(f)->line(l)->set_boundary_indicator(6); - } + else if (std::abs(dy - outer_radius) < eps) + cell->face(f)->set_boundary_indicator(3); - } - else - { - Point c = cell->face(f)->center(); - c(2) = 0; - double d = c.norm(); - if (d-inner_radius < 0) - { - cell->face(f)->set_boundary_indicator(1); - for (unsigned int l=0; l::lines_per_face; ++l) - cell->face(f)->line(l)->set_boundary_indicator(1); - } - else - cell->face(f)->set_boundary_indicator(0); - } - } - } -} + 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) + cell->face(f)->line(l)->set_boundary_indicator(6); + } + } + else + { + Point c = cell->face(f)->center(); + c(2) = 0; + double d = c.norm(); + if (d-inner_radius < 0) + { + cell->face(f)->set_boundary_indicator(1); + for (unsigned int l=0; l::lines_per_face; ++l) + cell->face(f)->line(l)->set_boundary_indicator(1); + } + else + cell->face(f)->set_boundary_indicator(0); + } + } + } + } +} // explicit instantiations #include "grid_generator.inst" diff --git a/deal.II/source/grid/grid_generator.inst.in b/deal.II/source/grid/grid_generator.inst.in index fadcea0bc7..39a7641c44 100644 --- a/deal.II/source/grid/grid_generator.inst.in +++ b/deal.II/source/grid/grid_generator.inst.in @@ -17,42 +17,45 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { +namespace GridGenerator +\{ #if deal_II_dimension <= deal_II_space_dimension template void - GridGenerator::hyper_rectangle ( + hyper_rectangle ( Triangulation &, const Point&, const Point&, const bool); template void - GridGenerator::colorize_hyper_rectangle ( + colorize_hyper_rectangle ( Triangulation &); template void - GridGenerator::hyper_cube ( + hyper_cube ( Triangulation &, const double, const double); template void - GridGenerator::merge_triangulations + merge_triangulations (const Triangulation &triangulation_1, const Triangulation &triangulation_2, Triangulation &result); #endif +\} } for (deal_II_dimension : DIMENSIONS) { - +namespace GridGenerator \{ template void - GridGenerator::subdivided_hyper_cube ( + subdivided_hyper_cube ( Triangulation &, const unsigned int, const double, const double); template void - GridGenerator::subdivided_hyper_rectangle + subdivided_hyper_rectangle (Triangulation &, const std::vector&, const Point&, @@ -61,7 +64,7 @@ for (deal_II_dimension : DIMENSIONS) template void - GridGenerator::subdivided_hyper_rectangle( + subdivided_hyper_rectangle( Triangulation &, const std::vector > &, const Point &, @@ -69,32 +72,31 @@ for (deal_II_dimension : DIMENSIONS) const bool ); template void - GridGenerator::parallelepiped ( + parallelepiped ( Triangulation&, const Point (&) [deal_II_dimension], const bool); template void - GridGenerator::subdivided_parallelepiped ( + subdivided_parallelepiped ( Triangulation&, - const size_type, + const unsigned int, const Point (&) [deal_II_dimension], const bool); template void - GridGenerator::subdivided_parallelepiped ( + subdivided_parallelepiped ( Triangulation&, - const size_type [deal_II_dimension], + const unsigned int [deal_II_dimension], const Point (&) [deal_II_dimension], const bool); #if deal_II_dimension > 1 template void - GridGenerator:: laplace_transformation (Triangulation &, const std::map > &); #endif - +\} } -- 2.39.5