From 23f2c4fc226f7b99c6dde7d72a6d9deb6a0c30c7 Mon Sep 17 00:00:00 2001 From: Olivier Guevremont Date: Fri, 15 Jan 2021 12:38:18 -0500 Subject: [PATCH] Add Step-04.cc converted to Simplex --- tests/simplex/step-04.cc | 315 ++++++++++++++++++ .../step-04.with_simplex_support=on.output | 15 + 2 files changed, 330 insertions(+) create mode 100644 tests/simplex/step-04.cc create mode 100644 tests/simplex/step-04.with_simplex_support=on.output diff --git a/tests/simplex/step-04.cc b/tests/simplex/step-04.cc new file mode 100644 index 0000000000..892c1efb04 --- /dev/null +++ b/tests/simplex/step-04.cc @@ -0,0 +1,315 @@ +/* --------------------------------------------------------------------- + * + * Copyright (C) 1999 - 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. + * + * --------------------------------------------------------------------- + + * + * Author: Wolfgang Bangerth, University of Heidelberg, 1999 + */ + + +// Step-04 on a simplex mesh. Following incompatible modifications had to be +// made: +// - Change the FE_Q to Simplex::FE_P +// - Put the MappingFE as a class member and use as an argument instead of +// default mapping +// - Change QGauss to Simplex::QGauss +// - Convert triangulation to a triangulation based on simplices + + +#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" + + +using namespace dealii; + +template +class Step4 +{ +public: + Step4(); + void + run(); + +private: + void + make_grid(); + void + setup_system(); + void + assemble_system(); + void + solve(); + void + output_results() const; + + Triangulation triangulation; + Simplex::FE_P fe; + DoFHandler dof_handler; + + MappingFE mapping; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; +}; + + +template +class RightHandSide : public Function +{ +public: + virtual double + value(const Point &p, const unsigned int component = 0) const override; +}; + + +template +class BoundaryValues : public Function +{ +public: + virtual double + value(const Point &p, const unsigned int component = 0) const override; +}; + + +template +double +RightHandSide::value(const Point &p, + const unsigned int /*component*/) const +{ + double return_value = 0.0; + for (unsigned int i = 0; i < dim; ++i) + return_value += 4.0 * std::pow(p(i), 4.0); + + return return_value; +} + + +template +double +BoundaryValues::value(const Point &p, + const unsigned int /*component*/) const +{ + return p.square(); +} + + + +template +Step4::Step4() + : fe(1) + , dof_handler(triangulation) + , mapping(fe) +{} + + +template +void +Step4::make_grid() +{ + Triangulation temp_tria; + GridGenerator::subdivided_hyper_cube(temp_tria, 8, -1., 1., false); + GridGenerator::convert_hypercube_to_simplex_mesh(temp_tria, + triangulation); + + deallog << " Number of active cells: " << triangulation.n_active_cells() + << std::endl + << " Total number of cells: " << triangulation.n_cells() + << std::endl; +} + + +template +void +Step4::setup_system() +{ + dof_handler.distribute_dofs(fe); + + deallog << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + DynamicSparsityPattern dsp(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 +Step4::assemble_system() +{ + Simplex::QGauss quadrature_formula(fe.degree + 1); + + RightHandSide right_hand_side; + + const UpdateFlags flag = update_JxW_values | update_values | + update_gradients | update_quadrature_points; + + FEValues fe_values(mapping, fe, quadrature_formula, flag); + + + const unsigned int dofs_per_cell = fe.n_dofs_per_cell(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + + std::vector local_dof_indices(dofs_per_cell); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; + + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + for (const unsigned int i : fe_values.dof_indices()) + { + for (const unsigned int j : fe_values.dof_indices()) + cell_matrix(i, j) += + (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q) + fe_values.shape_grad(j, q_index) * // grad phi_j(x_q) + fe_values.JxW(q_index)); // dx + + const auto x_q = fe_values.quadrature_point(q_index); + cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q) + right_hand_side.value(x_q) * // f(x_q) + fe_values.JxW(q_index)); // dx + } + + cell->get_dof_indices(local_dof_indices); + for (const unsigned int i : fe_values.dof_indices()) + { + for (const unsigned int j : fe_values.dof_indices()) + 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, BoundaryValues(), boundary_values); + MatrixTools::apply_boundary_values(boundary_values, + system_matrix, + solution, + system_rhs); +} + + + +template +void +Step4::solve() +{ + SolverControl solver_control(1000, 1e-12); + SolverCG> solver(solver_control); + solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity()); + + deallog << " " << solver_control.last_step() + << " CG iterations needed to obtain convergence." << std::endl; +} + + + +template +void +Step4::output_results() const +{ + DataOut data_out; + + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + + data_out.build_patches(mapping, 0); + + std::ofstream output(dim == 2 ? "solution-2d.vtk" : "solution-3d.vtk"); + // data_out.write_vtk(output); +} + + +template +void +Step4::run() +{ + deallog << "Solving problem in " << dim << " space dimensions." << std::endl; + + make_grid(); + setup_system(); + assemble_system(); + solve(); + output_results(); +} + + + +int +main() +{ + initlog(); + deallog.depth_console(2); + { + Step4<2> laplace_problem_2d; + laplace_problem_2d.run(); + } + + { + Step4<3> laplace_problem_3d; + laplace_problem_3d.run(); + } + + return 0; +} diff --git a/tests/simplex/step-04.with_simplex_support=on.output b/tests/simplex/step-04.with_simplex_support=on.output new file mode 100644 index 0000000000..f1a975bbd9 --- /dev/null +++ b/tests/simplex/step-04.with_simplex_support=on.output @@ -0,0 +1,15 @@ + +DEAL::Solving problem in 2 space dimensions. +DEAL:: Number of active cells: 512 +DEAL:: Total number of cells: 512 +DEAL:: Number of degrees of freedom: 289 +DEAL:cg::Starting value 11.6837 +DEAL:cg::Convergence step 33 value 2.19749e-16 +DEAL:: 33 CG iterations needed to obtain convergence. +DEAL::Solving problem in 3 space dimensions. +DEAL:: Number of active cells: 12288 +DEAL:: Total number of cells: 12288 +DEAL:: Number of degrees of freedom: 2457 +DEAL:cg::Starting value 13.9290 +DEAL:cg::Convergence step 46 value 8.50319e-13 +DEAL:: 46 CG iterations needed to obtain convergence. -- 2.39.5