From a9b2a255c9197f373af079ca4106d8b9d01a5df9 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 19 Dec 2020 15:47:41 +0100 Subject: [PATCH] Enable mixed meshes for MF and DG --- source/fe/mapping_fe.cc | 14 +- tests/simplex/matrix_free_face_integral_01.cc | 238 ++++ ...integral_01.with_simplex_support=on.output | 8 + tests/simplex/step-12b.cc | 1151 +++++++++++++++++ .../step-12b.with_simplex_support=on.output | 89 ++ tests/simplex/step-12c.cc | 1151 +++++++++++++++++ .../step-12c.with_simplex_support=on.output | 56 + 7 files changed, 2703 insertions(+), 4 deletions(-) create mode 100644 tests/simplex/matrix_free_face_integral_01.cc create mode 100644 tests/simplex/matrix_free_face_integral_01.with_simplex_support=on.output create mode 100644 tests/simplex/step-12b.cc create mode 100644 tests/simplex/step-12b.with_simplex_support=on.output create mode 100644 tests/simplex/step-12c.cc create mode 100644 tests/simplex/step-12c.with_simplex_support=on.output diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index 3a4cc76024..e0e6863333 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -1697,7 +1697,13 @@ namespace internal const typename Mapping::InternalDataBase &mapping_data, const ArrayView> & output) { - AssertDimension(input.size(), output.size()); + // In the case of wedges and pyramids, faces might have different + // numbers of quadrature points on each face with the result + // that input and output have different sizes, since input has + // the correct size but the size of output is the maximum of + // all possible sizes. + AssertIndexRange(input.size(), output.size() + 1); + Assert( (dynamic_cast< const typename dealii::MappingFE::InternalData *>( @@ -1717,7 +1723,7 @@ namespace internal typename FEValuesBase::ExcAccessToUninitializedField( "update_contravariant_transformation")); - for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int i = 0; i < input.size(); ++i) output[i] = apply_transformation(data.contravariant[i], input[i]); @@ -1738,7 +1744,7 @@ namespace internal if (rank != 1) return; - for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int i = 0; i < input.size(); ++i) { output[i] = apply_transformation(data.contravariant[i], input[i]); @@ -1756,7 +1762,7 @@ namespace internal typename FEValuesBase::ExcAccessToUninitializedField( "update_covariant_transformation")); - for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int i = 0; i < input.size(); ++i) output[i] = apply_transformation(data.covariant[i], input[i]); return; diff --git a/tests/simplex/matrix_free_face_integral_01.cc b/tests/simplex/matrix_free_face_integral_01.cc new file mode 100644 index 0000000000..418ce55b47 --- /dev/null +++ b/tests/simplex/matrix_free_face_integral_01.cc @@ -0,0 +1,238 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + +// Compute the right-hand side vector due to a Neumann boundary condition on +// different meshes (pure simplex, wedge, pyramid mesh). + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +#include "./simplex_grids.h" + +using namespace dealii; + + +template +void +test(const unsigned int v, const unsigned int degree) +{ + Triangulation tria; + + std::shared_ptr> fe; + std::shared_ptr> quad; + std::shared_ptr> fe_mapping; + hp::QCollection face_quad; + + if (v == 0) + { + GridGenerator::subdivided_hyper_cube_with_simplices(tria, + dim == 2 ? 16 : 8); + fe = std::make_shared>(degree); + quad = std::make_shared>(degree + 1); + fe_mapping = std::make_shared>(1); + face_quad = hp::QCollection{QGaussSimplex(degree + 1)}; + } + else if (v == 1) + { + GridGenerator::subdivided_hyper_cube_with_wedges(tria, dim == 2 ? 16 : 8); + fe = std::make_shared>(degree); + quad = std::make_shared>(degree + 1); + fe_mapping = std::make_shared>(1); + face_quad = hp::QCollection{QGaussSimplex(degree + 1), + QGaussSimplex(degree + 1), + QGauss(degree + 1), + QGauss(degree + 1), + QGauss(degree + 1)}; + } + else if (v == 2) + { + GridGenerator::subdivided_hyper_cube_with_pyramids(tria, + dim == 2 ? 16 : 8); + fe = std::make_shared>(degree); + quad = std::make_shared>(degree + 1); + fe_mapping = std::make_shared>(1); + face_quad = hp::QCollection{QGauss(degree + 1), + QGaussSimplex(degree + 1), + QGaussSimplex(degree + 1), + QGaussSimplex(degree + 1), + QGaussSimplex(degree + 1)}; + } + else + Assert(false, ExcNotImplemented()); + + MappingFE mapping(*fe_mapping); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(*fe); + + AffineConstraints constraints; + + using VectorType = LinearAlgebra::distributed::Vector; + using number = double; + + VectorType vec_mf, vec_mb; + + vec_mf.reinit(dof_handler.n_dofs()); + vec_mb.reinit(dof_handler.n_dofs()); + + // compute vector with FEFaceEvaluation + { + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags = update_gradients | update_values; + additional_data.mapping_update_flags_boundary_faces = + update_gradients | update_values; + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, *quad, additional_data); + + matrix_free.template loop( + [&](const auto &data, auto &dst, const auto &src, const auto cell_range) { + (void)data; + (void)dst; + (void)src; + (void)cell_range; + }, + [&](const auto &data, auto &dst, const auto &src, const auto face_range) { + (void)data; + (void)dst; + (void)src; + (void)face_range; + }, + [&](const auto &data, auto &dst, const auto &src, const auto face_range) { + FEFaceEvaluation fe_eval(data, face_range, true); + for (unsigned int face = face_range.first; face < face_range.second; + face++) + { + fe_eval.reinit(face); + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + fe_eval.submit_value(1.0, q); + + fe_eval.integrate_scatter(true, false, dst); + } + }, + vec_mf, + vec_mf); + } + + // compute vector with FEFaceValues + { + const UpdateFlags flag = update_JxW_values | update_values | + update_gradients | update_quadrature_points; + FEValues fe_values(mapping, *fe, *quad, flag); + + FEFaceValues fe_face_values(mapping, *fe, face_quad, flag); + + const unsigned int dofs_per_cell = fe->dofs_per_cell; + + std::vector dof_indices(dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + for (const auto &cell : dof_handler.cell_iterators()) + { + if (!cell->is_locally_owned()) + continue; + + fe_values.reinit(cell); + cell_rhs = 0; + + for (const auto &face : cell->face_iterators()) + if (face->at_boundary()) + { + fe_face_values.reinit(cell, face); + for (const auto q : fe_face_values.quadrature_point_indices()) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + cell_rhs(i) += + (1.0 * // 1.0 + fe_face_values.shape_value(i, q) * // phi_i(x_q) + fe_face_values.JxW(q)); // dx + } + + cell->get_dof_indices(dof_indices); + + constraints.distribute_local_to_global(cell_rhs, dof_indices, vec_mb); + } + } + +#if false + vec_mf.print(deallog.get_file_stream()); + vec_mb.print(deallog.get_file_stream()); +#endif + + for (unsigned int i = 0; i < vec_mf.size(); ++i) + Assert(std::abs(vec_mf[i] - vec_mb[i]) < 1e-8, ExcNotImplemented()); + + deallog << " dim=" << dim << " degree=" << degree << ": OK" << std::endl; +} + + +int +main(int argc, char **argv) +{ + initlog(); + + deallog.depth_file(2); + + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + + for (unsigned int i = 0; i <= 2; ++i) + { + if (i == 0) + deallog.push("SIMPLEX"); + else if (i == 1) + deallog.push("WEDGE "); + else if (i == 2) + deallog.push("PYRAMID"); + else + Assert(false, ExcNotImplemented()); + + if (i == 0) // 2D makes only sense for simplex + { + test<2>(i, /*degree=*/1); + test<2>(i, /*degree=*/2); + } + + test<3>(i, /*degree=*/1); + + if (i != 2) // for pyramids: no quadratic elements implemented + test<3>(i, /*degree=*/2); + + deallog.pop(); + } +} diff --git a/tests/simplex/matrix_free_face_integral_01.with_simplex_support=on.output b/tests/simplex/matrix_free_face_integral_01.with_simplex_support=on.output new file mode 100644 index 0000000000..e55e7637a4 --- /dev/null +++ b/tests/simplex/matrix_free_face_integral_01.with_simplex_support=on.output @@ -0,0 +1,8 @@ + +DEAL:SIMPLEX:: dim=2 degree=1: OK +DEAL:SIMPLEX:: dim=2 degree=2: OK +DEAL:SIMPLEX:: dim=3 degree=1: OK +DEAL:SIMPLEX:: dim=3 degree=2: OK +DEAL:WEDGE :: dim=3 degree=1: OK +DEAL:WEDGE :: dim=3 degree=2: OK +DEAL:PYRAMID:: dim=3 degree=1: OK diff --git a/tests/simplex/step-12b.cc b/tests/simplex/step-12b.cc new file mode 100644 index 0000000000..3673d19a13 --- /dev/null +++ b/tests/simplex/step-12b.cc @@ -0,0 +1,1151 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 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. +// +// --------------------------------------------------------------------- + + + +// a copy of tests/hp/step-12 testing different reference-cell types + + +#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" + +#include "./simplex_grids.h" + +template +struct ScratchData +{ + hp::MappingCollection mapping; + hp::FECollection fe; + hp::QCollection quadrature; + std::vector> face_quadrature; + + std::function &)> mesh_generator; +}; + + +template +class RHS : public Function +{ +public: + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const; +}; + + +template +class BoundaryValues : public Function +{ +public: + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const; +}; + + +template +class Beta +{ +public: + Beta() + {} + void + value_list(const std::vector> &points, + std::vector> & values) const; +}; + + +template +void +RHS::value_list(const std::vector> &points, + std::vector & values, + const unsigned int) const +{ + // Assert(values.size() == points.size(), + // ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < std::min(points.size(), values.size()); ++i) + values[i] = 0; +} + + +template +void +Beta::value_list(const std::vector> &points, + std::vector> & values) const +{ + // Assert(values.size() == points.size(), + // ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < std::min(points.size(), values.size()); ++i) + { + const Point &p = points[i]; + Point & beta = values[i]; + + beta(0) = -p(1); + beta(1) = p(0); + + if (beta.norm() > 1e-10) + beta /= std::sqrt(beta.square()); + } +} + + +template +void +BoundaryValues::value_list(const std::vector> &points, + std::vector & values, + const unsigned int) const +{ + // Assert(values.size() == points.size(), + // ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < std::min(values.size(), points.size()); ++i) + { + if (points[i](0) < 0.5) + values[i] = 1.; + else + values[i] = 0.; + } +} + + +template +class DGTransportEquation +{ +public: + DGTransportEquation(); + + void + assemble_cell_term(const hp::FEValues &fe_v, + FullMatrix & ui_vi_matrix, + Vector & cell_vector) const; + + void + assemble_boundary_term(const hp::FEFaceValues &fe_v, + FullMatrix & ui_vi_matrix, + Vector & cell_vector) const; + + template + void + assemble_face_term1(const X & fe_v, + const Y & fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix) const; + + template + void + assemble_face_term2(const X & fe_v, + const Y & fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix, + FullMatrix &ui_ve_matrix, + FullMatrix &ue_ve_matrix) const; + +private: + const Beta beta_function; + const RHS rhs_function; + const BoundaryValues boundary_function; +}; + + +template +DGTransportEquation::DGTransportEquation() + : beta_function() + , rhs_function() + , boundary_function() +{} + + +template +void +DGTransportEquation::assemble_cell_term(const hp::FEValues &fe_v, + FullMatrix &ui_vi_matrix, + Vector &cell_vector) const +{ + const std::vector &JxW = + fe_v.get_present_fe_values().get_JxW_values(); + + std::vector> beta( + fe_v.get_present_fe_values().n_quadrature_points); + std::vector rhs(fe_v.get_present_fe_values().n_quadrature_points); + + beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + beta); + rhs_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + rhs); + + for (unsigned int point = 0; + point < fe_v.get_present_fe_values().n_quadrature_points; + ++point) + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + { + for (unsigned int j = 0; j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_vi_matrix(i, j) -= + beta[point] * fe_v.get_present_fe_values().shape_grad(i, point) * + fe_v.get_present_fe_values().shape_value(j, point) * JxW[point]; + + cell_vector(i) += rhs[point] * + fe_v.get_present_fe_values().shape_value(i, point) * + JxW[point]; + } +} + + +template +void +DGTransportEquation::assemble_boundary_term( + const hp::FEFaceValues &fe_v, + FullMatrix & ui_vi_matrix, + Vector & cell_vector) const +{ + const std::vector &JxW = + fe_v.get_present_fe_values().get_JxW_values(); + const std::vector> &normals = + fe_v.get_present_fe_values().get_normal_vectors(); + + std::vector> beta( + fe_v.get_present_fe_values().n_quadrature_points); + std::vector g(fe_v.get_present_fe_values().n_quadrature_points); + + beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + beta); + boundary_function.value_list( + fe_v.get_present_fe_values().get_quadrature_points(), g); + + for (unsigned int point = 0; + point < fe_v.get_present_fe_values().n_quadrature_points; + ++point) + { + const double beta_n = beta[point] * normals[point]; + if (beta_n > 0) + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int j = 0; + j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_vi_matrix(i, j) += + beta_n * fe_v.get_present_fe_values().shape_value(j, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + else + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + cell_vector(i) -= beta_n * g[point] * + fe_v.get_present_fe_values().shape_value(i, point) * + JxW[point]; + } +} + + +template +template +void +DGTransportEquation::assemble_face_term1( + const X & fe_v, + const Y & fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix) const +{ + const std::vector &JxW = + fe_v.get_present_fe_values().get_JxW_values(); + const std::vector> &normals = + fe_v.get_present_fe_values().get_normal_vectors(); + + std::vector> beta( + fe_v.get_present_fe_values().n_quadrature_points); + beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + beta); + + for (unsigned int point = 0; + point < fe_v.get_present_fe_values().n_quadrature_points; + ++point) + { + const double beta_n = beta[point] * normals[point]; + if (beta_n > 0) + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int j = 0; + j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_vi_matrix(i, j) += + beta_n * fe_v.get_present_fe_values().shape_value(j, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + else + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int k = 0; + k < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++k) + ue_vi_matrix(i, k) += + beta_n * + fe_v_neighbor.get_present_fe_values().shape_value(k, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + } +} + + +template +template +void +DGTransportEquation::assemble_face_term2( + const X & fe_v, + const Y & fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix, + FullMatrix &ui_ve_matrix, + FullMatrix &ue_ve_matrix) const +{ + const std::vector &JxW = + fe_v.get_present_fe_values().get_JxW_values(); + const std::vector> &normals = + fe_v.get_present_fe_values().get_normal_vectors(); + + std::vector> beta( + fe_v.get_present_fe_values().n_quadrature_points); + + beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + beta); + + for (unsigned int point = 0; + point < fe_v.get_present_fe_values().n_quadrature_points; + ++point) + { + const double beta_n = beta[point] * normals[point]; + if (beta_n > 0) + { + for (unsigned int i = 0; + i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int j = 0; + j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_vi_matrix(i, j) += + beta_n * fe_v.get_present_fe_values().shape_value(j, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + + for (unsigned int k = 0; + k < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++k) + for (unsigned int j = 0; + j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_ve_matrix(k, j) -= + beta_n * fe_v.get_present_fe_values().shape_value(j, point) * + fe_v_neighbor.get_present_fe_values().shape_value(k, point) * + JxW[point]; + } + else + { + for (unsigned int i = 0; + i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int l = 0; + l < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++l) + ue_vi_matrix(i, l) += + beta_n * + fe_v_neighbor.get_present_fe_values().shape_value(l, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + + for (unsigned int k = 0; + k < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++k) + for (unsigned int l = 0; + l < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++l) + ue_ve_matrix(k, l) -= + beta_n * + fe_v_neighbor.get_present_fe_values().shape_value(l, point) * + fe_v_neighbor.get_present_fe_values().shape_value(k, point) * + JxW[point]; + } + } +} + + +template +class DGMethod +{ +public: + DGMethod(const ScratchData &scratch_data); + ~DGMethod(); + + void + run(); + +private: + void + setup_system(); + void + assemble_system1(); + void + assemble_system2(); + void + solve(Vector &solution); + void + refine_grid(); + void + output_results(const unsigned int cycle) const; + + Triangulation triangulation; + const hp::MappingCollection mapping; + + hp::FECollection fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + const hp::QCollection quadrature; + const std::vector> face_quadrature; + const hp::QCollection face_quadrature_dummy; + + Vector solution1; + Vector solution2; + Vector right_hand_side; + + const std::function &)> mesh_generator; + + const DGTransportEquation dg; +}; + + +template +DGMethod::DGMethod(const ScratchData &scratch_data) + : mapping(scratch_data.mapping) + , fe(scratch_data.fe) + , dof_handler(triangulation) + , quadrature(scratch_data.quadrature) + , face_quadrature(scratch_data.face_quadrature) + , mesh_generator(scratch_data.mesh_generator) + , dg() +{} + + +template +DGMethod::~DGMethod() +{ + dof_handler.clear(); +} + + +template +void +DGMethod::setup_system() +{ + dof_handler.distribute_dofs(fe); + + sparsity_pattern.reinit(dof_handler.n_dofs(), + dof_handler.n_dofs(), + (GeometryInfo::faces_per_cell * + GeometryInfo::max_children_per_face + + 1) * + fe[0].dofs_per_cell); + + DoFTools::make_flux_sparsity_pattern(dof_handler, sparsity_pattern); + + sparsity_pattern.compress(); + + system_matrix.reinit(sparsity_pattern); + + solution1.reinit(dof_handler.n_dofs()); + solution2.reinit(dof_handler.n_dofs()); + right_hand_side.reinit(dof_handler.n_dofs()); +} + + +template +void +DGMethod::assemble_system1() +{ + const unsigned int dofs_per_cell = dof_handler.get_fe(0).dofs_per_cell; + std::vector dofs(dofs_per_cell); + std::vector dofs_neighbor(dofs_per_cell); + + const UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + const UpdateFlags face_update_flags = + update_values | update_quadrature_points | update_JxW_values | + update_normal_vectors; + + const UpdateFlags neighbor_face_update_flags = update_values; + + hp::FEValues fe_v(mapping, fe, quadrature, update_flags); + + hp::FEFaceValues fe_v_face(mapping, + fe, + face_quadrature, + face_update_flags); + hp::FESubfaceValues fe_v_subface(mapping, + fe, + face_quadrature_dummy, + face_update_flags); + hp::FEFaceValues fe_v_face_neighbor(mapping, + fe, + face_quadrature, + neighbor_face_update_flags); + hp::FESubfaceValues fe_v_subface_neighbor(mapping, + fe, + face_quadrature_dummy, + neighbor_face_update_flags); + + FullMatrix ui_vi_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix ue_vi_matrix(dofs_per_cell, dofs_per_cell); + + Vector cell_vector(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + + for (; cell != endc; ++cell) + { + ui_vi_matrix = 0; + cell_vector = 0; + + fe_v.reinit(cell); + + dg.assemble_cell_term(fe_v, ui_vi_matrix, cell_vector); + + cell->get_dof_indices(dofs); + + for (const unsigned int face_no : cell->face_indices()) + { + typename DoFHandler::face_iterator face = cell->face(face_no); + + ue_vi_matrix = 0; + + if (face->at_boundary()) + { + fe_v_face.reinit(cell, face_no); + + dg.assemble_boundary_term(fe_v_face, ui_vi_matrix, cell_vector); + } + else + { + typename DoFHandler::cell_iterator neighbor = + cell->neighbor(face_no); + ; + + if (face->has_children()) + { + const unsigned int neighbor2 = + cell->neighbor_of_neighbor(face_no); + + + for (unsigned int subface_no = 0; + subface_no < face->n_children(); + ++subface_no) + { + typename DoFHandler::active_cell_iterator + neighbor_child = + cell->neighbor_child_on_subface(face_no, subface_no); + + Assert(neighbor_child->face(neighbor2) == + face->child(subface_no), + ExcInternalError()); + Assert(!neighbor_child->has_children(), + ExcInternalError()); + + ue_vi_matrix = 0; + + fe_v_subface.reinit(cell, face_no, subface_no); + fe_v_face_neighbor.reinit(neighbor_child, neighbor2); + + dg.assemble_face_term1(fe_v_subface, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_matrix); + + neighbor_child->get_dof_indices(dofs_neighbor); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int k = 0; k < dofs_per_cell; ++k) + system_matrix.add(dofs[i], + dofs_neighbor[k], + ue_vi_matrix(i, k)); + } + } + else + { + if (neighbor->level() == cell->level()) + { + const unsigned int neighbor2 = + cell->neighbor_of_neighbor(face_no); + + fe_v_face.reinit(cell, face_no); + fe_v_face_neighbor.reinit(neighbor, neighbor2); + + dg.assemble_face_term1(fe_v_face, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_matrix); + } + else + { + Assert(neighbor->level() < cell->level(), + ExcInternalError()); + + const std::pair + faceno_subfaceno = + cell->neighbor_of_coarser_neighbor(face_no); + const unsigned int neighbor_face_no = + faceno_subfaceno.first, + neighbor_subface_no = + faceno_subfaceno.second; + + Assert(neighbor->neighbor_child_on_subface( + neighbor_face_no, neighbor_subface_no) == cell, + ExcInternalError()); + + fe_v_face.reinit(cell, face_no); + fe_v_subface_neighbor.reinit(neighbor, + neighbor_face_no, + neighbor_subface_no); + + dg.assemble_face_term1(fe_v_face, + fe_v_subface_neighbor, + ui_vi_matrix, + ue_vi_matrix); + } + + neighbor->get_dof_indices(dofs_neighbor); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int k = 0; k < dofs_per_cell; ++k) + system_matrix.add(dofs[i], + dofs_neighbor[k], + ue_vi_matrix(i, k)); + } + } + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i, j)); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + right_hand_side(dofs[i]) += cell_vector(i); + } +} + + + +template +void +DGMethod::assemble_system2() +{ + const unsigned int dofs_per_cell = dof_handler.get_fe(0).dofs_per_cell; + std::vector dofs(dofs_per_cell); + std::vector dofs_neighbor(dofs_per_cell); + + const UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + const UpdateFlags face_update_flags = + update_values | update_quadrature_points | update_JxW_values | + update_normal_vectors; + + const UpdateFlags neighbor_face_update_flags = update_values; + + hp::FEValues fe_v(mapping, fe, quadrature, update_flags); + hp::FEFaceValues fe_v_face(mapping, + fe, + face_quadrature, + face_update_flags); + hp::FESubfaceValues fe_v_subface(mapping, + fe, + face_quadrature_dummy, + face_update_flags); + hp::FEFaceValues fe_v_face_neighbor(mapping, + fe, + face_quadrature, + neighbor_face_update_flags); + + + FullMatrix ui_vi_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix ue_vi_matrix(dofs_per_cell, dofs_per_cell); + + FullMatrix ui_ve_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix ue_ve_matrix(dofs_per_cell, dofs_per_cell); + + Vector cell_vector(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + ui_vi_matrix = 0; + cell_vector = 0; + + fe_v.reinit(cell); + + dg.assemble_cell_term(fe_v, ui_vi_matrix, cell_vector); + + cell->get_dof_indices(dofs); + + for (const unsigned int face_no : cell->face_indices()) + { + typename DoFHandler::face_iterator face = cell->face(face_no); + + if (face->at_boundary()) + { + fe_v_face.reinit(cell, face_no); + + dg.assemble_boundary_term(fe_v_face, ui_vi_matrix, cell_vector); + } + else + { + Assert(cell->neighbor(face_no).state() == IteratorState::valid, + ExcInternalError()); + typename DoFHandler::cell_iterator neighbor = + cell->neighbor(face_no); + if (face->has_children()) + { + const unsigned int neighbor2 = + cell->neighbor_of_neighbor(face_no); + + for (unsigned int subface_no = 0; + subface_no < face->n_children(); + ++subface_no) + { + typename DoFHandler::cell_iterator neighbor_child = + cell->neighbor_child_on_subface(face_no, subface_no); + Assert(neighbor_child->face(neighbor2) == + face->child(subface_no), + ExcInternalError()); + Assert(!neighbor_child->has_children(), + ExcInternalError()); + + ue_vi_matrix = 0; + ui_ve_matrix = 0; + ue_ve_matrix = 0; + + fe_v_subface.reinit(cell, face_no, subface_no); + fe_v_face_neighbor.reinit(neighbor_child, neighbor2); + + dg.assemble_face_term2(fe_v_subface, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_matrix, + ui_ve_matrix, + ue_ve_matrix); + + neighbor_child->get_dof_indices(dofs_neighbor); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + system_matrix.add(dofs[i], + dofs_neighbor[j], + ue_vi_matrix(i, j)); + system_matrix.add(dofs_neighbor[i], + dofs[j], + ui_ve_matrix(i, j)); + system_matrix.add(dofs_neighbor[i], + dofs_neighbor[j], + ue_ve_matrix(i, j)); + } + } + } + else + { + if (neighbor->level() == cell->level() && + neighbor->index() > cell->index()) + { + const unsigned int neighbor2 = + cell->neighbor_of_neighbor(face_no); + + ue_vi_matrix = 0; + ui_ve_matrix = 0; + ue_ve_matrix = 0; + + fe_v_face.reinit(cell, face_no); + fe_v_face_neighbor.reinit(neighbor, neighbor2); + + dg.assemble_face_term2(fe_v_face, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_matrix, + ui_ve_matrix, + ue_ve_matrix); + + neighbor->get_dof_indices(dofs_neighbor); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + system_matrix.add(dofs[i], + dofs_neighbor[j], + ue_vi_matrix(i, j)); + system_matrix.add(dofs_neighbor[i], + dofs[j], + ui_ve_matrix(i, j)); + system_matrix.add(dofs_neighbor[i], + dofs_neighbor[j], + ue_ve_matrix(i, j)); + } + } + } + } + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i, j)); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + right_hand_side(dofs[i]) += cell_vector(i); + } +} + + +template +void +DGMethod::solve(Vector &solution) +{ + ReductionControl solver_control(10000, 1e-4, 1e-3, false, false); + SolverRichardson<> solver(solver_control); + + PreconditionBlockSSOR> preconditioner; + preconditioner.initialize(system_matrix, fe[0].dofs_per_cell); + + deallog << " " << right_hand_side.l2_norm() << std::endl; + solver.solve(system_matrix, solution, right_hand_side, preconditioner); + deallog << " " << solution.l2_norm() << std::endl; +} + + +template +void +DGMethod::refine_grid() +{ + Vector gradient_indicator(triangulation.n_active_cells()); + + DerivativeApproximation::approximate_gradient(mapping[0], + dof_handler, + solution2, + gradient_indicator); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int cell_no = 0; cell != endc; ++cell, ++cell_no) + 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 +DGMethod::output_results(const unsigned int cycle) const +{ + return; + + std::string filename = "grid-"; + filename += ('0' + cycle); + Assert(cycle < 10, ExcInternalError()); + + filename += ".eps"; + deallog << "Writing grid to <" << filename << ">..." << std::endl; + + GridOut grid_out; + // grid_out.write_eps(triangulation, deallog.get_file_stream()); + + filename = "sol-"; + filename += ('0' + cycle); + Assert(cycle < 10, ExcInternalError()); + + filename += ".gnuplot"; + deallog << "Writing solution to <" << filename << ">..." << std::endl + << std::endl; + + DataOut> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution2, "u"); + + data_out.build_patches(mapping, 2); + +#if false + data_out.write_vtk(deallog.get_file_stream()); +#else + static unsigned int counter = 0; + + std::ofstream output("test." + std::to_string(dim) + "." + + std::to_string(counter++) + ".vtk"); + data_out.write_vtk(output); +#endif +} + + +template +void +DGMethod::run() +{ + for (unsigned int cycle = 0; cycle < 6; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + mesh_generator(triangulation); + + + 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; + + Timer assemble_timer; + assemble_system1(); + solve(solution1); + + system_matrix = 0; + right_hand_side = 0; + assemble_timer.reset(); + + assemble_timer.start(); + assemble_system2(); + solve(solution2); + + solution1 -= solution2; + const double difference = solution1.linfty_norm(); + if (difference > 1e-13) + deallog << "solution1 and solution2 differ!!" << std::endl; + else + deallog << "solution1 and solution2 coincide." << std::endl; + + output_results(cycle); + + break; + } +} + +int +main() +{ + try + { + initlog(); + deallog.get_file_stream().precision(2); + + // triangle + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<2> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<2>(MappingFE<2>(FE_SimplexP<2>(1))); + scratch_data.fe = hp::FECollection<2>(FE_SimplexDGP<2>(i)); + scratch_data.quadrature = hp::QCollection<2>(QGaussSimplex<2>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<1>(QGaussSimplex<1>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<2> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube_with_simplices(triangulation, + 16); + }; + + DGMethod<2> dgmethod(scratch_data); + dgmethod.run(); + } + + // quadrilateral + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<2> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<2>(MappingQGeneric<2>(1)); + scratch_data.fe = hp::FECollection<2>(FE_DGQ<2>(i)); + scratch_data.quadrature = hp::QCollection<2>(QGauss<2>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<1>(QGauss<1>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<2> &triangulation) -> void { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(3); + }; + + DGMethod<2> dgmethod(scratch_data); + dgmethod.run(); + } + + // tetrahedron + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<3> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<3>(MappingFE<3>(FE_SimplexP<3>(1))); + scratch_data.fe = hp::FECollection<3>(FE_SimplexDGP<3>(i)); + scratch_data.quadrature = hp::QCollection<3>(QGaussSimplex<3>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<2>(QGaussSimplex<2>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<3> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube_with_simplices(triangulation, + 8); + }; + + DGMethod<3> dgmethod(scratch_data); + dgmethod.run(); + } + + // pyramid + for (unsigned int i = 1; i <= 1; ++i) + { + ScratchData<3> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<3>(MappingFE<3>(FE_PyramidP<3>(1))); + scratch_data.fe = hp::FECollection<3>(FE_PyramidDGP<3>(i)); + scratch_data.quadrature = hp::QCollection<3>(QGaussPyramid<3>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<2>(QGauss<2>(i + 1), + QGaussSimplex<2>(i + 1), + QGaussSimplex<2>(i + 1), + QGaussSimplex<2>(i + 1), + QGaussSimplex<2>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<3> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube_with_pyramids(triangulation, + 8); + }; + + DGMethod<3> dgmethod(scratch_data); + dgmethod.run(); + } + + // wedge + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<3> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<3>(MappingFE<3>(FE_WedgeP<3>(1))); + scratch_data.fe = hp::FECollection<3>(FE_WedgeDGP<3>(i)); + scratch_data.quadrature = hp::QCollection<3>(QGaussWedge<3>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<2>(QGaussSimplex<2>(i + 1), + QGaussSimplex<2>(i + 1), + QGauss<2>(i + 1), + QGauss<2>(i + 1), + QGauss<2>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<3> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube_with_wedges(triangulation, 8); + }; + + DGMethod<3> dgmethod(scratch_data); + dgmethod.run(); + } + + // hexahedron + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<3> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<3>(MappingQGeneric<3>(1)); + scratch_data.fe = hp::FECollection<3>(FE_DGQ<3>(i)); + scratch_data.quadrature = hp::QCollection<3>(QGauss<3>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<2>(QGauss<2>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<3> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube(triangulation, 8); + }; + + DGMethod<3> dgmethod(scratch_data); + 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-12b.with_simplex_support=on.output b/tests/simplex/step-12b.with_simplex_support=on.output new file mode 100644 index 0000000000..89ed3cd545 --- /dev/null +++ b/tests/simplex/step-12b.with_simplex_support=on.output @@ -0,0 +1,89 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 512 +DEAL:: Number of degrees of freedom: 1536 +DEAL:: 0.125000 +DEAL:: 17.2926 +DEAL:: 0.125000 +DEAL:: 17.2926 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 512 +DEAL:: Number of degrees of freedom: 3072 +DEAL:: 0.125000 +DEAL:: 24.3239 +DEAL:: 0.125000 +DEAL:: 24.3239 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 64 +DEAL:: Number of degrees of freedom: 256 +DEAL:: 0.176777 +DEAL:: 7.04972 +DEAL:: 0.176777 +DEAL:: 7.04972 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 64 +DEAL:: Number of degrees of freedom: 576 +DEAL:: 0.176777 +DEAL:: 10.4560 +DEAL:: 0.176777 +DEAL:: 10.4560 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 2560 +DEAL:: Number of degrees of freedom: 10240 +DEAL:: 0.0360844 +DEAL:: 45.1474 +DEAL:: 0.0360844 +DEAL:: 45.1474 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 2560 +DEAL:: Number of degrees of freedom: 25600 +DEAL:: 0.0360844 +DEAL:: 69.8629 +DEAL:: 0.0360844 +DEAL:: 69.8629 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 3072 +DEAL:: Number of degrees of freedom: 15360 +DEAL:: 0.0441942 +DEAL:: 54.8297 +DEAL:: 0.0441942 +DEAL:: 54.8297 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 1024 +DEAL:: Number of degrees of freedom: 6144 +DEAL:: 0.0441942 +DEAL:: 34.6857 +DEAL:: 0.0441942 +DEAL:: 34.6857 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 1024 +DEAL:: Number of degrees of freedom: 18432 +DEAL:: 0.0441942 +DEAL:: 59.5766 +DEAL:: 0.0441942 +DEAL:: 59.5766 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 512 +DEAL:: Number of degrees of freedom: 4096 +DEAL:: 0.0441942 +DEAL:: 28.1989 +DEAL:: 0.0441942 +DEAL:: 28.1989 +DEAL::solution1 and solution2 coincide. +DEAL::Cycle 0: +DEAL:: Number of active cells: 512 +DEAL:: Number of degrees of freedom: 13824 +DEAL:: 0.0441942 +DEAL:: 51.2235 +DEAL:: 0.0441942 +DEAL:: 51.2235 +DEAL::solution1 and solution2 coincide. diff --git a/tests/simplex/step-12c.cc b/tests/simplex/step-12c.cc new file mode 100644 index 0000000000..086c091c54 --- /dev/null +++ b/tests/simplex/step-12c.cc @@ -0,0 +1,1151 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 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. +// +// --------------------------------------------------------------------- + + + +// a copy of tests/hp/step-12 testing different reference-cell types + + +#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" + +#include "./simplex_grids.h" + +template +struct ScratchData +{ + hp::MappingCollection mapping; + hp::FECollection fe; + hp::QCollection quadrature; + std::vector> face_quadrature; + + std::function &)> mesh_generator; +}; + + +template +class RHS : public Function +{ +public: + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const; +}; + + +template +class BoundaryValues : public Function +{ +public: + virtual void + value_list(const std::vector> &points, + std::vector & values, + const unsigned int component = 0) const; +}; + + +template +class Beta +{ +public: + Beta() + {} + void + value_list(const std::vector> &points, + std::vector> & values) const; +}; + + +template +void +RHS::value_list(const std::vector> &points, + std::vector & values, + const unsigned int) const +{ + // Assert(values.size() == points.size(), + // ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < std::min(points.size(), values.size()); ++i) + values[i] = 0; +} + + +template +void +Beta::value_list(const std::vector> &points, + std::vector> & values) const +{ + // Assert(values.size() == points.size(), + // ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < std::min(points.size(), values.size()); ++i) + { + const Point &p = points[i]; + Point & beta = values[i]; + + beta(0) = -p(1); + beta(1) = p(0); + + if (beta.norm() > 1e-10) + beta /= std::sqrt(beta.square()); + } +} + + +template +void +BoundaryValues::value_list(const std::vector> &points, + std::vector & values, + const unsigned int) const +{ + // Assert(values.size() == points.size(), + // ExcDimensionMismatch(values.size(), points.size())); + + for (unsigned int i = 0; i < std::min(values.size(), points.size()); ++i) + { + if (points[i](0) < 0.5) + values[i] = 1.; + else + values[i] = 0.; + } +} + + +template +class DGTransportEquation +{ +public: + DGTransportEquation(); + + void + assemble_cell_term(const hp::FEValues &fe_v, + FullMatrix & ui_vi_matrix, + Vector & cell_vector) const; + + void + assemble_boundary_term(const hp::FEFaceValues &fe_v, + FullMatrix & ui_vi_matrix, + Vector & cell_vector) const; + + template + void + assemble_face_term1(const X & fe_v, + const Y & fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix) const; + + template + void + assemble_face_term2(const X & fe_v, + const Y & fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix, + FullMatrix &ui_ve_matrix, + FullMatrix &ue_ve_matrix) const; + +private: + const Beta beta_function; + const RHS rhs_function; + const BoundaryValues boundary_function; +}; + + +template +DGTransportEquation::DGTransportEquation() + : beta_function() + , rhs_function() + , boundary_function() +{} + + +template +void +DGTransportEquation::assemble_cell_term(const hp::FEValues &fe_v, + FullMatrix &ui_vi_matrix, + Vector &cell_vector) const +{ + const std::vector &JxW = + fe_v.get_present_fe_values().get_JxW_values(); + + std::vector> beta( + fe_v.get_present_fe_values().n_quadrature_points); + std::vector rhs(fe_v.get_present_fe_values().n_quadrature_points); + + beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + beta); + rhs_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + rhs); + + for (unsigned int point = 0; + point < fe_v.get_present_fe_values().n_quadrature_points; + ++point) + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + { + for (unsigned int j = 0; j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_vi_matrix(i, j) -= + beta[point] * fe_v.get_present_fe_values().shape_grad(i, point) * + fe_v.get_present_fe_values().shape_value(j, point) * JxW[point]; + + cell_vector(i) += rhs[point] * + fe_v.get_present_fe_values().shape_value(i, point) * + JxW[point]; + } +} + + +template +void +DGTransportEquation::assemble_boundary_term( + const hp::FEFaceValues &fe_v, + FullMatrix & ui_vi_matrix, + Vector & cell_vector) const +{ + const std::vector &JxW = + fe_v.get_present_fe_values().get_JxW_values(); + const std::vector> &normals = + fe_v.get_present_fe_values().get_normal_vectors(); + + std::vector> beta( + fe_v.get_present_fe_values().n_quadrature_points); + std::vector g(fe_v.get_present_fe_values().n_quadrature_points); + + beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + beta); + boundary_function.value_list( + fe_v.get_present_fe_values().get_quadrature_points(), g); + + for (unsigned int point = 0; + point < fe_v.get_present_fe_values().n_quadrature_points; + ++point) + { + const double beta_n = beta[point] * normals[point]; + if (beta_n > 0) + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int j = 0; + j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_vi_matrix(i, j) += + beta_n * fe_v.get_present_fe_values().shape_value(j, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + else + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + cell_vector(i) -= beta_n * g[point] * + fe_v.get_present_fe_values().shape_value(i, point) * + JxW[point]; + } +} + + +template +template +void +DGTransportEquation::assemble_face_term1( + const X & fe_v, + const Y & fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix) const +{ + const std::vector &JxW = + fe_v.get_present_fe_values().get_JxW_values(); + const std::vector> &normals = + fe_v.get_present_fe_values().get_normal_vectors(); + + std::vector> beta( + fe_v.get_present_fe_values().n_quadrature_points); + beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + beta); + + for (unsigned int point = 0; + point < fe_v.get_present_fe_values().n_quadrature_points; + ++point) + { + const double beta_n = beta[point] * normals[point]; + if (beta_n > 0) + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int j = 0; + j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_vi_matrix(i, j) += + beta_n * fe_v.get_present_fe_values().shape_value(j, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + else + for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int k = 0; + k < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++k) + ue_vi_matrix(i, k) += + beta_n * + fe_v_neighbor.get_present_fe_values().shape_value(k, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + } +} + + +template +template +void +DGTransportEquation::assemble_face_term2( + const X & fe_v, + const Y & fe_v_neighbor, + FullMatrix &ui_vi_matrix, + FullMatrix &ue_vi_matrix, + FullMatrix &ui_ve_matrix, + FullMatrix &ue_ve_matrix) const +{ + const std::vector &JxW = + fe_v.get_present_fe_values().get_JxW_values(); + const std::vector> &normals = + fe_v.get_present_fe_values().get_normal_vectors(); + + std::vector> beta( + fe_v.get_present_fe_values().n_quadrature_points); + + beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(), + beta); + + for (unsigned int point = 0; + point < fe_v.get_present_fe_values().n_quadrature_points; + ++point) + { + const double beta_n = beta[point] * normals[point]; + if (beta_n > 0) + { + for (unsigned int i = 0; + i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int j = 0; + j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_vi_matrix(i, j) += + beta_n * fe_v.get_present_fe_values().shape_value(j, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + + for (unsigned int k = 0; + k < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++k) + for (unsigned int j = 0; + j < fe_v.get_present_fe_values().dofs_per_cell; + ++j) + ui_ve_matrix(k, j) -= + beta_n * fe_v.get_present_fe_values().shape_value(j, point) * + fe_v_neighbor.get_present_fe_values().shape_value(k, point) * + JxW[point]; + } + else + { + for (unsigned int i = 0; + i < fe_v.get_present_fe_values().dofs_per_cell; + ++i) + for (unsigned int l = 0; + l < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++l) + ue_vi_matrix(i, l) += + beta_n * + fe_v_neighbor.get_present_fe_values().shape_value(l, point) * + fe_v.get_present_fe_values().shape_value(i, point) * JxW[point]; + + for (unsigned int k = 0; + k < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++k) + for (unsigned int l = 0; + l < fe_v_neighbor.get_present_fe_values().dofs_per_cell; + ++l) + ue_ve_matrix(k, l) -= + beta_n * + fe_v_neighbor.get_present_fe_values().shape_value(l, point) * + fe_v_neighbor.get_present_fe_values().shape_value(k, point) * + JxW[point]; + } + } +} + + +template +class DGMethod +{ +public: + DGMethod(const ScratchData &scratch_data); + ~DGMethod(); + + void + run(); + +private: + void + setup_system(); + void + assemble_system2(); + void + solve(Vector &solution); + void + refine_grid(); + void + output_results(const unsigned int cycle) const; + + Triangulation triangulation; + const hp::MappingCollection mapping; + + hp::FECollection fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + const hp::QCollection quadrature; + const std::vector> face_quadrature; + const hp::QCollection face_quadrature_dummy; + + Vector solution1; + Vector solution2; + Vector right_hand_side; + + const std::function &)> mesh_generator; + + const DGTransportEquation dg; +}; + + +template +DGMethod::DGMethod(const ScratchData &scratch_data) + : mapping(scratch_data.mapping) + , fe(scratch_data.fe) + , dof_handler(triangulation) + , quadrature(scratch_data.quadrature) + , face_quadrature(scratch_data.face_quadrature) + , mesh_generator(scratch_data.mesh_generator) + , dg() +{} + + +template +DGMethod::~DGMethod() +{ + dof_handler.clear(); +} + + +template +void +DGMethod::setup_system() +{ + dof_handler.distribute_dofs(fe); + + sparsity_pattern.reinit(dof_handler.n_dofs(), + dof_handler.n_dofs(), + (GeometryInfo::faces_per_cell * + GeometryInfo::max_children_per_face + + 1) * + fe[0].dofs_per_cell); + + DoFTools::make_flux_sparsity_pattern(dof_handler, sparsity_pattern); + + sparsity_pattern.compress(); + + system_matrix.reinit(sparsity_pattern); + + solution1.reinit(dof_handler.n_dofs()); + solution2.reinit(dof_handler.n_dofs()); + right_hand_side.reinit(dof_handler.n_dofs()); +} + + + +template +void +DGMethod::assemble_system2() +{ + const unsigned int dofs_per_cell = dof_handler.get_fe(0).dofs_per_cell; + std::vector dofs(dofs_per_cell); + std::vector dofs_neighbor(dofs_per_cell); + + const UpdateFlags update_flags = update_values | update_gradients | + update_quadrature_points | update_JxW_values; + + const UpdateFlags face_update_flags = + update_values | update_quadrature_points | update_JxW_values | + update_normal_vectors; + + const UpdateFlags neighbor_face_update_flags = update_values; + + hp::FEValues fe_v(mapping, fe, quadrature, update_flags); + hp::FEFaceValues fe_v_face(mapping, + fe, + face_quadrature, + face_update_flags); + hp::FESubfaceValues fe_v_subface(mapping, + fe, + face_quadrature_dummy, + face_update_flags); + hp::FEFaceValues fe_v_face_neighbor(mapping, + fe, + face_quadrature, + neighbor_face_update_flags); + + + FullMatrix ui_vi_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix ue_vi_matrix(dofs_per_cell, dofs_per_cell); + + FullMatrix ui_ve_matrix(dofs_per_cell, dofs_per_cell); + FullMatrix ue_ve_matrix(dofs_per_cell, dofs_per_cell); + + Vector cell_vector(dofs_per_cell); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + ui_vi_matrix = 0; + cell_vector = 0; + + fe_v.reinit(cell); + + dg.assemble_cell_term(fe_v, ui_vi_matrix, cell_vector); + + cell->get_dof_indices(dofs); + + for (const unsigned int face_no : cell->face_indices()) + { + typename DoFHandler::face_iterator face = cell->face(face_no); + + if (face->at_boundary()) + { + fe_v_face.reinit(cell, face_no); + + dg.assemble_boundary_term(fe_v_face, ui_vi_matrix, cell_vector); + } + else + { + Assert(cell->neighbor(face_no).state() == IteratorState::valid, + ExcInternalError()); + typename DoFHandler::cell_iterator neighbor = + cell->neighbor(face_no); + if (face->has_children()) + { + const unsigned int neighbor2 = + cell->neighbor_of_neighbor(face_no); + + for (unsigned int subface_no = 0; + subface_no < face->n_children(); + ++subface_no) + { + typename DoFHandler::cell_iterator neighbor_child = + cell->neighbor_child_on_subface(face_no, subface_no); + Assert(neighbor_child->face(neighbor2) == + face->child(subface_no), + ExcInternalError()); + Assert(!neighbor_child->has_children(), + ExcInternalError()); + + ue_vi_matrix = 0; + ui_ve_matrix = 0; + ue_ve_matrix = 0; + + fe_v_subface.reinit(cell, face_no, subface_no); + fe_v_face_neighbor.reinit(neighbor_child, neighbor2); + + dg.assemble_face_term2(fe_v_subface, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_matrix, + ui_ve_matrix, + ue_ve_matrix); + + neighbor_child->get_dof_indices(dofs_neighbor); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + system_matrix.add(dofs[i], + dofs_neighbor[j], + ue_vi_matrix(i, j)); + system_matrix.add(dofs_neighbor[i], + dofs[j], + ui_ve_matrix(i, j)); + system_matrix.add(dofs_neighbor[i], + dofs_neighbor[j], + ue_ve_matrix(i, j)); + } + } + } + else + { + if (neighbor->level() == cell->level() && + neighbor->index() > cell->index()) + { + const unsigned int neighbor2 = + cell->neighbor_of_neighbor(face_no); + + ue_vi_matrix = 0; + ui_ve_matrix = 0; + ue_ve_matrix = 0; + + fe_v_face.reinit(cell, face_no); + fe_v_face_neighbor.reinit(neighbor, neighbor2); + + dg.assemble_face_term2(fe_v_face, + fe_v_face_neighbor, + ui_vi_matrix, + ue_vi_matrix, + ui_ve_matrix, + ue_ve_matrix); + + neighbor->get_dof_indices(dofs_neighbor); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + system_matrix.add(dofs[i], + dofs_neighbor[j], + ue_vi_matrix(i, j)); + system_matrix.add(dofs_neighbor[i], + dofs[j], + ui_ve_matrix(i, j)); + system_matrix.add(dofs_neighbor[i], + dofs_neighbor[j], + ue_ve_matrix(i, j)); + } + } + } + } + } + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i, j)); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + right_hand_side(dofs[i]) += cell_vector(i); + } +} + + + +template +class Operator +{ +public: + using VectorType = Vector; + using number = double; + + Operator(const MatrixFree &matrix_free) + : matrix_free(matrix_free) + {} + + void + initialize_dof_vector(VectorType &vec) + { + matrix_free.initialize_dof_vector(vec); + } + + void + rhs(VectorType &vec) const + { + const int dummy = 0; + + matrix_free.template loop( + [&](const auto &data, auto &dst, const auto &src, const auto cell_range) { + FEEvaluation phi(data, cell_range); + for (unsigned int cell = cell_range.first; cell < cell_range.second; + ++cell) + { + phi.reinit(cell); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_value(0.0, q); + + phi.integrate_scatter(true, false, dst); + } + }, + [&](const auto &data, auto &dst, const auto &src, const auto face_range) { + (void)data; + (void)dst; + (void)src; + (void)face_range; + }, + [&](const auto &data, auto &dst, const auto &src, const auto face_range) { + FEFaceEvaluation phi(data, face_range); + for (unsigned int face = face_range.first; face < face_range.second; + ++face) + { + phi.reinit(face); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + const auto beta_n = + beta(phi.quadrature_point(q)) * phi.get_normal_vector(q); + const auto beta_n_m = (-std::abs(beta_n) + beta_n) * 0.5; + phi.submit_value(-beta_n_m * + boundary_values(phi.quadrature_point(q)), + q); + } + + phi.integrate_scatter(true, false, dst); + } + }, + vec, + dummy, + true); + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + matrix_free.template loop( + [&](const auto &data, auto &dst, const auto &src, const auto cell_range) { + FEEvaluation phi(data, cell_range); + for (unsigned int cell = cell_range.first; cell < cell_range.second; + ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(src, true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(-beta(phi.quadrature_point(q)) * + phi.get_value(q), + q); + phi.integrate_scatter(false, true, dst); + } + }, + [&](const auto &data, auto &dst, const auto &src, const auto face_range) { + FEFaceEvaluation phi_m(data, face_range, true); + FEFaceEvaluation phi_p(data, face_range, false); + for (unsigned int cell = face_range.first; cell < face_range.second; + ++cell) + { + phi_m.reinit(cell); + phi_m.gather_evaluate(src, true, false); + phi_p.reinit(cell); + phi_p.gather_evaluate(src, true, false); + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + const auto beta_n = + beta(phi_m.quadrature_point(q)) * phi_m.get_normal_vector(q); + + const auto u_m = phi_m.get_value(q); + const auto u_p = phi_p.get_value(q); + const auto temp = + 0.5 * ((u_m + u_p) * beta_n + (u_m - u_p) * std::abs(beta_n)); + + phi_m.submit_value(+temp, q); + phi_p.submit_value(-temp, q); + } + + phi_m.integrate_scatter(true, false, dst); + phi_p.integrate_scatter(true, false, dst); + } + }, + [&](const auto &data, auto &dst, const auto &src, const auto face_range) { + FEFaceEvaluation phi(data, face_range); + for (unsigned int cell = face_range.first; cell < face_range.second; + ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(src, true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + const auto beta_n = + beta(phi.quadrature_point(q)) * phi.get_normal_vector(q); + const auto beta_n_p = (std::abs(beta_n) + beta_n) * 0.5; + phi.submit_value(beta_n_p * phi.get_value(q), q); + } + + phi.integrate_scatter(true, false, dst); + } + }, + dst, + src, + true); + } + + Tensor<1, dim, VectorizedArray> + beta(const Point> &points) const + { + Tensor<1, dim, VectorizedArray> betas; + + for (unsigned int v = 0; v < VectorizedArray::size(); ++v) + { + Point p; + Tensor<1, dim, double> beta; + + for (int d = 0; d < dim; ++d) + p[d] = points[d][v]; + + beta[0] = -p[1]; + beta[1] = +p[0]; + + if (beta.norm() > 1e-10) + beta /= beta.norm(); + + for (int d = 0; d < dim; ++d) + betas[d][v] = beta[d]; + } + + return betas; + } + + VectorizedArray + boundary_values(const Point> &points) const + { + VectorizedArray betas; + + for (unsigned int v = 0; v < VectorizedArray::size(); ++v) + betas[v] = (points[0][v] < 0.5) ? 1.0 : 0.0; + + return betas; + } + +private: + const MatrixFree &matrix_free; +}; + + +template +void +DGMethod::solve(Vector &solution) +{ + ReductionControl solver_control(10000, 1e-4, 1e-3, false, false); + SolverRichardson<> solver(solver_control); + + PreconditionBlockSSOR> preconditioner; + preconditioner.initialize(system_matrix, fe[0].dofs_per_cell); + + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags = + update_gradients | update_values | update_quadrature_points; + additional_data.mapping_update_flags_inner_faces = + update_gradients | update_values | update_quadrature_points; + additional_data.mapping_update_flags_boundary_faces = + update_gradients | update_values | update_quadrature_points; + additional_data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + + AffineConstraints constraints; + constraints.close(); + + MatrixFree matrix_free; + matrix_free.reinit( + mapping, dof_handler, constraints, quadrature, additional_data); + + Operator op(matrix_free); + + op.initialize_dof_vector(solution); + op.initialize_dof_vector(right_hand_side); + + op.rhs(right_hand_side); + + deallog << " " << right_hand_side.l2_norm() << std::endl; + solver.solve(op, solution, right_hand_side, preconditioner); + deallog << " " << solution.l2_norm() << std::endl; +} + + +template +void +DGMethod::refine_grid() +{ + Vector gradient_indicator(triangulation.n_active_cells()); + + DerivativeApproximation::approximate_gradient(mapping[0], + dof_handler, + solution2, + gradient_indicator); + + typename DoFHandler::active_cell_iterator cell = + dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int cell_no = 0; cell != endc; ++cell, ++cell_no) + 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 +DGMethod::output_results(const unsigned int cycle) const +{ + return; + + std::string filename = "grid-"; + filename += ('0' + cycle); + Assert(cycle < 10, ExcInternalError()); + + filename += ".eps"; + deallog << "Writing grid to <" << filename << ">..." << std::endl; + + GridOut grid_out; + // grid_out.write_eps(triangulation, deallog.get_file_stream()); + + filename = "sol-"; + filename += ('0' + cycle); + Assert(cycle < 10, ExcInternalError()); + + filename += ".gnuplot"; + deallog << "Writing solution to <" << filename << ">..." << std::endl + << std::endl; + + DataOut> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution1, "u"); + + data_out.build_patches(mapping, 2); + +#if false + data_out.write_vtk(deallog.get_file_stream()); +#else + static unsigned int counter = 0; + + std::ofstream output("test." + std::to_string(dim) + "." + + std::to_string(counter++) + ".vtk"); + data_out.write_vtk(output); +#endif +} + + +template +void +DGMethod::run() +{ + deallog << "Cycle " << 0 << ':' << std::endl; + mesh_generator(triangulation); + deallog << " Number of active cells: " + << triangulation.n_active_cells() << std::endl; + + setup_system(); + assemble_system2(); + + deallog << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + solve(solution1); + + output_results(0); +} + +int +main() +{ + try + { + initlog(); + deallog.get_file_stream().precision(2); + + // triangle + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<2> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<2>(MappingFE<2>(FE_SimplexP<2>(1))); + scratch_data.fe = hp::FECollection<2>(FE_SimplexDGP<2>(i)); + scratch_data.quadrature = hp::QCollection<2>(QGaussSimplex<2>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<1>(QGaussSimplex<1>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<2> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube_with_simplices(triangulation, + 16); + }; + + DGMethod<2> dgmethod(scratch_data); + dgmethod.run(); + } + + // quadrilateral + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<2> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<2>(MappingQGeneric<2>(1)); + scratch_data.fe = hp::FECollection<2>(FE_DGQ<2>(i)); + scratch_data.quadrature = hp::QCollection<2>(QGauss<2>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<1>(QGauss<1>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<2> &triangulation) -> void { + GridGenerator::hyper_cube(triangulation); + triangulation.refine_global(3); + }; + + DGMethod<2> dgmethod(scratch_data); + dgmethod.run(); + } + + // tetrahedron + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<3> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<3>(MappingFE<3>(FE_SimplexP<3>(1))); + scratch_data.fe = hp::FECollection<3>(FE_SimplexDGP<3>(i)); + scratch_data.quadrature = hp::QCollection<3>(QGaussSimplex<3>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<2>(QGaussSimplex<2>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<3> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube_with_simplices(triangulation, + 8); + }; + + DGMethod<3> dgmethod(scratch_data); + dgmethod.run(); + } + + // pyramid + for (unsigned int i = 1; i <= 1; ++i) + { + ScratchData<3> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<3>(MappingFE<3>(FE_PyramidP<3>(1))); + scratch_data.fe = hp::FECollection<3>(FE_PyramidDGP<3>(i)); + scratch_data.quadrature = hp::QCollection<3>(QGaussPyramid<3>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<2>(QGauss<2>(i + 1), + QGaussSimplex<2>(i + 1), + QGaussSimplex<2>(i + 1), + QGaussSimplex<2>(i + 1), + QGaussSimplex<2>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<3> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube_with_pyramids(triangulation, + 8); + }; + + DGMethod<3> dgmethod(scratch_data); + dgmethod.run(); + } + + // wedge + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<3> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<3>(MappingFE<3>(FE_WedgeP<3>(1))); + scratch_data.fe = hp::FECollection<3>(FE_WedgeDGP<3>(i)); + scratch_data.quadrature = hp::QCollection<3>(QGaussWedge<3>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<2>(QGaussSimplex<2>(i + 1), + QGaussSimplex<2>(i + 1), + QGauss<2>(i + 1), + QGauss<2>(i + 1), + QGauss<2>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<3> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube_with_wedges(triangulation, 8); + }; + + DGMethod<3> dgmethod(scratch_data); + dgmethod.run(); + } + + // hexahedron + for (unsigned int i = 1; i <= 2; ++i) + { + ScratchData<3> scratch_data; + + scratch_data.mapping = + hp::MappingCollection<3>(MappingQGeneric<3>(1)); + scratch_data.fe = hp::FECollection<3>(FE_DGQ<3>(i)); + scratch_data.quadrature = hp::QCollection<3>(QGauss<3>(i + 1)); + scratch_data.face_quadrature = std::vector>{ + hp::QCollection<2>(QGauss<2>(i + 1))}; + scratch_data.mesh_generator = + [](Triangulation<3> &triangulation) -> void { + GridGenerator::subdivided_hyper_cube(triangulation, 8); + }; + + DGMethod<3> dgmethod(scratch_data); + 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-12c.with_simplex_support=on.output b/tests/simplex/step-12c.with_simplex_support=on.output new file mode 100644 index 0000000000..b8e96aab02 --- /dev/null +++ b/tests/simplex/step-12c.with_simplex_support=on.output @@ -0,0 +1,56 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 512 +DEAL:: Number of degrees of freedom: 1536 +DEAL:: 0.125000 +DEAL:: 17.2926 +DEAL::Cycle 0: +DEAL:: Number of active cells: 512 +DEAL:: Number of degrees of freedom: 3072 +DEAL:: 0.125000 +DEAL:: 24.3239 +DEAL::Cycle 0: +DEAL:: Number of active cells: 64 +DEAL:: Number of degrees of freedom: 256 +DEAL:: 0.176777 +DEAL:: 7.04972 +DEAL::Cycle 0: +DEAL:: Number of active cells: 64 +DEAL:: Number of degrees of freedom: 576 +DEAL:: 0.176777 +DEAL:: 10.4560 +DEAL::Cycle 0: +DEAL:: Number of active cells: 2560 +DEAL:: Number of degrees of freedom: 10240 +DEAL:: 0.0360844 +DEAL:: 45.1474 +DEAL::Cycle 0: +DEAL:: Number of active cells: 2560 +DEAL:: Number of degrees of freedom: 25600 +DEAL:: 0.0360844 +DEAL:: 69.8629 +DEAL::Cycle 0: +DEAL:: Number of active cells: 3072 +DEAL:: Number of degrees of freedom: 15360 +DEAL:: 0.0441942 +DEAL:: 54.8297 +DEAL::Cycle 0: +DEAL:: Number of active cells: 1024 +DEAL:: Number of degrees of freedom: 6144 +DEAL:: 0.0441942 +DEAL:: 34.6857 +DEAL::Cycle 0: +DEAL:: Number of active cells: 1024 +DEAL:: Number of degrees of freedom: 18432 +DEAL:: 0.0441942 +DEAL:: 59.5766 +DEAL::Cycle 0: +DEAL:: Number of active cells: 512 +DEAL:: Number of degrees of freedom: 4096 +DEAL:: 0.0441942 +DEAL:: 28.1989 +DEAL::Cycle 0: +DEAL:: Number of active cells: 512 +DEAL:: Number of degrees of freedom: 13824 +DEAL:: 0.0441942 +DEAL:: 51.2235 -- 2.39.5