From c445c41b16ea4ee88a11edee1bb68a48f88641c2 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 29 Jun 2018 14:28:54 +0200 Subject: [PATCH] examples/step-39: Update indenting and modernize --- examples/step-39/doc/results.dox | 79 ++++++++++++------------ examples/step-39/output.reference.dat | 24 ++++---- examples/step-39/plot_errors.gpl | 11 ++++ examples/step-39/step-39.cc | 80 +++++++++++++------------ include/deal.II/meshworker/functional.h | 5 +- 5 files changed, 111 insertions(+), 88 deletions(-) create mode 100644 examples/step-39/plot_errors.gpl diff --git a/examples/step-39/doc/results.dox b/examples/step-39/doc/results.dox index 1928dbd24c..250ca43f67 100644 --- a/examples/step-39/doc/results.dox +++ b/examples/step-39/doc/results.dox @@ -6,88 +6,93 @@ First, the program produces the usual logfile here stored in deallog. I @code DEAL::Element: FE_DGQ<2>(3) DEAL::Step 0 -DEAL::Triangulation 4 cells, 1 levels -DEAL::DoFHandler 64 dofs, level dofs 64 +DEAL::Triangulation 16 cells, 2 levels +DEAL::DoFHandler 256 dofs, level dofs 64 256 DEAL::Assemble matrix DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve -DEAL:cg::Starting value 27.1275 -DEAL:cg::Convergence step 1 value 1.97998e-14 -DEAL::Error 0.161172 -DEAL::Estimate 1.35839 +DEAL:cg::Starting value 37.4071 +DEAL:cg::Convergence step 13 value 1.64974e-13 +DEAL::energy-error: 0.297419 +DEAL::L2-error: 0.00452447 +DEAL::Estimate 0.990460 DEAL::Writing solution to ... DEAL:: DEAL::Step 1 -DEAL::Triangulation 10 cells, 2 levels -DEAL::DoFHandler 160 dofs, level dofs 64 128 +DEAL::Triangulation 25 cells, 3 levels +DEAL::DoFHandler 400 dofs, level dofs 64 256 192 DEAL::Assemble matrix DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve -DEAL:cg::Starting value 35.5356 -DEAL:cg::Convergence step 14 value 3.21479e-13 -DEAL::Error 0.164760 -DEAL::Estimate 1.08528 +DEAL:cg::Starting value 37.4071 +DEAL:cg::Convergence step 14 value 3.72262e-13 +DEAL::energy-error: 0.258559 +DEAL::L2-error: 0.00288510 +DEAL::Estimate 0.738624 DEAL::Writing solution to ... DEAL:: DEAL::Step 2 -DEAL::Triangulation 16 cells, 2 levels -DEAL::DoFHandler 256 dofs, level dofs 64 256 +DEAL::Triangulation 34 cells, 4 levels +DEAL::DoFHandler 544 dofs, level dofs 64 256 256 128 DEAL::Assemble matrix DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve -DEAL:cg::Starting value 37.0552 -DEAL:cg::Convergence step 14 value 6.05416e-13 -DEAL::Error 0.113503 -DEAL::Estimate 0.990460 +DEAL:cg::Starting value 37.4071 +DEAL:cg::Convergence step 15 value 1.91610e-13 +DEAL::energy-error: 0.189234 +DEAL::L2-error: 0.00147954 +DEAL::Estimate 0.657507 DEAL::Writing solution to ... ... DEAL::Step 10 -DEAL::Triangulation 124 cells, 9 levels -DEAL::DoFHandler 1984 dofs, level dofs 64 256 512 512 256 256 256 256 256 +DEAL::Triangulation 232 cells, 11 levels +DEAL::DoFHandler 3712 dofs, level dofs 64 256 896 768 768 640 512 256 256 256 256 DEAL::Assemble matrix DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve -DEAL:cg::Starting value 38.5798 -DEAL:cg::Convergence step 17 value 2.64999e-13 -DEAL::Error 0.0101278 -DEAL::Estimate 0.0957571 +DEAL:cg::Starting value 51.1571 +DEAL:cg::Convergence step 15 value 7.19599e-13 +DEAL::energy-error: 0.0132475 +DEAL::L2-error: 1.00423e-05 +DEAL::Estimate 0.0470724 DEAL::Writing solution to ... DEAL:: DEAL::Step 11 -DEAL::Triangulation 163 cells, 10 levels -DEAL::DoFHandler 2608 dofs, level dofs 64 256 768 576 512 256 256 256 256 256 +DEAL::Triangulation 322 cells, 12 levels +DEAL::DoFHandler 5152 dofs, level dofs 64 256 1024 1024 896 768 768 640 448 320 320 320 DEAL::Assemble matrix DEAL::Assemble multilevel matrix DEAL::Assemble right hand side DEAL::Solve -DEAL:cg::Starting value 44.1721 -DEAL:cg::Convergence step 17 value 3.18657e-13 -DEAL::Error 0.00716962 -DEAL::Estimate 0.0681646 +DEAL:cg::Starting value 52.2226 +DEAL:cg::Convergence step 15 value 8.15195e-13 +DEAL::energy-error: 0.00934891 +DEAL::L2-error: 5.41095e-06 +DEAL::Estimate 0.0329102 DEAL::Writing solution to ... DEAL:: @endcode This log for instance shows that the number of conjugate gradient -iteration steps is constant at approximately 17. +iteration steps is constant at approximately 15.

Postprocessing of the logfile

Using the perl script postprocess.pl, we extract relevant data into output.dat, which can be used to plot graphs with -gnuplot. The graph above for instance was produced with +gnuplot. The graph above for instance was produced using the gnuplot +script plot_errors.gpl via @code -set style data linespoints -set logscale -set xrange [50:3000] -plot "output.dat" using 2:3 title "error", "" using 2:4 title "estimate", \ - "" using 2:(3000*$2**-1.5) title "3rd order" +perl postprocess.pl deallog &> output.dat +gnuplot plot_errors.gpl @endcode + +Reference data can be found in output.reference.dat. diff --git a/examples/step-39/output.reference.dat b/examples/step-39/output.reference.dat index 35de96d7f0..e544e6e087 100644 --- a/examples/step-39/output.reference.dat +++ b/examples/step-39/output.reference.dat @@ -1,13 +1,13 @@ #step dofs error estimate l2error iterations efficiency order l2order -0 256 2.974190e-01 9.904600e-01 4.524470e-03 14 0.300284 0.000000 0.000000 -1 400 2.585590e-01 7.386240e-01 2.885100e-03 16 0.350055 0.627480 2.016374 -2 544 1.892340e-01 6.575070e-01 1.479540e-03 16 0.287805 2.030277 4.343815 -3 688 1.503620e-01 4.547120e-01 9.005740e-04 16 0.330675 1.958261 4.228028 -4 832 1.053170e-01 3.548790e-01 4.750150e-04 17 0.296769 3.747254 6.731991 -5 1024 7.459510e-02 2.539050e-01 2.362060e-04 17 0.293791 3.322106 6.729380 -6 1216 5.291920e-02 1.836100e-01 1.224410e-04 17 0.288215 3.995439 7.647065 -7 1504 3.756820e-02 1.335280e-01 6.804430e-05 17 0.281351 3.223616 5.527534 -8 1984 2.657860e-02 9.576030e-02 3.998000e-05 17 0.277553 2.498686 3.839747 -9 2608 1.880990e-02 6.816520e-02 1.951950e-05 17 0.275946 2.528429 5.243492 -10 3472 1.328020e-02 4.780300e-02 1.006070e-05 17 0.277811 2.433078 4.632422 -11 4672 9.367220e-03 3.336360e-02 5.766390e-06 17 0.280762 2.351695 3.749897 +0 256 2.974190e-01 9.904600e-01 4.524470e-03 13 0.300284 0.000000 0.000000 +1 400 2.585590e-01 7.386240e-01 2.885100e-03 14 0.350055 0.627480 2.016374 +2 544 1.892340e-01 6.575070e-01 1.479540e-03 15 0.287805 2.030277 4.343815 +3 688 1.503620e-01 4.547120e-01 9.005740e-04 15 0.330675 1.958261 4.228028 +4 832 1.053170e-01 3.548790e-01 4.750150e-04 15 0.296769 3.747254 6.731991 +5 1024 7.459510e-02 2.539050e-01 2.362060e-04 15 0.293791 3.322106 6.729380 +6 1216 5.291920e-02 1.836100e-01 1.224410e-04 15 0.288215 3.995439 7.647065 +7 1504 3.756820e-02 1.335280e-01 6.804430e-05 15 0.281351 3.223616 5.527534 +8 2080 2.653150e-02 9.483870e-02 3.619830e-05 16 0.279754 2.145482 3.893088 +9 2752 1.877250e-02 6.708930e-02 1.797770e-05 16 0.279814 2.471350 4.999857 +10 3712 1.324750e-02 4.707240e-02 1.004230e-05 15 0.281428 2.329774 3.891995 +11 5152 9.348910e-03 3.291020e-02 5.410950e-06 15 0.284073 2.126504 3.772756 diff --git a/examples/step-39/plot_errors.gpl b/examples/step-39/plot_errors.gpl new file mode 100644 index 0000000000..6305e3b701 --- /dev/null +++ b/examples/step-39/plot_errors.gpl @@ -0,0 +1,11 @@ +set logscale +set style data linespoints +set xrange [200:6000] +set xtics (200,400,800,1600,3200) +set xlabel "#dofs" +set yrange [.004:2] +set terminal png +set output "step-39-convergence.png" +plot "output.dat" using 2:3 title "error", \ + "" using 2:4 title "estimate", \ + "" using 2:(5000*$2**-1.5) title "3rd order" diff --git a/examples/step-39/step-39.cc b/examples/step-39/step-39.cc index 910cd3856d..de84658ca7 100644 --- a/examples/step-39/step-39.cc +++ b/examples/step-39/step-39.cc @@ -131,11 +131,11 @@ namespace Step39 MeshWorker::DoFInfo & dinfo, typename MeshWorker::IntegrationInfo &info) const { - const unsigned int deg = info.fe_values(0).get_fe().tensor_degree(); + const unsigned int degree = info.fe_values(0).get_fe().tensor_degree(); LocalIntegrators::Laplace::nitsche_matrix( dinfo.matrix(0, false).matrix, info.fe_values(0), - LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, deg, deg)); + LocalIntegrators::Laplace::compute_penalty(dinfo, dinfo, degree, degree)); } // Interior faces use the interior penalty method @@ -146,7 +146,7 @@ namespace Step39 typename MeshWorker::IntegrationInfo &info1, typename MeshWorker::IntegrationInfo &info2) const { - const unsigned int deg = info1.fe_values(0).get_fe().tensor_degree(); + const unsigned int degree = info1.fe_values(0).get_fe().tensor_degree(); LocalIntegrators::Laplace::ip_matrix( dinfo1.matrix(0, false).matrix, dinfo1.matrix(0, true).matrix, @@ -154,7 +154,8 @@ namespace Step39 dinfo2.matrix(0, false).matrix, info1.fe_values(0), info2.fe_values(0), - LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); + LocalIntegrators::Laplace::compute_penalty( + dinfo1, dinfo2, degree, degree)); } // The second local integrator builds the right hand side. In our example, @@ -194,16 +195,16 @@ namespace Step39 std::vector boundary_values(fe.n_quadrature_points); exact_solution.value_list(fe.get_quadrature_points(), boundary_values); - const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = - 2. * deg * (deg + 1) * dinfo.face->measure() / dinfo.cell->measure(); + const unsigned int degree = fe.get_fe().tensor_degree(); + const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() / + dinfo.cell->measure(); for (unsigned k = 0; k < fe.n_quadrature_points; ++k) for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) local_vector(i) += - (-fe.shape_value(i, k) * penalty * boundary_values[k] + - (fe.normal_vector(k) * fe.shape_grad(i, k)) * boundary_values[k]) * - fe.JxW(k); + (-fe.shape_value(i, k) * penalty // -sigma * v_i(x_k) + + fe.normal_vector(k) * fe.shape_grad(i, k)) // n * grad v_i(x_k) + * boundary_values[k] * fe.JxW(k); // u^D(x_k) * dx } @@ -268,13 +269,15 @@ namespace Step39 const std::vector &uh = info.values[0][0]; - const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = - 2. * deg * (deg + 1) * dinfo.face->measure() / dinfo.cell->measure(); + const unsigned int degree = fe.get_fe().tensor_degree(); + const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() / + dinfo.cell->measure(); for (unsigned k = 0; k < fe.n_quadrature_points; ++k) - dinfo.value(0) += penalty * (boundary_values[k] - uh[k]) * - (boundary_values[k] - uh[k]) * fe.JxW(k); + { + const double diff = boundary_values[k] - uh[k]; + dinfo.value(0) += penalty * diff * diff * fe.JxW(k); + } dinfo.value(0) = std::sqrt(dinfo.value(0)); } @@ -294,18 +297,18 @@ namespace Step39 const std::vector> &Duh1 = info1.gradients[0][0]; const std::vector> &Duh2 = info2.gradients[0][0]; - const unsigned int deg = fe.get_fe().tensor_degree(); + const unsigned int degree = fe.get_fe().tensor_degree(); const double penalty1 = - deg * (deg + 1) * dinfo1.face->measure() / dinfo1.cell->measure(); + degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure(); const double penalty2 = - deg * (deg + 1) * dinfo2.face->measure() / dinfo2.cell->measure(); + degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure(); const double penalty = penalty1 + penalty2; const double h = dinfo1.face->measure(); for (unsigned k = 0; k < fe.n_quadrature_points; ++k) { - double diff1 = uh1[k] - uh2[k]; - double diff2 = + const double diff1 = uh1[k] - uh2[k]; + const double diff2 = fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k]; dinfo1.value(0) += (penalty * diff1 * diff1 + h * diff2 * diff2) * fe.JxW(k); @@ -398,9 +401,9 @@ namespace Step39 const std::vector &uh = info.values[0][0]; - const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = - 2. * deg * (deg + 1) * dinfo.face->measure() / dinfo.cell->measure(); + const unsigned int degree = fe.get_fe().tensor_degree(); + const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() / + dinfo.cell->measure(); for (unsigned k = 0; k < fe.n_quadrature_points; ++k) { @@ -422,16 +425,16 @@ namespace Step39 const std::vector &uh1 = info1.values[0][0]; const std::vector &uh2 = info2.values[0][0]; - const unsigned int deg = fe.get_fe().tensor_degree(); + const unsigned int degree = fe.get_fe().tensor_degree(); const double penalty1 = - deg * (deg + 1) * dinfo1.face->measure() / dinfo1.cell->measure(); + degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure(); const double penalty2 = - deg * (deg + 1) * dinfo2.face->measure() / dinfo2.cell->measure(); + degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure(); const double penalty = penalty1 + penalty2; for (unsigned k = 0; k < fe.n_quadrature_points; ++k) { - double diff = uh1[k] - uh2[k]; + const double diff = uh1[k] - uh2[k]; dinfo1.value(0) += (penalty * diff * diff) * fe.JxW(k); } dinfo1.value(0) = std::sqrt(dinfo1.value(0)); @@ -776,7 +779,10 @@ namespace Step39 { // The results of the estimator are stored in a vector with one entry per // cell. Since cells in deal.II are not numbered, we have to create our - // own numbering in order to use this vector. + // own numbering in order to use this vector. For the assembler used below + // the information in which component of a vector the result is stored is + // transmitted by the user_index variable for each cell. We need to set this + // numering up here. // // On the other hand, somebody might have used the user indices // already. So, let's be good citizens and save them before tampering with @@ -786,11 +792,8 @@ namespace Step39 estimates.block(0).reinit(triangulation.n_active_cells()); unsigned int i = 0; - for (typename Triangulation::active_cell_iterator cell = - triangulation.begin_active(); - cell != triangulation.end(); - ++cell, ++i) - cell->set_user_index(i); + for (const auto &cell : triangulation.active_cell_iterators()) + cell->set_user_index(i++); // This starts like before, MeshWorker::IntegrationInfoBox info_box; @@ -860,12 +863,12 @@ namespace Step39 BlockVector errors(2); errors.block(0).reinit(triangulation.n_active_cells()); errors.block(1).reinit(triangulation.n_active_cells()); + + std::vector old_user_indices; + triangulation.save_user_indices(old_user_indices); unsigned int i = 0; - for (typename Triangulation::active_cell_iterator cell = - triangulation.begin_active(); - cell != triangulation.end(); - ++cell, ++i) - cell->set_user_index(i); + for (const auto &cell : triangulation.active_cell_iterators()) + cell->set_user_index(i++); MeshWorker::IntegrationInfoBox info_box; const unsigned int n_gauss_points = @@ -899,6 +902,7 @@ namespace Step39 info_box, integrator, assembler); + triangulation.load_user_indices(old_user_indices); deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl; deallog << "L2-error: " << errors.block(1).l2_norm() << std::endl; diff --git a/include/deal.II/meshworker/functional.h b/include/deal.II/meshworker/functional.h index 69dab9eb79..208a7c8179 100644 --- a/include/deal.II/meshworker/functional.h +++ b/include/deal.II/meshworker/functional.h @@ -94,7 +94,10 @@ namespace MeshWorker /** * Compute cell and face contributions of one or several functionals, - * typically for error estimates. + * typically for error estimates. The information in which component the + * result is stored for a given cell or face is transmitted by its + * user_index variable. Hence, you need to make sure to set these variables + * appropriately before using this class. * * @ingroup MeshWorker * @author Guido Kanschat, 2009 -- 2.39.5