From 57a4d38e2bd4daf8ab4b4fd2679a5a9a2c85c162 Mon Sep 17 00:00:00 2001 From: Nils Much Date: Wed, 13 Jan 2021 13:39:44 +0100 Subject: [PATCH] Add step-8 to simplex test --- tests/simplex/step-08.cc | 384 ++++++++++++++++++ .../step-08.with_simplex_support=on.output | 15 + 2 files changed, 399 insertions(+) create mode 100644 tests/simplex/step-08.cc create mode 100644 tests/simplex/step-08.with_simplex_support=on.output diff --git a/tests/simplex/step-08.cc b/tests/simplex/step-08.cc new file mode 100644 index 0000000000..eececbf1b4 --- /dev/null +++ b/tests/simplex/step-08.cc @@ -0,0 +1,384 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Step-08 on a simplex mesh. Following incompatible modifications had to be +// made: +// - Change the FE_Q to Simplex::FE_P. +// - Change QGauss to Simplex::QGauss. +// - Use MappingFE (Do not use default mapping). +// - Convert triangulation to a triangulation based on simplices. +// - Use refine_global() instead of execute_coarsening_and_refinement(). +// - Reduce number of cycles from 7 to 4 to reduce runtime. + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +namespace Step8 +{ + using namespace dealii; + + template + class ElasticProblem + { + public: + ElasticProblem(); + void + run(); + + private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + void + output_results(const unsigned int cycle) const; + + Triangulation triangulation; + DoFHandler dof_handler; + MappingFE mapping; + + FESystem fe; + + AffineConstraints constraints; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + Vector solution; + Vector system_rhs; + }; + + template + void + right_hand_side(const std::vector> &points, + std::vector> & values) + { + Assert(values.size() == points.size(), + ExcDimensionMismatch(values.size(), points.size())); + Assert(dim >= 2, ExcNotImplemented()); + + Point point_1, point_2; + point_1(0) = 0.5; + point_2(0) = -0.5; + + for (unsigned int point_n = 0; point_n < points.size(); ++point_n) + { + if (((points[point_n] - point_1).norm_square() < 0.2 * 0.2) || + ((points[point_n] - point_2).norm_square() < 0.2 * 0.2)) + values[point_n][0] = 1.0; + else + values[point_n][0] = 0.0; + + if (points[point_n].norm_square() < 0.2 * 0.2) + values[point_n][1] = 1.0; + else + values[point_n][1] = 0.0; + } + } + + + template + ElasticProblem::ElasticProblem() + : dof_handler(triangulation) + , mapping(Simplex::FE_P(1)) + , fe(Simplex::FE_P(1), dim) + {} + + template + void + ElasticProblem::setup_system() + { + dof_handler.distribute_dofs(fe); + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); + + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + VectorTools::interpolate_boundary_values( + mapping, dof_handler, 0, Functions::ZeroFunction(dim), constraints); + constraints.close(); + + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + sparsity_pattern.copy_from(dsp); + + system_matrix.reinit(sparsity_pattern); + } + + template + void + ElasticProblem::assemble_system() + { + Simplex::QGauss quadrature_formula(fe.degree + 1); + 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 local_dof_indices(dofs_per_cell); + + std::vector lambda_values(n_q_points); + std::vector mu_values(n_q_points); + + Functions::ConstantFunction lambda(1.), mu(1.); + + std::vector> rhs_values(n_q_points); + + for (const auto &cell : dof_handler.active_cell_iterators()) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit(cell); + + lambda.value_list(fe_values.get_quadrature_points(), lambda_values); + mu.value_list(fe_values.get_quadrature_points(), mu_values); + right_hand_side(fe_values.get_quadrature_points(), rhs_values); + + for (const unsigned int i : fe_values.dof_indices()) + { + const unsigned int component_i = + fe.system_to_component_index(i).first; + + for (const unsigned int j : fe_values.dof_indices()) + { + const unsigned int component_j = + fe.system_to_component_index(j).first; + + for (const unsigned int q_point : + fe_values.quadrature_point_indices()) + { + cell_matrix(i, j) += + ((fe_values.shape_grad(i, q_point)[component_i] * + fe_values.shape_grad(j, q_point)[component_j] * + lambda_values[q_point]) + + (fe_values.shape_grad(i, q_point)[component_j] * + fe_values.shape_grad(j, q_point)[component_i] * + mu_values[q_point]) + + ((component_i == component_j) ? // + (fe_values.shape_grad(i, q_point) * // + fe_values.shape_grad(j, q_point) * // + mu_values[q_point]) : // + 0) // + ) * // + fe_values.JxW(q_point); // + } + } + } + + for (const unsigned int i : fe_values.dof_indices()) + { + const unsigned int component_i = + fe.system_to_component_index(i).first; + + for (const unsigned int q_point : + fe_values.quadrature_point_indices()) + cell_rhs(i) += fe_values.shape_value(i, q_point) * + rhs_values[q_point][component_i] * + fe_values.JxW(q_point); + } + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global( + cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs); + } + } + + template + void + ElasticProblem::solve() + { + SolverControl solver_control(1000, 1e-12); + SolverCG> cg(solver_control); + + PreconditionSSOR> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + cg.solve(system_matrix, solution, system_rhs, preconditioner); + + constraints.distribute(solution); + } + + template + void + ElasticProblem::refine_grid() + { + Vector estimated_error_per_cell(triangulation.n_active_cells()); + + Simplex::QGauss quadrature(fe.degree + 1); + KellyErrorEstimator::estimate( + mapping, dof_handler, quadrature, {}, solution, estimated_error_per_cell); + + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + estimated_error_per_cell, + 0.3, + 0.03); + + triangulation.refine_global(); + } + + template + void + ElasticProblem::output_results(const unsigned int cycle) const + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + + std::vector solution_names; + switch (dim) + { + case 1: + solution_names.emplace_back("displacement"); + break; + case 2: + solution_names.emplace_back("x_displacement"); + solution_names.emplace_back("y_displacement"); + break; + case 3: + solution_names.emplace_back("x_displacement"); + solution_names.emplace_back("y_displacement"); + solution_names.emplace_back("z_displacement"); + break; + default: + Assert(false, ExcNotImplemented()); + } + + data_out.add_data_vector(solution, solution_names); + data_out.build_patches(mapping); + + std::ofstream output("solution-" + std::to_string(cycle) + ".vtk"); + data_out.write_vtk(output); + } + + template + void + ElasticProblem::run() + { + for (unsigned int cycle = 0; cycle < 5; ++cycle) + { + std::cout << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { + Triangulation temp; + GridGenerator::hyper_cube(temp, -1, 1); + GridGenerator::convert_hypercube_to_simplex_mesh(temp, + triangulation); + } + else + refine_grid(); + + std::cout << " Number of active cells: " + << triangulation.n_active_cells() << std::endl; + + setup_system(); + + std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + assemble_system(); + solve(); + output_results(cycle); + } + } +} // namespace Step8 + +int +main() +{ + try + { + Step8::ElasticProblem<2> elastic_problem_2d; + elastic_problem_2d.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-08.with_simplex_support=on.output b/tests/simplex/step-08.with_simplex_support=on.output new file mode 100644 index 0000000000..a92534cf0b --- /dev/null +++ b/tests/simplex/step-08.with_simplex_support=on.output @@ -0,0 +1,15 @@ +Cycle 0: + Number of active cells: 8 + Number of degrees of freedom: 18 +Cycle 1: + Number of active cells: 32 + Number of degrees of freedom: 50 +Cycle 2: + Number of active cells: 128 + Number of degrees of freedom: 162 +Cycle 3: + Number of active cells: 512 + Number of degrees of freedom: 578 +Cycle 4: + Number of active cells: 2048 + Number of degrees of freedom: 2178 -- 2.39.5