From a08ce69b7e58da08003195c1dfa479db737beb57 Mon Sep 17 00:00:00 2001 From: Wenyu Lei Date: Sat, 16 Jan 2021 17:13:21 +0100 Subject: [PATCH] added a simplex test for step-38 indent and fix --- tests/simplex/step-38.cc | 445 ++++++++++++++++++ .../step-38.with_simplex_support=on.output | 3 + 2 files changed, 448 insertions(+) create mode 100644 tests/simplex/step-38.cc create mode 100644 tests/simplex/step-38.with_simplex_support=on.output diff --git a/tests/simplex/step-38.cc b/tests/simplex/step-38.cc new file mode 100644 index 0000000000..82f919a5c3 --- /dev/null +++ b/tests/simplex/step-38.cc @@ -0,0 +1,445 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 2010 - 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. + * + * --------------------------------------------------------------------- + */ + +// Test step-38 on a simplex mesh. +// Note that currently the quadratic mapping for simplices is not implemented. + +#include +#include + +#include +#include +#include + +#include +#include +#include // new header +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +// simplex +#include +#include +#include + +#include +#include + +#define USE_SIMPLEX + +namespace Step38 +{ + using namespace dealii; + + template + class LaplaceBeltramiProblem + { + public: + LaplaceBeltramiProblem(const unsigned degree = 1); + void + run(); + + private: + static constexpr unsigned int dim = spacedim - 1; + + void + make_grid_and_dofs(); + void + assemble_system(); + void + solve(); + void + output_results() const; + void + compute_error() const; + + Triangulation triangulation; + DoFHandler dof_handler; + +#ifdef USE_SIMPLEX + Simplex::FE_P fe; + MappingFE mapping; +#else + FE_Q fe; + MappingQ mapping; +#endif + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + }; + + template + class Solution : public Function + { + public: + virtual double + value(const Point &p, const unsigned int component = 0) const override; + + virtual Tensor<1, dim> + gradient(const Point & p, + const unsigned int component = 0) const override; + }; + + template <> + double + Solution<2>::value(const Point<2> &p, const unsigned int) const + { + return (-2. * p(0) * p(1)); + } + + template <> + Tensor<1, 2> + Solution<2>::gradient(const Point<2> &p, const unsigned int) const + { + Tensor<1, 2> return_value; + return_value[0] = -2. * p(1) * (1 - 2. * p(0) * p(0)); + return_value[1] = -2. * p(0) * (1 - 2. * p(1) * p(1)); + + return return_value; + } + + template <> + double + Solution<3>::value(const Point<3> &p, const unsigned int) const + { + return (std::sin(numbers::PI * p(0)) * std::cos(numbers::PI * p(1)) * + exp(p(2))); + } + + template <> + Tensor<1, 3> + Solution<3>::gradient(const Point<3> &p, const unsigned int) const + { + using numbers::PI; + + Tensor<1, 3> return_value; + + return_value[0] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + return_value[1] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); + return_value[2] = sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + + return return_value; + } + + template + class RightHandSide : public Function + { + public: + virtual double + value(const Point &p, const unsigned int component = 0) const override; + }; + + template <> + double + RightHandSide<2>::value(const Point<2> &p, + const unsigned int /*component*/) const + { + return (-8. * p(0) * p(1)); + } + + template <> + double + RightHandSide<3>::value(const Point<3> &p, + const unsigned int /*component*/) const + { + using numbers::PI; + + Tensor<2, 3> hessian; + + hessian[0][0] = -PI * PI * sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + hessian[1][1] = -PI * PI * sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + hessian[2][2] = sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + + hessian[0][1] = -PI * PI * cos(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); + hessian[1][0] = -PI * PI * cos(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); + + hessian[0][2] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + hessian[2][0] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + + hessian[1][2] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); + hessian[2][1] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); + + Tensor<1, 3> gradient; + gradient[0] = PI * cos(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + gradient[1] = -PI * sin(PI * p(0)) * sin(PI * p(1)) * exp(p(2)); + gradient[2] = sin(PI * p(0)) * cos(PI * p(1)) * exp(p(2)); + + Point<3> normal = p; + normal /= p.norm(); + + return (-trace(hessian) + 2 * (gradient * normal) + + (hessian * normal) * normal); + } + + template + LaplaceBeltramiProblem::LaplaceBeltramiProblem( + const unsigned degree) + : dof_handler(triangulation) + , fe(degree) +#ifdef USE_SIMPLEX + // Quadratic mapping for simplex is still not implemented. + , mapping(Simplex::FE_P(degree)) +#else + , mapping(degree) +#endif + {} + + template + void + LaplaceBeltramiProblem::make_grid_and_dofs() + { + { + Triangulation volume_mesh; + GridGenerator::half_hyper_ball(volume_mesh); + + std::set boundary_ids; + boundary_ids.insert(0); + +#ifdef USE_SIMPLEX + Triangulation tria_temp; + GridGenerator::extract_boundary_mesh(volume_mesh, + tria_temp, + boundary_ids); + tria_temp.set_all_manifold_ids(0); + tria_temp.set_manifold(0, SphericalManifold()); + GridGenerator::convert_hypercube_to_simplex_mesh(tria_temp, + triangulation); + for (const auto i : tria_temp.get_manifold_ids()) + if (i != numbers::flat_manifold_id) + triangulation.set_manifold(i, tria_temp.get_manifold(i)); + +#else + GridGenerator::extract_boundary_mesh(volume_mesh, + triangulation, + boundary_ids); + triangulation.set_all_manifold_ids(0); + triangulation.set_manifold(0, SphericalManifold()); +#endif + triangulation.refine_global(4); + } + + std::cout << "Surface mesh has " << triangulation.n_active_cells() + << " cells." << std::endl; + + dof_handler.distribute_dofs(fe); + + std::cout << "Surface mesh has " << dof_handler.n_dofs() + << " degrees of freedom." << std::endl; + + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + } + + template + void + LaplaceBeltramiProblem::assemble_system() + { + system_matrix = 0; + system_rhs = 0; +#ifdef USE_SIMPLEX + const Simplex::QGauss quadrature_formula(2 * fe.degree); +#else + const QGauss quadrature_formula(2 * fe.degree); +#endif + FEValues fe_values(mapping, + fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + const unsigned int n_q_points = quadrature_formula.size(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector rhs_values(n_q_points); + std::vector local_dof_indices(dofs_per_cell); + + RightHandSide rhs; + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit(cell); + + rhs.value_list(fe_values.get_quadrature_points(), rhs_values); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + cell_matrix(i, j) += fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * + fe_values.JxW(q_point); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + cell_rhs(i) += fe_values.shape_value(i, q_point) * + rhs_values[q_point] * fe_values.JxW(q_point); + + cell->get_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + system_matrix.add(local_dof_indices[i], + local_dof_indices[j], + cell_matrix(i, j)); + + system_rhs(local_dof_indices[i]) += cell_rhs(i); + } + } + + std::map boundary_values; + VectorTools::interpolate_boundary_values( + mapping, dof_handler, 0, Solution(), boundary_values); + + MatrixTools::apply_boundary_values( + boundary_values, system_matrix, solution, system_rhs, false); + } + + template + void + LaplaceBeltramiProblem::solve() + { + SolverControl solver_control(solution.size(), 1e-7 * system_rhs.l2_norm()); + SolverCG> cg(solver_control); + + PreconditionSSOR> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve(system_matrix, solution, system_rhs, preconditioner); + } + + template + void + LaplaceBeltramiProblem::output_results() const + { + DataOut> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector( + solution, + "solution", + DataOut>::type_dof_data); + data_out.build_patches(mapping, mapping.get_degree()); + + const std::string filename = + "solution-" + std::to_string(spacedim) + "d.vtk"; + std::ofstream output(filename); + data_out.write_vtk(output); + } + + template + void + LaplaceBeltramiProblem::compute_error() const + { + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + Solution(), + difference_per_cell, +#ifdef USE_SIMPLEX + Simplex::QGauss(2 * fe.degree + 1), +#else + QGauss(2 * fe.degree + + 1), // This also works on + // triangular meshes +#endif + VectorTools::H1_norm); + + double h1_error = VectorTools::compute_global_error(triangulation, + difference_per_cell, + VectorTools::H1_norm); + std::cout << "H1 error = " << h1_error << std::endl; + } + + template + void + LaplaceBeltramiProblem::run() + { + make_grid_and_dofs(); + assemble_system(); + solve(); + output_results(); + compute_error(); + } +} // namespace Step38 + +int +main() +{ + try + { + using namespace Step38; + + LaplaceBeltramiProblem<3> laplace_beltrami(1); + laplace_beltrami.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-38.with_simplex_support=on.output b/tests/simplex/step-38.with_simplex_support=on.output new file mode 100644 index 0000000000..9bd139ac6b --- /dev/null +++ b/tests/simplex/step-38.with_simplex_support=on.output @@ -0,0 +1,3 @@ +Surface mesh has 10240 cells. +Surface mesh has 5185 degrees of freedom. +H1 error = 0.384717 -- 2.39.5