]> https://gitweb.dealii.org/ - dealii.git/commitdiff
examples/step-39: Update indenting and modernize
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Fri, 29 Jun 2018 12:28:54 +0000 (14:28 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Fri, 29 Jun 2018 12:31:14 +0000 (14:31 +0200)
examples/step-39/doc/results.dox
examples/step-39/output.reference.dat
examples/step-39/plot_errors.gpl [new file with mode: 0644]
examples/step-39/step-39.cc
include/deal.II/meshworker/functional.h

index 1928dbd24c8c3917609c993676efdc88f272060e..250ca43f67216f51eef0e41c01293f6ac0b8db61 100644 (file)
@@ -6,88 +6,93 @@ First, the program produces the usual logfile here stored in <tt>deallog</tt>. 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 <sol-00.gnuplot>...
 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 <sol-01.gnuplot>...
 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 <sol-02.gnuplot>...
 
 ...
 
 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 <sol-10.gnuplot>...
 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 <sol-11.gnuplot>...
 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.
 
 <h2>Postprocessing of the logfile</h2>
 
 <img src="https://www.dealii.org/images/steps/developer/step-39-convergence.png" alt="">
 Using the perl script <tt>postprocess.pl</tt>, we extract relevant
 data into <tt>output.dat</tt>, which can be used to plot graphs with
-<tt>gnuplot</tt>. The graph above for instance was produced with
+<tt>gnuplot</tt>. The graph above for instance was produced using the gnuplot
+script <tt>plot_errors.gpl</tt> 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 <tt>output.reference.dat</tt>.
index 35de96d7f058602c4e36c4ee2d8991193c3aa90a..e544e6e087a3086ed0f083482a9b88983976a343 100644 (file)
@@ -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 (file)
index 0000000..6305e3b
--- /dev/null
@@ -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"
index 910cd3856df36ecfa0222228b21f0779c03fe293..de84658ca754a8c9299ed69a7d67908dcafaa815 100644 (file)
@@ -131,11 +131,11 @@ namespace Step39
     MeshWorker::DoFInfo<dim> &                 dinfo,
     typename MeshWorker::IntegrationInfo<dim> &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<dim> &info1,
     typename MeshWorker::IntegrationInfo<dim> &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<double> 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<double> &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<Tensor<1, dim>> &Duh1 = info1.gradients[0][0];
     const std::vector<Tensor<1, dim>> &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<double> &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<double> &uh1 = info1.values[0][0];
     const std::vector<double> &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<dim>::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<dim> info_box;
@@ -860,12 +863,12 @@ namespace Step39
     BlockVector<double> errors(2);
     errors.block(0).reinit(triangulation.n_active_cells());
     errors.block(1).reinit(triangulation.n_active_cells());
+
+    std::vector<unsigned int> old_user_indices;
+    triangulation.save_user_indices(old_user_indices);
     unsigned int i = 0;
-    for (typename Triangulation<dim>::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<dim> 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;
index 69dab9eb795844ac08019c996690d4ed38e6a0cb..208a7c81795f67488fb3a5b692455d4e8a415944 100644 (file)
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.