From 555051c48c48cb98bf2f060220812e09a09bcea4 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 3 Jan 2021 11:35:35 +0100 Subject: [PATCH] Convert step-74 to a simplex test --- tests/simplex/step-74.cc | 1052 +++++++++++++++++ .../step-74.with_simplex_support=on.output | 30 + 2 files changed, 1082 insertions(+) create mode 100644 tests/simplex/step-74.cc create mode 100644 tests/simplex/step-74.with_simplex_support=on.output diff --git a/tests/simplex/step-74.cc b/tests/simplex/step-74.cc new file mode 100644 index 0000000000..7c9ce0553c --- /dev/null +++ b/tests/simplex/step-74.cc @@ -0,0 +1,1052 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Step-74 on a simplex mesh. Following incompatible modifications had to be +// made: +// - The functions GridGenerator::hyper_cube() and GridGenerator::hyper_cube() +// are only working for hypercube meshes. Use these functions to create +// temporary meshes and convert them to simplex mesh. +// - Local refinement is replaced by global refinement. + + +#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 + +#include "../tests.h" + +#define USE_SIMPLEX + +namespace Step74 +{ + using namespace dealii; + + enum class Test_Case + { + convergence_rate, + l_singularity + }; + + template + class SmoothSolution : public Function + { + public: + SmoothSolution() + : Function() + {} + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const override; + virtual Tensor<1, dim> + gradient(const Point & point, + const unsigned int component = 0) const override; + }; + + template + void + SmoothSolution::value_list(const std::vector> &points, + std::vector & values, + const unsigned int /*component*/) const + { + using numbers::PI; + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = + std::sin(2. * PI * points[i][0]) * std::sin(2. * PI * points[i][1]); + } + + template + Tensor<1, dim> + SmoothSolution::gradient(const Point &point, + const unsigned int /*component*/) const + { + Tensor<1, dim> return_value; + using numbers::PI; + return_value[0] = + 2. * PI * std::cos(2. * PI * point[0]) * std::sin(2. * PI * point[1]); + return_value[1] = + 2. * PI * std::sin(2. * PI * point[0]) * std::cos(2. * PI * point[1]); + return return_value; + } + + template + class SmoothRightHandSide : public Function + { + public: + SmoothRightHandSide() + : Function() + {} + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int /*component*/) const override; + }; + + template + void + SmoothRightHandSide::value_list(const std::vector> &points, + std::vector & values, + const unsigned int /*component*/) const + { + using numbers::PI; + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = 8. * PI * PI * std::sin(2. * PI * points[i][0]) * + std::sin(2. * PI * points[i][1]); + } + + template + class SingularRightHandSide : public Function + { + public: + SingularRightHandSide() + : Function() + {} + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int /*component*/) const override; + + private: + Functions::LSingularityFunction ref; + }; + + template + void + SingularRightHandSide::value_list(const std::vector> &points, + std::vector & values, + const unsigned int /*component*/) const + { + for (unsigned int i = 0; i < values.size(); ++i) + values[i] = -ref.laplacian(points[i]); + } + + template + void + get_function_jump(const FEInterfaceValues &fe_iv, + const Vector & solution, + std::vector & jump) + { + const unsigned n_q = fe_iv.n_quadrature_points; + std::array, 2> face_values; + jump.resize(n_q); + for (unsigned i = 0; i < 2; ++i) + { + face_values[i].resize(n_q); + fe_iv.get_fe_face_values(i).get_function_values(solution, + face_values[i]); + } + for (unsigned int q = 0; q < n_q; ++q) + jump[q] = face_values[0][q] - face_values[1][q]; + } + + template + void + get_function_gradient_jump(const FEInterfaceValues &fe_iv, + const Vector & solution, + std::vector> & gradient_jump) + { + const unsigned n_q = fe_iv.n_quadrature_points; + std::vector> face_gradients[2]; + gradient_jump.resize(n_q); + for (unsigned i = 0; i < 2; ++i) + { + face_gradients[i].resize(n_q); + fe_iv.get_fe_face_values(i).get_function_gradients(solution, + face_gradients[i]); + } + for (unsigned int q = 0; q < n_q; ++q) + gradient_jump[q] = face_gradients[0][q] - face_gradients[1][q]; + } + + double + compute_penalty(const unsigned int fe_degree, + const double cell_extend_left, + const double cell_extend_right) + { + const double degree = std::max(1., static_cast(fe_degree)); + return degree * (degree + 1.) * 0.5 * + (1. / cell_extend_left + 1. / cell_extend_right); + } + + + struct CopyDataFace + { + FullMatrix cell_matrix; + std::vector joint_dof_indices; + std::array values; + std::array cell_indices; + }; + + struct CopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + std::vector face_data; + double value; + unsigned int cell_index; + template + void + reinit(const Iterator &cell, unsigned int dofs_per_cell) + { + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_rhs.reinit(dofs_per_cell); + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + } + }; + + template + class SIPGLaplace + { + public: + SIPGLaplace(const Test_Case &test_case); + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + void + output_results(const unsigned int cycle) const; + + void + compute_errors(); + void + compute_error_estimate(); + double + compute_energy_norm(); + + Triangulation triangulation; + const unsigned degree; + +#ifdef USE_SIMPLEX + const Simplex::QGauss quadrature; + const Simplex::QGauss face_quadrature; + const Simplex::QGauss quadrature_2; + const Simplex::QGauss face_quadrature_2; + const MappingFE mapping; + const Simplex::FE_DGP fe; +#else + const QGauss quadrature; + const QGauss face_quadrature; + const QGauss quadrature_2; + const QGauss face_quadrature_2; + const MappingQ1 mapping; + const FE_DGQ fe; +#endif + + using ScratchData = MeshWorker::ScratchData; + + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + + Vector estimated_error_square_per_cell; + Vector energy_norm_square_per_cell; + + ConvergenceTable convergence_table; + + const double diffusion_coefficient = 1.; + + const Test_Case test_case; + + std::unique_ptr> exact_solution; + std::unique_ptr> rhs_function; + }; + + template + SIPGLaplace::SIPGLaplace(const Test_Case &test_case) + : degree(2) + , quadrature(degree + 1) + , face_quadrature(degree + 1) + , quadrature_2(degree + 1) + , face_quadrature_2(degree + 1) +#ifdef USE_SIMPLEX + , mapping(Simplex::FE_P(1)) +#else + , mapping() +#endif + , fe(degree) + , dof_handler(triangulation) + , test_case(test_case) + { + if (test_case == Test_Case::convergence_rate) + { + exact_solution = std::make_unique>(); + rhs_function = std::make_unique>(); + } + + else if (test_case == Test_Case::l_singularity) + { + exact_solution = std::make_unique(); + rhs_function = std::make_unique>(); + } + else + AssertThrow(false, ExcNotImplemented()); + } + + template + void + SIPGLaplace::setup_system() + { + dof_handler.distribute_dofs(fe); + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_flux_sparsity_pattern(dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + } + + template + void + SIPGLaplace::assemble_system() + { + const auto cell_worker = + [&](const auto &cell, auto &scratch_data, auto ©_data) { + const FEValues &fe_v = scratch_data.reinit(cell); + const unsigned int dofs_per_cell = fe_v.dofs_per_cell; + copy_data.reinit(cell, dofs_per_cell); + + const auto & q_points = scratch_data.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const std::vector &JxW = scratch_data.get_JxW_values(); + + std::vector rhs(n_q_points); + rhs_function->value_list(q_points, rhs); + + for (unsigned int point = 0; point < n_q_points; ++point) + for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) += + diffusion_coefficient * // nu + fe_v.shape_grad(i, point) * // grad v_h + fe_v.shape_grad(j, point) * // grad u_h + JxW[point]; // dx + + copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h + rhs[point] * // f + JxW[point]; // dx + } + }; + + const auto boundary_worker = [&](const auto & cell, + const unsigned int &face_no, + auto & scratch_data, + auto & copy_data) { + const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); + + const auto & q_points = scratch_data.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const unsigned int dofs_per_cell = fe_fv.dofs_per_cell; + + const std::vector & JxW = scratch_data.get_JxW_values(); + const std::vector> &normals = + scratch_data.get_normal_vectors(); + + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); + + const double extent1 = cell->measure() / cell->face(face_no)->measure(); + const double penalty = compute_penalty(degree, extent1, extent1); + + for (unsigned int point = 0; point < n_q_points; ++point) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + copy_data.cell_matrix(i, j) += + (-diffusion_coefficient * // - nu + fe_fv.shape_value(i, point) * // v_h + (fe_fv.shape_grad(j, point) * // (grad u_h . + normals[point]) // n) + + - diffusion_coefficient * // - nu + (fe_fv.shape_grad(i, point) * // (grad v_h . + normals[point]) * // n) + fe_fv.shape_value(j, point) // u_h + + + diffusion_coefficient * penalty * // + nu sigma + fe_fv.shape_value(i, point) * // v_h + fe_fv.shape_value(j, point) // u_h + + ) * + JxW[point]; // dx + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + copy_data.cell_rhs(i) += + (-diffusion_coefficient * // - nu + (fe_fv.shape_grad(i, point) * // (grad v_h . + normals[point]) * // n) + g[point] // g + + + + diffusion_coefficient * penalty * // + nu sigma + fe_fv.shape_value(i, point) * g[point] // v_h g + + ) * + JxW[point]; // dx + } + }; + + const auto face_worker = [&](const auto & cell, + const unsigned int &f, + const unsigned int &sf, + const auto & ncell, + const unsigned int &nf, + const unsigned int &nsf, + auto & scratch_data, + auto & copy_data) { + const FEInterfaceValues &fe_iv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + + const auto & q_points = fe_iv.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + + copy_data.face_data.emplace_back(); + CopyDataFace & copy_data_face = copy_data.face_data.back(); + const unsigned int n_dofs_face = fe_iv.n_current_interface_dofs(); + copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices(); + copy_data_face.cell_matrix.reinit(n_dofs_face, n_dofs_face); + + const std::vector & JxW = fe_iv.get_JxW_values(); + const std::vector> &normals = fe_iv.get_normal_vectors(); + + const double extent1 = cell->measure() / cell->face(f)->measure(); + const double extent2 = ncell->measure() / ncell->face(nf)->measure(); + const double penalty = compute_penalty(degree, extent1, extent2); + + for (unsigned int point = 0; point < n_q_points; ++point) + { + for (unsigned int i = 0; i < n_dofs_face; ++i) + for (unsigned int j = 0; j < n_dofs_face; ++j) + copy_data_face.cell_matrix(i, j) += + (-diffusion_coefficient * // - nu + fe_iv.jump(i, point) * // [v_h] + (fe_iv.average_gradient(j, point) * // ({grad u_h} . + normals[point]) // n) + + - diffusion_coefficient * // - nu + (fe_iv.average_gradient(i, point) * // (grad v_h . + normals[point]) * // n) + fe_iv.jump(j, point) // [u_h] + + + diffusion_coefficient * penalty * // + nu sigma + fe_iv.jump(i, point) * // [v_h] + fe_iv.jump(j, point) // [u_h] + + ) * + JxW[point]; // dx + } + }; + + AffineConstraints constraints; + constraints.close(); + const auto copier = [&](const auto &c) { + constraints.distribute_local_to_global(c.cell_matrix, + c.cell_rhs, + c.local_dof_indices, + system_matrix, + system_rhs); + + for (auto &cdf : c.face_data) + { + constraints.distribute_local_to_global(cdf.cell_matrix, + cdf.joint_dof_indices, + system_matrix); + } + }; + + UpdateFlags cell_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + UpdateFlags face_flags = update_values | update_gradients | + update_quadrature_points | update_normal_vectors | + update_JxW_values; + + ScratchData scratch_data( + mapping, fe, quadrature, cell_flags, face_quadrature, face_flags); + CopyData cd; + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + cd, + MeshWorker::assemble_own_cells | + MeshWorker::assemble_boundary_faces | + MeshWorker::assemble_own_interior_faces_once, + boundary_worker, + face_worker); + } + + template + void + SIPGLaplace::solve() + { + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix); + A_direct.vmult(solution, system_rhs); + } + + template + void + SIPGLaplace::output_results(const unsigned int cycle) const + { + std::string filename = "sol_Q" + Utilities::int_to_string(degree, 1) + "-" + + Utilities::int_to_string(cycle, 2) + ".vtu"; + std::ofstream output(filename); + + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "u", DataOut::type_dof_data); + data_out.build_patches(mapping); + data_out.write_vtu(output); + } + + template + void + SIPGLaplace::compute_error_estimate() + { + estimated_error_square_per_cell.reinit(triangulation.n_active_cells()); + + const auto cell_worker = + [&](const auto &cell, auto &scratch_data, auto ©_data) { + const FEValues &fe_v = scratch_data.reinit(cell); + + copy_data.cell_index = cell->active_cell_index(); + + const auto & q_points = fe_v.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const std::vector &JxW = fe_v.get_JxW_values(); + + std::vector> hessians(n_q_points); + fe_v.get_function_hessians(solution, hessians); + + std::vector rhs(n_q_points); + rhs_function->value_list(q_points, rhs); + + const double hk = cell->diameter(); + double residual_norm_square = 0; + + for (unsigned int point = 0; point < n_q_points; ++point) + { + const double residual = + rhs[point] + diffusion_coefficient * trace(hessians[point]); + residual_norm_square += residual * residual * JxW[point]; + } + copy_data.value = hk * hk * residual_norm_square; + }; + + const auto boundary_worker = [&](const auto & cell, + const unsigned int &face_no, + auto & scratch_data, + auto & copy_data) { + const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); + + const auto & q_points = fe_fv.get_quadrature_points(); + const unsigned n_q_points = q_points.size(); + + const std::vector &JxW = fe_fv.get_JxW_values(); + + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); + + std::vector sol_u(n_q_points); + fe_fv.get_function_values(solution, sol_u); + + const double extent1 = cell->measure() / cell->face(face_no)->measure(); + const double penalty = compute_penalty(degree, extent1, extent1); + + double difference_norm_square = 0.; + for (unsigned int point = 0; point < q_points.size(); ++point) + { + const double diff = (g[point] - sol_u[point]); + difference_norm_square += diff * diff * JxW[point]; + } + copy_data.value += penalty * difference_norm_square; + }; + + const auto face_worker = [&](const auto & cell, + const unsigned int &f, + const unsigned int &sf, + const auto & ncell, + const unsigned int &nf, + const unsigned int &nsf, + auto & scratch_data, + auto & copy_data) { + const FEInterfaceValues &fe_iv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + + copy_data.face_data.emplace_back(); + CopyDataFace ©_data_face = copy_data.face_data.back(); + + copy_data_face.cell_indices[0] = cell->active_cell_index(); + copy_data_face.cell_indices[1] = ncell->active_cell_index(); + + const std::vector & JxW = fe_iv.get_JxW_values(); + const std::vector> &normals = fe_iv.get_normal_vectors(); + + const auto & q_points = fe_iv.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + + std::vector jump(n_q_points); + get_function_jump(fe_iv, solution, jump); + + std::vector> grad_jump(n_q_points); + get_function_gradient_jump(fe_iv, solution, grad_jump); + + const double h = cell->face(f)->diameter(); + + const double extent1 = cell->measure() / cell->face(f)->measure(); + const double extent2 = ncell->measure() / ncell->face(nf)->measure(); + const double penalty = compute_penalty(degree, extent1, extent2); + + double flux_jump_square = 0; + double u_jump_square = 0; + for (unsigned int point = 0; point < n_q_points; ++point) + { + u_jump_square += jump[point] * jump[point] * JxW[point]; + const double flux_jump = grad_jump[point] * normals[point]; + flux_jump_square += + diffusion_coefficient * flux_jump * flux_jump * JxW[point]; + } + copy_data_face.values[0] = + 0.5 * h * (flux_jump_square + penalty * u_jump_square); + copy_data_face.values[1] = copy_data_face.values[0]; + }; + + const auto copier = [&](const auto ©_data) { + if (copy_data.cell_index != numbers::invalid_unsigned_int) + estimated_error_square_per_cell[copy_data.cell_index] += + copy_data.value; + for (auto &cdf : copy_data.face_data) + for (unsigned int j = 0; j < 2; ++j) + estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j]; + }; + + UpdateFlags cell_flags = + update_hessians | update_quadrature_points | update_JxW_values; + UpdateFlags face_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values | + update_normal_vectors; + + ScratchData scratch_data( + mapping, fe, quadrature, cell_flags, face_quadrature, face_flags); + + CopyData cd; + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + cd, + MeshWorker::assemble_own_cells | + MeshWorker::assemble_own_interior_faces_once | + MeshWorker::assemble_boundary_faces, + boundary_worker, + face_worker); + } + + template + double + SIPGLaplace::compute_energy_norm() + { + energy_norm_square_per_cell.reinit(triangulation.n_active_cells()); + + const auto cell_worker = + [&](const auto &cell, auto &scratch_data, auto ©_data) { + const FEValues &fe_v = scratch_data.reinit(cell); + + copy_data.cell_index = cell->active_cell_index(); + + const auto & q_points = fe_v.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + const std::vector &JxW = fe_v.get_JxW_values(); + + std::vector> grad_u(n_q_points); + fe_v.get_function_gradients(solution, grad_u); + + std::vector> grad_exact(n_q_points); + exact_solution->gradient_list(q_points, grad_exact); + + double norm_square = 0; + for (unsigned int point = 0; point < n_q_points; ++point) + { + norm_square += + (grad_u[point] - grad_exact[point]).norm_square() * JxW[point]; + } + copy_data.value = norm_square; + }; + + const auto boundary_worker = [&](const auto & cell, + const unsigned int &face_no, + auto & scratch_data, + auto & copy_data) { + const FEFaceValuesBase &fe_fv = scratch_data.reinit(cell, face_no); + + const auto & q_points = fe_fv.get_quadrature_points(); + const unsigned n_q_points = q_points.size(); + + const std::vector &JxW = fe_fv.get_JxW_values(); + + std::vector g(n_q_points); + exact_solution->value_list(q_points, g); + + std::vector sol_u(n_q_points); + fe_fv.get_function_values(solution, sol_u); + + const double extent1 = cell->measure() / cell->face(face_no)->measure(); + const double penalty = compute_penalty(degree, extent1, extent1); + + double difference_norm_square = 0.; + for (unsigned int point = 0; point < q_points.size(); ++point) + { + const double diff = (g[point] - sol_u[point]); + difference_norm_square += diff * diff * JxW[point]; + } + copy_data.value += penalty * difference_norm_square; + }; + + const auto face_worker = [&](const auto & cell, + const unsigned int &f, + const unsigned int &sf, + const auto & ncell, + const unsigned int &nf, + const unsigned int &nsf, + auto & scratch_data, + auto & copy_data) { + const FEInterfaceValues &fe_iv = + scratch_data.reinit(cell, f, sf, ncell, nf, nsf); + + copy_data.face_data.emplace_back(); + CopyDataFace ©_data_face = copy_data.face_data.back(); + + copy_data_face.cell_indices[0] = cell->active_cell_index(); + copy_data_face.cell_indices[1] = ncell->active_cell_index(); + + const std::vector &JxW = fe_iv.get_JxW_values(); + + const auto & q_points = fe_iv.get_quadrature_points(); + const unsigned int n_q_points = q_points.size(); + + std::vector jump(n_q_points); + get_function_jump(fe_iv, solution, jump); + + const double extent1 = cell->measure() / cell->face(f)->measure(); + const double extent2 = ncell->measure() / ncell->face(nf)->measure(); + const double penalty = compute_penalty(degree, extent1, extent2); + + double u_jump_square = 0; + for (unsigned int point = 0; point < n_q_points; ++point) + { + u_jump_square += jump[point] * jump[point] * JxW[point]; + } + copy_data_face.values[0] = 0.5 * penalty * u_jump_square; + copy_data_face.values[1] = copy_data_face.values[0]; + }; + + const auto copier = [&](const auto ©_data) { + if (copy_data.cell_index != numbers::invalid_unsigned_int) + energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value; + for (auto &cdf : copy_data.face_data) + for (unsigned int j = 0; j < 2; ++j) + energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j]; + }; + + UpdateFlags cell_flags = + update_gradients | update_quadrature_points | update_JxW_values; + UpdateFlags face_flags = + update_values | update_quadrature_points | update_JxW_values; + + ScratchData scratch_data( + mapping, fe, quadrature_2, cell_flags, face_quadrature_2, face_flags); + + CopyData cd; + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + cd, + MeshWorker::assemble_own_cells | + MeshWorker::assemble_own_interior_faces_once | + MeshWorker::assemble_boundary_faces, + boundary_worker, + face_worker); + const double energy_error = + std::sqrt(energy_norm_square_per_cell.l1_norm()); + return energy_error; + } + + template + void + SIPGLaplace::refine_grid() + { + const double refinement_fraction = 0.1; + + GridRefinement::refine_and_coarsen_fixed_number( + triangulation, estimated_error_square_per_cell, refinement_fraction, 0.); + + triangulation.execute_coarsening_and_refinement(); + } + + template + void + SIPGLaplace::compute_errors() + { + double L2_error, H1_error; + + { + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + *(exact_solution.get()), + difference_per_cell, + quadrature_2, + VectorTools::L2_norm); + + L2_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::L2_norm); + } + + { + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + *(exact_solution.get()), + difference_per_cell, + quadrature_2, + VectorTools::H1_seminorm); + + H1_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::H1_seminorm); + } + + convergence_table.add_value("L2", L2_error); + convergence_table.add_value("H1", H1_error); + const double energy_error = compute_energy_norm(); + convergence_table.add_value("Energy", energy_error); + + deallog << " Error in the L2 norm : " << L2_error << std::endl + << " Error in the H1 seminorm : " << H1_error << std::endl + << " Error in the energy norm : " << energy_error << std::endl; + } + + template + void + SIPGLaplace::run() + { + unsigned int max_cycle = 3; + for (unsigned int cycle = 0; cycle < max_cycle; ++cycle) + { + deallog << "Cycle " << cycle << std::endl; + + switch (test_case) + { + case Test_Case::convergence_rate: + { + if (cycle == 0) + { +#ifdef USE_SIMPLEX + Triangulation temp; + GridGenerator::hyper_cube(temp); + GridGenerator::convert_hypercube_to_simplex_mesh( + temp, triangulation); +#else + GridGenerator::hyper_cube(triangulation); +#endif + triangulation.refine_global(2); + } + else + { + triangulation.refine_global(); + } + break; + } + case Test_Case::l_singularity: + { + if (cycle == 0) + { +#ifdef USE_SIMPLEX + Triangulation temp; + GridGenerator::hyper_L(temp); + GridGenerator::convert_hypercube_to_simplex_mesh( + temp, triangulation); +#else + GridGenerator::hyper_L(triangulation); +#endif + triangulation.refine_global(3); + } + else + { + triangulation.refine_global(); + } + break; + } + default: + { + Assert(false, ExcNotImplemented()); + } + } + deallog << " Number of active cells : " + << triangulation.n_active_cells() << std::endl; + setup_system(); + + deallog << " Number of degrees of freedom : " << dof_handler.n_dofs() + << std::endl; + + assemble_system(); + solve(); + output_results(cycle); + { + convergence_table.add_value("cycle", cycle); + convergence_table.add_value("cells", triangulation.n_active_cells()); + convergence_table.add_value("dofs", dof_handler.n_dofs()); + } + compute_errors(); + + if (test_case == Test_Case::l_singularity) + { + compute_error_estimate(); + deallog << " Estimated error : " + << std::sqrt(estimated_error_square_per_cell.l1_norm()) + << std::endl; + + convergence_table.add_value( + "Estimator", + std::sqrt(estimated_error_square_per_cell.l1_norm())); + } + deallog << std::endl; + } + { + convergence_table.set_precision("L2", 3); + convergence_table.set_precision("H1", 3); + convergence_table.set_precision("Energy", 3); + + convergence_table.set_scientific("L2", true); + convergence_table.set_scientific("H1", true); + convergence_table.set_scientific("Energy", true); + + if (test_case == Test_Case::l_singularity) + { + convergence_table.set_precision("Estimator", 3); + convergence_table.set_scientific("Estimator", true); + } + if (test_case == Test_Case::convergence_rate) + { + convergence_table.evaluate_convergence_rates( + "L2", ConvergenceTable::reduction_rate_log2); + convergence_table.evaluate_convergence_rates( + "H1", ConvergenceTable::reduction_rate_log2); + } + deallog << "degree = " << degree << std::endl; + convergence_table.write_text( + deallog.get_file_stream(), + TableHandler::TextOutputFormat::org_mode_table); + } + } +} // namespace Step74 + + +int +main() +{ + initlog(); + + try + { + using namespace dealii; + using namespace Step74; + Test_Case test_case = Test_Case::l_singularity; + SIPGLaplace<2> problem(test_case); + problem.run(); + } + 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/simplex/step-74.with_simplex_support=on.output b/tests/simplex/step-74.with_simplex_support=on.output new file mode 100644 index 0000000000..ccab956206 --- /dev/null +++ b/tests/simplex/step-74.with_simplex_support=on.output @@ -0,0 +1,30 @@ + +DEAL::Cycle 0 +DEAL:: Number of active cells : 1536 +DEAL:: Number of degrees of freedom : 9216 +DEAL:: Error in the L2 norm : 0.000101269 +DEAL:: Error in the H1 seminorm : 0.0222535 +DEAL:: Error in the energy norm : 0.0343502 +DEAL:: Estimated error : 0.0798085 +DEAL:: +DEAL::Cycle 1 +DEAL:: Number of active cells : 6144 +DEAL:: Number of degrees of freedom : 36864 +DEAL:: Error in the L2 norm : 3.19139e-05 +DEAL:: Error in the H1 seminorm : 0.0140189 +DEAL:: Error in the energy norm : 0.0216394 +DEAL:: Estimated error : 0.0502598 +DEAL:: +DEAL::Cycle 2 +DEAL:: Number of active cells : 24576 +DEAL:: Number of degrees of freedom : 147456 +DEAL:: Error in the L2 norm : 1.00600e-05 +DEAL:: Error in the H1 seminorm : 0.00883139 +DEAL:: Error in the energy norm : 0.0136320 +DEAL:: Estimated error : 0.0316564 +DEAL:: +DEAL::degree = 2 +| cycle | cells | dofs | L2 | H1 | Energy | Estimator | +| 0 | 1536 | 9216 | 1.013e-04 | 2.225e-02 | 3.435e-02 | 7.981e-02 | +| 1 | 6144 | 36864 | 3.191e-05 | 1.402e-02 | 2.164e-02 | 5.026e-02 | +| 2 | 24576 | 147456 | 1.006e-05 | 8.831e-03 | 1.363e-02 | 3.166e-02 | -- 2.39.5