From bde129fb865c68587f0a6217c7c532c3f93423bc Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 6 Aug 2013 13:30:49 +0000 Subject: [PATCH] First version of new tutorial example demonstrating hybridizable discontinuous Galerkin methods on the problem of step-7 git-svn-id: https://svn.dealii.org/trunk@30236 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-51/CMakeLists.txt | 31 + deal.II/examples/step-51/doc/builds-on | 1 + deal.II/examples/step-51/doc/intro.dox | 5 + deal.II/examples/step-51/doc/kind | 1 + deal.II/examples/step-51/doc/results.dox | 139 ++++ deal.II/examples/step-51/doc/tooltip | 1 + deal.II/examples/step-51/step-51.cc | 893 +++++++++++++++++++++++ 7 files changed, 1071 insertions(+) create mode 100644 deal.II/examples/step-51/CMakeLists.txt create mode 100644 deal.II/examples/step-51/doc/builds-on create mode 100644 deal.II/examples/step-51/doc/intro.dox create mode 100644 deal.II/examples/step-51/doc/kind create mode 100644 deal.II/examples/step-51/doc/results.dox create mode 100644 deal.II/examples/step-51/doc/tooltip create mode 100644 deal.II/examples/step-51/step-51.cc diff --git a/deal.II/examples/step-51/CMakeLists.txt b/deal.II/examples/step-51/CMakeLists.txt new file mode 100644 index 0000000000..6d7622ff1d --- /dev/null +++ b/deal.II/examples/step-51/CMakeLists.txt @@ -0,0 +1,31 @@ +## +# CMake script for the step-7 tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "step-51") + +# Declare all source files the target consists of: +SET(TARGET_SRC + ${TARGET}.cc + # You can specify additional files here! + ) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/deal.II/examples/step-51/doc/builds-on b/deal.II/examples/step-51/doc/builds-on new file mode 100644 index 0000000000..943473db87 --- /dev/null +++ b/deal.II/examples/step-51/doc/builds-on @@ -0,0 +1 @@ +step-7 diff --git a/deal.II/examples/step-51/doc/intro.dox b/deal.II/examples/step-51/doc/intro.dox new file mode 100644 index 0000000000..a5458f74da --- /dev/null +++ b/deal.II/examples/step-51/doc/intro.dox @@ -0,0 +1,5 @@ + +

Introduction

+ +This tutorial program presents the implementation of a hybridizable +discontinuous Galkerin method for the convection-diffusion equation. \ No newline at end of file diff --git a/deal.II/examples/step-51/doc/kind b/deal.II/examples/step-51/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/deal.II/examples/step-51/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/deal.II/examples/step-51/doc/results.dox b/deal.II/examples/step-51/doc/results.dox new file mode 100644 index 0000000000..ae4bf1cc53 --- /dev/null +++ b/deal.II/examples/step-51/doc/results.dox @@ -0,0 +1,139 @@ +

Results

+ + +The program writes convergence tables to the screen while running: + + +@code +examples/\step-51> make run +============================ Running \step-51 +Solving with Q1 elements, adaptive refinement +============================================= + +Cycle 0: + Number of active cells: 4 + Number of degrees of freedom: 9 +Cycle 1: + Number of active cells: 13 + Number of degrees of freedom: 22 +Cycle 2: + Number of active cells: 31 + Number of degrees of freedom: 46 +Cycle 3: + Number of active cells: 64 + Number of degrees of freedom: 87 +Cycle 4: + Number of active cells: 127 + Number of degrees of freedom: 160 +Cycle 5: + Number of active cells: 244 + Number of degrees of freedom: 297 +Cycle 6: + Number of active cells: 466 + Number of degrees of freedom: 543 + +cycle cells dofs L2 H1 Linfty + 0 4 9 1.198e+00 2.732e+00 1.383e+00 + 1 13 22 8.795e-02 1.193e+00 1.816e-01 + 2 31 46 8.147e-02 1.167e+00 1.654e-01 + 3 64 87 7.702e-02 1.077e+00 1.310e-01 + 4 127 160 4.643e-02 7.988e-01 6.745e-02 + 5 244 297 2.470e-02 5.568e-01 3.668e-02 + 6 466 543 1.622e-02 4.107e-01 2.966e-02 + +Solving with Q1 elements, global refinement +=========================================== + +Cycle 0: + Number of active cells: 4 + Number of degrees of freedom: 9 +Cycle 1: + Number of active cells: 16 + Number of degrees of freedom: 25 +Cycle 2: + Number of active cells: 64 + Number of degrees of freedom: 81 +Cycle 3: + Number of active cells: 256 + Number of degrees of freedom: 289 +Cycle 4: + Number of active cells: 1024 + Number of degrees of freedom: 1089 +Cycle 5: + Number of active cells: 4096 + Number of degrees of freedom: 4225 +Cycle 6: + Number of active cells: 16384 + Number of degrees of freedom: 16641 + +cycle cells dofs L2 H1 Linfty + 0 4 9 1.198e+00 2.732e+00 1.383e+00 + 1 16 25 8.281e-02 1.190e+00 1.808e-01 + 2 64 81 8.142e-02 1.129e+00 1.294e-01 + 3 256 289 2.113e-02 5.828e-01 4.917e-02 + 4 1024 1089 5.319e-03 2.934e-01 1.359e-02 + 5 4096 4225 1.332e-03 1.469e-01 3.482e-03 + 6 16384 16641 3.332e-04 7.350e-02 8.758e-04 + +n cells H1 L2 + 0 4 2.732e+00 - 1.198e+00 - - + 1 16 1.190e+00 1.20 8.281e-02 14.47 3.86 + 2 64 1.129e+00 0.08 8.142e-02 1.02 0.02 + 3 256 5.828e-01 0.95 2.113e-02 3.85 1.95 + 4 1024 2.934e-01 0.99 5.319e-03 3.97 1.99 + 5 4096 1.469e-01 1.00 1.332e-03 3.99 2.00 + 6 16384 7.350e-02 1.00 3.332e-04 4.00 2.00 + +Solving with Q2 elements, global refinement +=========================================== + +Cycle 0: + Number of active cells: 4 + Number of degrees of freedom: 25 +Cycle 1: + Number of active cells: 16 + Number of degrees of freedom: 81 +Cycle 2: + Number of active cells: 64 + Number of degrees of freedom: 289 +Cycle 3: + Number of active cells: 256 + Number of degrees of freedom: 1089 +Cycle 4: + Number of active cells: 1024 + Number of degrees of freedom: 4225 +Cycle 5: + Number of active cells: 4096 + Number of degrees of freedom: 16641 +Cycle 6: + Number of active cells: 16384 + Number of degrees of freedom: 66049 + +cycle cells dofs L2 H1 Linfty + 0 4 25 1.433e+00 2.445e+00 1.286e+00 + 1 16 81 7.912e-02 1.168e+00 1.728e-01 + 2 64 289 7.755e-03 2.511e-01 1.991e-02 + 3 256 1089 9.969e-04 6.235e-02 2.764e-03 + 4 1024 4225 1.265e-04 1.571e-02 3.527e-04 + 5 4096 16641 1.587e-05 3.937e-03 4.343e-05 + 6 16384 66049 1.986e-06 9.847e-04 5.402e-06 + +n cells H1 L2 + 0 4 2.445e+00 - 1.433e+00 - - + 1 16 1.168e+00 1.07 7.912e-02 18.11 4.18 + 2 64 2.511e-01 2.22 7.755e-03 10.20 3.35 + 3 256 6.235e-02 2.01 9.969e-04 7.78 2.96 + 4 1024 1.571e-02 1.99 1.265e-04 7.88 2.98 + 5 4096 3.937e-03 2.00 1.587e-05 7.97 2.99 + 6 16384 9.847e-04 2.00 1.986e-06 7.99 3.00 +@endcode + + +One can see the error reduction upon grid refinement, and for the +cases where global refinement was performed, also the convergence +rates can be seen. The linear and quadratic convergence rates of Q1 +and Q2 elements in the $H^1$ norm can clearly be seen, as +are the quadratic and cubic rates in the $L_2$ norm. + + + diff --git a/deal.II/examples/step-51/doc/tooltip b/deal.II/examples/step-51/doc/tooltip new file mode 100644 index 0000000000..784e257e2f --- /dev/null +++ b/deal.II/examples/step-51/doc/tooltip @@ -0,0 +1 @@ +Convection-diffusion equation. Hybridizable discontinuous Galerkin methods. Face elements. diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc new file mode 100644 index 0000000000..256e48dbba --- /dev/null +++ b/deal.II/examples/step-51/step-51.cc @@ -0,0 +1,893 @@ +/* --------------------------------------------------------------------- + * $Id$ + * + * Copyright (C) 2013 - 2013 by the deal.II authors + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + + * + * Author: Martin Kronbichler, TU Muenchen, + * Scott T. Miller, xxx, 2013 + */ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace Step51 +{ + using namespace dealii; + + // @sect3{Equation data} + + // The structure of the analytic solution is the same as in step-7. There + // are two exceptions. Firstly, we also create a solution for the 3d case, + // and secondly, we take into account the convection velocity in the right + // hand side that is variable in this case. + template + class SolutionBase + { + protected: + static const unsigned int n_source_centers = 3; + static const Point source_centers[n_source_centers]; + static const double width; + }; + + + template <> + const Point<1> + SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] + = { Point<1>(-1.0 / 3.0), + Point<1>(0.0), + Point<1>(+1.0 / 3.0) + }; + + + template <> + const Point<2> + SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] + = { Point<2>(-0.5, +0.5), + Point<2>(-0.5, -0.5), + Point<2>(+0.5, -0.5) + }; + + template <> + const Point<3> + SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers] + = { Point<3>(-0.5, +0.5, 0.25), + Point<3>(-0.6, -0.5, -0.125), + Point<3>(+0.5, -0.5, 0.5) }; + + template + const double SolutionBase::width = 1./5.; + + + + template + class ConvectionVelocity : public TensorFunction<1,dim> + { + public: + ConvectionVelocity() : TensorFunction<1,dim>() {} + + virtual Tensor<1,dim> value (const Point &p) const; + }; + + + + template + Tensor<1,dim> + ConvectionVelocity::value(const Point &p) const + { + Tensor<1,dim> convection; + switch (dim) + { + case 1: + convection[0] = 1; + break; + case 2: + convection[0] = p[1]; + convection[1] = -p[0]; + break; + case 3: + convection[0] = p[1]; + convection[1] = -p[0]; + convection[2] = 1; + break; + default: + Assert(false, ExcNotImplemented()); + } + return convection; + } + + + template + class Solution : public Function, + protected SolutionBase + { + public: + Solution () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + }; + + + + template + double Solution::value (const Point &p, + const unsigned int) const + { + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + return_value += std::exp(-x_minus_xi.square() / + (this->width * this->width)); + } + + return return_value / + Utilities::fixed_power(std::sqrt(2. * numbers::PI) * this->width); + } + + + + template + Tensor<1,dim> Solution::gradient (const Point &p, + const unsigned int) const + { + Tensor<1,dim> return_value; + + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + + return_value += (-2 / (this->width * this->width) * + std::exp(-x_minus_xi.square() / + (this->width * this->width)) * + x_minus_xi); + } + + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) * + this->width); + } + + + + template + class SolutionAndGradient : public Function, + protected SolutionBase + { + public: + SolutionAndGradient () : Function(dim) {} + + virtual void vector_value (const Point &p, + Vector &v) const + { + AssertDimension(v.size(), dim+1); + Solution solution; + Tensor<1,dim> grad = solution.gradient(p); + for (unsigned int d=0; d + class RightHandSide : public Function, + protected SolutionBase + { + public: + RightHandSide () : Function() {} + + virtual double value (const Point &p, + const unsigned int component = 0) const; + + private: + const ConvectionVelocity convection_velocity; + }; + + + template + double RightHandSide::value (const Point &p, + const unsigned int) const + { + Tensor<1,dim> convection = convection_velocity.value(p); + double return_value = 0; + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; + + return_value += + ((2*dim - 2*convection*x_minus_xi - 4*x_minus_xi.square()/ + (this->width * this->width)) / + (this->width * this->width) * + std::exp(-x_minus_xi.square() / + (this->width * this->width))); + } + + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) + * this->width); + } + + + + template + class Step51 + { + public: + enum RefinementMode + { + global_refinement, adaptive_refinement + }; + + Step51 (const unsigned int degree, + const RefinementMode refinement_mode); + void run (); + + private: + void setup_system (); + void assemble_system (const bool reconstruct_trace = false); + void solve (); + void postprocess (); + void refine_mesh (); + void output_results (const unsigned int cycle); + + Triangulation triangulation; + + const MappingQ mapping; + + FESystem fe_local; + DoFHandler dof_handler_local; + + FE_FaceQ fe; + DoFHandler dof_handler; + + FE_DGQ fe_u_post; + DoFHandler dof_handler_u_post; + + ConstraintMatrix constraints; + ChunkSparsityPattern sparsity_pattern; + ChunkSparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + Vector solution_local; + Vector solution_u_post; + + const RefinementMode refinement_mode; + + ConvergenceTable convergence_table; + }; + + + + template + Step51::Step51 (const unsigned int degree, + const RefinementMode refinement_mode) : + mapping (3), + fe_local (FE_DGQ(degree), dim, + FE_DGQ(degree), 1), + dof_handler_local (triangulation), + fe (degree), + dof_handler (triangulation), + fe_u_post (degree+1), + dof_handler_u_post (triangulation), + refinement_mode (refinement_mode) + {} + + + + template + void + Step51::setup_system () + { + dof_handler_local.distribute_dofs(fe_local); + dof_handler.distribute_dofs(fe); + dof_handler_u_post.distribute_dofs(fe_u_post); + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + solution_local.reinit (dof_handler_local.n_dofs()); + solution_u_post.reinit (dof_handler_u_post.n_dofs()); + + constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + std::map boundary_values; + typename FunctionMap::type boundary_functions; + Solution solution; + boundary_functions[0] = &solution; + VectorTools::project_boundary_values (mapping, dof_handler, + boundary_functions, + QGauss(fe.degree+1), + boundary_values); + for (std::map::iterator it = boundary_values.begin(); + it != boundary_values.end(); ++it) + if (constraints.is_constrained(it->first) == false) + { + constraints.add_line(it->first); + constraints.set_inhomogeneity(it->first, it->second); + } + constraints.close (); + + { + CompressedSimpleSparsityPattern csp (dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, false); + sparsity_pattern.copy_from(csp, fe.dofs_per_face); + } + system_matrix.reinit (sparsity_pattern); + } + + + + template + void + Step51::assemble_system (const bool trace_reconstruct) + { + QGauss quadrature_formula(fe.degree+1); + QGauss face_quadrature_formula(fe.degree+1); + + FEValues fe_values_local (mapping, fe_local, quadrature_formula, + update_values | update_gradients | + update_JxW_values | update_quadrature_points); + FEFaceValues fe_face_values (mapping, fe, face_quadrature_formula, + update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + FEFaceValues fe_face_values_local (mapping, fe_local, + face_quadrature_formula, + update_values); + + const unsigned int n_q_points = quadrature_formula.size(); + const unsigned int n_face_q_points = face_quadrature_formula.size(); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int loc_dofs_per_cell = fe_local.dofs_per_cell; + + FullMatrix ll_matrix (loc_dofs_per_cell, loc_dofs_per_cell); + FullMatrix lf_matrix (loc_dofs_per_cell, dofs_per_cell); + FullMatrix fl_matrix (dofs_per_cell, loc_dofs_per_cell); + FullMatrix tmp_matrix (dofs_per_cell, loc_dofs_per_cell); + FullMatrix ff_matrix (dofs_per_cell, dofs_per_cell); + Vector l_rhs (loc_dofs_per_cell); + Vector f_rhs (dofs_per_cell); + Vector tmp_rhs (loc_dofs_per_cell); + + std::vector dof_indices (dofs_per_cell); + std::vector loc_dof_indices (loc_dofs_per_cell); + + ConvectionVelocity convection; + std::vector > convection_values (n_q_points); + std::vector > convection_values_face (n_face_q_points); + + std::vector trace_values(n_face_q_points); + + // Choose stabilization parameter to be 5 * diffusion = 5 + const double tau_stab_diffusion = 5.; + std::vector tau_stab (n_q_points); + + RightHandSide right_hand_side; + std::vector rhs_values (n_q_points); + + const Solution exact_solution; + std::vector neumann_values (n_face_q_points); + + const FEValuesExtractors::Vector gradients (0); + const FEValuesExtractors::Scalar values (dim); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + loc_cell = dof_handler_local.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell, ++loc_cell) + { + if (!trace_reconstruct) + { + lf_matrix = 0; + fl_matrix = 0; + ff_matrix = 0; + f_rhs = 0; + } + fe_values_local.reinit (loc_cell); + right_hand_side.value_list (fe_values_local.get_quadrature_points(), + rhs_values); + convection.value_list(fe_values_local.get_quadrature_points(), + convection_values); + + for (unsigned int i=0; i::faces_per_cell; ++face) + { + fe_face_values_local.reinit(loc_cell, face); + fe_face_values.reinit(cell, face); + const std::vector &JxW = fe_face_values.get_JxW_values(); + const std::vector > &normals = + fe_face_values.get_normal_vectors(); + convection.value_list(fe_face_values.get_quadrature_points(), + convection_values_face); + for (unsigned int q=0; qface(face)->at_boundary() + && + (cell->face(face)->boundary_indicator() == 1)) + { + exact_solution.value_list(fe_face_values.get_quadrature_points(), + neumann_values); + for (unsigned int i=0; iget_dof_indices(dof_indices); + constraints.distribute_local_to_global (ff_matrix, f_rhs, + dof_indices, + system_matrix, system_rhs); + } + else + { + ll_matrix.vmult(tmp_rhs, l_rhs); + loc_cell->set_dof_values(tmp_rhs, solution_local); + } + } + } + + + + template + void Step51::solve () + { + SolverControl solver_control (system_matrix.m()*10, + 1e-10*system_rhs.l2_norm()); + SolverGMRES<> solver (solver_control, 50); + solver.solve (system_matrix, solution, system_rhs, + PreconditionIdentity()); + + std::cout << " Number of GMRES iterations: " << solver_control.last_step() + << std::endl; + + system_matrix.clear(); + sparsity_pattern.reinit(0,0,0,1); + constraints.distribute(solution); + + // update local values + assemble_system(true); + } + + + + template + void + Step51::postprocess() + { + const unsigned int n_active_cells=triangulation.n_active_cells(); + Vector difference_per_cell (triangulation.n_active_cells()); + + ComponentSelectFunction value_select (dim, dim+1); + VectorTools::integrate_difference (mapping, dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(fe.degree+2), + VectorTools::L2_norm, + &value_select); + const double L2_error = difference_per_cell.l2_norm(); + + ComponentSelectFunction gradient_select (std::pair(0, dim), + dim+1); + VectorTools::integrate_difference (mapping, dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(fe.degree+2), + VectorTools::L2_norm, + &gradient_select); + const double grad_error = difference_per_cell.l2_norm(); + + convergence_table.add_value("cells", n_active_cells); + convergence_table.add_value("dofs", dof_handler.n_dofs()); + convergence_table.add_value("val L2", L2_error); + convergence_table.add_value("grad L2", grad_error); + + // construct post-processed solution with (hopefully) higher order of + // accuracy + QGauss quadrature(fe_u_post.degree+1); + FEValues fe_values(mapping, fe_u_post, quadrature, + update_values | update_JxW_values | + update_gradients); + + const unsigned int n_q_points = quadrature.size(); + std::vector u_values(n_q_points); + std::vector > u_gradients(n_q_points); + FEValuesExtractors::Vector gradients(0); + FEValuesExtractors::Scalar values(dim); + FEValues fe_values_local(mapping, fe_local, quadrature, update_values); + FullMatrix cell_matrix(fe_u_post.dofs_per_cell, + fe_u_post.dofs_per_cell); + Vector cell_rhs(fe_u_post.dofs_per_cell); + Vector cell_sol(fe_u_post.dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell_loc = dof_handler_local.begin_active(), + cell = dof_handler_u_post.begin_active(), + endc = dof_handler_u_post.end(); + for ( ; cell != endc; ++cell, ++cell_loc) + { + fe_values.reinit(cell); + fe_values_local.reinit(cell_loc); + + fe_values_local[values].get_function_values(solution_local, u_values); + fe_values_local[gradients].get_function_values(solution_local, u_gradients); + for (unsigned int i=1; idistribute_local_to_global(cell_sol, solution_u_post); + } + + VectorTools::integrate_difference (mapping, dof_handler_u_post, + solution_u_post, + Solution(), + difference_per_cell, + QGauss(fe.degree+3), + VectorTools::L2_norm); + double post_error = difference_per_cell.l2_norm(); + convergence_table.add_value("val L2-post", post_error); + } + + + + template + void Step51::output_results (const unsigned int cycle) + { + std::string filename; + switch (refinement_mode) + { + case global_refinement: + filename = "solution-global"; + break; + case adaptive_refinement: + filename = "solution-adaptive"; + break; + default: + Assert (false, ExcNotImplemented()); + } + + filename += "-q" + Utilities::int_to_string(fe.degree,1); + filename += "-" + Utilities::int_to_string(cycle,2); + filename += ".vtk"; + std::ofstream output (filename.c_str()); + + DataOut data_out; + std::vector names (dim, "gradient"); + names.push_back ("solution"); + std::vector + component_interpretation + (dim+1, DataComponentInterpretation::component_is_part_of_vector); + component_interpretation[dim] + = DataComponentInterpretation::component_is_scalar; + data_out.add_data_vector (dof_handler_local, solution_local, + names, component_interpretation); + + data_out.build_patches (fe.degree); + data_out.write_vtk (output); + } + + + + + template + void Step51::run () + { + const bool do_cube = true; + if (!do_cube) + { + GridGenerator::hyper_ball (triangulation); + static const HyperBallBoundary boundary; + triangulation.set_boundary(0, boundary); + triangulation.refine_global(6-2*dim); + } + + for (unsigned int cycle=0; cycle<10; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + + if (do_cube) + { + triangulation.clear(); + GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); + triangulation.refine_global(3-dim+cycle/2); + } + else triangulation.refine_global(1); + + setup_system (); + assemble_system (false); + solve (); + postprocess(); + output_results (cycle); + } + + + + convergence_table.set_precision("val L2", 3); + convergence_table.set_scientific("val L2", true); + convergence_table.set_precision("grad L2", 3); + convergence_table.set_scientific("grad L2", true); + convergence_table.set_precision("val L2-post", 3); + convergence_table.set_scientific("val L2-post", true); + + convergence_table + .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim); + convergence_table + .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim); + convergence_table + .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim); + convergence_table.write_text(std::cout); + } +} + + +int main (int argc, char** argv) +{ + const unsigned int dim = 2; + + try + { + using namespace dealii; + using namespace Step51; + + deallog.depth_console (0); + + // Now for the three calls to the main class in complete analogy to + // step-7. + { + std::cout << "Solving with Q1 elements, adaptive refinement" << std::endl + << "=============================================" << std::endl + << std::endl; + + Step51 hdg_problem (1, Step51::adaptive_refinement); + hdg_problem.run (); + + std::cout << std::endl; + } + + { + std::cout << "Solving with Q1 elements, global refinement" << std::endl + << "===========================================" << std::endl + << std::endl; + + Step51 hdg_problem (1, Step51::global_refinement); + hdg_problem.run (); + + std::cout << std::endl; + } + + { + std::cout << "Solving with Q3 elements, global refinement" << std::endl + << "===========================================" << std::endl + << std::endl; + + Step51 hdg_problem (3, Step51::global_refinement); + hdg_problem.run (); + + std::cout << std::endl; + } + + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} -- 2.39.5