From: Marc Fehling Date: Thu, 21 Jan 2021 23:52:52 +0000 (-0700) Subject: simplex: step-06 X-Git-Tag: v9.3.0-rc1~510^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F11638%2Fhead;p=dealii.git simplex: step-06 --- diff --git a/tests/simplex/step-06.cc b/tests/simplex/step-06.cc new file mode 100644 index 0000000000..d355951591 --- /dev/null +++ b/tests/simplex/step-06.cc @@ -0,0 +1,394 @@ +// --------------------------------------------------------------------- +// +// 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-06 on a simplex mesh. +// +// The following changes had to be made due to incompatibilities: +// - Manifolds are disabled +// All manifold ids from the triangulation object generated with +// GridGenerator::hyper_ball() had to be removed +// - Error estimation has to be skipped +// KellyErrorEstimator requires QProjector::project_to_all_subfaces() +// to work with triangles +// - GridOut::write_gnuplot() has been skipped + + +#define USE_SIMPLEX +// #define OUTPUT_RESULTS + + +#ifdef USE_SIMPLEX +# include + +# include +# include +#else +# include + +# include +#endif + + +#ifdef OUTPUT_RESULTS +# include + +# include +#endif + + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include "../tests.h" + + + +using namespace dealii; + + +template +class Step6 +{ +public: + Step6(); + + void + run(); + +private: + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + void + output_results(const unsigned int cycle) const; + +#ifdef USE_SIMPLEX + Simplex::FE_P fe; +#else + FE_Q fe; +#endif + + Triangulation triangulation; + DoFHandler dof_handler; + + AffineConstraints constraints; + + SparseMatrix system_matrix; + SparsityPattern sparsity_pattern; + + Vector solution; + Vector system_rhs; +}; + + + +template +double +coefficient(const Point &p) +{ + if (p.square() < 0.5 * 0.5) + return 20; + else + return 1; +} + + + +template +Step6::Step6() + : fe(2) + , dof_handler(triangulation) +{} + + + +template +void +Step6::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(dof_handler, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, + dsp, + constraints, + /*keep_constrained_dofs = */ false); + sparsity_pattern.copy_from(dsp); + system_matrix.reinit(sparsity_pattern); +} + + + +template +void +Step6::assemble_system() +{ +#ifdef USE_SIMPLEX + const Simplex::QGauss quadrature_formula(fe.degree + 1); +#else + const QGauss quadrature_formula(fe.degree + 1); +#endif + + MappingFE mapping(fe); + 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(); + + 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()) + { + cell_matrix = 0; + cell_rhs = 0; + + fe_values.reinit(cell); + + for (const unsigned int q_index : fe_values.quadrature_point_indices()) + { + const double current_coefficient = + coefficient(fe_values.quadrature_point(q_index)); + for (const unsigned int i : fe_values.dof_indices()) + { + for (const unsigned int j : fe_values.dof_indices()) + cell_matrix(i, j) += + (current_coefficient * // a(x_q) + 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 + cell_rhs(i) += (1.0 * // f(x) + fe_values.shape_value(i, q_index) * // phi_i(x_q) + fe_values.JxW(q_index)); // dx + } + } + + 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 +Step6::solve() +{ + SolverControl solver_control(1000, 1e-12); + SolverCG> solver(solver_control); + + PreconditionSSOR> preconditioner; + preconditioner.initialize(system_matrix, 1.2); + + solver.solve(system_matrix, solution, system_rhs, preconditioner); + + constraints.distribute(solution); +} + + + +template +void +Step6::refine_grid() +{ + Vector estimated_error_per_cell(triangulation.n_active_cells()); + +#if false +# ifdef USE_SIMPLEX + KellyErrorEstimator::estimate(MappingFE(fe), + dof_handler, + Simplex::QGauss(fe.degree + 1), + {}, + solution, + estimated_error_per_cell); +# else + KellyErrorEstimator::estimate(MappingFE(fe), + dof_handler, + QGauss(fe.degree + 1), + {}, + solution, + estimated_error_per_cell); +# endif + + GridRefinement::refine_and_coarsen_fixed_number(triangulation, + estimated_error_per_cell, + 0.3, + 0.03); +#else + // incompatibility: KellyErrorEstimator does not yet work for simplices + // instead, we will mark cells for refinement by hand + for (const auto &cell : triangulation.active_cell_iterators()) + for (unsigned int v = 0; v < cell->n_vertices(); ++v) + { + const double vertex_square = cell->vertex(v).square(); + if (vertex_square > 0.48 * 0.48 && vertex_square < 0.52 * 0.52) + cell->set_refine_flag(); + } +#endif + + triangulation.execute_coarsening_and_refinement(); +} + + + +template +void +Step6::output_results(const unsigned int cycle) const +{ +#ifdef OUTPUT_RESULTS + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + data_out.build_patches(MappingFE(fe)); + + std::ofstream output("solution-" + std::to_string(cycle) + ".vtu"); + data_out.write_vtu(output); + } +#else + (void)cycle; +#endif +} + + + +template +void +Step6::run() +{ + for (unsigned int cycle = 0; cycle < 5; ++cycle) + { + deallog << "Cycle " << cycle << ':' << std::endl; + + if (cycle == 0) + { +#ifdef USE_SIMPLEX + Triangulation tria_quad; + GridGenerator::hyper_ball(tria_quad); + + // incompatibility: remove all manifold ids + for (const auto &cell : tria_quad.active_cell_iterators()) + { + cell->set_manifold_id(numbers::flat_manifold_id); + for (const auto &face : cell->face_iterators()) + face->set_manifold_id(numbers::flat_manifold_id); + } + + GridGenerator::convert_hypercube_to_simplex_mesh(tria_quad, + triangulation); +#else + GridGenerator::hyper_ball(triangulation); +#endif + + triangulation.refine_global(1); + } + else + refine_grid(); + + deallog << " Number of active cells: " + << triangulation.n_active_cells() << std::endl; + + setup_system(); + + deallog << " Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + assemble_system(); + solve(); + output_results(cycle); + } + + deallog << "OK" << std::endl; +} + + + +int +main() +{ + initlog(); + + try + { + Step6<2> laplace_problem_2d; + laplace_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-06.with_simplex_support=on.output b/tests/simplex/step-06.with_simplex_support=on.output new file mode 100644 index 0000000000..be395de651 --- /dev/null +++ b/tests/simplex/step-06.with_simplex_support=on.output @@ -0,0 +1,27 @@ + +DEAL::Cycle 0: +DEAL:: Number of active cells: 160 +DEAL:: Number of degrees of freedom: 337 +DEAL:cg::Starting value 0.129098 +DEAL:cg::Convergence step 63 value 3.90426e-13 +DEAL::Cycle 1: +DEAL:: Number of active cells: 208 +DEAL:: Number of degrees of freedom: 465 +DEAL:cg::Starting value 0.123383 +DEAL:cg::Convergence step 70 value 4.47635e-13 +DEAL::Cycle 2: +DEAL:: Number of active cells: 400 +DEAL:: Number of degrees of freedom: 945 +DEAL:cg::Starting value 0.116265 +DEAL:cg::Convergence step 95 value 9.94249e-13 +DEAL::Cycle 3: +DEAL:: Number of active cells: 904 +DEAL:: Number of degrees of freedom: 2161 +DEAL:cg::Starting value 0.109845 +DEAL:cg::Convergence step 144 value 8.27369e-13 +DEAL::Cycle 4: +DEAL:: Number of active cells: 2536 +DEAL:: Number of degrees of freedom: 5793 +DEAL:cg::Starting value 0.102828 +DEAL:cg::Convergence step 246 value 9.36235e-13 +DEAL::OK