From 9999086d6465d45e05cb4e506183452259741d0e Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Mon, 16 Sep 2013 20:08:25 +0000 Subject: [PATCH] Reindent file with astyle git-svn-id: https://svn.dealii.org/trunk@30743 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-51/doc/intro.dox | 20 +- deal.II/examples/step-51/step-51.cc | 1928 ++++++++++++------------ 2 files changed, 975 insertions(+), 973 deletions(-) diff --git a/deal.II/examples/step-51/doc/intro.dox b/deal.II/examples/step-51/doc/intro.dox index bdf0b8039d..dacb532c60 100644 --- a/deal.II/examples/step-51/doc/intro.dox +++ b/deal.II/examples/step-51/doc/intro.dox @@ -276,13 +276,13 @@ For this tutorial program, we consider almost the same test case as in step-7. The computational domain is $\Omega := [-1,1]^d$ and the exact solution corresponds to the one in step-7, except for a scaling. We use the following source centers xi for the exponentials @@ -298,13 +298,13 @@ equal to one and the convection as (y, -x, 1), & \textrm{dim}=3 \end{cases} \f] - +Note that the convection is divergence-free, $\nabla \cdot c = 0$.

Implementation

-Implementation notes: +Besides implementing the above equations, the implementation below provides the following features: -- WorkStream to parallelize local solvers. Workstream is already used in step-32, step-44. -- Reconstructing the trace +- WorkStream to parallelize local solvers. Workstream is already used in step-32, step-44. +- Reconstruct the local DG solution from the trace trace - Post-processing the solution for superconvergence -- DataOutFaces: direct output of the global solution +- DataOutFaces for direct output of the global skeleton solution diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index 2f578d4619..6d5399ce30 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -80,10 +80,12 @@ // the simulation. #include + +// We start by putting the class into its own namespace. namespace Step51 { -using namespace dealii; + using namespace dealii; // @sect3{Equation data} // @@ -91,98 +93,98 @@ using namespace dealii; // two exceptions. Firstly, we also create a solution for the 3d case, and // secondly, we scale the solution so its norm is of order unity for all // values of the solution width. -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 + 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<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<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<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 <> + 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 -class Solution : public Function, - protected SolutionBase -{ -public: - Solution () : Function() {} + template + const double SolutionBase::width = 1./5.; - 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 + 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 + 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]; + template + Tensor<1,dim> Solution::gradient (const Point &p, + const unsigned int) const + { + Tensor<1,dim> return_value; - return_value += (-2 / (this->width * this->width) * - std::exp(-x_minus_xi.square() / - (this->width * this->width)) * - x_minus_xi); - } + for (unsigned int i=0; in_source_centers; ++i) + { + const Point x_minus_xi = p - this->source_centers[i]; - return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) * - this->width); -} + 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); + } @@ -190,69 +192,69 @@ Tensor<1,dim> Solution::gradient (const Point &p, // gradient are collected together. This function is used when computing the // error of the HDG approximation and its implementation is to simply call // value and gradient function of the Solution class. -template -class SolutionAndGradient : public Function, - protected SolutionBase -{ -public: - SolutionAndGradient () : Function(dim) {} + template + class SolutionAndGradient : public Function, + protected SolutionBase + { + public: + SolutionAndGradient () : Function(dim) {} - virtual void vector_value (const Point &p, - Vector &v) const; -}; + virtual void vector_value (const Point &p, + Vector &v) const; + }; -template -void SolutionAndGradient::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 + void SolutionAndGradient::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 ConvectionVelocity : public TensorFunction<1,dim> -{ -public: - ConvectionVelocity() : TensorFunction<1,dim>() {} + template + class ConvectionVelocity : public TensorFunction<1,dim> + { + public: + ConvectionVelocity() : TensorFunction<1,dim>() {} - virtual Tensor<1,dim> value (const Point &p) const; -}; + 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 + 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; + } @@ -261,42 +263,42 @@ ConvectionVelocity::value(const Point &p) const // a convection term instead of the reaction term. Since the velocity field is // incompressible, i.e. $\nabla \cdot \mathbf{c} = 0$, this term simply reads // $\mathbf{c} \nabla \ve u$. -template -class RightHandSide : public Function, - protected SolutionBase -{ -public: - RightHandSide () : Function() {} + template + class RightHandSide : public Function, + protected SolutionBase + { + public: + RightHandSide () : Function() {} - virtual double value (const Point &p, - const unsigned int component = 0) const; + virtual double value (const Point &p, + const unsigned int component = 0) const; -private: - const ConvectionVelocity convection_velocity; -}; + 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))); - } + 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); -} + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) + * this->width); + } // @sect3{The HDG solver class} @@ -311,99 +313,99 @@ double RightHandSide::value (const Point &p, // generate the system matrix and once when we compute the element-interior // solutions from the skeleton values) and for the postprocessing where // we extract a solution that converges at higher order. -template -class HDG -{ -public: - enum RefinementMode + template + class HDG { - global_refinement, adaptive_refinement - }; + public: + enum RefinementMode + { + global_refinement, adaptive_refinement + }; - HDG (const unsigned int degree, - const RefinementMode refinement_mode); - void run (); + HDG (const unsigned int degree, + const RefinementMode refinement_mode); + void run (); -private: + private: // Data for the assembly and solution of the primal variables. - struct PerTaskData; - struct ScratchData; + struct PerTaskData; + struct ScratchData; // Post-processing the solution to obtain $u^*$ is an element-by-element // procedure; as such, we do not need to assemble any global data and do // not declare any 'task data' for WorkStream to use. - struct PostProcessScratchData; + struct PostProcessScratchData; - void setup_system (); - void assemble_system (const bool reconstruct_trace = false); - void solve (); - void postprocess (); + 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); + void refine_grid (const unsigned int cylce); + void output_results (const unsigned int cycle); // The following three functions are used by WorkStream to do the actual work of // the program. - void assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, - ScratchData &scratch, - PerTaskData &task_data); - - void copy_local_to_global(const PerTaskData &data); - - void postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, - PostProcessScratchData &scratch, - unsigned int &empty_data); - - - Triangulation triangulation; - - // The 'local' solutions are interior to each element. These - // represent the primal solution field $u$ as well as the auxiliary - // field $\mathbf{q} = -\nabla u$. - FESystem fe_local; - DoFHandler dof_handler_local; - Vector solution_local; - - // The new finite element type and corresponding DoFHandler are - // used for the global skeleton solution that couples the element-level local - // solutions. - FE_FaceQ fe; - DoFHandler dof_handler; - Vector solution; - Vector system_rhs; - - // As stated in the introduction, HDG solutions can be post-processed to - // attain superconvergence rates of $\mathcal{O}(h^{p+2})$. The - // post-processed solution is a discontinuous finite element solution - // representing the primal variable on the interior of each cell. We define - // a FE type of degree $p+1$ to represent this post-processed solution, - // which we only use for output after constructing it. - FE_DGQ fe_u_post; - DoFHandler dof_handler_u_post; - Vector solution_u_post; - - // The degrees of freedom corresponding to the skeleton strongly enforce - // Dirichlet boundary conditions, just as in a continuous Galerkin finite - // element method. We can enforce the boundary conditions in an analogous - // manner through the use of ConstrainMatrix constructs. In - // addition, hanging nodes where cells of different refinement levels meet - // are set as for continuous finite elements: For the face elements which - // only define degrees of freedom on the face, this process sets the - // solution on the refined to be the one from the coarse side. - ConstraintMatrix constraints; - - // The usage of the ChunkSparseMatrix class is similar to the usual sparse - // matrices: You need a sparsity pattern of type ChunkSparsityPattern and - // the actual matrix object. When creating the sparsity pattern, we just - // have to additionally pass the size of local blocks. - ChunkSparsityPattern sparsity_pattern; - ChunkSparseMatrix system_matrix; - - // Same as step-7: - const RefinementMode refinement_mode; - ConvergenceTable convergence_table; -}; + void assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data); + + void copy_local_to_global(const PerTaskData &data); + + void postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + unsigned int &empty_data); + + + Triangulation triangulation; + + // The 'local' solutions are interior to each element. These + // represent the primal solution field $u$ as well as the auxiliary + // field $\mathbf{q} = -\nabla u$. + FESystem fe_local; + DoFHandler dof_handler_local; + Vector solution_local; + + // The new finite element type and corresponding DoFHandler are + // used for the global skeleton solution that couples the element-level local + // solutions. + FE_FaceQ fe; + DoFHandler dof_handler; + Vector solution; + Vector system_rhs; + + // As stated in the introduction, HDG solutions can be post-processed to + // attain superconvergence rates of $\mathcal{O}(h^{p+2})$. The + // post-processed solution is a discontinuous finite element solution + // representing the primal variable on the interior of each cell. We define + // a FE type of degree $p+1$ to represent this post-processed solution, + // which we only use for output after constructing it. + FE_DGQ fe_u_post; + DoFHandler dof_handler_u_post; + Vector solution_u_post; + + // The degrees of freedom corresponding to the skeleton strongly enforce + // Dirichlet boundary conditions, just as in a continuous Galerkin finite + // element method. We can enforce the boundary conditions in an analogous + // manner through the use of ConstrainMatrix constructs. In + // addition, hanging nodes where cells of different refinement levels meet + // are set as for continuous finite elements: For the face elements which + // only define degrees of freedom on the face, this process sets the + // solution on the refined to be the one from the coarse side. + ConstraintMatrix constraints; + + // The usage of the ChunkSparseMatrix class is similar to the usual sparse + // matrices: You need a sparsity pattern of type ChunkSparsityPattern and + // the actual matrix object. When creating the sparsity pattern, we just + // have to additionally pass the size of local blocks. + ChunkSparsityPattern sparsity_pattern; + ChunkSparseMatrix system_matrix; + + // Same as step-7: + const RefinementMode refinement_mode; + ConvergenceTable convergence_table; + }; // @sect3{The HDG class implementation} @@ -412,18 +414,18 @@ private: // FiniteElement objects. Note that we create a system of finite // elements for the local DG part, including the gradient/flux part and the // scalar part. -template -HDG::HDG (const unsigned int degree, - const RefinementMode refinement_mode) : - 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 + HDG::HDG (const unsigned int degree, + const RefinementMode refinement_mode) : + 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) + {} @@ -432,43 +434,43 @@ HDG::HDG (const unsigned int degree, // of the other tutorial programs. We are careful to distribute dofs with // all of our DoFHandler objects. The @p solution and @p system_matrix // objects go with the global skeleton solution. -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); - - 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_function; - boundary_functions[0] = &solution_function; - VectorTools::project_boundary_values (dof_handler, - boundary_functions, - QGauss(fe.degree+1), - constraints); - constraints.close (); - + template + void + HDG::setup_system () { - CompressedSimpleSparsityPattern csp (dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern (dof_handler, csp, - constraints, false); - sparsity_pattern.copy_from(csp, fe.dofs_per_face); + 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_function; + boundary_functions[0] = &solution_function; + VectorTools::project_boundary_values (dof_handler, + boundary_functions, + QGauss(fe.degree+1), + constraints); + 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); } - system_matrix.reinit (sparsity_pattern); -} @@ -486,22 +488,22 @@ HDG::setup_system () // the same function but only switch between the two based on a flag that we // set when starting the assembly. Since we need to pass this information on // to the local worker routines, we store it once in the task data. -template -struct HDG::PerTaskData -{ - FullMatrix cell_matrix; - Vector cell_vector; - std::vector dof_indices; - - bool trace_reconstruct; - - PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct) - : cell_matrix(n_dofs, n_dofs), - cell_vector(n_dofs), - dof_indices(n_dofs), - trace_reconstruct(trace_reconstruct) - {} -}; + template + struct HDG::PerTaskData + { + FullMatrix cell_matrix; + Vector cell_vector; + std::vector dof_indices; + + bool trace_reconstruct; + + PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct) + : cell_matrix(n_dofs, n_dofs), + cell_vector(n_dofs), + dof_indices(n_dofs), + trace_reconstruct(trace_reconstruct) + {} + }; @@ -517,102 +519,102 @@ struct HDG::PerTaskData // constructor and store it once for all cells that we work on. Had we not // stored this information, we would be forced to assemble a large number of // zero terms on each cell, which would significantly slow the program. -template -struct HDG::ScratchData -{ - FEValues fe_values_local; - FEFaceValues fe_face_values_local; - FEFaceValues fe_face_values; - - FullMatrix ll_matrix; - FullMatrix lf_matrix; - FullMatrix fl_matrix; - FullMatrix tmp_matrix; - Vector l_rhs; - Vector tmp_rhs; - - std::vector > q_phi; - std::vector q_phi_div; - std::vector u_phi; - std::vector > u_phi_grad; - std::vector tr_phi; - std::vector trace_values; - - std::vector > fe_local_support_on_face; - std::vector > fe_support_on_face; - - ConvectionVelocity convection_velocity; - RightHandSide right_hand_side; - const Solution exact_solution; - - ScratchData(const FiniteElement &fe, - const FiniteElement &fe_local, - const QGauss &quadrature_formula, - const QGauss &face_quadrature_formula, - const UpdateFlags local_flags, - const UpdateFlags local_face_flags, - const UpdateFlags flags) - : - fe_values_local (fe_local, quadrature_formula, local_flags), - fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags), - fe_face_values (fe, face_quadrature_formula, flags), - ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell), - lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell), - fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), - tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), - l_rhs (fe_local.dofs_per_cell), - tmp_rhs (fe_local.dofs_per_cell), - q_phi (fe_local.dofs_per_cell), - q_phi_div (fe_local.dofs_per_cell), - u_phi (fe_local.dofs_per_cell), - u_phi_grad (fe_local.dofs_per_cell), - tr_phi (fe.dofs_per_cell), - trace_values(face_quadrature_formula.size()), - fe_local_support_on_face(GeometryInfo::faces_per_cell), - fe_support_on_face(GeometryInfo::faces_per_cell) + template + struct HDG::ScratchData { - for (unsigned int face=0; face::faces_per_cell; ++face) - for (unsigned int i=0; i::faces_per_cell; ++face) - for (unsigned int i=0; i fe_values_local; + FEFaceValues fe_face_values_local; + FEFaceValues fe_face_values; + + FullMatrix ll_matrix; + FullMatrix lf_matrix; + FullMatrix fl_matrix; + FullMatrix tmp_matrix; + Vector l_rhs; + Vector tmp_rhs; + + std::vector > q_phi; + std::vector q_phi_div; + std::vector u_phi; + std::vector > u_phi_grad; + std::vector tr_phi; + std::vector trace_values; + + std::vector > fe_local_support_on_face; + std::vector > fe_support_on_face; + + ConvectionVelocity convection_velocity; + RightHandSide right_hand_side; + const Solution exact_solution; + + ScratchData(const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &quadrature_formula, + const QGauss &face_quadrature_formula, + const UpdateFlags local_flags, + const UpdateFlags local_face_flags, + const UpdateFlags flags) + : + fe_values_local (fe_local, quadrature_formula, local_flags), + fe_face_values_local (fe_local, face_quadrature_formula, local_face_flags), + fe_face_values (fe, face_quadrature_formula, flags), + ll_matrix (fe_local.dofs_per_cell, fe_local.dofs_per_cell), + lf_matrix (fe_local.dofs_per_cell, fe.dofs_per_cell), + fl_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), + tmp_matrix (fe.dofs_per_cell, fe_local.dofs_per_cell), + l_rhs (fe_local.dofs_per_cell), + tmp_rhs (fe_local.dofs_per_cell), + q_phi (fe_local.dofs_per_cell), + q_phi_div (fe_local.dofs_per_cell), + u_phi (fe_local.dofs_per_cell), + u_phi_grad (fe_local.dofs_per_cell), + tr_phi (fe.dofs_per_cell), + trace_values(face_quadrature_formula.size()), + fe_local_support_on_face(GeometryInfo::faces_per_cell), + fe_support_on_face(GeometryInfo::faces_per_cell) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + for (unsigned int i=0; i::faces_per_cell; ++face) + for (unsigned int i=0; i::ScratchData // @p PostProcessScratchData contains the data used by WorkStream // when post-processing the local solution $u^*$. It is similar, but much // simpler, than @p ScratchData. -template -struct HDG::PostProcessScratchData -{ - FEValues fe_values_local; - FEValues fe_values; - - std::vector u_values; - std::vector > u_gradients; - FullMatrix cell_matrix; - - Vector cell_rhs; - Vector cell_sol; - - PostProcessScratchData(const FiniteElement &fe, - const FiniteElement &fe_local, - const QGauss &quadrature_formula, - const UpdateFlags local_flags, - const UpdateFlags flags) - : - fe_values_local (fe_local, quadrature_formula, local_flags), - fe_values (fe, quadrature_formula, flags), - u_values (quadrature_formula.size()), - u_gradients (quadrature_formula.size()), - cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell), - cell_rhs (fe.dofs_per_cell), - cell_sol (fe.dofs_per_cell) - {} - - PostProcessScratchData(const PostProcessScratchData &sd) - : - fe_values_local (sd.fe_values_local.get_fe(), - sd.fe_values_local.get_quadrature(), - sd.fe_values_local.get_update_flags()), - fe_values (sd.fe_values.get_fe(), - sd.fe_values.get_quadrature(), - sd.fe_values.get_update_flags()), - u_values (sd.u_values), - u_gradients (sd.u_gradients), - cell_matrix (sd.cell_matrix), - cell_rhs (sd.cell_rhs), - cell_sol (sd.cell_sol) - {} -}; + template + struct HDG::PostProcessScratchData + { + FEValues fe_values_local; + FEValues fe_values; + + std::vector u_values; + std::vector > u_gradients; + FullMatrix cell_matrix; + + Vector cell_rhs; + Vector cell_sol; + + PostProcessScratchData(const FiniteElement &fe, + const FiniteElement &fe_local, + const QGauss &quadrature_formula, + const UpdateFlags local_flags, + const UpdateFlags flags) + : + fe_values_local (fe_local, quadrature_formula, local_flags), + fe_values (fe, quadrature_formula, flags), + u_values (quadrature_formula.size()), + u_gradients (quadrature_formula.size()), + cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell), + cell_rhs (fe.dofs_per_cell), + cell_sol (fe.dofs_per_cell) + {} + + PostProcessScratchData(const PostProcessScratchData &sd) + : + fe_values_local (sd.fe_values_local.get_fe(), + sd.fe_values_local.get_quadrature(), + sd.fe_values_local.get_update_flags()), + fe_values (sd.fe_values.get_fe(), + sd.fe_values.get_quadrature(), + sd.fe_values.get_update_flags()), + u_values (sd.u_values), + u_gradients (sd.u_gradients), + cell_matrix (sd.cell_matrix), + cell_rhs (sd.cell_rhs), + cell_sol (sd.cell_sol) + {} + }; // @sect4{HDG::copy_local_to_global} // If we are in the first step of the solution, i.e. @p trace_reconstruct=false, // then we assemble the global system. -template -void HDG::copy_local_to_global(const PerTaskData &data) -{ - if (data.trace_reconstruct == false) - constraints.distribute_local_to_global (data.cell_matrix, - data.cell_vector, - data.dof_indices, - system_matrix, system_rhs); -} + template + void HDG::copy_local_to_global(const PerTaskData &data) + { + if (data.trace_reconstruct == false) + constraints.distribute_local_to_global (data.cell_matrix, + data.cell_vector, + data.dof_indices, + system_matrix, system_rhs); + } @@ -688,39 +690,39 @@ void HDG::copy_local_to_global(const PerTaskData &data) // The @p trace_reconstruct input parameter is used to decide whether we are // solving for the local solution (true) or the global skeleton solution // (false). -template -void -HDG::assemble_system (const bool trace_reconstruct) -{ - const QGauss quadrature_formula(fe.degree+1); - const QGauss face_quadrature_formula(fe.degree+1); - - const UpdateFlags local_flags (update_values | update_gradients | - update_JxW_values | update_quadrature_points); - - const UpdateFlags local_face_flags (update_values); - - const UpdateFlags flags ( update_values | update_normal_vectors | - update_quadrature_points | - update_JxW_values); - - PerTaskData task_data (fe.dofs_per_cell, - trace_reconstruct); - ScratchData scratch (fe, fe_local, - quadrature_formula, - face_quadrature_formula, - local_flags, - local_face_flags, - flags); - - WorkStream::run(dof_handler.begin_active(), - dof_handler.end(), - *this, - &HDG::assemble_system_one_cell, - &HDG::copy_local_to_global, - scratch, - task_data); -} + template + void + HDG::assemble_system (const bool trace_reconstruct) + { + const QGauss quadrature_formula(fe.degree+1); + const QGauss face_quadrature_formula(fe.degree+1); + + const UpdateFlags local_flags (update_values | update_gradients | + update_JxW_values | update_quadrature_points); + + const UpdateFlags local_face_flags (update_values); + + const UpdateFlags flags ( update_values | update_normal_vectors | + update_quadrature_points | + update_JxW_values); + + PerTaskData task_data (fe.dofs_per_cell, + trace_reconstruct); + ScratchData scratch (fe, fe_local, + quadrature_formula, + face_quadrature_formula, + local_flags, + local_face_flags, + flags); + + WorkStream::run(dof_handler.begin_active(), + dof_handler.end(), + *this, + &HDG::assemble_system_one_cell, + &HDG::copy_local_to_global, + scratch, + task_data); + } @@ -728,258 +730,258 @@ HDG::assemble_system (const bool trace_reconstruct) // The real work of the HDG program is done by @p assemble_system_one_cell. // Assembling the local matrices $A, B, C$ is done here, along with the // local contributions of the global matrix $D$. -template -void -HDG::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, - ScratchData &scratch, - PerTaskData &task_data) -{ + template + void + HDG::assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, + ScratchData &scratch, + PerTaskData &task_data) + { // Construct iterator for dof_handler_local for FEValues reinit function. - typename DoFHandler::active_cell_iterator - loc_cell (&triangulation, - cell->level(), - cell->index(), - &dof_handler_local); + typename DoFHandler::active_cell_iterator + loc_cell (&triangulation, + cell->level(), + cell->index(), + &dof_handler_local); - const unsigned int n_q_points = scratch.fe_values_local.get_quadrature().size(); - const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size(); + const unsigned int n_q_points = scratch.fe_values_local.get_quadrature().size(); + const unsigned int n_face_q_points = scratch.fe_face_values_local.get_quadrature().size(); - const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell; + const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell; - const FEValuesExtractors::Vector fluxes (0); - const FEValuesExtractors::Scalar scalar (dim); + const FEValuesExtractors::Vector fluxes (0); + const FEValuesExtractors::Scalar scalar (dim); - scratch.ll_matrix = 0; - scratch.l_rhs = 0; - if (!task_data.trace_reconstruct) - { - scratch.lf_matrix = 0; - scratch.fl_matrix = 0; - task_data.cell_matrix = 0; - task_data.cell_vector = 0; - } - scratch.fe_values_local.reinit (loc_cell); - - // We first compute the cell-interior contribution to @p ll_matrix matrix - // (referred to as matrix $A$ in the introduction) corresponding to - // local-local coupling, as well as the local right-hand-side vector. We - // store the values at each quadrature point for the basis functions, the - // right-hand-side value, and the convection velocity. - for (unsigned int q=0; q convection - = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q)); - const double JxW = scratch.fe_values_local.JxW(q); - for (unsigned int k=0; k::faces_per_cell; ++face) - { - scratch.fe_face_values_local.reinit(loc_cell, face); - scratch.fe_face_values.reinit(cell, face); + scratch.ll_matrix = 0; + scratch.l_rhs = 0; + if (!task_data.trace_reconstruct) + { + scratch.lf_matrix = 0; + scratch.fl_matrix = 0; + task_data.cell_matrix = 0; + task_data.cell_vector = 0; + } + scratch.fe_values_local.reinit (loc_cell); - // The already obtained $\hat{u}$ values are needed when solving for the - // local variables. - if (task_data.trace_reconstruct) - scratch.fe_face_values.get_function_values (solution, scratch.trace_values); + // We first compute the cell-interior contribution to @p ll_matrix matrix + // (referred to as matrix $A$ in the introduction) corresponding to + // local-local coupling, as well as the local right-hand-side vector. We + // store the values at each quadrature point for the basis functions, the + // right-hand-side value, and the convection velocity. + for (unsigned int q=0; q convection + = scratch.convection_velocity.value(scratch.fe_values_local.quadrature_point(q)); + const double JxW = scratch.fe_values_local.JxW(q); + for (unsigned int k=0; k quadrature_point = - scratch.fe_face_values.quadrature_point(q); - const Point normal = scratch.fe_face_values.normal_vector(q); - const Tensor<1,dim> convection - = scratch.convection_velocity.value(quadrature_point); - - // Here we compute the stabilization parameter discussed in the - // introduction: since the diffusion is one and the diffusion length - // scale is set to 1/5, it simply results in a contribution of 5 for - // the diffusion part and the magnitude of convection through the - // element boundary in a centered-like scheme for the convection - // part. - const double tau_stab = (5. + - std::abs(convection * normal)); - - // We store the non-zero flux and scalar values, making use of the - // support_on_face information we calculated in @p ScratchData. - for (unsigned int k=0; k::faces_per_cell; ++face) + { + scratch.fe_face_values_local.reinit(loc_cell, face); + scratch.fe_face_values.reinit(cell, face); + + // The already obtained $\hat{u}$ values are needed when solving for the + // local variables. + if (task_data.trace_reconstruct) + scratch.fe_face_values.get_function_values (solution, scratch.trace_values); + + for (unsigned int q=0; q quadrature_point = + scratch.fe_face_values.quadrature_point(q); + const Point normal = scratch.fe_face_values.normal_vector(q); + const Tensor<1,dim> convection + = scratch.convection_velocity.value(quadrature_point); + + // Here we compute the stabilization parameter discussed in the + // introduction: since the diffusion is one and the diffusion length + // scale is set to 1/5, it simply results in a contribution of 5 for + // the diffusion part and the magnitude of convection through the + // element boundary in a centered-like scheme for the convection + // part. + const double tau_stab = (5. + + std::abs(convection * normal)); + + // We store the non-zero flux and scalar values, making use of the + // support_on_face information we calculated in @p ScratchData. + for (unsigned int k=0; kface(face)->at_boundary() - && - (cell->face(face)->boundary_indicator() == 1)) - { - const double neumann_value = - - scratch.exact_solution.gradient (quadrature_point) * normal - + convection * normal * scratch.exact_solution.value(quadrature_point); - for (unsigned int i=0; i_{\partial \mathcal T}$ to the local matrix. As opposed - // to the face matrices above, we need it in both assembly stages. - for (unsigned int i=0; iface(face)->at_boundary() + && + (cell->face(face)->boundary_indicator() == 1)) + { + const double neumann_value = + - scratch.exact_solution.gradient (quadrature_point) * normal + + convection * normal * scratch.exact_solution.value(quadrature_point); + for (unsigned int i=0; i_{\partial \mathcal T}$ to the local matrix. As opposed + // to the face matrices above, we need it in both assembly stages. for (unsigned int i=0; iget_dof_indices(task_data.dof_indices); - } - // For (2), we are simply solving (ll_matrix).(solution_local) = (l_rhs). - // Hence, we multiply @p l_rhs by our already inverted local-local matrix - // and store the result using the set_dof_values function. - else - { - scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs); - loc_cell->set_dof_values(scratch.tmp_rhs, solution_local); - } -} + // Once assembly of all of the local contributions is complete, we must either: + // (1) assemble the global system, or (2) compute the local solution values and + // save them. + // In either case, the first step is to invert the local-local matrix. + scratch.ll_matrix.gauss_jordan(); + + // For (1), we compute the Schur complement and add it to the @p + // cell_matrix, matrix $D$ in the introduction. + if (task_data.trace_reconstruct == false) + { + scratch.fl_matrix.mmult(scratch.tmp_matrix, scratch.ll_matrix); + scratch.tmp_matrix.vmult_add(task_data.cell_vector, scratch.l_rhs); + scratch.tmp_matrix.mmult(task_data.cell_matrix, scratch.lf_matrix, true); + cell->get_dof_indices(task_data.dof_indices); + } + // For (2), we are simply solving (ll_matrix).(solution_local) = (l_rhs). + // Hence, we multiply @p l_rhs by our already inverted local-local matrix + // and store the result using the set_dof_values function. + else + { + scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs); + loc_cell->set_dof_values(scratch.tmp_rhs, solution_local); + } + } // @sect4{HDG::solve} // The skeleton solution is solved for by using a BiCGStab solver with // identity preconditioner. -template -void HDG::solve () -{ - SolverControl solver_control (system_matrix.m()*10, - 1e-11*system_rhs.l2_norm()); - SolverBicgstab<> solver (solver_control, false); - solver.solve (system_matrix, solution, system_rhs, - PreconditionIdentity()); + template + void HDG::solve () + { + SolverControl solver_control (system_matrix.m()*10, + 1e-11*system_rhs.l2_norm()); + SolverBicgstab<> solver (solver_control, false); + solver.solve (system_matrix, solution, system_rhs, + PreconditionIdentity()); - std::cout << " Number of BiCGStab iterations: " << solver_control.last_step() - << std::endl; + std::cout << " Number of BiCGStab iterations: " << solver_control.last_step() + << std::endl; - system_matrix.clear(); - sparsity_pattern.reinit(0,0,0,1); + system_matrix.clear(); + sparsity_pattern.reinit(0,0,0,1); - constraints.distribute(solution); + constraints.distribute(solution); - // Once we have solved for the skeleton solution, - // we can solve for the local solutions in an element-by-element - // fashion. We do this by re-using the same @p assemble_system function - // but switching @p trace_reconstruct to true. - assemble_system(true); -} + // Once we have solved for the skeleton solution, + // we can solve for the local solutions in an element-by-element + // fashion. We do this by re-using the same @p assemble_system function + // but switching @p trace_reconstruct to true. + assemble_system(true); + } @@ -1005,69 +1007,69 @@ void HDG::solve () // SolutionAndGradient class introduced above that contains the analytic // parts of either of them. Eventually, we also compute the L2-error of the // post-processed solution and add the results into the convergence table. -template -void -HDG::postprocess() -{ + template + void + HDG::postprocess() { - const QGauss quadrature_formula(fe_u_post.degree+1); - const UpdateFlags local_flags (update_values); - const UpdateFlags flags ( update_values | update_gradients | - update_JxW_values); + { + const QGauss quadrature_formula(fe_u_post.degree+1); + const UpdateFlags local_flags (update_values); + const UpdateFlags flags ( update_values | update_gradients | + update_JxW_values); + + PostProcessScratchData scratch (fe_u_post, fe_local, + quadrature_formula, + local_flags, + flags); + + WorkStream::run(dof_handler_u_post.begin_active(), + dof_handler_u_post.end(), + std_cxx1x::bind (&HDG::postprocess_one_cell, + std_cxx1x::ref(*this), + std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3), + std_cxx1x::function(), + scratch, + 0U); + } - PostProcessScratchData scratch (fe_u_post, fe_local, - quadrature_formula, - local_flags, - flags); - - WorkStream::run(dof_handler_u_post.begin_active(), - dof_handler_u_post.end(), - std_cxx1x::bind (&HDG::postprocess_one_cell, - std_cxx1x::ref(*this), - std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::_3), - std_cxx1x::function(), - scratch, - 0U); + Vector difference_per_cell (triangulation.n_active_cells()); + + ComponentSelectFunction value_select (dim, dim+1); + VectorTools::integrate_difference (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 (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(); + + VectorTools::integrate_difference (dof_handler_u_post, + solution_u_post, + Solution(), + difference_per_cell, + QGauss(fe.degree+3), + VectorTools::L2_norm); + const double post_error = difference_per_cell.l2_norm(); + + convergence_table.add_value("cells", triangulation.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); + convergence_table.add_value("val L2-post", post_error); } - Vector difference_per_cell (triangulation.n_active_cells()); - - ComponentSelectFunction value_select (dim, dim+1); - VectorTools::integrate_difference (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 (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(); - - VectorTools::integrate_difference (dof_handler_u_post, - solution_u_post, - Solution(), - difference_per_cell, - QGauss(fe.degree+3), - VectorTools::L2_norm); - const double post_error = difference_per_cell.l2_norm(); - - convergence_table.add_value("cells", triangulation.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); - convergence_table.add_value("val L2-post", post_error); -} - // @sect4{HDG::postprocess_one_cell} @@ -1089,71 +1091,71 @@ HDG::postprocess() // row would correspond to the constant part already and deleting e.g. the // last row would give us a singular system. This way, our program can also // be used for those elements. -template -void -HDG::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, - PostProcessScratchData &scratch, - unsigned int &) -{ - typename DoFHandler::active_cell_iterator - loc_cell (&triangulation, - cell->level(), - cell->index(), - &dof_handler_local); + template + void + HDG::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, + PostProcessScratchData &scratch, + unsigned int &) + { + typename DoFHandler::active_cell_iterator + loc_cell (&triangulation, + cell->level(), + cell->index(), + &dof_handler_local); - scratch.fe_values_local.reinit (loc_cell); - scratch.fe_values.reinit(cell); + scratch.fe_values_local.reinit (loc_cell); + scratch.fe_values.reinit(cell); - FEValuesExtractors::Vector fluxes(0); - FEValuesExtractors::Scalar scalar(dim); + FEValuesExtractors::Vector fluxes(0); + FEValuesExtractors::Scalar scalar(dim); - const unsigned int n_q_points = scratch.fe_values.get_quadrature().size(); - const unsigned int dofs_per_cell = scratch.fe_values.dofs_per_cell; + const unsigned int n_q_points = scratch.fe_values.get_quadrature().size(); + const unsigned int dofs_per_cell = scratch.fe_values.dofs_per_cell; - scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values); - scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients); + scratch.fe_values_local[scalar].get_function_values(solution_local, scratch.u_values); + scratch.fe_values_local[fluxes].get_function_values(solution_local, scratch.u_gradients); - double sum = 0; - for (unsigned int i=1; idistribute_local_to_global(scratch.cell_sol, solution_u_post); -} + // Having assembled all terms, we can again go on and solve the linear + // system. We again invert the matrix and then multiply the inverse by the + // right hand side. An alternative (and more numerically stable) would have + // been to only factorize the matrix and apply the factorization. + scratch.cell_matrix.gauss_jordan(); + scratch.cell_matrix.vmult(scratch.cell_sol, scratch.cell_rhs); + cell->distribute_local_to_global(scratch.cell_sol, solution_u_post); + } @@ -1165,78 +1167,78 @@ HDG::postprocess_one_cell (const typename DoFHandler::active_cell_iter // to the same vtk file, even though they correspond to different DoFHandler // objects. The graphical output for the skeleton variable is done through // use of the DataOutFaces class. -template -void HDG::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()); - } + template + void HDG::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()); + } - std::string face_out(filename); - face_out += "-face"; + 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()); + 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; + DataOut data_out; // We first define the names and types of the local solution, // and add the data to @p 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); + 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); // The second data item we add is the post-processed solution. // In this case, it is a single scalar variable belonging to // a different DoFHandler. - 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); + 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); - data_out.build_patches (fe.degree); - data_out.write_vtk (output); + data_out.build_patches (fe.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()); + 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()); // The DataOutFaces class works analagously to the DataOut // class when we have a DoFHandler that defines the solution on // the skeleton of the triangulation. We treat it as such here, and the code is // similar to that above. - 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); -} + 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); + } // @sect4{HDG::refine_grid} @@ -1250,117 +1252,117 @@ void HDG::output_results (const unsigned int cycle) // The adaptive_refinement mode uses the KellyErrorEstimator to // give a decent indication of the non-regular regions in the scalar local // solutions. -template -void HDG::refine_grid (const unsigned int cycle) -{ - if (cycle == 0) - { - GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1); - triangulation.refine_global(3-dim); - } - else - switch (refinement_mode) - { - case global_refinement: + template + void HDG::refine_grid (const unsigned int cycle) + { + if (cycle == 0) { - triangulation.clear(); - GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); - triangulation.refine_global(3-dim+cycle/2); - break; + GridGenerator::subdivided_hyper_cube (triangulation, 2, -1, 1); + triangulation.refine_global(3-dim); } + else + switch (refinement_mode) + { + case global_refinement: + { + triangulation.clear(); + GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); + triangulation.refine_global(3-dim+cycle/2); + break; + } - case adaptive_refinement: - { - Vector estimated_error_per_cell (triangulation.n_active_cells()); + 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)); + 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.); + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.); - triangulation.execute_coarsening_and_refinement (); + triangulation.execute_coarsening_and_refinement (); - break; - } + break; + } - default: - { - Assert (false, ExcNotImplemented()); - } - } + default: + { + Assert (false, ExcNotImplemented()); + } + } - // Just as in step-7, we set the boundary indicator of one of the faces to 1 - // where we want to specify Neumann boundary conditions instead of Dirichlet - // conditions. Since we re-create the triangulation every time for global - // refinement, the flags are set in every refinement step, not just at the - // beginning. - typename Triangulation::cell_iterator - cell = triangulation.begin (), - endc = triangulation.end(); - for (; cell!=endc; ++cell) - for (unsigned int face=0; face::faces_per_cell; ++face) - if (cell->face(face)->at_boundary()) - if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12) - || - (std::fabs(cell->face(face)->center()(1) - (-1)) < 1e-12)) - cell->face(face)->set_boundary_indicator (1); -} + // Just as in step-7, we set the boundary indicator of one of the faces to 1 + // where we want to specify Neumann boundary conditions instead of Dirichlet + // conditions. Since we re-create the triangulation every time for global + // refinement, the flags are set in every refinement step, not just at the + // beginning. + typename Triangulation::cell_iterator + cell = triangulation.begin (), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12) + || + (std::fabs(cell->face(face)->center()(1) - (-1)) < 1e-12)) + cell->face(face)->set_boundary_indicator (1); + } // @sect4{HDG::run} // The functionality here is basically the same as Step-7. // We loop over 10 cycles, refining the grid on each one. At the end, // convergence tables are created. -template -void HDG::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); - } + template + void HDG::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); + } - 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); - // There is one minor change for the convergence table compared to step-7: - // Since we did not refine our mesh by a factor two in each cycle (but - // rather used the sequence 2, 3, 4, 6, 8, 12, ...), we need to tell the - // convergence rate evaluation about this. We do this by setting the number - // of cells as a reference column and additionally specifying the dimension - // of the problem, which gives the computation the necessary information for - // how much the mesh was refinement given a certain increase in the number - // of cells. - if (refinement_mode == global_refinement) - { - 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); -} + // There is one minor change for the convergence table compared to step-7: + // Since we did not refine our mesh by a factor two in each cycle (but + // rather used the sequence 2, 3, 4, 6, 8, 12, ...), we need to tell the + // convergence rate evaluation about this. We do this by setting the number + // of cells as a reference column and additionally specifying the dimension + // of the problem, which gives the computation the necessary information for + // how much the mesh was refinement given a certain increase in the number + // of cells. + if (refinement_mode == global_refinement) + { + 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); + } } // end of namespace Step51 -- 2.39.5