From 959f368eabce6e6ee2895f244591db51a18a6548 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Tue, 6 Aug 2013 17:03:45 +0000 Subject: [PATCH] Remove namespace step-51 as the class name already contains step-51. Speed up assembly by caching the output of FEValuesExtractors and only working on dofs with support on faces git-svn-id: https://svn.dealii.org/trunk@30237 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-51/doc/results.dox | 159 +-- deal.II/examples/step-51/step-51.cc | 1299 +++++++++++----------- 2 files changed, 679 insertions(+), 779 deletions(-) diff --git a/deal.II/examples/step-51/doc/results.dox b/deal.II/examples/step-51/doc/results.dox index ae4bf1cc53..1b7193f9dc 100644 --- a/deal.II/examples/step-51/doc/results.dox +++ b/deal.II/examples/step-51/doc/results.dox @@ -5,135 +5,38 @@ The program writes convergence tables to the screen while running: @code -examples/\step-51> make run -============================ Running \step-51 -Solving with Q1 elements, adaptive refinement -============================================= - -Cycle 0: - Number of active cells: 4 - Number of degrees of freedom: 9 -Cycle 1: - Number of active cells: 13 - Number of degrees of freedom: 22 -Cycle 2: - Number of active cells: 31 - Number of degrees of freedom: 46 -Cycle 3: - Number of active cells: 64 - Number of degrees of freedom: 87 -Cycle 4: - Number of active cells: 127 - Number of degrees of freedom: 160 -Cycle 5: - Number of active cells: 244 - Number of degrees of freedom: 297 -Cycle 6: - Number of active cells: 466 - Number of degrees of freedom: 543 - -cycle cells dofs L2 H1 Linfty - 0 4 9 1.198e+00 2.732e+00 1.383e+00 - 1 13 22 8.795e-02 1.193e+00 1.816e-01 - 2 31 46 8.147e-02 1.167e+00 1.654e-01 - 3 64 87 7.702e-02 1.077e+00 1.310e-01 - 4 127 160 4.643e-02 7.988e-01 6.745e-02 - 5 244 297 2.470e-02 5.568e-01 3.668e-02 - 6 466 543 1.622e-02 4.107e-01 2.966e-02 - -Solving with Q1 elements, global refinement -=========================================== - -Cycle 0: - Number of active cells: 4 - Number of degrees of freedom: 9 -Cycle 1: - Number of active cells: 16 - Number of degrees of freedom: 25 -Cycle 2: - Number of active cells: 64 - Number of degrees of freedom: 81 -Cycle 3: - Number of active cells: 256 - Number of degrees of freedom: 289 -Cycle 4: - Number of active cells: 1024 - Number of degrees of freedom: 1089 -Cycle 5: - Number of active cells: 4096 - Number of degrees of freedom: 4225 -Cycle 6: - Number of active cells: 16384 - Number of degrees of freedom: 16641 - -cycle cells dofs L2 H1 Linfty - 0 4 9 1.198e+00 2.732e+00 1.383e+00 - 1 16 25 8.281e-02 1.190e+00 1.808e-01 - 2 64 81 8.142e-02 1.129e+00 1.294e-01 - 3 256 289 2.113e-02 5.828e-01 4.917e-02 - 4 1024 1089 5.319e-03 2.934e-01 1.359e-02 - 5 4096 4225 1.332e-03 1.469e-01 3.482e-03 - 6 16384 16641 3.332e-04 7.350e-02 8.758e-04 - -n cells H1 L2 - 0 4 2.732e+00 - 1.198e+00 - - - 1 16 1.190e+00 1.20 8.281e-02 14.47 3.86 - 2 64 1.129e+00 0.08 8.142e-02 1.02 0.02 - 3 256 5.828e-01 0.95 2.113e-02 3.85 1.95 - 4 1024 2.934e-01 0.99 5.319e-03 3.97 1.99 - 5 4096 1.469e-01 1.00 1.332e-03 3.99 2.00 - 6 16384 7.350e-02 1.00 3.332e-04 4.00 2.00 - -Solving with Q2 elements, global refinement -=========================================== - -Cycle 0: - Number of active cells: 4 - Number of degrees of freedom: 25 -Cycle 1: - Number of active cells: 16 - Number of degrees of freedom: 81 -Cycle 2: - Number of active cells: 64 - Number of degrees of freedom: 289 -Cycle 3: - Number of active cells: 256 - Number of degrees of freedom: 1089 -Cycle 4: - Number of active cells: 1024 - Number of degrees of freedom: 4225 -Cycle 5: - Number of active cells: 4096 - Number of degrees of freedom: 16641 -Cycle 6: - Number of active cells: 16384 - Number of degrees of freedom: 66049 - -cycle cells dofs L2 H1 Linfty - 0 4 25 1.433e+00 2.445e+00 1.286e+00 - 1 16 81 7.912e-02 1.168e+00 1.728e-01 - 2 64 289 7.755e-03 2.511e-01 1.991e-02 - 3 256 1089 9.969e-04 6.235e-02 2.764e-03 - 4 1024 4225 1.265e-04 1.571e-02 3.527e-04 - 5 4096 16641 1.587e-05 3.937e-03 4.343e-05 - 6 16384 66049 1.986e-06 9.847e-04 5.402e-06 - -n cells H1 L2 - 0 4 2.445e+00 - 1.433e+00 - - - 1 16 1.168e+00 1.07 7.912e-02 18.11 4.18 - 2 64 2.511e-01 2.22 7.755e-03 10.20 3.35 - 3 256 6.235e-02 2.01 9.969e-04 7.78 2.96 - 4 1024 1.571e-02 1.99 1.265e-04 7.88 2.98 - 5 4096 3.937e-03 2.00 1.587e-05 7.97 2.99 - 6 16384 9.847e-04 2.00 1.986e-06 7.99 3.00 +Q1 elements: + 16 80 4.570e+00 - 1.221e+01 - 4.333e+00 - + 36 168 1.869e+00 2.20 7.299e+00 1.27 1.734e+00 2.26 + 64 288 7.177e-01 3.33 4.218e+00 1.91 2.538e-01 6.68 + 144 624 2.729e-01 2.38 1.867e+00 2.01 6.110e-02 3.51 + 256 1088 1.493e-01 2.10 1.046e+00 2.01 2.878e-02 2.62 + 576 2400 6.964e-02 1.88 4.847e-01 1.90 9.202e-03 2.81 + 1024 4224 4.018e-02 1.91 2.785e-01 1.93 4.027e-03 2.87 + 2304 9408 1.831e-02 1.94 1.264e-01 1.95 1.236e-03 2.91 + 4096 16640 1.043e-02 1.96 7.185e-02 1.96 5.306e-04 2.94 + 9216 37248 4.690e-03 1.97 3.228e-02 1.97 1.600e-04 2.96 + +Q3 elements: + 16 160 2.398e-01 - 1.873e+00 - 1.354e-01 - + 36 336 5.843e-02 3.48 5.075e-01 3.22 1.882e-02 4.87 + 64 576 3.466e-02 1.82 2.534e-01 2.41 4.326e-03 5.11 + 144 1248 8.297e-03 3.53 5.925e-02 3.58 6.330e-04 4.74 + 256 2176 2.254e-03 4.53 1.636e-02 4.47 1.403e-04 5.24 + 576 4800 4.558e-04 3.94 3.278e-03 3.96 1.844e-05 5.01 + 1024 8448 1.471e-04 3.93 1.052e-03 3.95 4.378e-06 5.00 + 2304 18816 2.956e-05 3.96 2.104e-04 3.97 5.751e-07 5.01 + 4096 33280 9.428e-06 3.97 6.697e-05 3.98 1.362e-07 5.01 + 9216 74496 1.876e-06 3.98 1.330e-05 3.99 1.817e-08 4.97 @endcode -One can see the error reduction upon grid refinement, and for the -cases where global refinement was performed, also the convergence -rates can be seen. The linear and quadratic convergence rates of Q1 -and Q2 elements in the $H^1$ norm can clearly be seen, as -are the quadratic and cubic rates in the $L_2$ norm. - - +One can see the error reduction upon grid refinement, and for the cases where +global refinement was performed, also the convergence rates can be seen. The +quadratic convergence rates of Q1 elements in the $L_2$ norm for both the +scalar variable and the gradient variable can clearly be seen, as is the cubic +rate for the postprocessed scalar variable in the $L_2$ norm. Likewise, the +scalar variable and gradient for Q3 elements converge at fourth order and the +postprocessed scalar variable at fifth order. +The same convergence rates are observed in 3d. \ No newline at end of file diff --git a/deal.II/examples/step-51/step-51.cc b/deal.II/examples/step-51/step-51.cc index 256e48dbba..41d3975d36 100644 --- a/deal.II/examples/step-51/step-51.cc +++ b/deal.II/examples/step-51/step-51.cc @@ -52,769 +52,767 @@ #include #include -namespace Step51 -{ - using namespace dealii; +using namespace dealii; - // @sect3{Equation data} +// @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; - }; +// The structure of the analytic solution is the same as in step-7. There +// are two exceptions. Firstly, we also create a solution for the 3d case, +// and secondly, we take into account the convection velocity in the right +// hand side that is variable in this case. +template +class SolutionBase +{ +protected: + static const unsigned int n_source_centers = 3; + static const Point source_centers[n_source_centers]; + static const double width; +}; - template <> - const Point<1> - SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers] - = { Point<1>(-1.0 / 3.0), - Point<1>(0.0), - Point<1>(+1.0 / 3.0) - }; +template <> +const Point<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<2> +SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers] += { Point<2>(-0.5, +0.5), + Point<2>(-0.5, -0.5), + Point<2>(+0.5, -0.5) +}; - template <> - const Point<3> - SolutionBase<3>::source_centers[SolutionBase<3>::n_source_centers] - = { Point<3>(-0.5, +0.5, 0.25), - Point<3>(-0.6, -0.5, -0.125), - Point<3>(+0.5, -0.5, 0.5) }; +template <> +const 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 double SolutionBase::width = 1./5.; - template - 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; +} - template - class Solution : public Function, - protected SolutionBase - { - public: - Solution () : Function() {} +template +class Solution : public Function, + protected SolutionBase +{ +public: + Solution () : 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; - virtual Tensor<1,dim> gradient (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)); - } +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); - } + 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; +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]; + 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_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); - } + return return_value / Utilities::fixed_power(std::sqrt(2 * numbers::PI) * + this->width); +} - 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 - { - AssertDimension(v.size(), dim+1); - Solution solution; - Tensor<1,dim> grad = solution.gradient(p); - for (unsigned int d=0; d &p, + Vector &v) const + { + AssertDimension(v.size(), dim+1); + Solution solution; + Tensor<1,dim> grad = solution.gradient(p); + for (unsigned int d=0; d - class RightHandSide : public Function, - protected SolutionBase - { - public: - RightHandSide () : Function() {} +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); +} - template - class Step51 - { - public: - enum RefinementMode +template +class Step51 +{ +public: + enum RefinementMode { global_refinement, adaptive_refinement }; - Step51 (const unsigned int degree, - const RefinementMode refinement_mode); - void run (); + Step51 (const unsigned int degree, + const RefinementMode refinement_mode); + void run (); + +private: + void setup_system (); + void assemble_system (const bool reconstruct_trace = false); + void solve (); + void postprocess (); + void refine_mesh (); + void output_results (const unsigned int cycle); - 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; - Triangulation triangulation; + const MappingQ mapping; - const MappingQ mapping; + FESystem fe_local; + DoFHandler dof_handler_local; - FESystem fe_local; - DoFHandler dof_handler_local; + FE_FaceQ fe; + DoFHandler dof_handler; - FE_FaceQ fe; - DoFHandler dof_handler; + FE_DGQ fe_u_post; + DoFHandler dof_handler_u_post; - FE_DGQ fe_u_post; - DoFHandler dof_handler_u_post; + ConstraintMatrix constraints; + ChunkSparsityPattern sparsity_pattern; + ChunkSparseMatrix system_matrix; - ConstraintMatrix constraints; - ChunkSparsityPattern sparsity_pattern; - ChunkSparseMatrix system_matrix; + Vector solution; + Vector system_rhs; - Vector solution; - Vector system_rhs; + Vector solution_local; + Vector solution_u_post; - Vector solution_local; - Vector solution_u_post; + const RefinementMode refinement_mode; - const RefinementMode refinement_mode; + ConvergenceTable convergence_table; +}; - ConvergenceTable convergence_table; - }; +template +Step51::Step51 (const unsigned int degree, + const RefinementMode refinement_mode) : + mapping (3), + fe_local (FE_DGQ(degree), dim, + FE_DGQ(degree), 1), + dof_handler_local (triangulation), + fe (degree), + dof_handler (triangulation), + fe_u_post (degree+1), + dof_handler_u_post (triangulation), + refinement_mode (refinement_mode) +{} - template - Step51::Step51 (const unsigned int degree, - const RefinementMode refinement_mode) : - mapping (3), - fe_local (FE_DGQ(degree), dim, - FE_DGQ(degree), 1), - dof_handler_local (triangulation), - fe (degree), - dof_handler (triangulation), - fe_u_post (degree+1), - dof_handler_u_post (triangulation), - refinement_mode (refinement_mode) - {} +template +void +Step51::setup_system () +{ + dof_handler_local.distribute_dofs(fe_local); + dof_handler.distribute_dofs(fe); + dof_handler_u_post.distribute_dofs(fe_u_post); + + std::cout << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + solution.reinit (dof_handler.n_dofs()); + system_rhs.reinit (dof_handler.n_dofs()); + + solution_local.reinit (dof_handler_local.n_dofs()); + solution_u_post.reinit (dof_handler_u_post.n_dofs()); + + constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, constraints); + std::map boundary_values; + typename FunctionMap::type boundary_functions; + Solution solution; + boundary_functions[0] = &solution; + VectorTools::project_boundary_values (mapping, dof_handler, + boundary_functions, + QGauss(fe.degree+1), + boundary_values); + for (std::map::iterator it = boundary_values.begin(); + it != boundary_values.end(); ++it) + if (constraints.is_constrained(it->first) == false) + { + constraints.add_line(it->first); + constraints.set_inhomogeneity(it->first, it->second); + } + constraints.close (); - 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); + 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); +} - 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::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 + 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); - 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) + for (unsigned int q=0; qfirst); - constraints.set_inhomogeneity(it->first, it->second); + const double rhs_value + = right_hand_side.value(fe_values_local.quadrature_point(q)); + const Tensor<1,dim> 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; k - void - Step51::assemble_system (const bool trace_reconstruct) - { - QGauss quadrature_formula(fe.degree+1); - QGauss face_quadrature_formula(fe.degree+1); - - FEValues fe_values_local (mapping, fe_local, quadrature_formula, - update_values | update_gradients | - update_JxW_values | update_quadrature_points); - FEFaceValues fe_face_values (mapping, fe, face_quadrature_formula, - update_values | update_normal_vectors | - update_quadrature_points | - update_JxW_values); - FEFaceValues fe_face_values_local (mapping, fe_local, - face_quadrature_formula, - update_values); - - const unsigned int n_q_points = quadrature_formula.size(); - const unsigned int n_face_q_points = face_quadrature_formula.size(); - - const unsigned int dofs_per_cell = fe.dofs_per_cell; - const unsigned int loc_dofs_per_cell = fe_local.dofs_per_cell; - - FullMatrix ll_matrix (loc_dofs_per_cell, loc_dofs_per_cell); - FullMatrix lf_matrix (loc_dofs_per_cell, dofs_per_cell); - FullMatrix fl_matrix (dofs_per_cell, loc_dofs_per_cell); - FullMatrix tmp_matrix (dofs_per_cell, loc_dofs_per_cell); - FullMatrix ff_matrix (dofs_per_cell, dofs_per_cell); - Vector l_rhs (loc_dofs_per_cell); - Vector f_rhs (dofs_per_cell); - Vector tmp_rhs (loc_dofs_per_cell); - - std::vector dof_indices (dofs_per_cell); - std::vector loc_dof_indices (loc_dofs_per_cell); - - ConvectionVelocity convection; - std::vector > convection_values (n_q_points); - std::vector > convection_values_face (n_face_q_points); - - std::vector trace_values(n_face_q_points); - - // Choose stabilization parameter to be 5 * diffusion = 5 - const double tau_stab_diffusion = 5.; - std::vector tau_stab (n_q_points); - - RightHandSide right_hand_side; - std::vector rhs_values (n_q_points); - - const Solution exact_solution; - std::vector neumann_values (n_face_q_points); - - const FEValuesExtractors::Vector gradients (0); - const FEValuesExtractors::Scalar values (dim); - - typename DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - loc_cell = dof_handler_local.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell, ++loc_cell) - { - if (!trace_reconstruct) - { - lf_matrix = 0; - fl_matrix = 0; - ff_matrix = 0; - f_rhs = 0; - } - fe_values_local.reinit (loc_cell); - right_hand_side.value_list (fe_values_local.get_quadrature_points(), - rhs_values); - convection.value_list(fe_values_local.get_quadrature_points(), - convection_values); - - for (unsigned int i=0; i::faces_per_cell; ++face) - { - fe_face_values_local.reinit(loc_cell, face); - fe_face_values.reinit(cell, face); - const std::vector &JxW = fe_face_values.get_JxW_values(); - const std::vector > &normals = - fe_face_values.get_normal_vectors(); - convection.value_list(fe_face_values.get_quadrature_points(), - convection_values_face); - for (unsigned int q=0; qface(face)->at_boundary() + && + (cell->face(face)->boundary_indicator() == 1)) { - double sum_lf = 0., sum_fl = 0.; - for (unsigned int q=0; qface(face)->at_boundary() - && - (cell->face(face)->boundary_indicator() == 1)) + } + + 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); - } - } - } + } + } + ll_matrix.gauss_jordan(); + if (trace_reconstruct == false) + { + fl_matrix.mmult(tmp_matrix, ll_matrix); + tmp_matrix.vmult_add(f_rhs, l_rhs); + 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()); - std::cout << " Number of GMRES iterations: " << solver_control.last_step() - << std::endl; +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()); - system_matrix.clear(); - sparsity_pattern.reinit(0,0,0,1); - constraints.distribute(solution); + std::cout << " Number of GMRES iterations: " << solver_control.last_step() + << std::endl; - // update local values - assemble_system(true); - } + system_matrix.clear(); + sparsity_pattern.reinit(0,0,0,1); + constraints.distribute(solution); + // update local values + assemble_system(true); +} - template - void - Step51::postprocess() - { - const unsigned int n_active_cells=triangulation.n_active_cells(); - Vector difference_per_cell (triangulation.n_active_cells()); - - ComponentSelectFunction value_select (dim, dim+1); - VectorTools::integrate_difference (mapping, dof_handler_local, - solution_local, - SolutionAndGradient(), - difference_per_cell, - QGauss(fe.degree+2), - VectorTools::L2_norm, - &value_select); - const double L2_error = difference_per_cell.l2_norm(); - - ComponentSelectFunction gradient_select (std::pair(0, dim), - dim+1); - VectorTools::integrate_difference (mapping, dof_handler_local, - solution_local, - SolutionAndGradient(), - difference_per_cell, - QGauss(fe.degree+2), - VectorTools::L2_norm, - &gradient_select); - const double grad_error = difference_per_cell.l2_norm(); - - convergence_table.add_value("cells", n_active_cells); - convergence_table.add_value("dofs", dof_handler.n_dofs()); - convergence_table.add_value("val L2", L2_error); - convergence_table.add_value("grad L2", grad_error); - - // construct post-processed solution with (hopefully) higher order of - // accuracy - QGauss quadrature(fe_u_post.degree+1); - FEValues fe_values(mapping, fe_u_post, quadrature, - update_values | update_JxW_values | - update_gradients); - - const unsigned int n_q_points = quadrature.size(); - std::vector u_values(n_q_points); - std::vector > u_gradients(n_q_points); - FEValuesExtractors::Vector gradients(0); - FEValuesExtractors::Scalar values(dim); - FEValues fe_values_local(mapping, fe_local, quadrature, update_values); - FullMatrix cell_matrix(fe_u_post.dofs_per_cell, - fe_u_post.dofs_per_cell); - Vector cell_rhs(fe_u_post.dofs_per_cell); - Vector cell_sol(fe_u_post.dofs_per_cell); - - typename DoFHandler::active_cell_iterator - cell_loc = dof_handler_local.begin_active(), - cell = dof_handler_u_post.begin_active(), - endc = dof_handler_u_post.end(); - for ( ; cell != endc; ++cell, ++cell_loc) - { - fe_values.reinit(cell); - fe_values_local.reinit(cell_loc); - - fe_values_local[values].get_function_values(solution_local, u_values); - fe_values_local[gradients].get_function_values(solution_local, u_gradients); - for (unsigned int i=1; i +void +Step51::postprocess() +{ + const unsigned int n_active_cells=triangulation.n_active_cells(); + Vector difference_per_cell (triangulation.n_active_cells()); + + ComponentSelectFunction value_select (dim, dim+1); + VectorTools::integrate_difference (mapping, dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(fe.degree+2), + VectorTools::L2_norm, + &value_select); + const double L2_error = difference_per_cell.l2_norm(); + + ComponentSelectFunction gradient_select (std::pair(0, dim), + dim+1); + VectorTools::integrate_difference (mapping, dof_handler_local, + solution_local, + SolutionAndGradient(), + difference_per_cell, + QGauss(fe.degree+2), + VectorTools::L2_norm, + &gradient_select); + const double grad_error = difference_per_cell.l2_norm(); + + convergence_table.add_value("cells", n_active_cells); + convergence_table.add_value("dofs", dof_handler.n_dofs()); + convergence_table.add_value("val L2", L2_error); + convergence_table.add_value("grad L2", grad_error); + + // construct post-processed solution with (hopefully) higher order of + // accuracy + QGauss quadrature(fe_u_post.degree+1); + FEValues fe_values(mapping, fe_u_post, quadrature, + update_values | update_JxW_values | + update_gradients); + + const unsigned int n_q_points = quadrature.size(); + std::vector u_values(n_q_points); + std::vector > u_gradients(n_q_points); + FEValuesExtractors::Vector gradients(0); + FEValuesExtractors::Scalar values(dim); + FEValues fe_values_local(mapping, fe_local, quadrature, update_values); + FullMatrix cell_matrix(fe_u_post.dofs_per_cell, + fe_u_post.dofs_per_cell); + Vector cell_rhs(fe_u_post.dofs_per_cell); + Vector cell_sol(fe_u_post.dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell_loc = dof_handler_local.begin_active(), + cell = dof_handler_u_post.begin_active(), + endc = dof_handler_u_post.end(); + for ( ; cell != endc; ++cell, ++cell_loc) + { + fe_values.reinit(cell); + fe_values_local.reinit(cell_loc); + + fe_values_local[values].get_function_values(solution_local, u_values); + fe_values_local[gradients].get_function_values(solution_local, u_gradients); + for (unsigned int i=1; idistribute_local_to_global(cell_sol, solution_u_post); + { + double sum = 0; + for (unsigned int q=0; q(), - 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); + } + VectorTools::integrate_difference (mapping, dof_handler_u_post, + solution_u_post, + Solution(), + difference_per_cell, + QGauss(fe.degree+3), + VectorTools::L2_norm); + double post_error = difference_per_cell.l2_norm(); + convergence_table.add_value("val L2-post", post_error); +} - template - void Step51::output_results (const unsigned int cycle) - { - std::string filename; - switch (refinement_mode) - { - case global_refinement: - filename = "solution-global"; - break; - case adaptive_refinement: - filename = "solution-adaptive"; - break; - default: - Assert (false, ExcNotImplemented()); - } - filename += "-q" + Utilities::int_to_string(fe.degree,1); - filename += "-" + Utilities::int_to_string(cycle,2); - filename += ".vtk"; - std::ofstream output (filename.c_str()); - - DataOut data_out; - std::vector names (dim, "gradient"); - names.push_back ("solution"); - std::vector - component_interpretation - (dim+1, DataComponentInterpretation::component_is_part_of_vector); - component_interpretation[dim] - = DataComponentInterpretation::component_is_scalar; - data_out.add_data_vector (dof_handler_local, solution_local, - names, component_interpretation); - - data_out.build_patches (fe.degree); - data_out.write_vtk (output); - } +template +void Step51::output_results (const unsigned int cycle) +{ + std::string filename; + switch (refinement_mode) + { + case global_refinement: + filename = "solution-global"; + break; + case adaptive_refinement: + filename = "solution-adaptive"; + break; + default: + Assert (false, ExcNotImplemented()); + } + filename += "-q" + Utilities::int_to_string(fe.degree,1); + filename += "-" + Utilities::int_to_string(cycle,2); + filename += ".vtk"; + std::ofstream output (filename.c_str()); + + DataOut data_out; + std::vector names (dim, "gradient"); + names.push_back ("solution"); + std::vector + component_interpretation + (dim+1, DataComponentInterpretation::component_is_part_of_vector); + component_interpretation[dim] + = DataComponentInterpretation::component_is_scalar; + data_out.add_data_vector (dof_handler_local, solution_local, + names, component_interpretation); + + data_out.build_patches (fe.degree); + data_out.write_vtk (output); +} - template - void Step51::run () - { - const bool do_cube = true; - if (!do_cube) - { - GridGenerator::hyper_ball (triangulation); - static const HyperBallBoundary boundary; - triangulation.set_boundary(0, boundary); - triangulation.refine_global(6-2*dim); - } - for (unsigned int cycle=0; cycle<10; ++cycle) - { - std::cout << "Cycle " << cycle << ':' << std::endl; - - if (do_cube) - { - triangulation.clear(); - GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); - triangulation.refine_global(3-dim+cycle/2); - } - else triangulation.refine_global(1); - - setup_system (); - assemble_system (false); - solve (); - postprocess(); - output_results (cycle); - } +template +void Step51::run () +{ + const bool do_cube = true; + if (!do_cube) + { + GridGenerator::hyper_ball (triangulation); + static const HyperBallBoundary boundary; + triangulation.set_boundary(0, boundary); + triangulation.refine_global(6-2*dim); + } + for (unsigned int cycle=0; cycle<10; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + if (do_cube) + { + triangulation.clear(); + GridGenerator::subdivided_hyper_cube (triangulation, 2+(cycle%2), -1, 1); + triangulation.refine_global(3-dim+cycle/2); + } + else triangulation.refine_global(1); - 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); + setup_system (); + assemble_system (false); + solve (); + postprocess(); + output_results (cycle); + } - 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.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); } @@ -825,7 +823,6 @@ int main (int argc, char** argv) try { using namespace dealii; - using namespace Step51; deallog.depth_console (0); -- 2.39.5