From: David Wells Date: Sat, 13 Apr 2024 11:43:08 +0000 (-0400) Subject: Add a test. X-Git-Tag: v9.6.0-rc1~377^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=362ca484b8e8537081dd7fcb6d0276a185197c32;p=dealii.git Add a test. This fails on master due to twisted elements. --- diff --git a/tests/simplex/orientation_06.cc b/tests/simplex/orientation_06.cc new file mode 100644 index 0000000000..b71172c78b --- /dev/null +++ b/tests/simplex/orientation_06.cc @@ -0,0 +1,223 @@ +/* ------------------------------------------------------------------------ + * + * SPDX-License-Identifier: LGPL-2.1-or-later + * Copyright (C) 2023 by the deal.II authors + * + * This file is part of the deal.II library. + * + * Part of the source code is dual licensed under Apache-2.0 WITH + * LLVM-exception OR LGPL-2.1-or-later. Detailed license information + * governing the source code and code contributions can be found in + * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. + * + * ------------------------------------------------------------------------ + */ + +// Make sure that we get the correct convergence order for a combination of P3 +// and Q3 elements and and two reversed quadrilateral line orientations. + +#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" + +// u +class Exact : public Function<2> +{ +public: + double + value(const Point<2> &p, unsigned int /*component*/ = 0) const + { + return std::sin(p[0]) * std::sin(p[1]); + } +}; + +// -Laplacian(u) + u +class Forcing : public Function<2> +{ +public: + double + value(const Point<2> &p, unsigned int /*component*/ = 0) const + { + return 3.0 * std::sin(p[0]) * std::sin(p[1]); + } +}; + +int +main() +{ + initlog(); + + const std::vector> vertices{{0.0, 0.0}, + {1.0, 0.0}, + {0.0, 2.0}, + {3.0, 1.0}, + {3.0, 3.0}, + {2.0, 3.0}, + {0.0, 3.0}, + {3.0, 0.0}}; + + // faces 1 and 3 of the quadrilateral will have reversed orientations + std::vector> cells(5); + cells[0].vertices = {0u, 1u, 2u}; + cells[1].vertices = {3u, 4u, 5u}; + cells[2].vertices = {1u, 3u, 7u}; + cells[3].vertices = {5u, 2u, 6u}; + cells[4].vertices = {1u, 3u, 2u, 5u}; + + const hp::MappingCollection<2> mapping(MappingFE<2>(FE_SimplexP<2>(1)), + MappingFE<2>(FE_Q<2>(1))); + + const hp::FECollection<2> fe(FE_SimplexP<2>(3), + FE_Q<2>(QIterated<1>(QTrapezoid<1>(), 3))); + const hp::QCollection<2> quadrature_formula(QGaussSimplex<2>(4), + QGauss<2>(4)); + + Triangulation<2> tria; + GridTools::invert_cells_with_negative_measure(vertices, cells); + tria.create_triangulation(vertices, cells, SubCellData()); + + for (unsigned int r = 0; r < 5; ++r) + { + tria.refine_global(1); + DoFHandler<2> dof_handler(tria); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + if (cell->reference_cell() == ReferenceCells::Triangle) + cell->set_active_fe_index(0); + else if (cell->reference_cell() == ReferenceCells::Quadrilateral) + cell->set_active_fe_index(1); + } + dof_handler.distribute_dofs(fe); + + // set up matrices and vectors: + AffineConstraints constraints(dof_handler.locally_owned_dofs(), + dof_handler.locally_owned_dofs()); + VectorTools::interpolate_boundary_values( + mapping, dof_handler, 0, Exact(), constraints); + constraints.close(); + SparsityPattern sparsity_pattern; + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + sparsity_pattern.copy_from(dsp); + + SparseMatrix system_matrix(sparsity_pattern); + Vector system_rhs(dof_handler.n_dofs()); + Vector solution(dof_handler.n_dofs()); + + hp::FEValues<2> hp_fe_values(mapping, + fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + + Forcing forcing; + for (const auto &cell : dof_handler.active_cell_iterators()) + { + hp_fe_values.reinit(cell); + const auto &fe_values = hp_fe_values.get_present_fe_values(); + + const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell(); + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_rhs.reinit(dofs_per_cell); + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + 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) * + fe_values.shape_grad(j, q_index) + + fe_values.shape_value(i, q_index) * + fe_values.shape_value(j, q_index)) * + fe_values.JxW(q_index); + + for (const unsigned int i : fe_values.dof_indices()) + cell_rhs(i) += + fe_values.shape_value(i, q_index) * + forcing.value(fe_values.quadrature_point(q_index)) * + fe_values.JxW(q_index); + } + + constraints.distribute_local_to_global(cell_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + + SolverControl solver_control(1000, 1e-12); + SolverCG> solver(solver_control); + + // Cut down on a lot of solver iterations by interpolating the solution as + // the initial guess: + VectorTools::interpolate(mapping, dof_handler, Exact(), solution); + + PreconditionSSOR<> precondition; + precondition.initialize(system_matrix, + PreconditionSSOR<>::AdditionalData()); + deallog << "rhs l2 = " << system_rhs.l2_norm() << std::endl; + solver.solve(system_matrix, solution, system_rhs, precondition); + constraints.distribute(solution); + + Vector cell_errors(tria.n_active_cells()); + VectorTools::integrate_difference(mapping, + dof_handler, + solution, + Exact(), + cell_errors, + quadrature_formula, + VectorTools::L2_norm); + + deallog << "global error = " + << VectorTools::compute_global_error(tria, + cell_errors, + VectorTools::L2_norm) + << std::endl; + + std::ofstream out("out-" + std::to_string(r) + ".vtu"); + DataOut<2> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + data_out.add_data_vector(cell_errors, "cell_errors"); + data_out.build_patches(mapping); + data_out.write_vtu(out); + } +} diff --git a/tests/simplex/orientation_06.output b/tests/simplex/orientation_06.output new file mode 100644 index 0000000000..c921f99fba --- /dev/null +++ b/tests/simplex/orientation_06.output @@ -0,0 +1,21 @@ + +DEAL::rhs l2 = 2.27397 +DEAL:cg::Starting value 1.73734 +DEAL:cg::Convergence step 23 value 8.54649e-13 +DEAL::global error = 0.00166288 +DEAL::rhs l2 = 2.08466 +DEAL:cg::Starting value 2.45428 +DEAL:cg::Convergence step 39 value 8.26837e-13 +DEAL::global error = 0.000105423 +DEAL::rhs l2 = 2.69354 +DEAL:cg::Starting value 3.47449 +DEAL:cg::Convergence step 67 value 7.08130e-13 +DEAL::global error = 6.61142e-06 +DEAL::rhs l2 = 3.77050 +DEAL:cg::Starting value 4.91796 +DEAL:cg::Convergence step 105 value 9.99712e-13 +DEAL::global error = 4.13368e-07 +DEAL::rhs l2 = 5.33304 +DEAL:cg::Starting value 6.95866 +DEAL:cg::Convergence step 169 value 9.77791e-13 +DEAL::global error = 2.58314e-08