From d50e6556347b6b2e98ded4bf3f6e0448f56b82d5 Mon Sep 17 00:00:00 2001 From: Roland Date: Thu, 25 Apr 2019 17:06:20 +0200 Subject: [PATCH] Made degree of quadrature depending on the degree of the finite element used in the program --- doc/news/changes/major/20190425RolandRichter | 3 ++ examples/step-13/step-13.cc | 2 +- examples/step-14/step-14.cc | 4 +- examples/step-15/step-15.cc | 6 +-- examples/step-16b/doc/results.dox | 56 ++++++++++---------- examples/step-16b/step-16b.cc | 2 +- examples/step-17/step-17.cc | 4 +- examples/step-18/step-18.cc | 4 +- examples/step-23/step-23.cc | 10 ++-- examples/step-24/step-24.cc | 10 ++-- examples/step-25/step-25.cc | 12 +++-- examples/step-29/step-29.cc | 4 +- examples/step-3/doc/results.dox | 2 +- examples/step-3/step-3.cc | 2 +- examples/step-33/step-33.cc | 4 +- examples/step-36/step-36.cc | 2 +- examples/step-4/step-4.cc | 2 +- examples/step-40/step-40.cc | 4 +- examples/step-5/step-5.cc | 2 +- examples/step-51/step-51.cc | 2 +- examples/step-55/step-55.cc | 2 +- examples/step-6/step-6.cc | 4 +- examples/step-7/doc/results.dox | 51 +++++++++--------- examples/step-7/step-7.cc | 17 +++--- examples/step-8/step-8.cc | 4 +- 25 files changed, 111 insertions(+), 104 deletions(-) create mode 100644 doc/news/changes/major/20190425RolandRichter diff --git a/doc/news/changes/major/20190425RolandRichter b/doc/news/changes/major/20190425RolandRichter new file mode 100644 index 0000000000..7477577a75 --- /dev/null +++ b/doc/news/changes/major/20190425RolandRichter @@ -0,0 +1,3 @@ +New: The tutorial examples now all use an interpolation formula with a degree depending on the degree of the finite element, which can be changed dynamically. +
+(Roland Richter, 2019/04/25) diff --git a/examples/step-13/step-13.cc b/examples/step-13/step-13.cc index b37aac8107..71613bf13e 100644 --- a/examples/step-13/step-13.cc +++ b/examples/step-13/step-13.cc @@ -1221,7 +1221,7 @@ namespace Step13 this->triangulation->n_active_cells()); KellyErrorEstimator::estimate( this->dof_handler, - QGauss(3), + QGauss(this->fe->degree + 1), std::map *>(), this->solution, estimated_error_per_cell); diff --git a/examples/step-14/step-14.cc b/examples/step-14/step-14.cc index b79fef2d2f..6e1c319a56 100644 --- a/examples/step-14/step-14.cc +++ b/examples/step-14/step-14.cc @@ -878,7 +878,7 @@ namespace Step14 this->triangulation->n_active_cells()); KellyErrorEstimator::estimate( this->dof_handler, - QGauss(3), + QGauss(this->fe->degree + 1), std::map *>(), this->solution, estimated_error_per_cell); @@ -1530,7 +1530,7 @@ namespace Step14 // Initialize a FEValues object with a quadrature formula, // have abbreviations for the number of quadrature points and shape // functions... - QGauss quadrature(4); + QGauss quadrature(dof_handler.get_fe().degree + 1); FEValues fe_values(dof_handler.get_fe(), quadrature, update_gradients | update_quadrature_points | diff --git a/examples/step-15/step-15.cc b/examples/step-15/step-15.cc index 0e980c9e7e..8f5c01b9ef 100644 --- a/examples/step-15/step-15.cc +++ b/examples/step-15/step-15.cc @@ -234,7 +234,7 @@ namespace Step15 template void MinimalSurfaceProblem::assemble_system() { - const QGauss quadrature_formula(3); + const QGauss quadrature_formula(fe.degree + 1); system_matrix = 0; system_rhs = 0; @@ -376,7 +376,7 @@ namespace Step15 KellyErrorEstimator::estimate( dof_handler, - QGauss(3), + QGauss(fe.degree + 1), std::map *>(), present_solution, estimated_error_per_cell); @@ -504,7 +504,7 @@ namespace Step15 evaluation_point = present_solution; evaluation_point.add(alpha, newton_update); - const QGauss quadrature_formula(3); + const QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, update_gradients | update_quadrature_points | diff --git a/examples/step-16b/doc/results.dox b/examples/step-16b/doc/results.dox index 4f5baa31f4..cd4951e1da 100644 --- a/examples/step-16b/doc/results.dox +++ b/examples/step-16b/doc/results.dox @@ -15,38 +15,38 @@ DEAL:: Number of degrees of freedom: 25 (by level: 8, 25) DEAL:cg::Starting value 0.510691 DEAL:cg::Convergence step 6 value 4.59193e-14 DEAL::Cycle 1 -DEAL:: Number of active cells: 41 -DEAL:: Number of degrees of freedom: 52 (by level: 8, 25, 41) -DEAL:cg::Starting value 0.455356 -DEAL:cg::Convergence step 8 value 3.09682e-13 +DEAL:: Number of active cells: 44 +DEAL:: Number of degrees of freedom: 55 (by level: 8, 25, 45) +DEAL:cg::Starting value 0.440678 +DEAL:cg::Convergence step 8 value 1.99419e-13 DEAL::Cycle 2 -DEAL:: Number of active cells: 80 -DEAL:: Number of degrees of freedom: 100 (by level: 8, 25, 61, 52) -DEAL:cg::Starting value 0.394469 -DEAL:cg::Convergence step 9 value 1.96993e-13 +DEAL:: Number of active cells: 86 +DEAL:: Number of degrees of freedom: 105 (by level: 8, 25, 69, 49) +DEAL:cg::Starting value 0.371855 +DEAL:cg::Convergence step 9 value 1.13984e-13 DEAL::Cycle 3 -DEAL:: Number of active cells: 161 -DEAL:: Number of degrees of freedom: 190 (by level: 8, 25, 77, 160) -DEAL:cg::Starting value 0.322156 -DEAL:cg::Convergence step 9 value 2.94418e-13 +DEAL:: Number of active cells: 170 +DEAL:: Number of degrees of freedom: 200 (by level: 8, 25, 77, 174) +DEAL:cg::Starting value 0.318967 +DEAL:cg::Convergence step 9 value 2.62112e-13 DEAL::Cycle 4 -DEAL:: Number of active cells: 311 -DEAL:: Number of degrees of freedom: 364 (by level: 8, 25, 86, 227, 174) -DEAL:cg::Starting value 0.279667 -DEAL:cg::Convergence step 10 value 3.45746e-13 +DEAL:: Number of active cells: 332 +DEAL:: Number of degrees of freedom: 388 (by level: 8, 25, 86, 231, 204) +DEAL:cg::Starting value 0.276534 +DEAL:cg::Convergence step 10 value 1.69562e-13 DEAL::Cycle 5 -DEAL:: Number of active cells: 593 -DEAL:: Number of degrees of freedom: 667 (by level: 8, 25, 89, 231, 490, 96) -DEAL:cg::Starting value 0.215917 -DEAL:cg::Convergence step 10 value 1.03758e-13 +DEAL:: Number of active cells: 632 +DEAL:: Number of degrees of freedom: 714 (by level: 8, 25, 89, 231, 514, 141) +DEAL:cg::Starting value 0.215300 +DEAL:cg::Convergence step 10 value 6.47463e-13 DEAL::Cycle 6 -DEAL:: Number of active cells: 1127 -DEAL:: Number of degrees of freedom: 1251 (by level: 8, 25, 89, 274, 760, 417, 178) -DEAL:cg::Starting value 0.185906 -DEAL:cg::Convergence step 10 value 3.40351e-13 +DEAL:: Number of active cells: 1202 +DEAL:: Number of degrees of freedom: 1332 (by level: 8, 25, 89, 282, 771, 435, 257) +DEAL:cg::Starting value 0.175848 +DEAL:cg::Convergence step 10 value 1.80664e-13 DEAL::Cycle 7 -DEAL:: Number of active cells: 2144 -DEAL:: Number of degrees of freedom: 2359 (by level: 8, 25, 89, 308, 779, 1262, 817) -DEAL:cg::Starting value 0.141519 -DEAL:cg::Convergence step 10 value 5.74965e-13 +DEAL:: Number of active cells: 2288 +DEAL:: Number of degrees of freedom: 2511 (by level: 8, 25, 89, 318, 779, 1420, 829, 30) +DEAL:cg::Starting value 0.136724 +DEAL:cg::Convergence step 11 value 9.73331e-14 diff --git a/examples/step-16b/step-16b.cc b/examples/step-16b/step-16b.cc index b35c6d4f55..101a570c46 100644 --- a/examples/step-16b/step-16b.cc +++ b/examples/step-16b/step-16b.cc @@ -561,7 +561,7 @@ namespace Step16 KellyErrorEstimator::estimate( dof_handler, - QGauss(3), + QGauss(fe.degree + 1), std::map *>(), solution, estimated_error_per_cell); diff --git a/examples/step-17/step-17.cc b/examples/step-17/step-17.cc index c837d09506..14cd33023b 100644 --- a/examples/step-17/step-17.cc +++ b/examples/step-17/step-17.cc @@ -462,7 +462,7 @@ namespace Step17 template void ElasticProblem::assemble_system() { - QGauss quadrature_formula(2); + QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, update_values | update_gradients | @@ -765,7 +765,7 @@ namespace Step17 Vector local_error_per_cell(triangulation.n_active_cells()); KellyErrorEstimator::estimate( dof_handler, - QGauss(2), + QGauss(fe.degree + 1), std::map *>(), localized_solution, local_error_per_cell, diff --git a/examples/step-18/step-18.cc b/examples/step-18/step-18.cc index 5796801b57..17a231225a 100644 --- a/examples/step-18/step-18.cc +++ b/examples/step-18/step-18.cc @@ -716,7 +716,7 @@ namespace Step18 : triangulation(MPI_COMM_WORLD) , fe(FE_Q(1), dim) , dof_handler(triangulation) - , quadrature_formula(2) + , quadrature_formula(fe.degree + 1) , present_time(0.0) , present_timestep(1.0) , end_time(10.0) @@ -1441,7 +1441,7 @@ namespace Step18 Vector error_per_cell(triangulation.n_active_cells()); KellyErrorEstimator::estimate( dof_handler, - QGauss(2), + QGauss(fe.degree + 1), std::map *>(), incremental_displacement, error_per_cell, diff --git a/examples/step-23/step-23.cc b/examples/step-23/step-23.cc index a120d0991f..5b81444550 100644 --- a/examples/step-23/step-23.cc +++ b/examples/step-23/step-23.cc @@ -373,9 +373,11 @@ namespace Step23 matrix_u.reinit(sparsity_pattern); matrix_v.reinit(sparsity_pattern); - MatrixCreator::create_mass_matrix(dof_handler, QGauss(3), mass_matrix); + MatrixCreator::create_mass_matrix(dof_handler, + QGauss(fe.degree + 1), + mass_matrix); MatrixCreator::create_laplace_matrix(dof_handler, - QGauss(3), + QGauss(fe.degree + 1), laplace_matrix); // The rest of the function is spent on setting vector sizes to the @@ -478,12 +480,12 @@ namespace Step23 VectorTools::project(dof_handler, constraints, - QGauss(3), + QGauss(fe.degree + 1), InitialValuesU(), old_solution_u); VectorTools::project(dof_handler, constraints, - QGauss(3), + QGauss(fe.degree + 1), InitialValuesV(), old_solution_v); diff --git a/examples/step-24/step-24.cc b/examples/step-24/step-24.cc index 737d5cd553..ed0622a618 100644 --- a/examples/step-24/step-24.cc +++ b/examples/step-24/step-24.cc @@ -289,9 +289,11 @@ namespace Step24 mass_matrix.reinit(sparsity_pattern); laplace_matrix.reinit(sparsity_pattern); - MatrixCreator::create_mass_matrix(dof_handler, QGauss(3), mass_matrix); + MatrixCreator::create_mass_matrix(dof_handler, + QGauss(fe.degree + 1), + mass_matrix); MatrixCreator::create_laplace_matrix(dof_handler, - QGauss(3), + QGauss(fe.degree + 1), laplace_matrix); // The second difference, as mentioned, to step-23 is that we need to @@ -336,7 +338,7 @@ namespace Step24 // do something only if that particular face is at the boundary of the // domain. Like this: { - const QGauss quadrature_formula(3); + const QGauss quadrature_formula(fe.degree + 1); FEFaceValues fe_values(fe, quadrature_formula, update_values | update_JxW_values); @@ -472,7 +474,7 @@ namespace Step24 VectorTools::project(dof_handler, constraints, - QGauss(3), + QGauss(fe.degree + 1), InitialValuesP(), old_solution_p); old_solution_v = 0; diff --git a/examples/step-25/step-25.cc b/examples/step-25/step-25.cc index 8f8d4d7b21..7372a47980 100644 --- a/examples/step-25/step-25.cc +++ b/examples/step-25/step-25.cc @@ -308,9 +308,11 @@ namespace Step25 mass_matrix.reinit(sparsity_pattern); laplace_matrix.reinit(sparsity_pattern); - MatrixCreator::create_mass_matrix(dof_handler, QGauss(3), mass_matrix); + MatrixCreator::create_mass_matrix(dof_handler, + QGauss(fe.degree + 1), + mass_matrix); MatrixCreator::create_laplace_matrix(dof_handler, - QGauss(3), + QGauss(fe.degree + 1), laplace_matrix); solution.reinit(dof_handler.n_dofs()); @@ -398,7 +400,7 @@ namespace Step25 Vector &nl_term) const { nl_term = 0; - const QGauss quadrature_formula(3); + const QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, update_values | update_JxW_values | @@ -462,7 +464,7 @@ namespace Step25 const Vector &new_data, SparseMatrix &nl_matrix) const { - QGauss quadrature_formula(3); + QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, update_values | update_JxW_values | @@ -601,7 +603,7 @@ namespace Step25 constraints.close(); VectorTools::project(dof_handler, constraints, - QGauss(3), + QGauss(fe.degree + 1), InitialValues(1, time), solution); } diff --git a/examples/step-29/step-29.cc b/examples/step-29/step-29.cc index cc2d2b4590..805f55801f 100644 --- a/examples/step-29/step-29.cc +++ b/examples/step-29/step-29.cc @@ -554,8 +554,8 @@ namespace Step29 // used. Since our bilinear form involves boundary integrals on // $\Gamma_2$, we also need a quadrature rule for surface integration on // the faces, which are $dim-1$ dimensional: - QGauss quadrature_formula(2); - QGauss face_quadrature_formula(2); + QGauss quadrature_formula(fe.degree + 1); + QGauss face_quadrature_formula(fe.degree + 1); const unsigned int n_q_points = quadrature_formula.size(), n_face_q_points = face_quadrature_formula.size(), diff --git a/examples/step-3/doc/results.dox b/examples/step-3/doc/results.dox index a39a00addc..6e61813cf7 100644 --- a/examples/step-3/doc/results.dox +++ b/examples/step-3/doc/results.dox @@ -192,7 +192,7 @@ suggestions: @code std::cout << "Mean value: " << VectorTools::compute_mean_value (dof_handler, - QGauss<2>(3), + QGauss<2>(fe.degree + 1), solution, 0) << std::endl; diff --git a/examples/step-3/step-3.cc b/examples/step-3/step-3.cc index 1e50eef4db..9f35b73279 100644 --- a/examples/step-3/step-3.cc +++ b/examples/step-3/step-3.cc @@ -272,7 +272,7 @@ void Step3::assemble_system() // 2D. This quadrature formula integrates polynomials of degrees up to three // exactly (in 1D). It is easy to check that this is sufficient for the // present problem: - QGauss<2> quadrature_formula(2); + QGauss<2> quadrature_formula(fe.degree + 1); // And we initialize the object which we have briefly talked about above. It // needs to be told which finite element we want to use, and the quadrature // points and their weights (jointly described by a Quadrature object). As diff --git a/examples/step-33/step-33.cc b/examples/step-33/step-33.cc index 917b51fcd3..a24a5665fc 100644 --- a/examples/step-33/step-33.cc +++ b/examples/step-33/step-33.cc @@ -1381,8 +1381,8 @@ namespace Step33 : mapping() , fe(FE_Q(1), EulerEquations::n_components) , dof_handler(triangulation) - , quadrature(2) - , face_quadrature(2) + , quadrature(fe.degree + 1) + , face_quadrature(fe.degree + 1) , verbose_cout(std::cout, false) { ParameterHandler prm; diff --git a/examples/step-36/step-36.cc b/examples/step-36/step-36.cc index fe6e593e15..6213053568 100644 --- a/examples/step-36/step-36.cc +++ b/examples/step-36/step-36.cc @@ -237,7 +237,7 @@ namespace Step36 template void EigenvalueProblem::assemble_system() { - QGauss quadrature_formula(2); + QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, diff --git a/examples/step-4/step-4.cc b/examples/step-4/step-4.cc index a4b0aa5390..e3d1447b22 100644 --- a/examples/step-4/step-4.cc +++ b/examples/step-4/step-4.cc @@ -308,7 +308,7 @@ void Step4::setup_system() template void Step4::assemble_system() { - QGauss quadrature_formula(2); + QGauss quadrature_formula(fe.degree + 1); // We wanted to have a non-constant right hand side, so we use an object of // the class declared above to generate the necessary data. Since this right diff --git a/examples/step-40/step-40.cc b/examples/step-40/step-40.cc index 0837dd65f6..25837b4726 100644 --- a/examples/step-40/step-40.cc +++ b/examples/step-40/step-40.cc @@ -366,7 +366,7 @@ namespace Step40 { TimerOutput::Scope t(computing_timer, "assembly"); - const QGauss quadrature_formula(3); + const QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, @@ -523,7 +523,7 @@ namespace Step40 Vector estimated_error_per_cell(triangulation.n_active_cells()); KellyErrorEstimator::estimate( dof_handler, - QGauss(3), + QGauss(fe.degree + 1), std::map *>(), locally_relevant_solution, estimated_error_per_cell); diff --git a/examples/step-5/step-5.cc b/examples/step-5/step-5.cc index 8b694d173a..2c0af87460 100644 --- a/examples/step-5/step-5.cc +++ b/examples/step-5/step-5.cc @@ -162,7 +162,7 @@ void Step5::setup_system() template void Step5::assemble_system() { - QGauss quadrature_formula(2); + QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, diff --git a/examples/step-51/step-51.cc b/examples/step-51/step-51.cc index c44dbcda51..36a1d91145 100644 --- a/examples/step-51/step-51.cc +++ b/examples/step-51/step-51.cc @@ -1333,7 +1333,7 @@ namespace Step51 std::map *> neumann_boundary; KellyErrorEstimator::estimate(dof_handler_local, - QGauss(3), + QGauss(fe.degree + 1), neumann_boundary, solution_local, estimated_error_per_cell, diff --git a/examples/step-55/step-55.cc b/examples/step-55/step-55.cc index 2d42ef85f4..10b6ac513f 100644 --- a/examples/step-55/step-55.cc +++ b/examples/step-55/step-55.cc @@ -691,7 +691,7 @@ namespace Step55 FEValuesExtractors::Vector velocities(0); KellyErrorEstimator::estimate( dof_handler, - QGauss(3), + QGauss(fe.degree + 1), std::map *>(), locally_relevant_solution, estimated_error_per_cell, diff --git a/examples/step-6/step-6.cc b/examples/step-6/step-6.cc index b45caaf7bd..18bb3b4768 100644 --- a/examples/step-6/step-6.cc +++ b/examples/step-6/step-6.cc @@ -269,7 +269,7 @@ void Step6::setup_system() template void Step6::assemble_system() { - const QGauss quadrature_formula(3); + const QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, @@ -422,7 +422,7 @@ void Step6::refine_grid() KellyErrorEstimator::estimate( dof_handler, - QGauss(3), + QGauss(fe.degree + 1), std::map *>(), solution, estimated_error_per_cell); diff --git a/examples/step-7/doc/results.dox b/examples/step-7/doc/results.dox index 6671886260..fdc5e9a67b 100644 --- a/examples/step-7/doc/results.dox +++ b/examples/step-7/doc/results.dox @@ -53,15 +53,15 @@ Cycle 8: Number of degrees of freedom: 15912 cycle cells dofs L2 H1 Linfty - 0 64 81 1.576e-01 1.418e+00 2.707e-01 - 1 124 157 4.285e-02 1.285e+00 1.469e-01 - 2 280 341 1.593e-02 7.909e-01 8.034e-02 - 3 577 690 9.359e-03 5.096e-01 2.784e-02 - 4 1099 1264 2.865e-03 3.038e-01 9.822e-03 - 5 2191 2452 1.480e-03 2.106e-01 5.679e-03 - 6 4165 4510 6.907e-04 1.462e-01 2.338e-03 - 7 7915 8440 4.743e-04 1.055e-01 1.442e-03 - 8 15196 15912 1.920e-04 7.468e-02 7.259e-04 + 0 64 81 1.840e+00 2.858e+00 1.835e+00 + 1 124 157 5.190e-02 1.200e+00 1.344e-01 + 2 280 341 1.439e-02 7.892e-01 7.554e-02 + 3 577 690 8.627e-03 5.061e-01 2.805e-02 + 4 1099 1264 3.217e-03 3.030e-01 1.073e-02 + 5 2191 2452 1.445e-03 2.097e-01 5.073e-03 + 6 4165 4510 8.387e-04 1.460e-01 2.013e-03 + 7 7915 8440 7.051e-04 1.053e-01 1.804e-03 + 8 15196 15912 2.774e-04 7.463e-02 6.911e-04 Solving with Q1 elements, global refinement =========================================== @@ -83,18 +83,18 @@ Cycle 4: Number of degrees of freedom: 16641 cycle cells dofs L2 H1 Linfty - 0 64 81 1.576e-01 1.418e+00 2.707e-01 - 1 256 289 4.280e-02 1.285e+00 1.444e-01 - 2 1024 1089 1.352e-02 7.556e-01 7.772e-02 - 3 4096 4225 3.423e-03 3.822e-01 2.332e-02 - 4 16384 16641 8.586e-04 1.917e-01 6.097e-03 - -n cells H1 L2 -0 64 1.418e+00 - - 1.576e-01 - - -1 256 1.285e+00 1.10 0.14 4.280e-02 3.68 1.88 -2 1024 7.556e-01 1.70 0.77 1.352e-02 3.17 1.66 -3 4096 3.822e-01 1.98 0.98 3.423e-03 3.95 1.98 -4 16384 1.917e-01 1.99 1.00 8.586e-04 3.99 2.00 + 0 64 81 1.840e+00 2.858e+00 1.835e+00 + 1 256 289 3.570e-02 1.199e+00 1.307e-01 + 2 1024 1089 1.192e-02 7.565e-01 7.168e-02 + 3 4096 4225 3.047e-03 3.823e-01 2.128e-02 + 4 16384 16641 7.660e-04 1.917e-01 5.554e-03 + +n cells H1 L2 +0 64 2.858e+00 - - 1.840e+00 - - +1 256 1.199e+00 2.38 1.25 3.570e-02 51.54 5.69 +2 1024 7.565e-01 1.58 0.66 1.192e-02 2.99 1.58 +3 4096 3.823e-01 1.98 0.98 3.047e-03 3.91 1.97 +4 16384 1.917e-01 1.99 1.00 7.660e-04 3.98 1.99 Solving with Q2 elements, global refinement =========================================== @@ -169,7 +169,7 @@ cycle cells dofs L2 H1 Linfty 5 2059 9223 7.738e-05 1.974e-02 7.270e-04 6 3913 17887 2.925e-05 8.772e-03 1.463e-04 7 7441 33807 1.024e-05 4.121e-03 8.567e-05 - 8 14212 64731 3.761e-06 2.108e-03 2.167e-05 + 8 14212 64731 3.761e-06 2.108e-03 2.167e-05 @endcode @@ -192,11 +192,8 @@ here). Go ahead and run the program with higher order elements (Q3, Q4, ...). You will notice that assertions in several parts of the code will trigger (for -example in the generation of the filename for the data output). After fixing -these you will not see the correct convergence orders that the theory -predicts. This is because the orders for the quadrature formulas are -hard-coded in this program and this order is not enough for higher order -discretizations. What is a good way to pick the orders dynamically? +example in the generation of the filename for the data output). You might have to address these, +but it should not be very hard to get the program to work!

Convergence Comparison

diff --git a/examples/step-7/step-7.cc b/examples/step-7/step-7.cc index 97c0e7daa9..223460cf7d 100644 --- a/examples/step-7/step-7.cc +++ b/examples/step-7/step-7.cc @@ -541,8 +541,8 @@ namespace Step7 template void HelmholtzProblem::assemble_system() { - QGauss quadrature_formula(3); - QGauss face_quadrature_formula(3); + QGauss quadrature_formula(fe->degree + 1); + QGauss face_quadrature_formula(fe->degree + 1); const unsigned int n_q_points = quadrature_formula.size(); const unsigned int n_face_q_points = face_quadrature_formula.size(); @@ -797,7 +797,7 @@ namespace Step7 KellyErrorEstimator::estimate( dof_handler, - QGauss(3), + QGauss(fe->degree + 1), std::map *>(), solution, estimated_error_per_cell); @@ -853,7 +853,7 @@ namespace Step7 solution, Solution(), difference_per_cell, - QGauss(3), + QGauss(fe->degree + 1), VectorTools::L2_norm); const double L2_error = VectorTools::compute_global_error(triangulation, @@ -871,7 +871,7 @@ namespace Step7 solution, Solution(), difference_per_cell, - QGauss(3), + QGauss(fe->degree + 1), VectorTools::H1_seminorm); const double H1_error = VectorTools::compute_global_error(triangulation, @@ -883,8 +883,9 @@ namespace Step7 // points. Since this depends quite sensitively on the quadrature rule // being used, and since we would like to avoid false results due to // super-convergence effects at some points, we use a special quadrature - // rule that is obtained by iterating the trapezoidal rule five times in - // each space direction. Note that the constructor of the QIterated class + // rule that is obtained by iterating the trapezoidal rule by the degree of + // of the finite element times two plus one in each space direction. + // Note that the constructor of the QIterated class // takes a one-dimensional quadrature rule and a number that tells it how // often it shall use this rule in each space direction. // @@ -893,7 +894,7 @@ namespace Step7 // from the L infinity errors on each cell with a call to // VectorTools::compute_global_error. const QTrapez<1> q_trapez; - const QIterated q_iterated(q_trapez, 5); + const QIterated q_iterated(q_trapez, fe->degree * 2 + 1); VectorTools::integrate_difference(dof_handler, solution, Solution(), diff --git a/examples/step-8/step-8.cc b/examples/step-8/step-8.cc index 4927b72858..3f7318fd21 100644 --- a/examples/step-8/step-8.cc +++ b/examples/step-8/step-8.cc @@ -269,7 +269,7 @@ namespace Step8 template void ElasticProblem::assemble_system() { - QGauss quadrature_formula(2); + QGauss quadrature_formula(fe.degree + 1); FEValues fe_values(fe, quadrature_formula, @@ -490,7 +490,7 @@ namespace Step8 KellyErrorEstimator::estimate( dof_handler, - QGauss(2), + QGauss(fe.degree + 1), std::map *>(), solution, estimated_error_per_cell); -- 2.39.5