From: bangerth Date: Sun, 1 Jul 2007 16:48:31 +0000 (+0000) Subject: Remove some of the stuff we won't need any more for the tutorial program X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6a9daecff988525f5a67c22c253107f8b7915f17;p=dealii-svn.git Remove some of the stuff we won't need any more for the tutorial program git-svn-id: https://svn.dealii.org/trunk@14841 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index 03c210f5c6..4c4d2d2319 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -1,10 +1,10 @@ /* $Id$ */ -/* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */ +/* Author: Wolfgang Bangerth, Texas A&M University, 2006, 2007 */ /* $Id$ */ /* Version: $Name$ */ /* */ -/* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007 by the deal.II authors */ +/* Copyright (C) 2006, 2007 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -25,7 +25,7 @@ #include #include #include -#include +#include #include #include #include @@ -42,15 +42,15 @@ #include #include -#include -#include - #include #include #include #include #include +#include +#include + #include @@ -62,7 +62,7 @@ template class LaplaceProblem { public: - LaplaceProblem (const bool condense_glob = true); + LaplaceProblem (); ~LaplaceProblem (); void run (); @@ -91,10 +91,7 @@ class LaplaceProblem Vector solution; Vector system_rhs; - Timer distr, condense, hang, assemble, solver; - const unsigned int max_degree; - const bool condense_global; }; @@ -138,10 +135,10 @@ RightHandSide::value (const Point &p, template -LaplaceProblem::LaplaceProblem (bool condense_glob) : - dof_handler (triangulation), - max_degree (dim == 2 ? 7 : 5), - condense_global (condense_glob) +LaplaceProblem::LaplaceProblem () + : + dof_handler (triangulation), + max_degree (dim == 2 ? 7 : 5) { for (unsigned int degree=2; degree<=max_degree; ++degree) { @@ -161,61 +158,22 @@ LaplaceProblem::~LaplaceProblem () template void LaplaceProblem::setup_system () { - distr.reset(); - distr.start(); dof_handler.distribute_dofs (fe_collection); - distr.stop(); solution.reinit (dof_handler.n_dofs()); system_rhs.reinit (dof_handler.n_dofs()); hanging_node_constraints.clear (); - hang.reset(); - hang.start(); DoFTools::make_hanging_node_constraints (dof_handler, hanging_node_constraints); hanging_node_constraints.close (); - hang.stop(); - if (dim < 3) - { - sparsity_pattern.reinit (dof_handler.n_dofs(), - dof_handler.n_dofs(), - dof_handler.max_couplings_between_dofs()); - condense.reset(); - if (condense_global) - { - DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); - condense.start(); - hanging_node_constraints.condense (sparsity_pattern); - condense.stop(); - } - else - DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern, - hanging_node_constraints); - - sparsity_pattern.compress(); - } - else - { - CompressedSparsityPattern csp (dof_handler.n_dofs(), - dof_handler.n_dofs()); - condense.reset(); - if (condense_global) - { - DoFTools::make_sparsity_pattern (dof_handler, csp); - - condense.start(); - hanging_node_constraints.condense (csp); - condense.stop(); - } - else - DoFTools::make_sparsity_pattern (dof_handler, csp, - hanging_node_constraints); - - sparsity_pattern.copy_from (csp); - } + CompressedSetSparsityPattern csp (dof_handler.n_dofs(), + dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, csp, + hanging_node_constraints); + sparsity_pattern.copy_from (csp); system_matrix.reinit (sparsity_pattern); } @@ -223,9 +181,6 @@ void LaplaceProblem::setup_system () template void LaplaceProblem::assemble_system () { - assemble.reset (); - assemble.start (); - hp::FEValues hp_fe_values (fe_collection, quadrature_collection, update_values | update_gradients | @@ -270,40 +225,15 @@ void LaplaceProblem::assemble_system () cell->get_dof_indices (local_dof_indices); - if (condense_global) - { - for (unsigned int i=0; i boundary_values; @@ -590,22 +520,13 @@ void LaplaceProblem::output_results (const unsigned int cycle) const Vector smoothness_indicators (triangulation.n_active_cells()); estimate_smoothness (smoothness_indicators); - Vector smoothness_field (dof_handler.n_dofs()); - DoFTools::distribute_cell_to_dof_vector (dof_handler, - smoothness_indicators, - smoothness_field); - Vector fe_indices (triangulation.n_active_cells()); { typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (unsigned int index=0; cell!=endc; ++cell, ++index) - { - - fe_indices(index) = cell->active_fe_index(); -// smoothness_indicators(index) *= std::sqrt(cell->diameter()); - } + fe_indices(index) = cell->active_fe_index(); } const std::string filename = "solution-" + @@ -615,8 +536,7 @@ void LaplaceProblem::output_results (const unsigned int cycle) const data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); - data_out.add_data_vector (smoothness_indicators, "smoothness1"); - data_out.add_data_vector (smoothness_field, "smoothness2"); + data_out.add_data_vector (smoothness_indicators, "smoothness"); data_out.add_data_vector (fe_indices, "fe_index"); data_out.build_patches (); @@ -698,586 +618,6 @@ void LaplaceProblem<2>::create_coarse_grid () -namespace BreastPhantom -{ - - - // 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, - 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) - { - CellData<3> cell; - for (unsigned int j=0; j<4; ++j) - { - cell.vertices[j] = connectivity[i][j]+first_upper_vertex; - cell.vertices[j+4] = connectivity[i][j]; - } - cell.material_id = 0; - cells.push_back (cell); - } - - // mark the 8 vertical edges with - // the correct boundary indicator - for (unsigned int i=0; i<8; ++i) - { - 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); - } - - // 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 - { - // 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]); - } - } - } - - - // 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); -} - - template void LaplaceProblem::run () @@ -1296,9 +636,6 @@ void LaplaceProblem::run () << triangulation.n_active_cells() << std::endl; - Timer all; - all.reset(); - all.start(); setup_system (); std::cout << " Number of degrees of freedom: " @@ -1311,20 +648,7 @@ void LaplaceProblem::run () assemble_system (); - solver.reset(); - solver.start(); solve (); - solver.stop(); - - all.stop(); - - std::cout << " All: " << all() - << ", distr: " << distr() - << ", hang: " << hang() - << ", condense: " << condense() - << ", assemble: " << assemble() - << ", solver: " << solver() - << std::endl; output_results (cycle); } @@ -1336,7 +660,7 @@ int main () { deallog.depth_console (0); - LaplaceProblem<3> laplace_problem (false); + LaplaceProblem<2> laplace_problem; laplace_problem.run (); } catch (std::exception &exc)