From 33c18670a4761144a1e55075ccd5f376f548b2b8 Mon Sep 17 00:00:00 2001 From: kanschat Date: Mon, 25 Feb 2013 15:25:17 +0000 Subject: [PATCH] add test for adaptive dg multigrid git-svn-id: https://svn.dealii.org/trunk@28553 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/multigrid/step-39.cc | 741 ++++++++++++++++++++++++++++ tests/multigrid/step-39/cmp/generic | 74 +++ 2 files changed, 815 insertions(+) create mode 100644 tests/multigrid/step-39.cc create mode 100644 tests/multigrid/step-39/cmp/generic diff --git a/tests/multigrid/step-39.cc b/tests/multigrid/step-39.cc new file mode 100644 index 0000000000..f611ca7289 --- /dev/null +++ b/tests/multigrid/step-39.cc @@ -0,0 +1,741 @@ +/* Author: Guido Kanschat, Texas A&M University, 2009 */ + +/* $Id: step-39.cc 28552 2013-02-25 15:18:15Z kanschat $ */ +/* */ +/* Copyright (C) 2010-2013 by the deal.II authors */ +/* */ +/* This file is subject to QPL and may not be distributed */ +/* without copyright and license information. Please refer */ +/* to the file deal.II/doc/license.html for the text and */ +/* further information on this license. */ + +#include "../tests.h" +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include + +namespace Step39 +{ + using namespace dealii; + + Functions::SlitSingularityFunction<2> exact_solution; + + + + + template + class MatrixIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void MatrixIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0,false).matrix, info.fe_values()); + } + + + template + void MatrixIntegrator::boundary( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + const unsigned int deg = 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)); + } + + template + void MatrixIntegrator::face( + MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const unsigned int deg = 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,true).matrix, dinfo2.matrix(0,false).matrix, + info1.fe_values(0), info2.fe_values(0), + LocalIntegrators::Laplace::compute_penalty(dinfo1, dinfo2, deg, deg)); + } + + template + class RHSIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void RHSIntegrator::cell(MeshWorker::DoFInfo &, typename MeshWorker::IntegrationInfo &) const + {} + + + template + void RHSIntegrator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + Vector &local_vector = dinfo.vector(0).block(0); + + 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(); + + for (unsigned k=0; k + void RHSIntegrator::face(MeshWorker::DoFInfo &, + MeshWorker::DoFInfo &, + typename MeshWorker::IntegrationInfo &, + typename MeshWorker::IntegrationInfo &) const + {} + + + template + class Estimator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void Estimator::cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + const std::vector > &DDuh = info.hessians[0][0]; + for (unsigned k=0; kdiameter() * trace(DDuh[k]); + dinfo.value(0) += t*t * fe.JxW(k); + } + dinfo.value(0) = std::sqrt(dinfo.value(0)); + } + + template + void Estimator::boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const + { + const FEValuesBase &fe = info.fe_values(); + + std::vector boundary_values(fe.n_quadrature_points); + exact_solution.value_list(fe.get_quadrature_points(), boundary_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 Estimator::face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const + { + const FEValuesBase &fe = info1.fe_values(); + const std::vector &uh1 = info1.values[0][0]; + const std::vector &uh2 = info2.values[0][0]; + 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 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; + const double h = dinfo1.face->measure(); + + for (unsigned k=0; k + class ErrorIntegrator : public MeshWorker::LocalIntegrator + { + public: + void cell(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void boundary(MeshWorker::DoFInfo &dinfo, typename MeshWorker::IntegrationInfo &info) const; + void face(MeshWorker::DoFInfo &dinfo1, + MeshWorker::DoFInfo &dinfo2, + typename MeshWorker::IntegrationInfo &info1, + typename MeshWorker::IntegrationInfo &info2) const; + }; + + + template + void ErrorIntegrator::cell( + MeshWorker::DoFInfo &dinfo, + typename MeshWorker::IntegrationInfo &info) const + { + 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 + { + 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 + { + 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 + class InteriorPenaltyProblem + { + public: + typedef MeshWorker::IntegrationInfo CellInfo; + + InteriorPenaltyProblem(const FiniteElement &fe); + + void run(unsigned int n_steps); + + private: + void setup_system (); + void assemble_matrix (); + void assemble_mg_matrix (); + void assemble_right_hand_side (); + void error (); + double estimate (); + void solve (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + const MappingQ1 mapping; + const FiniteElement &fe; + MGDoFHandler mg_dof_handler; + DoFHandler &dof_handler; + + SparsityPattern sparsity; + SparseMatrix matrix; + Vector solution; + Vector right_hand_side; + BlockVector estimates; + + MGLevelObject mg_sparsity; + MGLevelObject > mg_matrix; + + MGLevelObject mg_sparsity_dg_interface; + MGLevelObject > mg_matrix_dg_down; + MGLevelObject > mg_matrix_dg_up; + }; + + + template + InteriorPenaltyProblem::InteriorPenaltyProblem(const FiniteElement &fe) + : + mapping(), + fe(fe), + mg_dof_handler(triangulation), + dof_handler(mg_dof_handler), + estimates(1) + { + GridGenerator::hyper_cube_slit(triangulation, -1, 1); + } + + + template + void + InteriorPenaltyProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + unsigned int n_dofs = dof_handler.n_dofs(); + solution.reinit(n_dofs); + right_hand_side.reinit(n_dofs); + + CompressedSparsityPattern c_sparsity(n_dofs); + DoFTools::make_flux_sparsity_pattern(dof_handler, c_sparsity); + sparsity.copy_from(c_sparsity); + matrix.reinit(sparsity); + + const unsigned int n_levels = triangulation.n_levels(); + mg_matrix.resize(0, n_levels-1); + mg_matrix.clear(); + mg_matrix_dg_up.resize(0, n_levels-1); + mg_matrix_dg_up.clear(); + mg_matrix_dg_down.resize(0, n_levels-1); + mg_matrix_dg_down.clear(); + mg_sparsity.resize(0, n_levels-1); + mg_sparsity_dg_interface.resize(0, n_levels-1); + + for (unsigned int level=mg_sparsity.min_level(); + level<=mg_sparsity.max_level(); ++level) + { + CompressedSparsityPattern c_sparsity(mg_dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern(mg_dof_handler, c_sparsity, level); + mg_sparsity[level].copy_from(c_sparsity); + mg_matrix[level].reinit(mg_sparsity[level]); + + if (level>0) + { + CompressedSparsityPattern ci_sparsity; + ci_sparsity.reinit(mg_dof_handler.n_dofs(level-1), mg_dof_handler.n_dofs(level)); + MGTools::make_flux_sparsity_pattern_edge(mg_dof_handler, ci_sparsity, level); + mg_sparsity_dg_interface[level].copy_from(ci_sparsity); + mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]); + mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]); + } + } + } + + + template + void + InteriorPenaltyProblem::assemble_matrix() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::MatrixSimple > assembler; + assembler.initialize(matrix); + + MatrixIntegrator integrator; + MeshWorker::integration_loop( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + } + + + template + void + InteriorPenaltyProblem::assemble_mg_matrix() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(mg_dof_handler); + + MeshWorker::Assembler::MGMatrixSimple > assembler; + assembler.initialize(mg_matrix); + assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); + + MatrixIntegrator integrator; + MeshWorker::integration_loop ( + mg_dof_handler.begin(), mg_dof_handler.end(), + dof_info, info_box, + integrator, assembler); + } + + + template + void + InteriorPenaltyProblem::assemble_right_hand_side() + { + MeshWorker::IntegrationInfoBox info_box; + UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; + info_box.add_update_flags_all(update_flags); + info_box.initialize(fe, mapping); + + MeshWorker::DoFInfo dof_info(dof_handler); + + MeshWorker::Assembler::ResidualSimple > assembler; + NamedData* > data; + Vector *rhs = &right_hand_side; + data.add(rhs, "RHS"); + assembler.initialize(data); + + RHSIntegrator integrator; + MeshWorker::integration_loop( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + right_hand_side *= -1.; + } + + + template + void + InteriorPenaltyProblem::solve() + { + SolverControl control(1000, 1.e-12); + SolverCG > solver(control); + + MGTransferPrebuilt > mg_transfer; + mg_transfer.build_matrices(mg_dof_handler); + + FullMatrix coarse_matrix; + coarse_matrix.copy_from (mg_matrix[0]); + MGCoarseGridHouseholder > mg_coarse; + mg_coarse.initialize(coarse_matrix); + + GrowingVectorMemory > mem; + typedef PreconditionSOR > RELAXATION; + mg::SmootherRelaxation > + mg_smoother; + RELAXATION::AdditionalData smoother_data(1.); + mg_smoother.initialize(mg_matrix, smoother_data); + + mg_smoother.set_steps(2); + mg_smoother.set_symmetric(true); + mg_smoother.set_variable(false); + + MGMatrix, Vector > mgmatrix(&mg_matrix); + MGMatrix, Vector > mgdown(&mg_matrix_dg_down); + MGMatrix, Vector > mgup(&mg_matrix_dg_up); + + Multigrid > mg(mg_dof_handler, mgmatrix, + mg_coarse, mg_transfer, + mg_smoother, mg_smoother); + mg.set_edge_flux_matrices(mgdown, mgup); + + PreconditionMG, + MGTransferPrebuilt > > + preconditioner(mg_dof_handler, mg, mg_transfer); + solver.solve(matrix, solution, right_hand_side, preconditioner); + } + + + template + double + InteriorPenaltyProblem::estimate() + { + std::vector old_user_indices; + triangulation.save_user_indices(old_user_indices); + + 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); + + 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", false, false, true); + info_box.boundary_selector.add("solution", true, true, false); + info_box.face_selector.add("solution", true, true, false); + + 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 = &estimates; + out_data.add(est, "cells"); + assembler.initialize(out_data, false); + + Estimator integrator; + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + triangulation.load_user_indices(old_user_indices); + return estimates.block(0).l2_norm(); + } + + + template + void + InteriorPenaltyProblem::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); + + ErrorIntegrator integrator; + MeshWorker::integration_loop ( + dof_handler.begin_active(), dof_handler.end(), + dof_info, info_box, + integrator, assembler); + + deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl; + deallog << "L2-error: " << errors.block(1).l2_norm() << std::endl; + } + + + template + void InteriorPenaltyProblem::output_results (const unsigned int cycle) const + { + char *fn = new char[100]; + sprintf(fn, "sol-%02d", cycle); + + std::string filename(fn); + filename += ".gnuplot"; + deallog << "Writing solution to <" << filename << ">..." + << std::endl << std::endl; + std::ofstream gnuplot_output (filename.c_str()); + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (solution, "u"); + data_out.add_data_vector (estimates.block(0), "est"); + + data_out.build_patches (); + + data_out.write_gnuplot(gnuplot_output); + } + + template + void + InteriorPenaltyProblem::run(unsigned int n_steps) + { + deallog << "Element: " << fe.get_name() << std::endl; + for (unsigned int s=0; s fe1(2); + InteriorPenaltyProblem<2> test1(fe1); + test1.run(6); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/tests/multigrid/step-39/cmp/generic b/tests/multigrid/step-39/cmp/generic new file mode 100644 index 0000000000..169e36cf07 --- /dev/null +++ b/tests/multigrid/step-39/cmp/generic @@ -0,0 +1,74 @@ +JobId unknown Mon Feb 25 16:24:28 2013 +DEAL::Element: FE_DGQ<2>(2) +DEAL::Step 0 +DEAL::Triangulation 16 cells, 2 levels +DEAL::DoFHandler 144 dofs, level dofs 36 144 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 23.0730 +DEAL:cg::Convergence step 11 value 5.78162e-14 +DEAL::energy-error: 0.439211 +DEAL::L2-error: 0.0109342 +DEAL::Estimate 0.979555 +DEAL::Step 1 +DEAL::Triangulation 25 cells, 3 levels +DEAL::DoFHandler 225 dofs, level dofs 36 144 108 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 23.0730 +DEAL:cg::Convergence step 11 value 9.37499e-13 +DEAL::energy-error: 0.332582 +DEAL::L2-error: 0.00548996 +DEAL::Estimate 0.838060 +DEAL::Step 2 +DEAL::Triangulation 37 cells, 4 levels +DEAL::DoFHandler 333 dofs, level dofs 36 144 144 108 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 23.0730 +DEAL:cg::Convergence step 12 value 3.41649e-13 +DEAL::energy-error: 0.237705 +DEAL::L2-error: 0.00250869 +DEAL::Estimate 0.611449 +DEAL::Step 3 +DEAL::Triangulation 58 cells, 5 levels +DEAL::DoFHandler 522 dofs, level dofs 36 144 180 180 144 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 23.4898 +DEAL:cg::Convergence step 12 value 5.58516e-13 +DEAL::energy-error: 0.170860 +DEAL::L2-error: 0.00120805 +DEAL::Estimate 0.452418 +DEAL::Step 4 +DEAL::Triangulation 85 cells, 6 levels +DEAL::DoFHandler 765 dofs, level dofs 36 144 360 180 180 108 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 25.9490 +DEAL:cg::Convergence step 13 value 2.71720e-13 +DEAL::energy-error: 0.122755 +DEAL::L2-error: 0.000602816 +DEAL::Estimate 0.332525 +DEAL::Step 5 +DEAL::Triangulation 130 cells, 7 levels +DEAL::DoFHandler 1170 dofs, level dofs 36 144 432 432 180 180 144 +DEAL::Assemble matrix +DEAL::Assemble multilevel matrix +DEAL::Assemble right hand side +DEAL::Solve +DEAL:cg::Starting value 27.1421 +DEAL:cg::Convergence step 13 value 2.50045e-13 +DEAL::energy-error: 0.0869445 +DEAL::L2-error: 0.000292783 +DEAL::Estimate 0.236647 -- 2.39.5