From 9328696b16f6922dfa4f16ca112d7307fad421d6 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Sat, 15 Sep 2018 18:21:39 +0200 Subject: [PATCH] Add test --- tests/numerics/project_bv_div_conf_03.cc | 263 +++++++++++++++++++ tests/numerics/project_bv_div_conf_03.output | 31 +++ 2 files changed, 294 insertions(+) create mode 100644 tests/numerics/project_bv_div_conf_03.cc create mode 100644 tests/numerics/project_bv_div_conf_03.output diff --git a/tests/numerics/project_bv_div_conf_03.cc b/tests/numerics/project_bv_div_conf_03.cc new file mode 100644 index 0000000000..4b204244ba --- /dev/null +++ b/tests/numerics/project_bv_div_conf_03.cc @@ -0,0 +1,263 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test VectorTools::project_boundary_values_div_conforming convergence rates +// for the case that the DoFHandler constains more than one +// FE_RaviartThomas element. + +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +template +class BoundaryFunctionDisp : public Function +{ +public: + BoundaryFunctionDisp() + : Function(dim) + {} + + virtual double + value(const Point & point, + const unsigned int component = 0) const override + { + Assert(component < dim, ExcNotImplemented()); + return std::sin(10 * point[component] + + point[component == 0 ? 1 : component - 1]); + } +}; + +template +class BoundaryFunctionVelo : public Function +{ +public: + BoundaryFunctionVelo() + : Function(dim) + {} + + virtual double + value(const Point & point, + const unsigned int component = 0) const override + { + Assert(component < dim, ExcNotImplemented()); + return std::cos(10 * point[component] + + point[component == 0 ? 1 : component - 1]); + } +}; + +template +class BoundaryFunctionPres : public Function +{ +public: + BoundaryFunctionPres() + : Function(2 * dim + 1) + {} + + // To make interpolate_boundary_values happy, the pressure is assigned to + // the 2*dim + 1th component + virtual double + value(const Point & point, + const unsigned int component = 0) const override + { + if (component != 2 * dim) + return std::numeric_limits::quiet_NaN(); + else + return std::sin(7.0 * point[0]) + std::cos(5.0 * point[dim - 1]); + } +}; + + +template +void +test_boundary_values(const FiniteElement &fe) +{ + Triangulation triangulation; + + GridGenerator::hyper_cube(triangulation); + MappingQ mapping(2); + QIterated face_quadrature(QTrapez<1>(), 5); + + double old_max_disp_error = 0.0; + double old_max_velo_error = 0.0; + double old_max_pres_error = 0.0; + for (unsigned int refinement_n = 0; refinement_n < 7 - dim; ++refinement_n) + { + triangulation.refine_global(1); + + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + const FEValuesExtractors::Vector displacements(0); + const FEValuesExtractors::Vector velocities(dim); + const FEValuesExtractors::Scalar pressure(2 * dim); + + BoundaryFunctionDisp boundary_function_disp; + BoundaryFunctionVelo boundary_function_velo; + BoundaryFunctionPres boundary_function_pres; + + AffineConstraints constraints; + + std::map *> boundary_functions; + boundary_functions[0] = &boundary_function_pres; + VectorTools::interpolate_boundary_values(mapping, + dof_handler, + boundary_functions, + constraints, + fe.component_mask(pressure)); + VectorTools::project_boundary_values_div_conforming( + dof_handler, + 0, /*first_vector_component*/ + boundary_function_disp, + 0, /*bdry_id*/ + constraints, + mapping); + VectorTools::project_boundary_values_div_conforming( + dof_handler, + dim, /*first_vector_component*/ + boundary_function_velo, + 0, /*bdry_id*/ + constraints, + mapping); + constraints.close(); + Vector solution(dof_handler.n_dofs()); + constraints.distribute(solution); + + FEFaceValues face_values(mapping, + fe, + face_quadrature, + update_values | update_normal_vectors | + update_quadrature_points); + double max_disp_error = 0.0; + double max_velo_error = 0.0; + double max_pres_error = 0.0; + + std::vector> cell_disp_values(face_quadrature.size()); + std::vector> cell_velo_values(face_quadrature.size()); + std::vector cell_pres_values(face_quadrature.size()); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + for (unsigned int face_n = 0; + face_n < GeometryInfo::faces_per_cell; + ++face_n) + { + auto face = cell->face(face_n); + if (face->at_boundary()) + { + face_values.reinit(cell, face_n); + face_values[displacements].get_function_values( + solution, cell_disp_values); + face_values[velocities].get_function_values(solution, + cell_velo_values); + face_values[pressure].get_function_values(solution, + cell_pres_values); + + for (unsigned int q_point_n = 0; + q_point_n < face_quadrature.size(); + ++q_point_n) + { + const Point q_point = + face_values.quadrature_point(q_point_n); + const Tensor<1, dim> n_vector = + face_values.normal_vector(q_point_n); + Tensor<1, dim> exact_disp; + Tensor<1, dim> exact_velo; + for (unsigned int dim_n = 0; dim_n < dim; ++dim_n) + { + exact_disp[dim_n] = + boundary_function_disp.value(q_point, dim_n); + exact_velo[dim_n] = + boundary_function_velo.value(q_point, dim_n); + } + const double exact_pres = + boundary_function_pres.value(q_point, 2 * dim); + + max_disp_error = + std::max(std::abs(cell_disp_values[q_point_n] * + n_vector - + exact_disp * n_vector), + max_disp_error); + max_velo_error = + std::max(std::abs(cell_velo_values[q_point_n] * + n_vector - + exact_velo * n_vector), + max_velo_error); + max_pres_error = + std::max(std::abs(cell_pres_values[q_point_n] - + exact_pres), + max_pres_error); + } + } + } + } + + if (refinement_n != 0) + { + deallog << std::endl; + deallog << "disp ratio: " << old_max_disp_error / max_disp_error + << std::endl; + deallog << "velo ratio: " << old_max_velo_error / max_velo_error + << std::endl; + deallog << "pres ratio: " << old_max_pres_error / max_pres_error + << std::endl; + } + + old_max_disp_error = max_disp_error; + old_max_velo_error = max_velo_error; + old_max_pres_error = max_pres_error; + } +} + +int +main() +{ + initlog(); + { + deallog << "2d:" << std::endl; + constexpr unsigned int dim = 2; + FE_RaviartThomas u(1); + FE_Q p(1); + FESystem fesys(u, 2, p, 1); + test_boundary_values(fesys); + } + + { + deallog << "3d:" << std::endl; + constexpr unsigned int dim = 3; + FE_RaviartThomas u(1); + FE_Q p(1); + FESystem fesys(u, 2, p, 1); + test_boundary_values(fesys); + } +} diff --git a/tests/numerics/project_bv_div_conf_03.output b/tests/numerics/project_bv_div_conf_03.output new file mode 100644 index 0000000000..8892464857 --- /dev/null +++ b/tests/numerics/project_bv_div_conf_03.output @@ -0,0 +1,31 @@ + +DEAL::2d: +DEAL:: +DEAL::disp ratio: 3.92698 +DEAL::velo ratio: 3.92520 +DEAL::pres ratio: 3.39935 +DEAL:: +DEAL::disp ratio: 3.98215 +DEAL::velo ratio: 3.98126 +DEAL::pres ratio: 3.67259 +DEAL:: +DEAL::disp ratio: 3.99576 +DEAL::velo ratio: 3.99531 +DEAL::pres ratio: 3.94304 +DEAL:: +DEAL::disp ratio: 3.99905 +DEAL::velo ratio: 3.99883 +DEAL::pres ratio: 3.98666 +DEAL::3d: +DEAL:: +DEAL::disp ratio: 3.92698 +DEAL::velo ratio: 3.92520 +DEAL::pres ratio: 3.12974 +DEAL:: +DEAL::disp ratio: 3.98215 +DEAL::velo ratio: 3.98126 +DEAL::pres ratio: 3.79870 +DEAL:: +DEAL::disp ratio: 3.99576 +DEAL::velo ratio: 3.99531 +DEAL::pres ratio: 3.91749 -- 2.39.5