From: kronbichler Date: Sun, 25 Aug 2013 18:43:27 +0000 (+0000) Subject: Write comments, cleanup second WorkStream call by using an empty function instead... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=38ca9f70419efb35757682ece19c5ed67d9baf6a;p=dealii-svn.git Write comments, cleanup second WorkStream call by using an empty function instead of a separate object, reinident file with astyle. git-svn-id: https://svn.dealii.org/trunk@30483 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index fcf883d205..1f838c2afd 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -62,12 +62,14 @@ // on codimension-2 surfaces. #include -// The second new file we include defines a new type of sparse matrix. -// The regular SparseMatrix type stores indices to all non-zero entries. -// The ChunkSparseMatrix takes advantage of the coupled nature of -// DG solutions. It stores an index to a matrix sub-block of a specified -// size. In the HDG context, this sub-block-size is actually the number -// of degrees of freedom per face defined by the skeleton solution field. +// The second new file we include defines a new type of sparse matrix. The +// regular SparseMatrix type stores indices to all non-zero +// entries. The ChunkSparseMatrix takes advantage of the coupled +// nature of DG solutions. It stores an index to a matrix sub-block of a +// specified size. In the HDG context, this sub-block-size is actually the +// number of degrees of freedom per face defined by the skeleton solution +// field. This reduces the memory consumption of the matrix by up to one third +// and results in similar speedups when using the matrix in solvers. #include // The final new include for this example deals with data output. Since @@ -82,10 +84,10 @@ 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. +// 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 scale the solution so its norm is of order unity for all +// values of the solution width. template class SolutionBase { @@ -102,7 +104,7 @@ 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 <> @@ -111,60 +113,23 @@ 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) }; + Point<3>(+0.5, -0.5, 0.5) + }; template const double SolutionBase::width = 1./5.; - -template -class ConvectionVelocity : public TensorFunction<1,dim> -{ -public: - ConvectionVelocity() : TensorFunction<1,dim>() {} - - virtual Tensor<1,dim> value (const Point &p) const; -}; - - - -template -Tensor<1,dim> -ConvectionVelocity::value(const Point &p) const -{ - Tensor<1,dim> convection; - switch (dim) - { - case 1: - convection[0] = 1; - break; - case 2: - convection[0] = p[1]; - convection[1] = -p[0]; - break; - case 3: - convection[0] = p[1]; - convection[1] = -p[0]; - convection[2] = 1; - break; - default: - Assert(false, ExcNotImplemented()); - } - return convection; -} - - template class Solution : public Function, - protected SolutionBase + protected SolutionBase { public: Solution () : Function() {} @@ -191,7 +156,7 @@ double Solution::value (const Point &p, } return return_value / - Utilities::fixed_power(std::sqrt(2. * numbers::PI) * this->width); + Utilities::fixed_power(std::sqrt(2. * numbers::PI) * this->width); } @@ -218,30 +183,84 @@ Tensor<1,dim> Solution::gradient (const Point &p, +// This class implements a function where the scalar solution and its negative +// 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 + protected SolutionBase { public: SolutionAndGradient () : Function(dim) {} virtual void vector_value (const Point &p, - Vector &v) const - { - AssertDimension(v.size(), dim+1); - Solution solution; - Tensor<1,dim> grad = solution.gradient(p); - for (unsigned int d=0; d &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 +class ConvectionVelocity : public TensorFunction<1,dim> +{ +public: + ConvectionVelocity() : TensorFunction<1,dim>() {} + + virtual Tensor<1,dim> value (const Point &p) const; +}; + + + +template +Tensor<1,dim> +ConvectionVelocity::value(const Point &p) const +{ + Tensor<1,dim> convection; + switch (dim) + { + case 1: + convection[0] = 1; + break; + case 2: + convection[0] = p[1]; + convection[1] = -p[0]; + break; + case 3: + convection[0] = p[1]; + convection[1] = -p[0]; + convection[2] = 1; + break; + default: + Assert(false, ExcNotImplemented()); + } + return convection; +} + +// The last function we implement is the right hand side for the manufactured +// solution. It is very similar to step-7, with the exception that we now have +// 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 + protected SolutionBase { public: RightHandSide () : Function() {} @@ -278,99 +297,97 @@ double RightHandSide::value (const Point &p, // @sect3{The Step51 HDG solver class} -// The HDG solution procedure follows closely that of step-7. The major -// difference is the use of 3 different sets of DoFHandler and FE objects, -// along with the ChunkSparseMatrix and the corresponding solutions vectors. -// We also use WorkStream to enable a multi-threaded local solution process which exploits -// the embarrassingly parallel nature of the local solver. - +// The HDG solution procedure follows closely that of step-7. The major +// difference is the use of 3 different sets of DoFHandler and FE +// objects, along with the ChunkSparseMatrix and the +// corresponding solutions vectors. We also use WorkStream to enable a +// multi-threaded local solution process which exploits the embarrassingly +// parallel nature of the local solver. For WorkStream, we define the local +// operations on a cell and a copy function into the global matrix and +// vector. We do this once for the assembly (which is run twice, once when we +// generate the system matrix and once when we compute the element-interior +// solutions from the skeleton values) and once for the postprocessing where +// we extract a solution that converges at higher order. template class Step51 { public: enum RefinementMode - { - global_refinement, adaptive_refinement - }; + { + global_refinement, adaptive_refinement + }; Step51 (const unsigned int degree, const RefinementMode refinement_mode); void run (); private: - + struct PerTaskData; struct ScratchData; - - struct EmptyData; + struct PostProcessScratchData; - void NullFunction(const EmptyData &data); - + void setup_system (); void assemble_system (const bool reconstruct_trace = false); void assemble_system_one_cell (const typename DoFHandler::active_cell_iterator &cell, - ScratchData &scratch, - PerTaskData &task_data); + ScratchData &scratch, + PerTaskData &task_data); void copy_local_to_global(const PerTaskData &data); void solve (); void postprocess (); void postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, - PostProcessScratchData &scratch, - EmptyData &task_data); + PostProcessScratchData &scratch, + unsigned int &empty_data); void refine_grid (const unsigned int cylce); void output_results (const unsigned int cycle); 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$. + + // 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 solution that couples the element-level local -// solution. + // The new finite element type and corresponding DoFHandler are + // used for the global 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 superconvegence 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 to represent this post-processed solution, which we -// only use for output after constructing it. + // As stated in the introduction, HDG solutions can be post-processed to + // attain superconvegence 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. + // 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; - - // Comment on chunk. + + // 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; - // Global/skeleton solution/rhs - Vector solution; - Vector system_rhs; - - // Local elementwise solution - 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; }; @@ -389,137 +406,141 @@ Step51::Step51 (const unsigned int degree, refinement_mode (refinement_mode) {} + + +// First come the definition of the local data structures for the parallel +// assembly. The first structure @p PerTaskData contains the local vector and +// matrix that are written into the global matrix, whereas the ScratchData +// contains all data that we need for the local assembly. template struct Step51::PerTaskData { - FullMatrix cell_matrix; - Vector cell_vector; - std::vector dof_indices; - - bool trace_reconstruct; - - PerTaskData(const unsigned int n_dofs, const bool trace_reconstruct) + 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) - {} - - void reset(){ - cell_matrix = 0.0; - cell_vector = 0.0; - } + {} + + void reset() + { + cell_matrix = 0.0; + cell_vector = 0.0; + } }; + + template struct Step51::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; - - // Full constructor - 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; ireset(), but -// need to have it for the WorkStream interface. - void reset() {} - -}; + 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; + + // Full constructor + 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; ireset(), but + // need to have it for the WorkStream interface. + void reset() {} -template -struct Step51::EmptyData -{ - EmptyData(){} - void reset(){} }; template @@ -527,65 +548,63 @@ struct Step51::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; - - // Full constructor + + // Full constructor 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) - {} - - // Copy constructor + 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) + {} + + // Copy constructor 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) - {} - - void reset() { - cell_matrix = 0.; - cell_rhs = 0.; - cell_sol = 0.; - } - -}; + : + 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) + {} + + void reset() + { + cell_matrix = 0.; + cell_rhs = 0.; + cell_sol = 0.; + } -template -void Step51::NullFunction(const EmptyData &data){} +}; template void Step51::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); + if (data.trace_reconstruct == false) + constraints.distribute_local_to_global (data.cell_matrix, + data.cell_vector, + data.dof_indices, + system_matrix, system_rhs); } template @@ -632,34 +651,34 @@ template void Step51::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, - &Step51::assemble_system_one_cell, - &Step51::copy_local_to_global, - scratch, - task_data); + 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, + &Step51::assemble_system_one_cell, + &Step51::copy_local_to_global, + scratch, + task_data); } @@ -669,17 +688,17 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce ScratchData &scratch, PerTaskData &task_data) { - // Construct iterator for dof_handler_local + // Construct iterator for dof_handler_local typename DoFHandler::active_cell_iterator - loc_cell (&triangulation, - cell->level(), - cell->index(), - &dof_handler_local); + 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 dofs_per_cell = scratch.fe_face_values.get_fe().dofs_per_cell; + // const unsigned int dofs_per_cell = scratch.fe_face_values.get_fe().dofs_per_cell; const unsigned int loc_dofs_per_cell = scratch.fe_values_local.get_fe().dofs_per_cell; // Choose stabilization parameter to be 5 * diffusion = 5 @@ -688,156 +707,156 @@ Step51::assemble_system_one_cell (const typename DoFHandler::active_ce const FEValuesExtractors::Vector fluxes (0); const FEValuesExtractors::Scalar scalar (dim); - scratch.ll_matrix = 0; - scratch.l_rhs = 0; - if (!task_data.trace_reconstruct) + scratch.ll_matrix = 0; + scratch.l_rhs = 0; + if (!task_data.trace_reconstruct) + { + scratch.lf_matrix = 0; + scratch.fl_matrix = 0; + task_data.reset(); + } + scratch.fe_values_local.reinit (loc_cell); + + 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 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) + for (unsigned int face=0; face::faces_per_cell; ++face) + { + scratch.fe_face_values_local.reinit(loc_cell, face); + scratch.fe_face_values.reinit(cell, face); + 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); + const double tau_stab = (tau_stab_diffusion + + std::abs(convection * normal)); - 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); - 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)) - { - const double neumann_value = - - scratch.exact_solution.gradient (quadrature_point) * normal - + convection * normal * scratch.exact_solution.value(quadrature_point); - 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; iget_dof_indices(task_data.dof_indices); - } - else - { - scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs); - loc_cell->set_dof_values(scratch.tmp_rhs, solution_local); + for (unsigned int i=0; iget_dof_indices(task_data.dof_indices); + } + else + { + scratch.ll_matrix.vmult(scratch.tmp_rhs, scratch.l_rhs); + loc_cell->set_dof_values(scratch.tmp_rhs, solution_local); + } } @@ -856,7 +875,7 @@ void Step51::solve () system_matrix.clear(); sparsity_pattern.reinit(0,0,0,1); - + constraints.distribute(solution); // update local values @@ -875,25 +894,25 @@ Step51::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); - + update_JxW_values); + EmptyData task_data; - + PostProcessScratchData scratch (fe_u_post, fe_local, - quadrature_formula, - local_flags, - flags); - + quadrature_formula, + local_flags, + flags); + WorkStream::run(dof_handler_u_post.begin_active(), dof_handler_u_post.end(), *this, &Step51::postprocess_one_cell, - &Step51::NullFunction, + std_cxx1x::function(), scratch, task_data); - } + } -// Compute some convergence rates, etc., and add to a table + // Compute some convergence rates, etc., and add to a table Vector difference_per_cell (triangulation.n_active_cells()); ComponentSelectFunction value_select (dim, dim+1); @@ -924,31 +943,31 @@ Step51::postprocess() 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); + + 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); } template void Step51::postprocess_one_cell (const typename DoFHandler::active_cell_iterator &cell, - PostProcessScratchData &scratch, - EmptyData &task_data) + PostProcessScratchData &scratch, + unsigned int &) { scratch.reset(); typename DoFHandler::active_cell_iterator - loc_cell (&triangulation, - cell->level(), - cell->index(), - &dof_handler_local); - + loc_cell (&triangulation, + cell->level(), + cell->index(), + &dof_handler_local); + scratch.fe_values_local.reinit (loc_cell); scratch.fe_values.reinit(cell); - + FEValuesExtractors::Vector fluxes(0); FEValuesExtractors::Scalar scalar(dim); @@ -957,38 +976,38 @@ Step51::postprocess_one_cell (const typename DoFHandler::active_cell_i 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; i::output_results (const unsigned int cycle) default: Assert (false, ExcNotImplemented()); } - + std::string face_out(filename); face_out += "-face"; @@ -1025,24 +1044,24 @@ void Step51::output_results (const unsigned int cycle) std::vector names (dim, "gradient"); names.push_back ("solution"); std::vector - component_interpretation - (dim+1, DataComponentInterpretation::component_is_part_of_vector); + component_interpretation + (dim+1, DataComponentInterpretation::component_is_part_of_vector); component_interpretation[dim] - = DataComponentInterpretation::component_is_scalar; + = DataComponentInterpretation::component_is_scalar; data_out.add_data_vector (dof_handler_local, solution_local, names, component_interpretation); - - // Post-processed solution: can now add more than 1 dof_handler to + + // Post-processed solution: can now add more than 1 dof_handler to // the DataOut object! std::vector post_name(1,"u_post"); - std::vector - post_comp_type(1, DataComponentInterpretation::component_is_scalar); + 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); + post_name, post_comp_type); 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"; @@ -1050,14 +1069,14 @@ void Step51::output_results (const unsigned int cycle) 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); - + 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); } @@ -1076,39 +1095,39 @@ void Step51::refine_grid (const unsigned int cycle) 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; - } + { + 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()); + { + 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()); - } + { + Assert (false, ExcNotImplemented()); + } } // Just as in step-7, we set the boundary indicator of one of the faces to 1 @@ -1117,8 +1136,8 @@ void Step51::refine_grid (const unsigned int cycle) // refinement, the flags are set in every refinement step, not just at the // beginning. typename Triangulation::cell_iterator - cell = triangulation.begin (), - endc = triangulation.end(); + cell = triangulation.begin (), + endc = triangulation.end(); for (; cell!=endc; ++cell) for (unsigned int face=0; face::faces_per_cell; ++face) if ((std::fabs(cell->face(face)->center()(0) - (-1)) < 1e-12) @@ -1137,7 +1156,7 @@ 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); @@ -1158,17 +1177,17 @@ void Step51::run () if (refinement_mode == global_refinement) { convergence_table - .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim); + .evaluate_convergence_rates("val L2", "cells", ConvergenceTable::reduction_rate_log2, dim); convergence_table - .evaluate_convergence_rates("grad L2", "cells", ConvergenceTable::reduction_rate_log2, dim); + .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); + .evaluate_convergence_rates("val L2-post", "cells", ConvergenceTable::reduction_rate_log2, dim); } convergence_table.write_text(std::cout); } -int main (int argc, char** argv) +int main (int argc, char **argv) { const unsigned int dim = 2;