@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>.
#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
--- /dev/null
+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"
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
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,
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,
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
}
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));
}
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);
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)
{
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));
{
// 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
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;
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 =
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;
/**
* 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
projects / dealii.git / commitdiff