From 4c443c0feb1c68650f6f9820bf372ff0fb64477b Mon Sep 17 00:00:00 2001 From: Scott Miller Date: Thu, 8 Aug 2013 02:07:07 +0000 Subject: [PATCH] Add output for step-51 git-svn-id: https://svn.dealii.org/trunk@30254 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-51/step-51.cc | 1427 ++++++++++++++------------- 1 file changed, 727 insertions(+), 700 deletions(-) diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index d67bbfb248..ab4bf33dae 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -16,7 +16,7 @@ * * Author: Martin Kronbichler, TU Muenchen, - * Scott T. Miller, xxx, 2013 + * Scott T. Miller, The Pennsylvania State University, 2013 */ @@ -32,829 +32,855 @@ #include #include #include +#include #include #include -#include #include #include -#include #include #include #include #include #include #include -#include #include #include +//! New include: output data on faces of a +// triangulation +#include + + +//! New include: fe_face.h +// Explain that it implements fe on +// codim=1 surfaces of a geometric discretization #include #include #include #include -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 +namespace Step51 { -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) }; + using namespace dealii; -template -const double SolutionBase::width = 1./5.; + // @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 -class ConvectionVelocity : public TensorFunction<1,dim> -{ -public: - ConvectionVelocity() : TensorFunction<1,dim>() {} + 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) + }; - virtual Tensor<1,dim> value (const Point &p) const; -}; + 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 -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 + const double SolutionBase::width = 1./5.; -template -class Solution : public Function, - protected SolutionBase -{ -public: - Solution () : Function() {} - virtual double value (const Point &p, - const unsigned int component = 0) const; + template + class ConvectionVelocity : public TensorFunction<1,dim> + { + public: + ConvectionVelocity() : TensorFunction<1,dim>() {} - virtual Tensor<1,dim> gradient (const Point &p, - const unsigned int component = 0) const; -}; + virtual Tensor<1,dim> value (const Point &p) 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)); - } + 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; + } - return return_value / - Utilities::fixed_power(std::sqrt(2. * numbers::PI) * this->width); -} + template + class Solution : public Function, + protected SolutionBase + { + public: + Solution () : Function() {} + virtual double value (const Point &p, + const unsigned int component = 0) const; -template -Tensor<1,dim> Solution::gradient (const Point &p, - const unsigned int) const -{ - Tensor<1,dim> return_value; + virtual Tensor<1,dim> gradient (const Point &p, + const unsigned int component = 0) const; + }; - 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 + 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 -class SolutionAndGradient : public Function, - protected SolutionBase -{ -public: - SolutionAndGradient () : Function(dim) {} - virtual void vector_value (const Point &p, - Vector &v) const + template + Tensor<1,dim> Solution::gradient (const Point &p, + const unsigned int) const { - AssertDimension(v.size(), dim+1); - Solution solution; - Tensor<1,dim> grad = solution.gradient(p); - for (unsigned int d=0; d 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); + } -template -class RightHandSide : public Function, - protected SolutionBase -{ -public: - RightHandSide () : Function() {} + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) * + this->width); + } - virtual double value (const Point &p, - const unsigned int component = 0) const; -private: - const ConvectionVelocity convection_velocity; -}; + template + class SolutionAndGradient : public Function, + protected SolutionBase + { + public: + SolutionAndGradient () : Function(dim) {} -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) + virtual void vector_value (const Point &p, + Vector &v) const { - 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))); + AssertDimension(v.size(), dim+1); + Solution solution; + Tensor<1,dim> grad = solution.gradient(p); + for (unsigned int d=0; d(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_grid (const unsigned int cylce); - void output_results (const unsigned int cycle); - Triangulation triangulation; - - const MappingQ mapping; - - FESystem fe_local; - DoFHandler dof_handler_local; + template + class RightHandSide : public Function, + protected SolutionBase + { + public: + RightHandSide () : Function() {} - FE_FaceQ fe; - DoFHandler dof_handler; + virtual double value (const Point &p, + const unsigned int component = 0) const; - FE_DGQ fe_u_post; - DoFHandler dof_handler_u_post; + private: + const ConvectionVelocity convection_velocity; + }; - ConstraintMatrix constraints; - ChunkSparsityPattern sparsity_pattern; - ChunkSparseMatrix system_matrix; - Vector solution; - Vector system_rhs; + 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))); + } - Vector solution_local; - Vector solution_u_post; + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) + * this->width); + } - const RefinementMode refinement_mode; - ConvergenceTable convergence_table; -}; + template + class HDG + { + public: + enum RefinementMode + { + global_refinement, adaptive_refinement + }; + HDG (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; + + // local (element interior) solutions + FESystem fe_local; + DoFHandler dof_handler_local; + + // global (trace/skeleton) solution + // Note that FE_FaceQ represents + // finite element data on the faces/edges + // of our triangulation + FE_FaceQ fe; + DoFHandler dof_handler; + + // post-processed solution + FE_DGQ fe_u_post; + DoFHandler dof_handler_u_post; + + // Dirichlet BCs are strongly enforced + // on the "skeleton" solution + ConstraintMatrix constraints; + + + // Comment on "chunk" here. + // First, set up objects for the global + // solution + ChunkSparsityPattern sparsity_pattern; + ChunkSparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + + // Local solution values + Vector solution_local; + + // HDG solutions can be post-processed + // to gain one order of accuracy. + // solution_u_post will be + // our post-processed DG solution on the + // interior of cells represented by a + // DG solution of order (degree+1) + Vector solution_u_post; + + // Same as step-7: + const RefinementMode refinement_mode; + + ConvergenceTable convergence_table; + }; + + + + template + HDG::HDG (const unsigned int degree, + const RefinementMode refinement_mode) : + mapping (1), + 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 + HDG::setup_system () + { + dof_handler_local.distribute_dofs(fe_local); + dof_handler.distribute_dofs(fe); + dof_handler_u_post.distribute_dofs(fe_u_post); -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) -{} + 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()); -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); - typename FunctionMap::type boundary_functions; - Solution solution; - boundary_functions[0] = &solution; - VectorTools::project_boundary_values (mapping, dof_handler, - boundary_functions, - QGauss(fe.degree+1), - constraints); - constraints.close (); + 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); + { + 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); } - 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); - - std::vector > q_phi (loc_dofs_per_cell); - std::vector q_phi_div (loc_dofs_per_cell); - std::vector u_phi (loc_dofs_per_cell); - std::vector > u_phi_grad (loc_dofs_per_cell); - std::vector tr_phi (dofs_per_cell); - - std::vector trace_values(n_face_q_points); - - // Choose stabilization parameter to be 5 * diffusion = 5 - const double tau_stab_diffusion = 5.; - - ConvectionVelocity convection_velocity; - RightHandSide right_hand_side; - const Solution exact_solution; - - const FEValuesExtractors::Vector fluxes (0); - const FEValuesExtractors::Scalar scalar (dim); - - std::vector > - fe_local_support_on_face(GeometryInfo::faces_per_cell); - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int i=0; i > - fe_support_on_face(GeometryInfo::faces_per_cell); - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int i=0; i::active_cell_iterator + template + void + HDG::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) - { - ll_matrix = 0; - l_rhs = 0; - if (!trace_reconstruct) - { - lf_matrix = 0; - fl_matrix = 0; - ff_matrix = 0; - f_rhs = 0; - } - fe_values_local.reinit (loc_cell); - - for (unsigned int q=0; q convection - = convection_velocity.value(fe_values_local.quadrature_point(q)); - const double JxW = fe_values_local.JxW(q); - for (unsigned int k=0; k::faces_per_cell; ++face) - { - fe_face_values_local.reinit(loc_cell, face); - fe_face_values.reinit(cell, face); - if (trace_reconstruct) - fe_face_values.get_function_values (solution, trace_values); - - for (unsigned int q=0; q normal = fe_face_values.normal_vector(q); - const Tensor<1,dim> convection - = convection_velocity.value(fe_face_values.quadrature_point(q)); - const double tau_stab = (tau_stab_diffusion + - std::abs(convection * normal)); - - for (unsigned int k=0; kface(face)->at_boundary() - && - (cell->face(face)->boundary_indicator() == 1)) + 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)) { - const unsigned int ii=fe_local_support_on_face[face][i]; - const unsigned int jj=fe_local_support_on_face[face][j]; - ll_matrix(ii,jj) += tau_stab * u_phi[i] * u_phi[j] * JxW; + 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); - } - } -} + } + } + + // invert ll_matrix and overwrite + ll_matrix.gauss_jordan(); + if (!trace_reconstruct) + { + // tmp_matrix = fl_matrix * ll_matrix + fl_matrix.mmult(tmp_matrix, ll_matrix); + + // f_rhs = tmp_matrix * l_rhs + tmp_matrix.vmult_add(f_rhs, l_rhs); + + // ff_matrix = ff_matrix + tmp_matrix * lf_matrix + tmp_matrix.mmult(ff_matrix, lf_matrix, true); + cell->get_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()); + template + void HDG::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; + 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); + system_matrix.clear(); + sparsity_pattern.reinit(0,0,0,1); + constraints.distribute(solution); - // update local values - assemble_system(true); -} + // 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 fluxes(0); - FEValuesExtractors::Scalar scalar(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[scalar].get_function_values(solution_local, u_values); - fe_values_local[fluxes].get_function_values(solution_local, u_gradients); - for (unsigned int i=1; i + void + HDG::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); -} + cell_matrix.gauss_jordan(); + cell_matrix.vmult(cell_sol, cell_rhs); + cell->distribute_local_to_global(cell_sol, solution_u_post); + } - -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); -} + 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::refine_grid (const unsigned int cycle) -{ - const bool do_cube = true; - if (cycle == 0) - { - if (!do_cube) - { - GridGenerator::hyper_ball (triangulation); - static const HyperBallBoundary boundary; - triangulation.set_boundary(0, boundary); - triangulation.refine_global(6-2*dim); - } - else - GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1); - } - else + template + void HDG::output_results (const unsigned int cycle) + { + std::string filename; switch (refinement_mode) { case global_refinement: - { - 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); - break; - } - + filename = "solution-global"; + break; case adaptive_refinement: - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - FEValuesExtractors::Scalar scalar(dim); - typename FunctionMap::type neumann_boundary; - KellyErrorEstimator::estimate (dof_handler_local, - QGauss(3), - neumann_boundary, - solution_local, - estimated_error_per_cell, - fe_local.component_mask(scalar)); - - GridRefinement::refine_and_coarsen_fixed_number (triangulation, - estimated_error_per_cell, - 0.3, 0.); - - triangulation.execute_coarsening_and_refinement (); - + filename = "solution-adaptive"; break; - } - default: - { Assert (false, ExcNotImplemented()); } - } + + std::string face_out(filename); + face_out += "-face"; + + 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); + + // Post-processed solution + std::vector post_name(1,"u_post"); + std::vector + post_comp_type(1, DataComponentInterpretation::component_is_scalar); + data_out.add_data_vector (dof_handler_u_post, solution_u_post, + post_name, post_comp_type); + + // build patches based on the highest degree, i.e. the post-proc'd soln + data_out.build_patches (fe_u_post.degree); + data_out.write_vtk (output); + + face_out += "-q" + Utilities::int_to_string(fe.degree,1); + face_out += "-" + Utilities::int_to_string(cycle,2); + face_out += ".vtk"; + std::ofstream face_output (face_out.c_str()); + + DataOutFaces data_out_face(false); + std::vector face_name(1,"lambda"); + std::vector + face_component_type(1, DataComponentInterpretation::component_is_scalar); + + data_out_face.add_data_vector (dof_handler, + solution, + face_name, + face_component_type); + + data_out_face.build_patches (fe.degree); + data_out_face.write_vtk (face_output); + } + template + void HDG::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); + } -template -void Step51::run () -{ - for (unsigned int cycle=0; cycle<10; ++cycle) - { - std::cout << "Cycle " << cycle << ':' << std::endl; - - refine_grid (cycle); - setup_system (); - assemble_system (false); - solve (); - postprocess(); - output_results (cycle); - } + 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.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); + 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); + } } @@ -865,6 +891,7 @@ int main (int argc, char** argv) try { using namespace dealii; + using namespace Step51; deallog.depth_console (0); @@ -875,7 +902,7 @@ int main (int argc, char** argv) << "=============================================" << std::endl << std::endl; - Step51 hdg_problem (1, Step51::adaptive_refinement); + HDG hdg_problem (1, HDG::adaptive_refinement); hdg_problem.run (); std::cout << std::endl; @@ -886,7 +913,7 @@ int main (int argc, char** argv) << "===========================================" << std::endl << std::endl; - Step51 hdg_problem (1, Step51::global_refinement); + HDG hdg_problem (1, HDG::global_refinement); hdg_problem.run (); std::cout << std::endl; @@ -897,7 +924,7 @@ int main (int argc, char** argv) << "===========================================" << std::endl << std::endl; - Step51 hdg_problem (3, Step51::global_refinement); + HDG hdg_problem (3, HDG::global_refinement); hdg_problem.run (); std::cout << std::endl; -- 2.39.5