From 09306d9b6a8a4f2a9638bed01e7f50dbf62c103f Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 3 Jan 2021 00:26:40 +0100 Subject: [PATCH] Convert step-12 to a simplex test --- tests/simplex/step-12.cc | 243 ++++---- .../step-12.with_simplex_support=on.output | 39 +- tests/simplex/step-12a.cc | 561 ++++++++++++++++++ .../step-12a.with_simplex_support=on.output | 15 + 4 files changed, 702 insertions(+), 156 deletions(-) create mode 100644 tests/simplex/step-12a.cc create mode 100644 tests/simplex/step-12a.with_simplex_support=on.output diff --git a/tests/simplex/step-12.cc b/tests/simplex/step-12.cc index 29484a3aa2..94d64a0856 100644 --- a/tests/simplex/step-12.cc +++ b/tests/simplex/step-12.cc @@ -14,7 +14,13 @@ // --------------------------------------------------------------------- -// Step-12 with tetrahedron mesh. +// Step-02 on a simplex mesh. Following incompatible modifications had to be +// made: +// - The function GridGenerator::hyper_cube() is only working for hypercube +// meshes. Use this function to create a temporary mesh and convert it to +// simplex mesh. +// - Local refinement is replaced by global refinement. + #include #include @@ -26,6 +32,7 @@ #include #include #include +#include #include #include @@ -41,20 +48,12 @@ #include #include +#include + #include #include #include -#include "../tests.h" - -// Finally, the new include file for using the mesh_loop from the MeshWorker -// framework -#include - -#include - -#include - #include #include #include @@ -62,8 +61,6 @@ #include #include -//#define HEX - namespace Step12 { using namespace dealii; @@ -110,20 +107,21 @@ namespace Step12 wind_field(0) = -p(1); wind_field(1) = p(0); - if (wind_field.norm() > 1e-6) + if (wind_field.norm() > 1e-10) wind_field /= wind_field.norm(); return wind_field; } + template struct ScratchData { ScratchData(const Mapping & mapping, const FiniteElement & fe, - const Quadrature & quad, - const Quadrature &quad_face, + const Quadrature & quadrature, + const Quadrature &quadrature_face, const UpdateFlags update_flags = update_values | update_gradients | update_quadrature_points | @@ -131,8 +129,11 @@ namespace Step12 const UpdateFlags interface_update_flags = update_values | update_gradients | update_quadrature_points | update_JxW_values | update_normal_vectors) - : fe_values(mapping, fe, quad, update_flags) - , fe_interface_values(mapping, fe, quad_face, interface_update_flags) + : fe_values(mapping, fe, quadrature, update_flags) + , fe_interface_values(mapping, + fe, + quadrature_face, + interface_update_flags) {} @@ -141,12 +142,10 @@ namespace Step12 scratch_data.fe_values.get_fe(), scratch_data.fe_values.get_quadrature(), scratch_data.fe_values.get_update_flags()) - , fe_interface_values( - scratch_data.fe_values - .get_mapping(), // TODO: implement for fe_interface_values - scratch_data.fe_values.get_fe(), - scratch_data.fe_interface_values.get_quadrature(), - scratch_data.fe_interface_values.get_update_flags()) + , fe_interface_values(scratch_data.fe_interface_values.get_mapping(), + scratch_data.fe_interface_values.get_fe(), + scratch_data.fe_interface_values.get_quadrature(), + scratch_data.fe_interface_values.get_update_flags()) {} FEValues fe_values; @@ -203,21 +202,14 @@ namespace Step12 void output_results(const unsigned int cycle) const; - Triangulation triangulation; -#ifdef HEX - const MappingQ1 mapping; -#else - Simplex::FE_P fe_mapping; + Triangulation triangulation; const MappingFE mapping; -#endif - // Furthermore we want to use DG elements. -#ifdef HEX - FE_DGQ fe; -#else - Simplex::FE_DGP fe; -#endif - DoFHandler dof_handler; + const Simplex::FE_DGP fe; + DoFHandler dof_handler; + + const Simplex::QGauss quadrature; + const Simplex::QGauss quadrature_face; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; @@ -229,15 +221,11 @@ namespace Step12 template AdvectionProblem::AdvectionProblem() -#ifdef HEX - : mapping() - , fe(2) -#else - : fe_mapping(1) - , mapping(fe_mapping) - , fe(2) -#endif + : mapping(Simplex::FE_P(1)) + , fe(1) , dof_handler(triangulation) + , quadrature(fe.tensor_degree() + 1) + , quadrature_face(fe.tensor_degree() + 1) {} @@ -263,10 +251,11 @@ namespace Step12 using Iterator = typename DoFHandler::active_cell_iterator; const BoundaryValues boundary_function; - auto cell_worker = [&](const Iterator & cell, - ScratchData &scratch_data, - CopyData & copy_data) { - const unsigned int n_dofs = scratch_data.fe_values.get_fe().dofs_per_cell; + const auto cell_worker = [&](const Iterator & cell, + ScratchData &scratch_data, + CopyData & copy_data) { + const unsigned int n_dofs = + scratch_data.fe_values.get_fe().n_dofs_per_cell(); copy_data.reinit(cell, n_dofs); scratch_data.fe_values.reinit(cell); @@ -290,10 +279,10 @@ namespace Step12 } }; - auto boundary_worker = [&](const Iterator & cell, - const unsigned int &face_no, - ScratchData & scratch_data, - CopyData & copy_data) { + const auto boundary_worker = [&](const Iterator & cell, + const unsigned int &face_no, + ScratchData & scratch_data, + CopyData & copy_data) { scratch_data.fe_interface_values.reinit(cell, face_no); const FEFaceValuesBase &fe_face = scratch_data.fe_interface_values.get_fe_face_values(0); @@ -330,14 +319,14 @@ namespace Step12 } }; - auto face_worker = [&](const Iterator & cell, - const unsigned int &f, - const unsigned int &sf, - const Iterator & ncell, - const unsigned int &nf, - const unsigned int &nsf, - ScratchData & scratch_data, - CopyData & copy_data) { + const auto face_worker = [&](const Iterator & cell, + const unsigned int &f, + const unsigned int &sf, + const Iterator & ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData & scratch_data, + CopyData & copy_data) { FEInterfaceValues &fe_iv = scratch_data.fe_interface_values; fe_iv.reinit(cell, f, sf, ncell, nf, nsf); const auto &q_points = fe_iv.get_quadrature_points(); @@ -367,9 +356,9 @@ namespace Step12 } }; - AffineConstraints constraints; + const AffineConstraints constraints; - auto copier = [&](const CopyData &c) { + const auto copier = [&](const CopyData &c) { constraints.distribute_local_to_global(c.cell_matrix, c.cell_rhs, c.local_dof_indices, @@ -384,19 +373,7 @@ namespace Step12 } }; - const unsigned int degree = dof_handler.get_fe().degree; - -#ifdef HEX - QGauss quad(degree + 1); - - QGauss face_quad(degree + 1); -#else - Simplex::QGauss quad(degree + 1); - - Simplex::QGauss face_quad(degree + 1); -#endif - - ScratchData scratch_data(mapping, fe, quad, face_quad); + ScratchData scratch_data(mapping, fe, quadrature, quadrature_face); CopyData copy_data; MeshWorker::mesh_loop(dof_handler.begin_active(), @@ -421,12 +398,12 @@ namespace Step12 PreconditionBlockSSOR> preconditioner; - preconditioner.initialize(system_matrix, fe.dofs_per_cell); + preconditioner.initialize(system_matrix, fe.n_dofs_per_cell()); solver.solve(system_matrix, solution, right_hand_side, preconditioner); - deallog << " Solver converged in " << solver_control.last_step() - << " iterations." << std::endl; + std::cout << " Solver converged in " << solver_control.last_step() + << " iterations." << std::endl; } @@ -434,24 +411,7 @@ namespace Step12 void AdvectionProblem::refine_grid() { - Vector gradient_indicator(triangulation.n_active_cells()); - - DerivativeApproximation::approximate_gradient(mapping, - dof_handler, - solution, - gradient_indicator); - - unsigned int cell_no = 0; - for (const auto &cell : dof_handler.active_cell_iterators()) - gradient_indicator(cell_no++) *= - std::pow(cell->diameter(), 1 + 1.0 * dim / 2); - - GridRefinement::refine_and_coarsen_fixed_number(triangulation, - gradient_indicator, - 0.3, - 0.1); - - triangulation.execute_coarsening_and_refinement(); + triangulation.refine_global(); } @@ -459,27 +419,17 @@ namespace Step12 void AdvectionProblem::output_results(const unsigned int cycle) const { -#if false -# ifdef HEX - const std::string filename = - dim == 2 ? ("step12-quad-" + std::to_string(cycle) + ".vtk") : - ("step12-hex-" + std::to_string(cycle) + ".vtk"); -# else - const std::string filename = - dim == 2 ? ("step12-tri-" + std::to_string(cycle) + ".vtk") : - ("step12-tet-" + std::to_string(cycle) + ".vtk"); -# endif - deallog << " Writing solution to <" << filename << ">" << std::endl; + const std::string filename = "solution-" + std::to_string(cycle) + ".vtk"; + std::cout << " Writing solution to <" << filename << ">" << std::endl; 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, 2); + data_out.build_patches(mapping); data_out.write_vtk(output); -#endif { Vector values(triangulation.n_active_cells()); @@ -488,13 +438,13 @@ namespace Step12 solution, Functions::ZeroFunction(), values, - QGauss(fe.degree + 1), + quadrature, VectorTools::Linfty_norm); const double l_infty = VectorTools::compute_global_error(triangulation, values, VectorTools::Linfty_norm); - deallog << " L-infinity norm: " << l_infty << std::endl; + std::cout << " L-infinity norm: " << l_infty << std::endl; } } @@ -503,36 +453,27 @@ namespace Step12 void AdvectionProblem::run() { - //#ifdef HEX - // for (unsigned int cycle = 0; cycle < 6; ++cycle) - //#else - for (unsigned int cycle = 0; cycle < 1; ++cycle) - //#endif + for (unsigned int cycle = 0; cycle < 4; ++cycle) { - deallog << "Cycle " << cycle << std::endl; + std::cout << "Cycle " << cycle << std::endl; if (cycle == 0) { -#ifdef HEX - // GridGenerator::hyper_cube(triangulation); - // triangulation.refine_global(3); - GridGenerator::subdivided_hyper_cube(triangulation, 16); -#else - GridGenerator::subdivided_hyper_cube_with_simplices(triangulation, - dim == 2 ? 32 : - 8); -#endif + Triangulation temp; + GridGenerator::hyper_cube(temp); + GridGenerator::convert_hypercube_to_simplex_mesh(temp, + triangulation); } else refine_grid(); - deallog << " Number of active cells: " - << triangulation.n_active_cells() << std::endl; + std::cout << " Number of active cells: " + << triangulation.n_active_cells() << std::endl; setup_system(); - deallog << " Number of degrees of freedom: " << dof_handler.n_dofs() - << std::endl; + std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; assemble_system(); solve(); @@ -546,16 +487,36 @@ namespace Step12 int main() { - initlog(); - - { - Step12::AdvectionProblem<2> dgmethod; - dgmethod.run(); - } - { - Step12::AdvectionProblem<3> dgmethod; - dgmethod.run(); - } + try + { + Step12::AdvectionProblem<2> dgmethod; + dgmethod.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-12.with_simplex_support=on.output b/tests/simplex/step-12.with_simplex_support=on.output index 7fbe5eb0f7..f6149cf944 100644 --- a/tests/simplex/step-12.with_simplex_support=on.output +++ b/tests/simplex/step-12.with_simplex_support=on.output @@ -1,15 +1,24 @@ - -DEAL::Cycle 0 -DEAL:: Number of active cells: 2048 -DEAL:: Number of degrees of freedom: 12288 -DEAL:Richardson::Starting value 0.0883883 -DEAL:Richardson::Convergence step 16 value 9.71637e-17 -DEAL:: Solver converged in 16 iterations. -DEAL:: L-infinity norm: 2.33280 -DEAL::Cycle 0 -DEAL:: Number of active cells: 2560 -DEAL:: Number of degrees of freedom: 25600 -DEAL:Richardson::Starting value 0.0360844 -DEAL:Richardson::Convergence step 13 value 2.52800e-14 -DEAL:: Solver converged in 13 iterations. -DEAL:: L-infinity norm: 13.9260 +Cycle 0 + Number of active cells: 8 + Number of degrees of freedom: 24 + Solver converged in 7 iterations. + Writing solution to + L-infinity norm: 1.00799 +Cycle 1 + Number of active cells: 32 + Number of degrees of freedom: 96 + Solver converged in 8 iterations. + Writing solution to + L-infinity norm: 1.14371 +Cycle 2 + Number of active cells: 128 + Number of degrees of freedom: 384 + Solver converged in 11 iterations. + Writing solution to + L-infinity norm: 1.09928 +Cycle 3 + Number of active cells: 512 + Number of degrees of freedom: 1536 + Solver converged in 21 iterations. + Writing solution to + L-infinity norm: 1.11521 diff --git a/tests/simplex/step-12a.cc b/tests/simplex/step-12a.cc new file mode 100644 index 0000000000..29484a3aa2 --- /dev/null +++ b/tests/simplex/step-12a.cc @@ -0,0 +1,561 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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-12 with tetrahedron mesh. + +#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" + +// Finally, the new include file for using the mesh_loop from the MeshWorker +// framework +#include + +#include + +#include + +#include +#include +#include + +#include +#include + +//#define HEX + +namespace Step12 +{ + using namespace dealii; + + template + class BoundaryValues : public Function + { + public: + BoundaryValues() = default; + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const override; + }; + + template + void + BoundaryValues::value_list(const std::vector> &points, + std::vector & values, + const unsigned int component) const + { + (void)component; + AssertIndexRange(component, 1); + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < values.size(); ++i) + { + if (points[i](0) < 0.5) + values[i] = 1.; + else + values[i] = 0.; + } + } + + + template + Tensor<1, dim> + beta(const Point &p) + { + Assert(dim >= 2, ExcNotImplemented()); + + Point wind_field; + wind_field(0) = -p(1); + wind_field(1) = p(0); + + if (wind_field.norm() > 1e-6) + wind_field /= wind_field.norm(); + + return wind_field; + } + + + template + struct ScratchData + { + ScratchData(const Mapping & mapping, + const FiniteElement & fe, + const Quadrature & quad, + const Quadrature &quad_face, + const UpdateFlags update_flags = update_values | + update_gradients | + update_quadrature_points | + update_JxW_values, + const UpdateFlags interface_update_flags = + update_values | update_gradients | update_quadrature_points | + update_JxW_values | update_normal_vectors) + : fe_values(mapping, fe, quad, update_flags) + , fe_interface_values(mapping, fe, quad_face, interface_update_flags) + {} + + + ScratchData(const ScratchData &scratch_data) + : fe_values(scratch_data.fe_values.get_mapping(), + scratch_data.fe_values.get_fe(), + scratch_data.fe_values.get_quadrature(), + scratch_data.fe_values.get_update_flags()) + , fe_interface_values( + scratch_data.fe_values + .get_mapping(), // TODO: implement for fe_interface_values + scratch_data.fe_values.get_fe(), + scratch_data.fe_interface_values.get_quadrature(), + scratch_data.fe_interface_values.get_update_flags()) + {} + + FEValues fe_values; + FEInterfaceValues fe_interface_values; + }; + + + + struct CopyDataFace + { + FullMatrix cell_matrix; + std::vector joint_dof_indices; + }; + + + + struct CopyData + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + std::vector face_data; + + 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 AdvectionProblem + { + public: + AdvectionProblem(); + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + void + output_results(const unsigned int cycle) const; + + Triangulation triangulation; +#ifdef HEX + const MappingQ1 mapping; +#else + Simplex::FE_P fe_mapping; + const MappingFE mapping; +#endif + + // Furthermore we want to use DG elements. +#ifdef HEX + FE_DGQ fe; +#else + Simplex::FE_DGP fe; +#endif + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector right_hand_side; + }; + + + template + AdvectionProblem::AdvectionProblem() +#ifdef HEX + : mapping() + , fe(2) +#else + : fe_mapping(1) + , mapping(fe_mapping) + , fe(2) +#endif + , dof_handler(triangulation) + {} + + + template + void + AdvectionProblem::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()); + right_hand_side.reinit(dof_handler.n_dofs()); + } + + template + void + AdvectionProblem::assemble_system() + { + using Iterator = typename DoFHandler::active_cell_iterator; + const BoundaryValues boundary_function; + + auto cell_worker = [&](const Iterator & cell, + ScratchData &scratch_data, + CopyData & copy_data) { + const unsigned int n_dofs = scratch_data.fe_values.get_fe().dofs_per_cell; + copy_data.reinit(cell, n_dofs); + scratch_data.fe_values.reinit(cell); + + const auto &q_points = scratch_data.fe_values.get_quadrature_points(); + + const FEValues & fe_v = scratch_data.fe_values; + const std::vector &JxW = fe_v.get_JxW_values(); + + for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point) + { + auto beta_q = beta(q_points[point]); + for (unsigned int i = 0; i < n_dofs; ++i) + for (unsigned int j = 0; j < n_dofs; ++j) + { + copy_data.cell_matrix(i, j) += + -beta_q // -\beta + * fe_v.shape_grad(i, point) // \nabla \phi_i + * fe_v.shape_value(j, point) // \phi_j + * JxW[point]; // dx + } + } + }; + + auto boundary_worker = [&](const Iterator & cell, + const unsigned int &face_no, + ScratchData & scratch_data, + CopyData & copy_data) { + scratch_data.fe_interface_values.reinit(cell, face_no); + const FEFaceValuesBase &fe_face = + scratch_data.fe_interface_values.get_fe_face_values(0); + + const auto &q_points = fe_face.get_quadrature_points(); + + const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell(); + const std::vector & JxW = fe_face.get_JxW_values(); + const std::vector> &normals = fe_face.get_normal_vectors(); + + std::vector g(q_points.size()); + boundary_function.value_list(q_points, g); + + for (unsigned int point = 0; point < q_points.size(); ++point) + { + const double beta_dot_n = beta(q_points[point]) * normals[point]; + + if (beta_dot_n > 0) + { + for (unsigned int i = 0; i < n_facet_dofs; ++i) + for (unsigned int j = 0; j < n_facet_dofs; ++j) + copy_data.cell_matrix(i, j) += + fe_face.shape_value(i, point) // \phi_i + * fe_face.shape_value(j, point) // \phi_j + * beta_dot_n // \beta . n + * JxW[point]; // dx + } + else + for (unsigned int i = 0; i < n_facet_dofs; ++i) + copy_data.cell_rhs(i) += -fe_face.shape_value(i, point) // \phi_i + * g[point] // g + * beta_dot_n // \beta . n + * JxW[point]; // dx + } + }; + + auto face_worker = [&](const Iterator & cell, + const unsigned int &f, + const unsigned int &sf, + const Iterator & ncell, + const unsigned int &nf, + const unsigned int &nsf, + ScratchData & scratch_data, + CopyData & copy_data) { + FEInterfaceValues &fe_iv = scratch_data.fe_interface_values; + fe_iv.reinit(cell, f, sf, ncell, nf, nsf); + const auto &q_points = fe_iv.get_quadrature_points(); + + copy_data.face_data.emplace_back(); + CopyDataFace ©_data_face = copy_data.face_data.back(); + + const unsigned int n_dofs = 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, n_dofs); + + const std::vector & JxW = fe_iv.get_JxW_values(); + const std::vector> &normals = fe_iv.get_normal_vectors(); + + for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) + { + const double beta_dot_n = beta(q_points[qpoint]) * normals[qpoint]; + for (unsigned int i = 0; i < n_dofs; ++i) + for (unsigned int j = 0; j < n_dofs; ++j) + copy_data_face.cell_matrix(i, j) += + fe_iv.jump(i, qpoint) // [\phi_i] + * + fe_iv.shape_value((beta_dot_n > 0), j, qpoint) // phi_j^{upwind} + * beta_dot_n // (\beta . n) + * JxW[qpoint]; // dx + } + }; + + AffineConstraints constraints; + + auto copier = [&](const CopyData &c) { + constraints.distribute_local_to_global(c.cell_matrix, + c.cell_rhs, + c.local_dof_indices, + system_matrix, + right_hand_side); + + for (auto &cdf : c.face_data) + { + constraints.distribute_local_to_global(cdf.cell_matrix, + cdf.joint_dof_indices, + system_matrix); + } + }; + + const unsigned int degree = dof_handler.get_fe().degree; + +#ifdef HEX + QGauss quad(degree + 1); + + QGauss face_quad(degree + 1); +#else + Simplex::QGauss quad(degree + 1); + + Simplex::QGauss face_quad(degree + 1); +#endif + + ScratchData scratch_data(mapping, fe, quad, face_quad); + CopyData copy_data; + + MeshWorker::mesh_loop(dof_handler.begin_active(), + dof_handler.end(), + cell_worker, + copier, + scratch_data, + copy_data, + MeshWorker::assemble_own_cells | + MeshWorker::assemble_boundary_faces | + MeshWorker::assemble_own_interior_faces_once, + boundary_worker, + face_worker); + } + + template + void + AdvectionProblem::solve() + { + SolverControl solver_control(1000, 1e-12); + SolverRichardson> solver(solver_control); + + PreconditionBlockSSOR> preconditioner; + + preconditioner.initialize(system_matrix, fe.dofs_per_cell); + + solver.solve(system_matrix, solution, right_hand_side, preconditioner); + + deallog << " Solver converged in " << solver_control.last_step() + << " iterations." << std::endl; + } + + + template + void + AdvectionProblem::refine_grid() + { + Vector gradient_indicator(triangulation.n_active_cells()); + + DerivativeApproximation::approximate_gradient(mapping, + dof_handler, + solution, + gradient_indicator); + + unsigned int cell_no = 0; + for (const auto &cell : dof_handler.active_cell_iterators()) + gradient_indicator(cell_no++) *= + std::pow(cell->diameter(), 1 + 1.0 * dim / 2); + + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + gradient_indicator, + 0.3, + 0.1); + + triangulation.execute_coarsening_and_refinement(); + } + + + template + void + AdvectionProblem::output_results(const unsigned int cycle) const + { +#if false +# ifdef HEX + const std::string filename = + dim == 2 ? ("step12-quad-" + std::to_string(cycle) + ".vtk") : + ("step12-hex-" + std::to_string(cycle) + ".vtk"); +# else + const std::string filename = + dim == 2 ? ("step12-tri-" + std::to_string(cycle) + ".vtk") : + ("step12-tet-" + std::to_string(cycle) + ".vtk"); +# endif + deallog << " Writing solution to <" << filename << ">" << std::endl; + 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, 2); + + data_out.write_vtk(output); +#endif + + { + Vector values(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + Functions::ZeroFunction(), + values, + QGauss(fe.degree + 1), + VectorTools::Linfty_norm); + const double l_infty = + VectorTools::compute_global_error(triangulation, + values, + VectorTools::Linfty_norm); + deallog << " L-infinity norm: " << l_infty << std::endl; + } + } + + + template + void + AdvectionProblem::run() + { + //#ifdef HEX + // for (unsigned int cycle = 0; cycle < 6; ++cycle) + //#else + for (unsigned int cycle = 0; cycle < 1; ++cycle) + //#endif + { + deallog << "Cycle " << cycle << std::endl; + + if (cycle == 0) + { +#ifdef HEX + // GridGenerator::hyper_cube(triangulation); + // triangulation.refine_global(3); + GridGenerator::subdivided_hyper_cube(triangulation, 16); +#else + GridGenerator::subdivided_hyper_cube_with_simplices(triangulation, + dim == 2 ? 32 : + 8); +#endif + } + else + refine_grid(); + + 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); + } + } +} // namespace Step12 + + +int +main() +{ + initlog(); + + { + Step12::AdvectionProblem<2> dgmethod; + dgmethod.run(); + } + { + Step12::AdvectionProblem<3> dgmethod; + dgmethod.run(); + } + + return 0; +} diff --git a/tests/simplex/step-12a.with_simplex_support=on.output b/tests/simplex/step-12a.with_simplex_support=on.output new file mode 100644 index 0000000000..7fbe5eb0f7 --- /dev/null +++ b/tests/simplex/step-12a.with_simplex_support=on.output @@ -0,0 +1,15 @@ + +DEAL::Cycle 0 +DEAL:: Number of active cells: 2048 +DEAL:: Number of degrees of freedom: 12288 +DEAL:Richardson::Starting value 0.0883883 +DEAL:Richardson::Convergence step 16 value 9.71637e-17 +DEAL:: Solver converged in 16 iterations. +DEAL:: L-infinity norm: 2.33280 +DEAL::Cycle 0 +DEAL:: Number of active cells: 2560 +DEAL:: Number of degrees of freedom: 25600 +DEAL:Richardson::Starting value 0.0360844 +DEAL:Richardson::Convergence step 13 value 2.52800e-14 +DEAL:: Solver converged in 13 iterations. +DEAL:: L-infinity norm: 13.9260 -- 2.39.5