From: Wolfgang Bangerth Date: Tue, 9 Sep 2003 14:11:48 +0000 (+0000) Subject: Include Joergs subdivided_hyper_rectangle. X-Git-Tag: v8.0.0~16244 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=362ca1777d16ff09322918a7d59b5204e049ef22;p=dealii.git Include Joergs subdivided_hyper_rectangle. git-svn-id: https://svn.dealii.org/trunk@7966 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/grid/grid_generator.h b/deal.II/deal.II/include/grid/grid_generator.h index 526716da81..2e2e7233eb 100644 --- a/deal.II/deal.II/include/grid/grid_generator.h +++ b/deal.II/deal.II/include/grid/grid_generator.h @@ -49,6 +49,12 @@ template class SparseMatrix; * @ref{GridGenerator}@p{::hyper_rectangle (tria, p1, p2)}, with two * opposite corner points @p{p1} and @p{p2}. * + * @item Rectangular coordinate-parallel domains with different numbers of + * cells in each coordinate direction are generated by + * @ref{GridGenerator}@p{::subdivided_hyper_rectangle + * (tria, int[dim] repetitions, p1, p2, colorize)}, with two + * opposite corner points @p{p1} and @p{p2}. + * * @item Generalized L-shape domain: * using the @ref{GridGenerator}@p{::hyper_L (tria, a,b)} function produces * the hypercube with the interval $[a,b]$ without the hypercube @@ -84,7 +90,7 @@ template class SparseMatrix; * the triangulation to be refined at the outer boundaries. The class * @ref{HalfHyperBallBoundary} will provide a boundary object. * -* @item Hyper shell: A hyper shell is the region between two hyper + * @item Hyper shell: A hyper shell is the region between two hyper * sphere with the same origin. Therefore, it is a ring in two * spatial dimensions. To triangulation it, call the function * @ref{GridGenerator}@p{::hyper_shell (tria, origin, inner_radius, outer_radius, N)}, @@ -131,7 +137,7 @@ template class SparseMatrix; * transform (simple-shaped) grids to a more complicated ones, like a * shell onto a grid of an airfoil, for example. * - * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat, Stefan Nauber, 1998, 1999, 2000, 2001, 2002, 2003. + * @author Wolfgang Bangerth, Ralf Hartmann, Guido Kanschat, Stefan Nauber, Joerg Weimar, 1998, 1999, 2000, 2001, 2002, 2003. */ class GridGenerator { @@ -205,6 +211,49 @@ class GridGenerator const Point &p2, const bool colorize = false); + /** + * Create a coordinate-parallel + * parallelepiped from the two + * diagonally opposite corner + * points @p{p1} and @p{p2}. In + * dimension @p{i}, + * @p{repetitions[i]} cells are + * generated. + * + * To get cells with an aspect + * ratio different from that of + * the recatangle, use different + * numbers of subdivisions in + * different coordinate + * directions. The minimum number + * of subdivisions in each + * direction is + * 1. @p{repetitions} is a list + * of integers denoting the + * number of subdivisions in each + * coordinate direction. + * + * If the @p{colorize} flag is + * set, the + * @p{boundary_indicator}s of the + * surfaces are assigned, such + * that the lower one in + * @p{x}-direction is 0, the + * upper one is 1. The indicators + * for the surfaces in + * @p{y}-direction are 2 and 3, + * the ones for @p{z} are 4 and + * 5. + */ + template + static + void + subdivided_hyper_rectangle (Triangulation &tria, + const std::vector &repetitions, + const Point &p1, + const Point &p2, + const bool colorize=false); + /** * Hypercube with a layer of * hypercubes around it. The @@ -591,6 +640,13 @@ class GridGenerator int, << "The number of repetitions " << arg1 << " must be >=1."); + /** + * Exception + */ + DeclException1 (ExcInvalidRepetitionsDimension, + int, + << "The vector of repetitions must have " + << arg1 <<" elements."); private: /** @@ -608,6 +664,35 @@ class GridGenerator */ static void colorize_hyper_rectangle (Triangulation<1> &tria); + /** + * Perform the action specified + * by the @p{colorize} flag of + * the + * @ref{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); + + /** + * Declaration of the above + * function for 1d. + */ + static + void + colorize_subdivided_hyper_rectangle (Triangulation<1> &tria, + const Point<1> &p1, + const Point<1> &p2, + const double epsilon); + /** * Solve the Laplace equation for * @p{laplace_transformation} @@ -617,9 +702,11 @@ class GridGenerator * a function of its own in order * to allow parallel execution. */ - static void laplace_solve (const SparseMatrix &S, - const std::map &m, - Vector &u); + static + void + laplace_solve (const SparseMatrix &S, + const std::map &m, + Vector &u); }; diff --git a/deal.II/deal.II/source/grid/grid_generator.cc b/deal.II/deal.II/source/grid/grid_generator.cc index 8645257c57..f55c3224d5 100644 --- a/deal.II/deal.II/source/grid/grid_generator.cc +++ b/deal.II/deal.II/source/grid/grid_generator.cc @@ -243,6 +243,203 @@ GridGenerator::subdivided_hyper_cube (Triangulation &tria, } + +template +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. + double delta[dim]; + + for (unsigned int i=0; i= 1, ExcInvalidRepetitions(repetitions[i])); + + delta[i] = (p2[i]-p1[i])/repetitions[i]; + } + + // then generate the necessary + // points + std::vector > points; + switch (dim) + { + case 1: + for (unsigned int x=0; x<=repetitions[0]; ++x) + points.push_back (Point (p1[0]+x*delta[0])); + break; + + 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])); + break; + + case 3: + for (unsigned int z=0; z<=repetitions[2]; ++z) + 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 1) + epsilon = std::min(epsilon,0.01*delta[1]); + + if (dim > 2) + epsilon = std::min(epsilon,0.01*delta[2]); + + // actual code is external since + // 1-D is different from 2/3D. + colorize_subdivided_hyper_rectangle (tria,p1,p2,epsilon); + } +} + + +#if deal_II_dimension == 1 + +void +GridGenerator::colorize_subdivided_hyper_rectangle (Triangulation<1> &, + const Point<1> &, + const Point<1> &, + const double ) +{ + // nothing to do in 1d + // boundary indicators are set to + // 0 (left) and 1 (right) by default. +} + +#else + +template +void +GridGenerator::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::raw_face_iterator face = tria.begin_raw_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()); + + } + } +} + +#endif + + + #if deal_II_dimension == 1 void GridGenerator::hyper_cube_slit (Triangulation<1> &, @@ -348,15 +545,15 @@ GridGenerator::hyper_cube_slit (Triangulation<2> &tria, { 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 ) }; + 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,2,3 }, { 9,4,5,2 }, { 3,2,6,7 }, @@ -390,13 +587,13 @@ GridGenerator::hyper_L (Triangulation<2> &tria, { const unsigned int dim=2; const Point vertices[8] = { Point (a,a), - Point ((a+b)/2,a), - Point (b,a), - Point (a,(a+b)/2), - Point ((a+b)/2,(a+b)/2), - Point (b,(a+b)/2), - Point (a,b), - Point ((a+b)/2,b) }; + Point ((a+b)/2,a), + Point (b,a), + Point (a,(a+b)/2), + Point ((a+b)/2,(a+b)/2), + Point (b,(a+b)/2), + Point (a,b), + Point ((a+b)/2,b) }; const int cell_vertices[3][4] = {{0, 1, 4, 3}, {1, 2, 5, 4}, {3, 4, 7, 6}}; @@ -427,13 +624,13 @@ GridGenerator::hyper_ball (Triangulation<2> &tria, // 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)) }; + 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, 3, 2}, {0, 2, 4, 6}, @@ -534,14 +731,14 @@ GridGenerator::cylinder (Triangulation<2> &tria, switch (f->boundary_indicator()) { case 0: - f->set_boundary_indicator(1); - break; + f->set_boundary_indicator(1); + break; case 1: - f->set_boundary_indicator(2); - break; + f->set_boundary_indicator(2); + break; default: - f->set_boundary_indicator(0); - break; + f->set_boundary_indicator(0); + break; } ++f; } @@ -559,13 +756,13 @@ GridGenerator::half_hyper_ball (Triangulation<2> &tria, // 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)) }; + 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, 3, 2}, {2, 3, 5, 4}, @@ -739,37 +936,37 @@ GridGenerator::hyper_L (Triangulation<3> &tria, // part of the cube const Point vertices[26] = { - // front face of the big cube - Point (a, a,a), - Point ((a+b)/2,a,a), - Point (b, a,a), - Point (a, a,(a+b)/2), - Point ((a+b)/2,a,(a+b)/2), - Point (b, a,(a+b)/2), - Point (a, a,b), - Point ((a+b)/2,a,b), - Point (b, a,b), - // middle face of the big cube - Point (a, (a+b)/2,a), - Point ((a+b)/2,(a+b)/2,a), - Point (b, (a+b)/2,a), - Point (a, (a+b)/2,(a+b)/2), - Point ((a+b)/2,(a+b)/2,(a+b)/2), - Point (b, (a+b)/2,(a+b)/2), - Point (a, (a+b)/2,b), - Point ((a+b)/2,(a+b)/2,b), - Point (b, (a+b)/2,b), - // back face of the big cube - // last (top right) point is missing - Point (a, b,a), - Point ((a+b)/2,b,a), - Point (b, b,a), - Point (a, b,(a+b)/2), - Point ((a+b)/2,b,(a+b)/2), - Point (b, b,(a+b)/2), - Point (a, b,b), - Point ((a+b)/2,b,b) - }; + // front face of the big cube + Point (a, a,a), + Point ((a+b)/2,a,a), + Point (b, a,a), + Point (a, a,(a+b)/2), + Point ((a+b)/2,a,(a+b)/2), + Point (b, a,(a+b)/2), + Point (a, a,b), + Point ((a+b)/2,a,b), + Point (b, a,b), + // middle face of the big cube + Point (a, (a+b)/2,a), + Point ((a+b)/2,(a+b)/2,a), + Point (b, (a+b)/2,a), + Point (a, (a+b)/2,(a+b)/2), + Point ((a+b)/2,(a+b)/2,(a+b)/2), + Point (b, (a+b)/2,(a+b)/2), + Point (a, (a+b)/2,b), + Point ((a+b)/2,(a+b)/2,b), + Point (b, (a+b)/2,b), + // back face of the big cube + // last (top right) point is missing + Point (a, b,a), + Point ((a+b)/2,b,a), + Point (b, b,a), + Point (a, b,(a+b)/2), + Point ((a+b)/2,b,(a+b)/2), + Point (b, b,(a+b)/2), + Point (a, b,b), + Point ((a+b)/2,b,b) + }; const int cell_vertices[7][8] = {{0, 1, 4, 3, 9, 10, 13, 12}, {1, 2, 5, 4, 10, 11, 14, 13}, {3, 4, 7, 6, 12, 13, 16, 15}, @@ -805,27 +1002,27 @@ GridGenerator::hyper_ball (Triangulation<3> &tria, 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)), - }; + // 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 @@ -866,31 +1063,31 @@ GridGenerator::cylinder (Triangulation<3> &tria, 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), - }; + 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) { @@ -995,7 +1192,7 @@ void GridGenerator::laplace_transformation (Triangulation &tria, // fill these maps using the data // given by new_points typename DoFHandler::cell_iterator cell=dof_handler.begin_active(), - endc=dof_handler.end(); + endc=dof_handler.end(); typename DoFHandler::face_iterator face; for (; cell!=endc; ++cell) { @@ -1013,7 +1210,7 @@ void GridGenerator::laplace_transformation (Triangulation &tria, for (unsigned int i=0; i ( - face->vertex_dof_index(vertex_no, 0), map_iter->second(i))); + face->vertex_dof_index(vertex_no, 0), map_iter->second(i))); } } } @@ -1088,6 +1285,14 @@ GridGenerator::subdivided_hyper_cube (Triangulation +(Triangulation &, + const std::vector&, + const Point&, + const Point&, + const bool); + #if deal_II_dimension != 1 template void