From dbada94f8a71f59cab6cc0339b6ff55460d2bd57 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Thu, 13 Sep 2018 22:53:23 +0200 Subject: [PATCH] Fix VectorTools::project_boundary_values_div_conforming for more than one FE_RT element --- .../deal.II/numerics/vector_tools.templates.h | 38 +++-- tests/numerics/project_bv_div_conf_02.cc | 132 ++++++++++++++++++ tests/numerics/project_bv_div_conf_02.output | 41 ++++++ 3 files changed, 203 insertions(+), 8 deletions(-) create mode 100644 tests/numerics/project_bv_div_conf_02.cc create mode 100644 tests/numerics/project_bv_div_conf_02.output diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index 755f61f47f..2a06dd02d0 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -6411,9 +6411,14 @@ namespace VectorTools for (unsigned int i = 0; i < fe.dofs_per_face; ++i) dof_values(i) += - tmp * - (normals[q_point] * - fe_values[vec].value(fe.face_to_cell_index(i, face), q_point)); + tmp * (normals[q_point] * + fe_values[vec].value( + fe.face_to_cell_index(i, + face, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face)), + q_point)); } std::vector face_dof_indices(fe.dofs_per_face); @@ -6424,7 +6429,13 @@ namespace VectorTools // Copy the computed values in the AffineConstraints only, if the degree // of freedom is not already constrained. for (unsigned int i = 0; i < fe.dofs_per_face; ++i) - if (!(constraints.is_constrained(face_dof_indices[i]))) + if (!(constraints.is_constrained(face_dof_indices[i])) && + fe.get_nonzero_components(fe.face_to_cell_index( + i, + face, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face)))[first_vector_component]) { constraints.add_line(face_dof_indices[i]); @@ -6509,9 +6520,14 @@ namespace VectorTools for (unsigned int i = 0; i < fe.dofs_per_face; ++i) dof_values_local(i) += - tmp * - (normals[q_point] * - fe_values[vec].value(fe.face_to_cell_index(i, face), q_point)); + tmp * (normals[q_point] * + fe_values[vec].value( + fe.face_to_cell_index(i, + face, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face)), + q_point)); } std::vector face_dof_indices(fe.dofs_per_face); @@ -6520,7 +6536,13 @@ namespace VectorTools cell->active_fe_index()); for (unsigned int i = 0; i < fe.dofs_per_face; ++i) - if (projected_dofs[face_dof_indices[i]] < fe.degree) + if (projected_dofs[face_dof_indices[i]] < fe.degree && + fe.get_nonzero_components(fe.face_to_cell_index( + i, + face, + cell->face_orientation(face), + cell->face_flip(face), + cell->face_rotation(face)))[first_vector_component]) { dof_values[face_dof_indices[i]] = dof_values_local(i); projected_dofs[face_dof_indices[i]] = fe.degree; diff --git a/tests/numerics/project_bv_div_conf_02.cc b/tests/numerics/project_bv_div_conf_02.cc new file mode 100644 index 0000000000..83e61bd60f --- /dev/null +++ b/tests/numerics/project_bv_div_conf_02.cc @@ -0,0 +1,132 @@ +// --------------------------------------------------------------------- +// +// 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 for the +// case that the DoFHandler constains is more than one FE_RaviartThomas element. + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + + +template +class BoundaryFunctionDisp : public dealii::Function +{ +public: + BoundaryFunctionDisp() + : Function(dim) + {} + + virtual double + value(const Point & /*point*/, + const unsigned int /*component*/) const override + { + return 1.; + } +}; + +template +class BoundaryFunctionVelo : public dealii::Function +{ +public: + BoundaryFunctionVelo() + : Function(dim) + {} + + virtual double + value(const Point & /*point*/, + const unsigned int /* component */) const override + { + return -1.; + } +}; + +template +void +test_boundary_values(const FiniteElement &fe) +{ + Triangulation triangulation; + + GridGenerator::hyper_cube(triangulation, -1., 1., true); + triangulation.refine_global(1); + + DoFHandler dof_handler(triangulation); + + dof_handler.distribute_dofs(fe); + + BoundaryFunctionDisp boundary_function_disp; + BoundaryFunctionVelo boundary_function_velo; + + AffineConstraints constraints; + + { + constraints.clear(); + VectorTools::project_boundary_values_div_conforming( + dof_handler, + 0, /*first_vector_component*/ + boundary_function_disp, + 0, /*bdry_id*/ + constraints); + VectorTools::project_boundary_values_div_conforming( + dof_handler, + dim, /*first_vector_component*/ + boundary_function_velo, + 0, /*bdry_id*/ + constraints); + constraints.close(); + } + + constraints.print(deallog.get_file_stream()); +} + +int +main() +{ + initlog(); + { + constexpr unsigned int dim = 2; + FE_RaviartThomas u(1); + FE_DGQ p(1); + FESystem fesys(u, 2, p, 1); + test_boundary_values(fesys); + } + + { + constexpr unsigned int dim = 3; + FE_RaviartThomas u(1); + FE_DGQ p(1); + FESystem fesys(u, 2, p, 1); + test_boundary_values(fesys); + } +} diff --git a/tests/numerics/project_bv_div_conf_02.output b/tests/numerics/project_bv_div_conf_02.output new file mode 100644 index 0000000000..96f2a107b8 --- /dev/null +++ b/tests/numerics/project_bv_div_conf_02.output @@ -0,0 +1,41 @@ + + 0 = 1.00000 + 1 = 0 + 2 = -1.00000 + 3 = 0 + 52 = 1.00000 + 53 = 0 + 54 = -1.00000 + 55 = 0 + 0 = 1.00000 + 1 = 0 + 2 = 0 + 3 = 0 + 4 = -1.00000 + 5 = 0 + 6 = 0 + 7 = 0 + 152 = 1.00000 + 153 = 0 + 154 = 0 + 155 = 0 + 156 = -1.00000 + 157 = 0 + 158 = 0 + 159 = 0 + 288 = 1.00000 + 289 = 0 + 290 = 0 + 291 = 0 + 292 = -1.00000 + 293 = 0 + 294 = 0 + 295 = 0 + 424 = 1.00000 + 425 = 0 + 426 = 0 + 427 = 0 + 428 = -1.00000 + 429 = 0 + 430 = 0 + 431 = 0 -- 2.39.5