From faeded4b15f3bca73802c019d6cf1aa5ddb33e94 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 23 Jan 2018 01:27:01 +0100 Subject: [PATCH] Fix periodicity constraints for multiply constrained DoFs --- source/dofs/dof_tools_constraints.cc | 343 ++++++++++++++++----------- tests/bits/periodicity_06.cc | 218 +++++++++++++++++ tests/bits/periodicity_06.output | 9 + tests/bits/periodicity_07.cc | 249 +++++++++++++++++++ tests/bits/periodicity_07.output | 5 + 5 files changed, 682 insertions(+), 142 deletions(-) create mode 100644 tests/bits/periodicity_06.cc create mode 100644 tests/bits/periodicity_06.output create mode 100644 tests/bits/periodicity_07.cc create mode 100644 tests/bits/periodicity_07.output diff --git a/source/dofs/dof_tools_constraints.cc b/source/dofs/dof_tools_constraints.cc index 00718d0672..4b2fe3cce1 100644 --- a/source/dofs/dof_tools_constraints.cc +++ b/source/dofs/dof_tools_constraints.cc @@ -1873,167 +1873,226 @@ namespace DoFTools // loop over all dofs on face 2 and constrain them against the ones on face 1 for (unsigned int i=0; i > *constraint_entries + = constraint_matrix.get_constraint_entries(new_dof); + if (constraint_entries->size()==1) + new_dof = (*constraint_entries)[0].first; + else + { + enter_constraint = true; + break; + } + } + else + { + enter_constraint = true; + break; + } } - else + + if (enter_constraint) { - is_identity_constrained = false; - identity_constraint_target = numbers::invalid_unsigned_int; - break; + constraint_matrix.add_line(dofs_1[j]); + constraint_matrix.add_entry(dofs_1[j], dofs_2[i], is_identity_constrained?1.0:-1.0); } } - } - - bool is_inverse_constrained = !is_identity_constrained; - unsigned int inverse_constraint_target = numbers::invalid_unsigned_int; - if (is_inverse_constrained) - for (unsigned int jj=0; jj > *constraint_entries + = constraint_matrix.get_constraint_entries(dofs_1[j]); + if (constraint_entries->size()==1 && (*constraint_entries)[0].first == dofs_2[i]) + { + if ((is_identity_constrained && std::abs((*constraint_entries)[0].second-1) > eps) || + (is_inverse_constrained && std::abs((*constraint_entries)[0].second+1) > eps)) + { + //this pair of constraints means that both dofs have to be constrained to 0. + constraint_matrix.add_line(dofs_2[i]); + } + } + else + { + // see if this would add an identity constraint cycle + types::global_dof_index new_dof = dofs_1[j]; + while (new_dof != dofs_2[i]) + if (constraint_matrix.is_constrained(new_dof)) + { + const std::vector > *constraint_entries + = constraint_matrix.get_constraint_entries(new_dof); + if (constraint_entries->size()==1) + new_dof = (*constraint_entries)[0].first; + else + { + enter_constraint = true; + break; + } + } + else + { + enter_constraint = true; + break; + } + } } - } - } - const unsigned int target = is_identity_constrained - ? identity_constraint_target - : inverse_constraint_target; - - // find out whether this dof also exists on face 1 - // if this is true and the constraint is no identity - // constraint to itself, set it to zero - bool constrained_set = false; - for (unsigned int j=0; j > *constraint_entries - = constraint_matrix.get_constraint_entries(dofs_1[j]); - if (constraint_entries->size()==1 && (*constraint_entries)[0].first == dofs_2[i]) - { - if ((is_identity_constrained && std::abs((*constraint_entries)[0].second-1) > eps) || - (is_inverse_constrained && std::abs((*constraint_entries)[0].second+1) > eps)) - { - //this pair of constraints means that both dofs have to be constrained to 0. - constraint_matrix.add_line(dofs_2[i]); - } - } - else - enter_constraint = true; - } - - if (enter_constraint) - { - constraint_matrix.add_line(dofs_2[i]); - constraint_matrix.add_entry(dofs_2[i], dofs_1[j], is_identity_constrained?1.0:-1.0); - } - } - else - { - // this is just a regular constraint. enter it piece by piece - constraint_matrix.add_line(dofs_2[i]); - for (unsigned int jj=0; jj +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +ConstraintMatrix make_constraint_matrix(const DoFHandler<2> &dof_handler, int version) +{ + constexpr int dim = 2; + ConstraintMatrix constraints; + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + std::vector::cell_iterator> > periodicity_vectorDof; + switch (version) + { + case 0: + GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof); + case 1: + GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof); + case 2: + GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof); + case 3: + GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof); + case 4: + GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof); + case 5: + GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof); + case 6: + GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof); + case 7: + GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof); + } + + DoFTools::make_periodicity_constraints >(periodicity_vectorDof, constraints); + + constraints.close(); + /* std::map > support_points; + DoFTools::map_dofs_to_support_points (MappingQ(1), dof_handler, support_points); + for (const auto &line: constraints.get_lines()) + for (const auto &entry: line.entries) + std::cout << "DoF " << line.index << " at " << support_points[line.index] + << " is constrained to " << " DoF " << entry.first << " at " + << support_points[entry.first] + << " with value " << entry.second << std::endl;*/ + return constraints; +} + +template +class PeriodicReference : public Function +{ +public: + PeriodicReference () : Function() {} + virtual double value (const Point &p, + const unsigned int component = 0) const override + { + if (dim==3) + return std::sin(p(0)+1.)*std::sin(p(1)+2.)*std::sin(p(2)+3.); + return std::sin(p(0)+1.)*std::sin(p(1)+2.); + } +}; + + +template +void get_point_value +(const DoFHandler &dof_handler, + const Point &point, + const Vector &solution, + Vector &value) +{ + VectorTools::point_value (dof_handler, solution, + point, value); +} + + +void check_periodicity(const DoFHandler<2> &dof_handler, Vector &solution, const unsigned int cycle) +{ + unsigned int n_points = 2; + for (unsigned int i = 0; i value1(1); + Vector value2(1); + + Point <2> point1; + point1(0)=-numbers::PI+2.*i/n_points+eps; + point1(1)=-numbers::PI; + Point <2> point2; + point2(0)=-numbers::PI+2.*i/n_points+eps; + point2(1)=numbers::PI; + + VectorTools::point_value (dof_handler, solution, point1, value1); + VectorTools::point_value (dof_handler, solution, point2, value2); + + if (std::abs(value2[0]-value1[0])>1e-8) + { + std::cout << point1 << "\t" << "fail" << std::endl; + std::cout< value1(1); + Vector value2(1); + + Point <2> point1; + point1(1)=-numbers::PI+2.*i/n_points+eps; + point1(0)=-numbers::PI; + Point <2> point2; + point2(1)=-numbers::PI+2.*i/n_points+eps; + point2(0)=numbers::PI; + + VectorTools::point_value (dof_handler, solution, point1, value1); + VectorTools::point_value (dof_handler, solution, point2, value2); + + if (std::abs(value2[0]-value1[0])>1e-8) + { + std::cout << point1 << "\t" << "fail" << std::endl; + std::cout< triangulation; + GridGenerator::hyper_cube (triangulation, -L, L, true); + + triangulation.refine_global(1); + typename Triangulation::active_cell_iterator cellBegin = triangulation.begin_active(); + cellBegin->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + + FE_Q fe(1); + DoFHandler dof_handler (triangulation); + dof_handler.distribute_dofs(fe); + + std::vector constraints(8); + + PeriodicReference periodic_function; + + std::vector > projection(8, Vector (dof_handler.n_dofs())); + + for (unsigned int i=0; i<8; ++i) + { + deallog << "Testing version " << i << std::endl; + constraints[i] = make_constraint_matrix (dof_handler, i); + VectorTools::project (dof_handler, constraints[i], QGauss(3), periodic_function, projection[i]); + check_periodicity(dof_handler, projection[i], i); + } +} + diff --git a/tests/bits/periodicity_06.output b/tests/bits/periodicity_06.output new file mode 100644 index 0000000000..d6ba424002 --- /dev/null +++ b/tests/bits/periodicity_06.output @@ -0,0 +1,9 @@ + +DEAL::Testing version 0 +DEAL::Testing version 1 +DEAL::Testing version 2 +DEAL::Testing version 3 +DEAL::Testing version 4 +DEAL::Testing version 5 +DEAL::Testing version 6 +DEAL::Testing version 7 diff --git a/tests/bits/periodicity_07.cc b/tests/bits/periodicity_07.cc new file mode 100644 index 0000000000..0f93f9eff1 --- /dev/null +++ b/tests/bits/periodicity_07.cc @@ -0,0 +1,249 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +// Make sure that periodic boundary conditions also work correctly if we have +// multiple periodic boundary pairs that meet at an edge. +// This tests the 3D case. + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +ConstraintMatrix make_constraint_matrix(const DoFHandler<3> &dof_handler, int version) +{ + constexpr int dim = 3; + + ConstraintMatrix constraints; + constraints.clear(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + std::vector::cell_iterator> > periodicity_vectorDof; + switch (version) + { + case 0: + GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 4, 5, 2, periodicity_vectorDof); + case 1: + GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 5, 4, 2, periodicity_vectorDof); + case 2: + GridTools::collect_periodic_faces(dof_handler, 4, 5, 2, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 2, 3, 1, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 0, 1, 0, periodicity_vectorDof); + case 3: + GridTools::collect_periodic_faces(dof_handler, 5, 4, 2, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 3, 2, 1, periodicity_vectorDof); + GridTools::collect_periodic_faces(dof_handler, 1, 0, 0, periodicity_vectorDof); + } + + DoFTools::make_periodicity_constraints >(periodicity_vectorDof, constraints); + + constraints.close(); + /*std::map > support_points; + DoFTools::map_dofs_to_support_points (MappingQ(1), dof_handler, support_points); + for (const auto &line: constraints.get_lines()) + for (const auto &entry: line.entries) + std::cout << "DoF " << line.index << " at " << support_points[line.index] + << " is constrained to " << " DoF " << entry.first << " at " + << support_points[entry.first] + << " with value " << entry.second << std::endl;*/ + return constraints; +} + +template +class PeriodicReference : public Function +{ +public: + PeriodicReference () : Function() {} + virtual double value (const Point &p, + const unsigned int component = 0) const override + { + if (dim==3) + return std::sin(p(0)+1.)*std::sin(p(1)+2.)*std::sin(p(2)+3.); + return std::sin(p(0)+1.)*std::sin(p(1)+2.); + } +}; + + +template +void get_point_value +(const DoFHandler &dof_handler, + const Point &point, + const Vector &solution, + Vector &value) +{ + VectorTools::point_value (dof_handler, solution, + point, value); +} + + +void check_periodicity(const DoFHandler<3> &dof_handler, Vector &solution, const unsigned int cycle) +{ + unsigned int n_points = 2; + for (unsigned int i = 0; i value1(1); + Vector value2(1); + + Point<3> point1; + point1(0)=-numbers::PI+2.*i/n_points+eps; + point1(1)=-numbers::PI; + point1(2)=-numbers::PI+2.*j/n_points+eps; + Point<3> point2; + point2(0)=-numbers::PI+2.*i/n_points+eps; + point2(1)=numbers::PI; + point2(2)=-numbers::PI+2.*j/n_points+eps; + + VectorTools::point_value (dof_handler, solution, point1, value1); + VectorTools::point_value (dof_handler, solution, point2, value2); + + if (std::abs(value2[0]-value1[0])>1e-8) + { + std::cout << point1 << "\t" << "fail" << std::endl; + std::cout< value1(1); + Vector value2(1); + + Point <3> point1; + point1(2)=-numbers::PI+2.*j/n_points+eps; + point1(1)=-numbers::PI+2.*i/n_points+eps; + point1(0)=-numbers::PI; + Point <3> point2; + point2(2)=-numbers::PI+2.*j/n_points+eps; + point2(1)=-numbers::PI+2.*i/n_points+eps; + point2(0)=numbers::PI; + + VectorTools::point_value (dof_handler, solution, point1, value1); + VectorTools::point_value (dof_handler, solution, point2, value2); + + if (std::abs(value2[0]-value1[0])>1e-8) + { + std::cout << point1 << "\t" << "fail" << std::endl; + std::cout< value1(1); + Vector value2(1); + + Point <3> point1; + point1(0)=-numbers::PI+2.*j/n_points+eps; + point1(1)=-numbers::PI+2.*i/n_points+eps; + point1(2)=-numbers::PI; + Point <3> point2; + point2(0)=-numbers::PI+2.*j/n_points+eps; + point2(1)=-numbers::PI+2.*i/n_points+eps; + point2(2)=numbers::PI; + + VectorTools::point_value (dof_handler, solution, point1, value1); + VectorTools::point_value (dof_handler, solution, point2, value2); + + if (std::abs(value2[0]-value1[0])>1e-8) + { + std::cout << point1 << "\t" << "fail" << std::endl; + std::cout< triangulation; + GridGenerator::hyper_cube (triangulation, -L, L, true); + + triangulation.refine_global(1); + typename Triangulation::active_cell_iterator cellBegin = triangulation.begin_active(); + cellBegin->set_refine_flag(); + triangulation.execute_coarsening_and_refinement(); + + FE_Q fe(1); + DoFHandler dof_handler (triangulation); + dof_handler.distribute_dofs(fe); + + std::vector constraints(4); + + PeriodicReference periodic_function; + + std::vector > projection(4, Vector (dof_handler.n_dofs())); + + for (unsigned int i=0; i<4; ++i) + { + deallog << "Testing version " << i << std::endl; + constraints[i] = make_constraint_matrix (dof_handler, i); + VectorTools::project (dof_handler, constraints[i], QGauss(3), periodic_function, projection[i]); + check_periodicity(dof_handler, projection[i], i); + } +} + diff --git a/tests/bits/periodicity_07.output b/tests/bits/periodicity_07.output new file mode 100644 index 0000000000..0132ef2ce9 --- /dev/null +++ b/tests/bits/periodicity_07.output @@ -0,0 +1,5 @@ + +DEAL::Testing version 0 +DEAL::Testing version 1 +DEAL::Testing version 2 +DEAL::Testing version 3 -- 2.39.5