From: bangerth Date: Wed, 14 Jul 2010 13:47:45 +0000 (+0000) Subject: Take over patch r21495 from mainline. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=783cc95333631ed2d8670ebaf15253ae5d11268d;p=dealii-svn.git Take over patch r21495 from mainline. git-svn-id: https://svn.dealii.org/branches/releases/Branch-6-3@21496 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index c390193784..635ed0b1fe 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -47,8 +47,9 @@ inconvenience this causes.

  • Fixed: GCC version 3.3.x failed to compile the files - lac/include/lac/precondition_block.h and - deal.II/source/multigrid/mg_dof_handler.cc. This problem has + lac/include/lac/precondition_block.h, + deal.II/source/multigrid/mg_dof_handler.cc and + examples/step-34/step-34.cc. These problems have now been worked around.
    (WB 2010/07/12) diff --git a/deal.II/examples/step-34/step-34.cc b/deal.II/examples/step-34/step-34.cc index d3e6840db0..10baaa1a99 100644 --- a/deal.II/examples/step-34/step-34.cc +++ b/deal.II/examples/step-34/step-34.cc @@ -1,8 +1,8 @@ //---------------------------- step-34.cc --------------------------- // $Id$ -// Version: $Name$ +// Version: $Name$ // -// Copyright (C) 2009 by the deal.II authors +// Copyright (C) 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -89,16 +89,16 @@ namespace LaplaceKernel { case 2: return (-std::log(R.norm()) / (2*numbers::PI) ); - + case 3: return (1./( R.norm()*4*numbers::PI ) ); - + default: Assert(false, ExcInternalError()); return 0.; } } - + template @@ -110,7 +110,7 @@ namespace LaplaceKernel return R / ( -2*numbers::PI * R.square()); case 3: return R / ( -4*numbers::PI * R.square() * R.norm() ); - + default: Assert(false, ExcInternalError()); return Point(); @@ -137,8 +137,8 @@ namespace LaplaceKernel // methods, and we won't comment too // much on them, except on the // differences. -template -class BEMProblem +template +class BEMProblem { public: BEMProblem(); @@ -146,13 +146,13 @@ class BEMProblem void run(); private: - + void read_parameters (const std::string &filename); - + void read_domain(); void refine_and_resize(); - + // The only really different // function that we find here is // the assembly routine. We wrote @@ -260,7 +260,7 @@ class BEMProblem // edges and $\frac 78$ on the 8 // nodes of the vertices. void compute_errors(const unsigned int cycle); - + // Once we obtained a solution on // the codimension one domain, we // want to interpolate it to the @@ -293,9 +293,9 @@ class BEMProblem // output_results() function, of // course. void compute_exterior_solution(); - + void output_results(const unsigned int cycle); - + // The usual deal.II classes can // be used for boundary element // methods by specifying the @@ -336,7 +336,7 @@ class BEMProblem // is non trivial, and we don't // treat this subject here. - FullMatrix system_matrix; + FullMatrix system_matrix; Vector system_rhs; // The next two variables will @@ -350,13 +350,13 @@ class BEMProblem // shape functions. Vector phi; Vector alpha; - + // The convergence table is used // to output errors in the exact // solution and in the computed // alphas. ConvergenceTable convergence_table; - + // The following variables are // the ones that we fill through // a parameter file. The new @@ -411,7 +411,7 @@ class BEMProblem std_cxx1x::shared_ptr > quadrature; unsigned int singular_quadrature_order; - + SolverControl solver_control; unsigned int n_cycles; @@ -459,14 +459,14 @@ BEMProblem::BEMProblem() {} -template +template void BEMProblem::read_parameters (const std::string &filename) { deallog << std::endl << "Parsing parameter file " << filename << std::endl << "for a " << dim << " dimensional simulation. " << std::endl; - + ParameterHandler prm; - + prm.declare_entry("Number of cycles", "4", Patterns::Integer()); prm.declare_entry("External refinement", "5", @@ -477,16 +477,16 @@ void BEMProblem::read_parameters (const std::string &filename) Patterns::Bool()); prm.declare_entry("Run 3d simulation", "true", Patterns::Bool()); - + prm.enter_subsection("Quadrature rules"); { - prm.declare_entry("Quadrature type", "gauss", + prm.declare_entry("Quadrature type", "gauss", Patterns::Selection(QuadratureSelector<(dim-1)>::get_quadrature_names())); prm.declare_entry("Quadrature order", "4", Patterns::Integer()); prm.declare_entry("Singular quadrature order", "5", Patterns::Integer()); } prm.leave_subsection(); - + // For both two and three // dimensions, we set the default // input data to be such that the @@ -503,7 +503,7 @@ void BEMProblem::read_parameters (const std::string &filename) // solution we pass to the program // needs to have the same value at // infinity for the error to be - // computed correctly. + // computed correctly. // // The use of the // Functions::ParsedFunction object @@ -555,7 +555,7 @@ void BEMProblem::read_parameters (const std::string &filename) prm.set("Function expression", "1; 1; 1"); } prm.leave_subsection(); - + prm.enter_subsection("Exact solution 2d"); { Functions::ParsedFunction<2>::declare_parameters(prm); @@ -590,10 +590,10 @@ void BEMProblem::read_parameters (const std::string &filename) // ParameterHandler object: prm.read_input(filename); - n_cycles = prm.get_integer("Number of cycles"); + n_cycles = prm.get_integer("Number of cycles"); external_refinement = prm.get_integer("External refinement"); extend_solution = prm.get_bool("Extend solution on the -2,2 box"); - + prm.enter_subsection("Quadrature rules"); { quadrature = @@ -603,7 +603,7 @@ void BEMProblem::read_parameters (const std::string &filename) singular_quadrature_order = prm.get_integer("Singular quadrature order"); } prm.leave_subsection(); - + prm.enter_subsection(std::string("Wind function ")+ Utilities::int_to_string(dim)+std::string("d")); { @@ -632,14 +632,14 @@ void BEMProblem::read_parameters (const std::string &filename) // setting the corresponding "Run // 2d simulation" or "Run 3d // simulation" flag to false: - run_in_this_dimension = prm.get_bool("Run " + + run_in_this_dimension = prm.get_bool("Run " + Utilities::int_to_string(dim) + "d simulation"); } // @sect4{BEMProblem::read_domain} - + // A boundary element method // triangulation is basically the // same as a (dim-1) dimensional @@ -691,11 +691,12 @@ void BEMProblem::read_parameters (const std::string &filename) // still has to be static to live at // least as long as the triangulation // object to which it is attached. - + template void BEMProblem::read_domain() { - static HyperBallBoundary boundary(Point(),1.); + static const Point center = Point(); + static const HyperBallBoundary boundary(center,1.); std::ifstream in; switch (dim) @@ -703,7 +704,7 @@ void BEMProblem::read_domain() case 2: in.open ("coarse_circle.inp"); break; - + case 3: in.open ("coarse_sphere.inp"); break; @@ -730,17 +731,17 @@ template void BEMProblem::refine_and_resize() { tria.refine_global(1); - + dh.distribute_dofs(fe); - + const unsigned int n_dofs = dh.n_dofs(); - + system_matrix.reinit(n_dofs, n_dofs); - + system_rhs.reinit(n_dofs); phi.reinit(n_dofs); alpha.reinit(n_dofs); -} +} // @sect4{BEMProblem::assemble_system} @@ -773,13 +774,13 @@ void BEMProblem::refine_and_resize() // dimensional case. template void BEMProblem::assemble_system() -{ - std::vector > sing_quadratures_3d; +{ + std::vector > sing_quadratures_3d; for(unsigned int i=0; i<4; ++i) sing_quadratures_3d.push_back (QGaussOneOverR<2>(singular_quadrature_order, i, true)); - - + + // Next, we initialize an FEValues // object with the quadrature // formula for the integration of @@ -795,14 +796,14 @@ void BEMProblem::assemble_system() update_cell_normal_vectors | update_quadrature_points | update_JxW_values); - + const unsigned int n_q_points = fe_v.n_quadrature_points; - + std::vector local_dof_indices(fe.dofs_per_cell); std::vector > cell_wind(n_q_points, Vector(dim) ); double normal_wind; - + // Unlike in finite element // methods, if we use a collocation // boundary element method, then in @@ -820,7 +821,7 @@ void BEMProblem::assemble_system() // object will hold this // information: Vector local_matrix_row_i(fe.dofs_per_cell); - + // The index $i$ runs on the // collocation points, which are // the support points of the $i$th @@ -836,7 +837,7 @@ void BEMProblem::assemble_system() AssertThrow(fe.dofs_per_cell == GeometryInfo::vertices_per_cell, ExcMessage("The code in this function can only be used for " "the usual Q1 elements.")); - + // Now that we have checked that // the number of vertices is equal // to the number of degrees of @@ -859,17 +860,17 @@ void BEMProblem::assemble_system() typename DoFHandler::active_cell_iterator cell = dh.begin_active(), endc = dh.end(); - + for(cell = dh.begin_active(); cell != endc; ++cell) { fe_v.reinit(cell); cell->get_dof_indices(local_dof_indices); - + const std::vector > &q_points = fe_v.get_quadrature_points(); const std::vector > &normals = fe_v.get_cell_normal_vectors(); wind.vector_value_list(q_points, cell_wind); - - + + // We then form the integral over // the current cell for all // degrees of freedom (note that @@ -888,13 +889,13 @@ void BEMProblem::assemble_system() // one is the singular index: for(unsigned int i=0; i::assemble_system() for(unsigned int q=0; q R = q_points[q] - support_points[i]; - - system_rhs(i) += ( LaplaceKernel::single_layer(R) * + + system_rhs(i) += ( LaplaceKernel::single_layer(R) * normal_wind * fe_v.JxW(q) ); - + for(unsigned int j=0; j::assemble_system() // integral alltogether // that needs to be // evaluated: - // + // // \f[ // \int_0^1 f(x)\ln(x/\alpha) dx = // \int_0^1 f(x)\ln(x) dx - \int_0^1 f(x) \ln(\alpha) dx. @@ -1067,7 +1068,7 @@ void BEMProblem::assemble_system() // quadrature points and weights to // take into consideration also the // second part of the integral. - // + // // A similar reasoning // should be done in the // three dimensional @@ -1164,23 +1165,23 @@ void BEMProblem::assemble_system() : 0)); Assert(singular_quadrature, ExcInternalError()); - - FEValues fe_v_singular (fe, *singular_quadrature, + + FEValues fe_v_singular (fe, *singular_quadrature, update_jacobians | update_values | update_cell_normal_vectors | update_quadrature_points ); fe_v_singular.reinit(cell); - - std::vector > singular_cell_wind( (*singular_quadrature).size(), + + std::vector > singular_cell_wind( (*singular_quadrature).size(), Vector(dim) ); - + const std::vector > &singular_normals = fe_v_singular.get_cell_normal_vectors(); const std::vector > &singular_q_points = fe_v_singular.get_quadrature_points(); - + wind.vector_value_list(singular_q_points, singular_cell_wind); - + for(unsigned int q=0; qsize(); ++q) { const Point R = singular_q_points[q] - support_points[i]; @@ -1188,11 +1189,11 @@ void BEMProblem::assemble_system() for(unsigned int d=0; d::assemble_system() fe_v_singular.JxW(q) ); } } - if(dim==2) + if(dim==2) delete singular_quadrature; } - + // Finally, we need to add // the contributions of the // current cell to the // global matrix. - for(unsigned int j=0; j::assemble_system() // and the corresponding matrix is // a diagonal one with entries // equal to $\alpha(\mathbf{x}_i)$. - + // One quick way to compute this // diagonal matrix of the solid // angles, is to use the Neumann @@ -1240,7 +1241,7 @@ void BEMProblem::assemble_system() // matrix: Vector ones(dh.n_dofs()); ones.add(-1.); - + system_matrix.vmult(alpha, ones); alpha.add(1); for(unsigned int i = 0; i::compute_errors(const unsigned int cycle) VectorTools::L2_norm); const double L2_error = difference_per_cell.l2_norm(); - + // The error in the alpha vector // can be computed directly using // the Vector::linfty_norm() @@ -1293,12 +1294,12 @@ void BEMProblem::compute_errors(const unsigned int cycle) // rates: Vector difference_per_node(alpha); difference_per_node.add(-.5); - + const double alpha_error = difference_per_node.linfty_norm(); const unsigned int n_active_cells=tria.n_active_cells(); const unsigned int n_dofs=dh.n_dofs(); - - deallog << "Cycle " << cycle << ':' + + deallog << "Cycle " << cycle << ':' << std::endl << " Number of active cells: " << n_active_cells @@ -1306,7 +1307,7 @@ void BEMProblem::compute_errors(const unsigned int cycle) << " Number of degrees of freedom: " << n_dofs << std::endl; - + convergence_table.add_value("cycle", cycle); convergence_table.add_value("cells", n_active_cells); convergence_table.add_value("dofs", n_dofs); @@ -1353,12 +1354,12 @@ void BEMProblem::compute_exterior_solution() FE_Q external_fe(1); DoFHandler external_dh (external_tria); - Vector external_phi; - + Vector external_phi; + external_tria.refine_global(external_refinement); external_dh.distribute_dofs(external_fe); external_phi.reinit(external_dh.n_dofs()); - + typename DoFHandler::active_cell_iterator cell = dh.begin_active(), endc = dh.end(); @@ -1369,15 +1370,15 @@ void BEMProblem::compute_exterior_solution() update_cell_normal_vectors | update_quadrature_points | update_JxW_values); - + const unsigned int n_q_points = fe_v.n_quadrature_points; - + std::vector dofs(fe.dofs_per_cell); - + std::vector local_phi(n_q_points); std::vector normal_wind(n_q_points); std::vector > local_wind(n_q_points, Vector(dim) ); - + typename DoFHandler::active_cell_iterator external_cell = external_dh.begin_active(), external_endc = external_dh.end(); @@ -1385,47 +1386,47 @@ void BEMProblem::compute_exterior_solution() std::vector > external_support_points(external_dh.n_dofs()); DoFTools::map_dofs_to_support_points( StaticMappingQ1::mapping, external_dh, external_support_points); - + for(cell = dh.begin_active(); cell != endc; ++cell) { fe_v.reinit(cell); - + const std::vector > &q_points = fe_v.get_quadrature_points(); const std::vector > &normals = fe_v.get_cell_normal_vectors(); - + cell->get_dof_indices(dofs); fe_v.get_function_values(phi, local_phi); - + wind.vector_value_list(q_points, local_wind); - + for(unsigned int q=0; q R = q_points[q] - external_support_points[i]; - - external_phi(i) += ( ( LaplaceKernel::single_layer(R) * + + external_phi(i) += ( ( LaplaceKernel::single_layer(R) * normal_wind[q] + - (LaplaceKernel::double_layer(R) * + (LaplaceKernel::double_layer(R) * normals[q] ) * local_phi[q] ) * fe_v.JxW(q) ); } } - + DataOut data_out; - + data_out.attach_dof_handler(external_dh); data_out.add_data_vector(external_phi, "external_phi"); data_out.build_patches(); - + const std::string filename = Utilities::int_to_string(dim) + "d_external.vtk"; std::ofstream file(filename.c_str()); @@ -1445,28 +1446,28 @@ template void BEMProblem::output_results(const unsigned int cycle) { DataOut > dataout; - + dataout.attach_dof_handler(dh); dataout.add_data_vector(phi, "phi"); dataout.add_data_vector(alpha, "alpha"); dataout.build_patches(); - - std::string filename = ( Utilities::int_to_string(dim) + + + std::string filename = ( Utilities::int_to_string(dim) + "d_boundary_solution_" + Utilities::int_to_string(cycle) + ".vtk" ); std::ofstream file(filename.c_str()); - + dataout.write_vtk(file); - + if(cycle == n_cycles-1) { convergence_table.set_precision("L2(phi)", 3); convergence_table.set_precision("Linfty(alpha)", 3); - + convergence_table.set_scientific("L2(phi)", true); convergence_table.set_scientific("Linfty(alpha)", true); - + convergence_table .evaluate_convergence_rates("L2(phi)", ConvergenceTable::reduction_rate_log2); convergence_table @@ -1485,19 +1486,19 @@ void BEMProblem::output_results(const unsigned int cycle) template void BEMProblem::run() { - + read_parameters("parameters.prm"); if(run_in_this_dimension == false) { - deallog << "Run in dimension " << dim - << " explicitly disabled in parameter file. " + deallog << "Run in dimension " << dim + << " explicitly disabled in parameter file. " << std::endl; return; } - + read_domain(); - + for(unsigned int cycle=0; cycle::run() compute_errors(cycle); output_results(cycle); } - + if(extend_solution == true) compute_exterior_solution(); } @@ -1525,7 +1526,7 @@ int main () BEMProblem<2> laplace_problem_2d; laplace_problem_2d.run(); - BEMProblem<3> laplace_problem_3d; + BEMProblem<3> laplace_problem_3d; laplace_problem_3d.run(); } catch (std::exception &exc) @@ -1538,10 +1539,10 @@ int main () << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; - + return 1; } - catch (...) + catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------"