From 3069e7bbd541fa1511c2417f443c6333e293b1fe Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 25 May 2007 20:03:21 +0000 Subject: [PATCH] Check in a more interesting 3d geometry git-svn-id: https://svn.dealii.org/trunk@14708 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-27/step-27.cc | 677 ++++++++++++++++++++++++---- 1 file changed, 588 insertions(+), 89 deletions(-) diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index 05db6bac68..7e28967828 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -34,6 +34,7 @@ #include #include #include +#include #include #include #include @@ -112,10 +113,23 @@ double RightHandSide::value (const Point &p, const unsigned int /*component*/) const { - double product = 1; - for (unsigned int d=0; dstd::fabs(p[1]) ? 1 : 0); + + default: + Assert (false, ExcNotImplemented()); + } + return 0.; } @@ -124,7 +138,7 @@ RightHandSide::value (const Point &p, template LaplaceProblem::LaplaceProblem () : dof_handler (triangulation), - max_degree (dim == 2 ? 7 : 4) + max_degree (dim == 2 ? 7 : 5) { for (unsigned int degree=2; degree<=max_degree; ++degree) { @@ -555,7 +569,7 @@ void LaplaceProblem::output_results (const unsigned int cycle) const const std::string filename = "solution-" + Utilities::int_to_string (cycle, 2) + - ".vtk"; + ".gmv"; DataOut > data_out; data_out.attach_dof_handler (dof_handler); @@ -566,7 +580,7 @@ void LaplaceProblem::output_results (const unsigned int cycle) const data_out.build_patches (); std::ofstream output (filename.c_str()); - data_out.write_vtk (output); + data_out.write_gmv (output); } } @@ -643,99 +657,584 @@ void LaplaceProblem<2>::create_coarse_grid () -template <> -void LaplaceProblem<3>::create_coarse_grid () +namespace BreastPhantom { - const unsigned int dim = 3; - - // GridGenerator::hyper_cube (triangulation); - // triangulation.refine_global (1); - - // Create a hollow cube, in analogy to the 2D example. - // The grid generation is done in two steps. First the - // cell data is created on a uniform grid. In the - // second step, the unused vertices are removed. - const unsigned char hollow [4][4] = {{1,1,1,1}, - {1,0,0,1}, - {1,0,0,1}, - {1,1,1,1}}; - const unsigned char solid [4][4] = {{1,1,1,1}, - {1,1,1,1}, - {1,1,1,1}, - {1,1,1,1}}; - const unsigned char (*layers[4])[4][4] = {&solid, &hollow, &hollow, &solid}; - - std::vector > cells; - std::vector vertex_used (5*5*5, false); - - for (unsigned int zc = 0; zc < 4; ++zc) - for (unsigned int yc = 0; yc < 4; ++yc) - for (unsigned int xc = 0; xc < 4; ++xc) + + + // Radius of the sphere of the breast + // phantom geometry + static const double hemisphere_radius = 5; + + // Radius of the disk underneath + static const double bottom_disk_radius = 10; + + // Bottom z-coordinate of the disk + // underneath + const double bottom_disk_floor = -3; + // Top z-coordinate of the disk + // underneath + const double bottom_disk_ceil = -.5; + + // radius of the inner set of cells of + // the sphere geometry + const double interior_hemisphere_radius + = hemisphere_radius/(1.+std::sqrt(2.0)); + + template + class SphereBoundary : public HyperBallBoundary + { + public: + SphereBoundary () + : + HyperBallBoundary (Point(), hemisphere_radius) + {} + }; + + + template + class CylinderBoundary : public StraightBoundary + { + public: + typedef + typename Triangulation::line_iterator + line_iterator; + + typedef + typename Triangulation::quad_iterator + quad_iterator; + + typedef + typename Triangulation::face_iterator + face_iterator; + + /** + * Constructor. + */ + CylinderBoundary (const double radius) + : + radius (radius) + {} + + + virtual Point + get_new_point_on_line (const line_iterator &line) const; + + virtual Point + get_new_point_on_quad (const quad_iterator &quad) const; + + virtual void + get_intermediate_points_on_line (const line_iterator &line, + std::vector > &points) const; + + virtual void + get_intermediate_points_on_quad (const quad_iterator &quad, + std::vector > &points) const; + + virtual void + get_normals_at_vertices (const face_iterator &face, + typename Boundary::FaceVertexNormals &face_vertex_normals) const; + + private: + const double radius; + + void + get_intermediate_points_between_points (const Point &p0, const Point &p1, + std::vector > &points) const; + }; + + + template <> + Point<3> + CylinderBoundary<3>:: + get_new_point_on_line (const line_iterator &line) const + { + const Point<3> middle = StraightBoundary<3>::get_new_point_on_line (line); + // project to boundary + Point<3> p(middle[0], middle[1], 0); + p *= radius/std::sqrt(p.square()); + + return Point<3> (p[0], p[1], middle[2]); + } + + + template<> + Point<3> + CylinderBoundary<3>:: + get_new_point_on_quad (const quad_iterator &quad) const + { + Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad); + + // project to boundary + Point<3> p(middle[0], middle[1], 0); + p *= radius/std::sqrt(p.square()); + + return Point<3> (p[0], p[1], middle[2]); + } + + + template + void + CylinderBoundary:: + get_intermediate_points_on_line (const line_iterator &line, + std::vector > &points) const + { + if (points.size()==1) + points[0] = get_new_point_on_line(line); + else + get_intermediate_points_between_points(line->vertex(0), line->vertex(1), points); + } + + + template + void + CylinderBoundary:: + get_intermediate_points_between_points (const Point &, + const Point &, + std::vector > &) const + { + Assert (false, ExcNotImplemented()); + } + + + template <> + void + CylinderBoundary<3>:: + get_intermediate_points_on_quad (const Triangulation<3>::quad_iterator &, + std::vector > &) const + { + Assert (false, ExcNotImplemented()); + } + + + template + void + CylinderBoundary:: + get_normals_at_vertices (const typename Triangulation::face_iterator &, + typename Boundary::FaceVertexNormals &) const + { + Assert (false, ExcNotImplemented()); + } + + + void + create_coarse_grid (Triangulation<3> &coarse_grid) + { + const unsigned int dim = 3; + + std::vector > vertices; + std::vector > cells; + SubCellData sub_cell_data; + + const unsigned char + bottom_cylinder_boundary_id = 10, + middle_cylinder_boundary_id = 11, + spherical_boundary_id = 12, + all_other_boundary_id = 13, + straight_nondirichlet_boundary = 14; + + + // first build up the cells of the + // cylinder + { + // the vertices in each plane of + // the cylinder are located on + // three concentric rings of radii + // interior_hemisphere_radius, + // hemisphere_radius, and + // bottom_disk_radius, + // respectively. first generate + // these three rings + const Point<3> ring_points[8] = { Point<3>(-1,0,0), + Point<3>(-1,-1,0) / std::sqrt(2.), + Point<3>(0,-1,0), + Point<3>(+1,-1,0) / std::sqrt(2.), + Point<3>(+1,0,0), + Point<3>(+1,+1,0) / std::sqrt(2.), + Point<3>(0,+1,0), + Point<3>(-1,+1,0) / std::sqrt(2.) }; + + // first the point in the middle + // and the rest of those on the + // upper surface + vertices.push_back (Point<3>(0,0,bottom_disk_ceil)); + for (unsigned int ring=0; ring<3; ++ring) + for (unsigned int i=0; i<8; ++i) + vertices.push_back (ring_points[i] * (ring == 0 ? + interior_hemisphere_radius : + (ring == 1 ? hemisphere_radius : + bottom_disk_radius)) + + + Point<3>(0,0,bottom_disk_ceil)); + + // then points on lower surface + vertices.push_back (Point<3>(0,0,bottom_disk_floor)); + for (unsigned int ring=0; ring<3; ++ring) + for (unsigned int i=0; i<8; ++i) + vertices.push_back (ring_points[i] * (ring == 0 ? + interior_hemisphere_radius : + (ring == 1 ? + hemisphere_radius : + bottom_disk_radius)) + + + Point<3>(0,0,bottom_disk_floor)); + + const unsigned int n_vertices_per_surface = 25; + Assert (vertices.size() == n_vertices_per_surface*2, + ExcInternalError()); + + // next create cells from these + // vertices. only store the + // vertices of the upper surface, + // the lower ones are the same + // +12 + { + const unsigned int connectivity[20][4] + = { { 1, 2, 3, 0 }, // four cells in the center + { 3, 4, 5, 0 }, + { 0, 5, 6, 7 }, + { 1, 0, 7, 8 }, + + { 9, 10, 2, 1 }, // eight cells of inner ring + { 10, 11, 3, 2 }, + { 11, 12, 4, 3 }, + { 4, 12, 13, 5 }, + { 5, 13, 14, 6 }, + { 6, 14, 15, 7 }, + { 8, 7, 15, 16 }, + { 9, 1, 8, 16 }, + + { 17, 18, 10, 9 }, // eight cells of outer ring + { 18, 19, 11, 10 }, + { 19, 20, 12, 11 }, + { 12, 20, 21, 13 }, + { 13, 21, 22, 14 }, + { 14, 22, 23, 15 }, + { 16, 15, 23, 24 }, + { 17, 9, 16, 24 } }; + + // now create cells out of this + for (unsigned int i=0; i<20; ++i) + { + CellData<3> cell; + for (unsigned int j=0; j<4; ++j) + { + cell.vertices[j] = connectivity[i][j]; + cell.vertices[j+4] = connectivity[i][j]+n_vertices_per_surface; + } + cell.material_id = 0; + cells.push_back (cell); + } + } + + // associate edges and faces on the + // outer boundary with boundary + // indicator of the cylinder + // boundary indicator. do this the + // same way as above, just this + // time with faces (edges follow + // from this immediately. some + // edges are duplicated since they + // belong to more than one cell, + // but that doesn't harm us here) + { + const unsigned int connectivity[8][2] + = { { 17,18 }, { 18, 19 }, { 19, 20 }, { 20, 21 }, + { 21,22 }, { 22, 23 }, { 23, 24 }, { 24, 17 }}; + + for (unsigned int i=0; i<8; ++i) + { + const CellData<2> face = + { { connectivity[i][0]+n_vertices_per_surface, + connectivity[i][1]+n_vertices_per_surface, + connectivity[i][1], + connectivity[i][0] }, + bottom_cylinder_boundary_id }; + sub_cell_data.boundary_quads.push_back (face); + + const CellData<1> edges[4] = + { { { connectivity[i][0], connectivity[i][1] }, + bottom_cylinder_boundary_id }, + { { connectivity[i][0]+n_vertices_per_surface, + connectivity[i][1]+n_vertices_per_surface }, + bottom_cylinder_boundary_id }, + { { connectivity[i][0]+n_vertices_per_surface, + connectivity[i][0] }, + bottom_cylinder_boundary_id }, + { { connectivity[i][1]+n_vertices_per_surface, + connectivity[i][1] }, + bottom_cylinder_boundary_id } }; + for (unsigned int i=0; i<4; ++i) + sub_cell_data.boundary_lines.push_back (edges[i]); + } + } + } + + // next build up the middle ring. for + // this, copy the first 17 vertices + // up to z=0 + { + const unsigned int first_upper_vertex = vertices.size(); + + for (unsigned int i=0; i<17; ++i) + vertices.push_back (Point<3>(vertices[i][0], vertices[i][1], 0)); + + // next create cells from these + // vertices. only store the + // vertices of the lower surface, + // the lower ones are the same + // +first_upper_vertex + const unsigned int connectivity[12][4] + = { { 1, 2, 3, 0 }, // four cells in the center + { 3, 4, 5, 0 }, + { 0, 5, 6, 7 }, + { 1, 0, 7, 8 }, + + { 9, 10, 2, 1 }, // eight cells of ring + { 10, 11, 3, 2 }, + { 11, 12, 4, 3 }, + { 4, 12, 13, 5 }, + { 5, 13, 14, 6 }, + { 6, 14, 15, 7 }, + { 8, 7, 15, 16 }, + { 9, 1, 8, 16 }}; + // now create cells out of this + for (unsigned int i=0; i<12; ++i) { - // Check if we have to create a cell - if ((*layers[zc])[xc][yc] == 1) + CellData<3> cell; + for (unsigned int j=0; j<4; ++j) { - const unsigned int z_vert = 25; - const unsigned int y_vert = 5; - unsigned int zoffs = zc * z_vert; - unsigned int yoffs = yc * y_vert; - unsigned int base_vert = zoffs + yoffs + xc; - - CellData cell; - cell.vertices[0] = base_vert; - cell.vertices[1] = cell.vertices[0] + 1; - cell.vertices[2] = cell.vertices[0] + y_vert; - cell.vertices[3] = cell.vertices[1] + y_vert; - cell.vertices[4] = cell.vertices[0] + z_vert; - cell.vertices[5] = cell.vertices[1] + z_vert; - cell.vertices[6] = cell.vertices[2] + z_vert; - cell.vertices[7] = cell.vertices[3] + z_vert; - cell.material_id = 0; - cells.push_back (cell); - - // Now add entries to the list of used - // vertices. - for (unsigned int i = 0; i < 8; ++i) - vertex_used[cell.vertices[i]] = true; + cell.vertices[j] = connectivity[i][j]+first_upper_vertex; + cell.vertices[j+4] = connectivity[i][j]; } + cell.material_id = 0; + cells.push_back (cell); } - // Now create vertices and renumber stuff; - std::vector > vertices; - std::vector vert_renumber (5*5*5, 0); - const double scale = 0.5; - unsigned int v_indx = 0; - - for (int zv = 0; zv < 5; ++zv) - for (int yv = 0; yv < 5; ++yv) - for (int xv = 0; xv < 5; ++xv) + // mark the 8 vertical edges with + // the correct boundary indicator + for (unsigned int i=0; i<8; ++i) { - Point p_new ((double)(xv-2) * scale, - (double)(yv-2) * scale, - (double)(zv-2) * scale); - - if (vertex_used[v_indx]) - { - vert_renumber[v_indx] = vertices.size (); - vertices.push_back (p_new); - } - v_indx++; + const CellData<1> edge = { { 9, 9+first_upper_vertex }, + middle_cylinder_boundary_id }; + sub_cell_data.boundary_lines.push_back (edge); + } + // likewise with the 8 tangential + // edges on the lower disk. the + // edges at the interface between + // the middle disk and the + // hemisphere are handled by the + // hemisphere boundary + for (unsigned int i=0; i<8; ++i) + { + const CellData<1> edge = { { 9+i, 9+(i+1)%8}, + middle_cylinder_boundary_id }; + sub_cell_data.boundary_lines.push_back (edge); } - // Finally renumber the vertex indices in the cells - std::vector >::iterator cell_iterator; - for (cell_iterator = cells.begin (); cell_iterator != cells.end (); - ++cell_iterator) + // then assign face indicators + for (unsigned int i=0; i<8; ++i) + { + const CellData<2> face = { { 9+i, + 9+(i+1)%8, + 9+(i+1)%8+first_upper_vertex, + 9+i+first_upper_vertex}, + middle_cylinder_boundary_id }; + sub_cell_data.boundary_quads.push_back (face); + } + } + + // the final part is setting the + // half-sphere on top of this { - for (unsigned int i = 0; i < 8; ++i) - cell_iterator->vertices[i] = - vert_renumber[cell_iterator->vertices[i]]; + // add four cubes to the top of + // the inner four cells, as well + // as 8 to their outside + { + // mirror the first nine vertices + // above the surface, and scale + // them to a certain distance + // outward + const double rx = hemisphere_radius / (1+std::sqrt(3.0)); + for (unsigned int i=0; i<9; ++i) + { + const Point<3> p (vertices[i][0], + vertices[i][1], + i == 0 ? + 1 + : + std::max(std::fabs(vertices[i][0]), + std::fabs(vertices[i][1]))); + vertices.push_back (p / std::sqrt(p.square()) * rx); + } + Assert (vertices.size() == 76, ExcInternalError()); + + // same with the next ring of + // vertices, except that they + // go to hemisphere_radius + for (unsigned int i=9; i<17; ++i) + { + Point<3> p (vertices[i][0], + vertices[i][1], + std::max(std::fabs(vertices[i][0]), + std::fabs(vertices[i][1]))); + vertices.push_back (p / std::sqrt(p.square()) * + hemisphere_radius); + } + Assert (vertices.size() == 84, ExcInternalError()); + + // make 12 cells out of this + const unsigned int connectivity[12][4] + = { { 1, 2, 3, 0 }, // four cells in the center + { 3, 4, 5, 0 }, + { 0, 5, 6, 7 }, + { 1, 0, 7, 8 }, + + { 9, 10, 2, 1 }, // eight cells of inner ring + { 10, 11, 3, 2 }, + { 11, 12, 4, 3 }, + { 4, 12, 13, 5 }, + { 5, 13, 14, 6 }, + { 6, 14, 15, 7 }, + { 8, 7, 15, 16 }, + { 9, 1, 8, 16 }, + }; + + for (unsigned int i=0; i<12; ++i) + { + CellData<3> cell; + for (unsigned int j=0; j<4; ++j) + { + cell.vertices[j] = connectivity[i][j]+67; + cell.vertices[j+4] = connectivity[i][j]+50; + } + cell.material_id = 0; + cells.push_back (cell); + } + } + + // assign boundary indicators to + // the faces and edges of these + // cells + { + // these are the numbers of the + // vertices on the top surface + // of the cylinder, with one + // "wrap-around": + const unsigned int vertices[9] = + { 9, 10, 11, 12, 13, 14, 15, 16, 9 }; + // their counter-parts are the + // same +67 + for (unsigned int i=0; i<8; ++i) + { + // generate a face + const CellData<2> face = + { { vertices[i]+50, vertices[i+1]+50, + vertices[i+1]+67, vertices[i]+67 }, + spherical_boundary_id }; + sub_cell_data.boundary_quads.push_back (face); + + // same for the faces + const CellData<1> edges[4] = + { { { vertices[i]+50, vertices[i+1]+50 }, + spherical_boundary_id }, + { { vertices[i]+67, vertices[i+1]+67 }, + spherical_boundary_id }, + { { vertices[i]+50, vertices[i]+67 }, + spherical_boundary_id }, + { { vertices[i+1]+50, vertices[i+1]+67 }, + spherical_boundary_id } }; + for (unsigned int j=0; j<4; ++j) + sub_cell_data.boundary_lines.push_back (edges[j]); + } + } + + + // finally top the building + // with four closing cells and + // the vertex at the top + { + vertices.push_back (Point<3> (0,0,hemisphere_radius)); + + const unsigned int connectivity[4][8] + = { { 59, 60, 61, 67, 51, 52, 53, 50 }, + { 61, 62, 63, 67, 53, 54, 55, 50 }, + { 67, 63, 64, 65, 50, 55, 56, 57 }, + { 59, 67, 65, 66, 51, 50, 57, 58 }}; + + for (unsigned int i=0; i<4; ++i) + { + CellData<3> cell; + for (unsigned int j=0; j<8; ++j) + cell.vertices[j] = connectivity[i][j]+17; + cell.material_id = 0; + cells.push_back (cell); + } + + // generate boundary + // information for these cells, + // too + for (unsigned int i=0; i<4; ++i) + { + const CellData<2> face = + { { connectivity[i][0]+17, connectivity[i][1]+17, + connectivity[i][2]+17, connectivity[i][3]+17 }, + spherical_boundary_id }; + sub_cell_data.boundary_quads.push_back (face); + + const CellData<1> edges[4] = + { { { connectivity[i][0]+17, connectivity[i][1]+17 }, + spherical_boundary_id }, + { { connectivity[i][1]+17, connectivity[i][2]+17 }, + spherical_boundary_id }, + { { connectivity[i][2]+17, connectivity[i][3]+17 }, + spherical_boundary_id }, + { { connectivity[i][3]+17, connectivity[i][0]+17 }, + spherical_boundary_id } }; + for (unsigned int j=0; j<4; ++j) + sub_cell_data.boundary_lines.push_back (edges[j]); + } + } } + - // Now create triangulation - triangulation.create_triangulation (vertices, - cells, - SubCellData()); + // finally generate a triangulation + // out of this + GridReordering<3>::reorder_cells (cells); + coarse_grid.create_triangulation_compatibility (vertices, cells, + sub_cell_data); + + // then associate boundary objects + // with the different boundary + // indicators + static const CylinderBoundary<3> + bottom_cylinder_boundary (bottom_disk_radius); + static const CylinderBoundary<3> + middle_cylinder_boundary (hemisphere_radius); + static const SphereBoundary<3> sphere_boundary; + + coarse_grid.set_boundary (bottom_cylinder_boundary_id, + bottom_cylinder_boundary); + coarse_grid.set_boundary (middle_cylinder_boundary_id, + middle_cylinder_boundary); + coarse_grid.set_boundary (spherical_boundary_id, + sphere_boundary); + + for (Triangulation::active_cell_iterator cell=coarse_grid.begin_active(); + cell != coarse_grid.end(); ++cell) + for (unsigned int f=0; f::faces_per_cell; ++f) + if ((cell->face(f)->boundary_indicator() == 0) + && + (cell->face(f)->center()[2] >= (bottom_disk_floor+bottom_disk_ceil)/2)) + cell->face(f)->set_boundary_indicator(straight_nondirichlet_boundary); + } +} + + + +template <> +void LaplaceProblem<3>::create_coarse_grid () +{ + BreastPhantom::create_coarse_grid (triangulation); } @@ -799,7 +1298,7 @@ int main () { deallog.depth_console (0); - LaplaceProblem<2> laplace_problem; + LaplaceProblem<3> laplace_problem; laplace_problem.run (); } catch (std::exception &exc) -- 2.39.5