From 9ba819ba716e693e4f192c8247509f2aaa7730cc Mon Sep 17 00:00:00 2001 From: kanschat Date: Thu, 28 Jul 2011 16:43:36 +0000 Subject: [PATCH] compute correct energy norm and L2 norm git-svn-id: https://svn.dealii.org/trunk@23980 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-39/postprocess.pl | 8 +- deal.II/examples/step-39/step-39.cc | 224 ++++++++++++++++++++++-- 2 files changed, 212 insertions(+), 20 deletions(-) diff --git a/deal.II/examples/step-39/postprocess.pl b/deal.II/examples/step-39/postprocess.pl index 5d60b44f60..567b14e685 100644 --- a/deal.II/examples/step-39/postprocess.pl +++ b/deal.II/examples/step-39/postprocess.pl @@ -8,19 +8,21 @@ use strict; my $step; # The iteration step in the adaptive loop my @dofs; # The number of degrees of freedom in each step -my @error; # The error of the solution +my @error; # The energy error of the solution +my @l2error; # The L2-error of the solution my @estimate; # The a posteriori error estimate my @steps; # The number of multigrid iteration steps while(<>) { $step = $1 if m/DEAL::Step\s*(\d+)/; $dofs[$step] = $1 if m/DEAL::DoFHandler\s*(\d+)/; - $error[$step] = $1 if m/DEAL::Error\s*(\S+)/; + $error[$step] = $1 if m/DEAL::energy-error:\s*(\S+)/; + $l2error[$step] = $1 if m/DEAL::L2-error:\s*(\S+)/; $estimate[$step] = $1 if m/DEAL::Estimate\s*(\S+)/; $steps[$step] = $1 if m/DEAL:\w+::Convergence step\s*(\S+)/; } for (my $i=0;$i<=$step;++$i) { - printf "%-3d\t%-7d\t%g\t%g\t%d\n", $i, $dofs[$i], $error[$i], $estimate[$i], $steps[$i]; + printf "%-3d\t%-7d\t%g\t%g\t%g\t%d\n", $i, $dofs[$i], $error[$i], $estimate[$i], $l2error[$i], $steps[$i]; } diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index 67a51756c0..5d3e16a537 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -321,6 +321,153 @@ void Estimator::face(MeshWorker::DoFInfo& dinfo1, dinfo2.value(0) = dinfo1.value(0); } + // Finally we have an integrator for + // the error. Since the energy norm + // for discontinuous Galerkin + // problems not only involves the + // difference of the gradient inside + // the cells, but also the jump terms + // across faces and at the boundary, + // we cannot just use + // VectorTools::integrate_difference(). + // Instead, we use the MeshWorker + // interface to compute the error + // ourselves. + + // There are several different ways + // to define this energy norm, but + // all of them are equivalent to each + // other uniformly with mesh size + // (some not uniformly with + // polynomial degree). Here, we + // choose + // @f[ + // \|u\|_{1,h} = \sum_{K\in \mathbb + // T_h} \|\nabla u\|_K^2 + // + \sum_{F \in F_h^i} + // 4\sigma_F\|\{\!\{ u \mathbf + // n\}\!\}\|^2_F + // + \sum_{F \in F_h^b} 2\sigma_F\|u\|^2_F + // @f] + +template +class ErrorIntegrator : public Subscriptor +{ + public: + static void cell(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationInfo& info); + static void boundary(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationInfo& info); + static void face(MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + typename MeshWorker::IntegrationInfo& info1, + typename MeshWorker::IntegrationInfo& info2); +}; + + // Here we have the integration on + // cells. There is currently no good + // interfce in MeshWorker that would + // allow us to access values of + // regular functions in the + // quadrature points. Thus, we have + // to create the vectors for the + // exact function's values and + // gradients inside the cell + // integrator. After that, everything + // is as before and we just add up + // the squares of the differences. + + // Additionally to computing the error + // in the energy norm, we use the + // capability of the mesh worker to + // compute two functionals at the + // same time and compute the + // L2-error in the + // same loop. Obviously, this one + // does not have any jump terms and + // only appears in the integration on + // cells. +template +void ErrorIntegrator::cell( + MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationInfo& info) +{ + const FEValuesBase& fe = info.fe_values(); + std::vector > exact_gradients(fe.n_quadrature_points); + std::vector exact_values(fe.n_quadrature_points); + + exact_solution.gradient_list(fe.get_quadrature_points(), exact_gradients); + exact_solution.value_list(fe.get_quadrature_points(), exact_values); + + const std::vector >& Duh = info.gradients[0][0]; + const std::vector& uh = info.values[0][0]; + + for (unsigned k=0;k +void ErrorIntegrator::boundary( + MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationInfo& info) +{ + const FEValuesBase& fe = info.fe_values(); + + std::vector exact_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), exact_values); + + 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(); + + for (unsigned k=0;k +void ErrorIntegrator::face( + MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + typename MeshWorker::IntegrationInfo& info1, + typename MeshWorker::IntegrationInfo& info2) +{ + const FEValuesBase& fe = info1.fe_values(); + 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 double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); + const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); + const double penalty = penalty1 + penalty2; + + for (unsigned k=0;k::solve() solver.solve(matrix, solution, right_hand_side, preconditioner); } - // Here we compare our finite element - // solution with the (known) exact - // solution and compute the mean - // quadratic error of the gradient. -template -void -Step39::error() -{ - const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+2; - Vector cell_errors(triangulation.n_active_cells()); - - QGauss quadrature(n_gauss_points); - VectorTools::integrate_difference(mapping, dof_handler, solution, exact_solution, - cell_errors, quadrature, VectorTools::H1_seminorm); - deallog << "Error " << cell_errors.l2_norm() << std::endl; -} - // Another clone of the assemble // function. The big difference to @@ -867,6 +997,66 @@ Step39::estimate() return estimates.block(0).l2_norm(); } + // Here we compare our finite element + // solution with the (known) exact + // solution and compute the mean + // quadratic error of the gradient + // and the function itself. This + // function is a clone of the + // estimation function right above. + + // Since we compute the error in the + // energy and the + // L2-norm, + // respectively, our block vector + // needs two blocks here. +template +void +Step39::error() +{ + BlockVector errors(2); + errors.block(0).reinit(triangulation.n_active_cells()); + errors.block(1).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); + + MeshWorker::IntegrationInfoBox info_box; + const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; + info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); + + NamedData* > solution_data; + solution_data.add(&solution, "solution"); + + info_box.cell_selector.add("solution", true, true, false); + info_box.boundary_selector.add("solution", true, false, false); + info_box.face_selector.add("solution", true, false, false); + + info_box.add_update_flags_cell(update_quadrature_points); + info_box.add_update_flags_boundary(update_quadrature_points); + info_box.initialize(fe, mapping, solution_data); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::CellsAndFaces assembler; + NamedData* > out_data; + BlockVector* est = &errors; + out_data.add(est, "cells"); + assembler.initialize(out_data, false); + + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + &ErrorIntegrator::cell, + &ErrorIntegrator::boundary, + &ErrorIntegrator::face, + assembler); + + deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl; + deallog << "L2-error: " << errors.block(1).l2_norm() << std::endl; +} + // Some graphical output template -- 2.39.5