// loop over all dofs on face 2 and constrain them against the ones on face 1
for (unsigned int i=0; i<dofs_per_face; ++i)
- if (!constraint_matrix.is_constrained(dofs_2[i]))
- if ((component_mask.n_selected_components(fe.n_components())
- == fe.n_components())
- ||
- component_mask[fe.face_system_to_component_index(i).first])
- {
- // as mentioned in the comment above this function, we need
- // to be careful about treating identity constraints differently.
- // consequently, find out whether this dof 'i' will be
- // identity constrained
- //
- // to check whether this is the case, first see whether there are
- // any weights other than 0 and 1, then in a first stage make sure
- // that if so there is only one weight equal to 1
- //
- // afterwards do the same for constraints of type dof1=-dof2
- bool is_identity_constrained = true;
- const double eps = 1.e-13;
+ if ((component_mask.n_selected_components(fe.n_components())
+ == fe.n_components())
+ ||
+ component_mask[fe.face_system_to_component_index(i).first])
+ {
+ // as mentioned in the comment above this function, we need
+ // to be careful about treating identity constraints differently.
+ // consequently, find out whether this dof 'i' will be
+ // identity constrained
+ //
+ // to check whether this is the case, first see whether there are
+ // any weights other than 0 and 1, then in a first stage make sure
+ // that if so there is only one weight equal to 1
+ //
+ // afterwards do the same for constraints of type dof1=-dof2
+ bool is_identity_constrained = true;
+ const double eps = 1.e-13;
+ for (unsigned int jj=0; jj<dofs_per_face; ++jj)
+ if (((std::abs(transformation(i,jj)) < eps) ||
+ (std::abs(transformation(i,jj)-1) < eps)) == false)
+ {
+ is_identity_constrained = false;
+ break;
+ }
+ unsigned int identity_constraint_target = numbers::invalid_unsigned_int;
+ if (is_identity_constrained == true)
+ {
+ bool one_identity_found = false;
+ for (unsigned int jj=0; jj<dofs_per_face; ++jj)
+ if (std::abs(transformation(i,jj)-1.) < eps)
+ {
+ if (one_identity_found == false)
+ {
+ one_identity_found = true;
+ identity_constraint_target = jj;
+ }
+ else
+ {
+ is_identity_constrained = false;
+ identity_constraint_target = numbers::invalid_unsigned_int;
+ break;
+ }
+ }
+ }
+
+ 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<dofs_per_face; ++jj)
if (((std::abs(transformation(i,jj)) < eps) ||
- (std::abs(transformation(i,jj)-1) < eps)) == false)
+ (std::abs(transformation(i,jj)+1) < eps)) == false)
{
- is_identity_constrained = false;
+ is_inverse_constrained = false;
break;
}
- unsigned int identity_constraint_target = numbers::invalid_unsigned_int;
- if (is_identity_constrained == true)
- {
- bool one_identity_found = false;
- for (unsigned int jj=0; jj<dofs_per_face; ++jj)
- if (std::abs(transformation(i,jj)-1.) < eps)
+ if (is_inverse_constrained)
+ {
+ bool one_identity_found = false;
+ for (unsigned int jj=0; jj<dofs_per_face; ++jj)
+ if (std::abs(transformation(i,jj)+1) < eps)
+ {
+ if (one_identity_found == false)
+ {
+ one_identity_found = true;
+ inverse_constraint_target = jj;
+ }
+ else
+ {
+ is_inverse_constrained = false;
+ inverse_constraint_target = numbers::invalid_unsigned_int;
+ 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 constraint_set = false;
+ for (unsigned int j=0; j<dofs_per_face; ++j)
+ {
+ if (dofs_2[i] == dofs_1[j])
+ if (!(is_identity_constrained && target==i))
+ {
+ constraint_matrix.add_line(dofs_2[i]);
+ constraint_set = true;
+ }
+ }
+
+ if (!constraint_set)
+ {
+ // now treat constraints, either as an equality constraint or
+ // as a sequence of constraints
+ if (is_identity_constrained == true || is_inverse_constrained == true)
+ {
+ // Query the correct face_index on face_1 respecting the given
+ // orientation:
+ const unsigned int j
+ = cell_to_rotated_face_index[fe.face_to_cell_index(target,
+ 0, /* It doesn't really matter, just assume
+ * we're on the first face...
+ */
+ face_orientation, face_flip, face_rotation)];
+
+ if (constraint_matrix.is_constrained(dofs_2[i]))
{
- if (one_identity_found == false)
+ // if the two aren't already identity constrained (whichever way
+ // around) or already identical (in case of rotated periodicity constraints),
+ // then enter the constraint. otherwise there is nothing for us still to do
+ bool enter_constraint = false;
+ // see if this would add an identity constraint cycle
+ if (!constraint_matrix.is_constrained(dofs_1[j]))
{
- one_identity_found = true;
- identity_constraint_target = jj;
+ types::global_dof_index new_dof = dofs_2[i];
+ while (new_dof != dofs_1[j])
+ if (constraint_matrix.is_constrained(new_dof))
+ {
+ const std::vector<std::pair<types::global_dof_index, double > > *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<dofs_per_face; ++jj)
- if (((std::abs(transformation(i,jj)) < eps) ||
- (std::abs(transformation(i,jj)+1) < eps)) == false)
- {
- is_inverse_constrained = false;
- break;
- }
- if (is_inverse_constrained)
- {
- bool one_identity_found = false;
- for (unsigned int jj=0; jj<dofs_per_face; ++jj)
- if (std::abs(transformation(i,jj)+1) < eps)
+ else
{
- if (one_identity_found == false)
+ // if the two aren't already identity constrained (whichever way
+ // around) or already identical (in case of rotated periodicity constraints),
+ // then enter the constraint. otherwise there is nothing for us still to do
+ bool enter_constraint = false;
+ if (!constraint_matrix.is_constrained(dofs_1[j]))
{
- one_identity_found = true;
- inverse_constraint_target = jj;
+ if (dofs_2[i] != dofs_1[j])
+ enter_constraint = true;
}
- else
+ else //dofs_1[j] is constrained, is it identity or inverse constrained?
{
- is_inverse_constrained = false;
- inverse_constraint_target = numbers::invalid_unsigned_int;
- break;
+ const std::vector<std::pair<types::global_dof_index, double > > *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<std::pair<types::global_dof_index, double > > *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<dofs_per_face; ++j)
- {
- if (dofs_2[i] == dofs_1[j])
- if (!(is_identity_constrained && target==i))
+ 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 if (!constraint_matrix.is_constrained(dofs_2[i]))
+ {
+ // this is just a regular constraint. enter it piece by piece
+ constraint_matrix.add_line(dofs_2[i]);
+ for (unsigned int jj=0; jj<dofs_per_face; ++jj)
{
- constraint_matrix.add_line(dofs_2[i]);
- constrained_set = true;
+ // Query the correct face_index on face_1 respecting the given
+ // orientation:
+ const unsigned int j =
+ cell_to_rotated_face_index[fe.face_to_cell_index
+ (jj, 0, face_orientation, face_flip, face_rotation)];
+
+ // And finally constrain the two DoFs respecting component_mask:
+ if (transformation(i,jj) != 0)
+ constraint_matrix.add_entry(dofs_2[i], dofs_1[j],
+ transformation(i,jj));
}
- }
-
- if (!constrained_set)
- {
- // now treat constraints, either as an equality constraint or
- // as a sequence of constraints
- if (is_identity_constrained == true || is_inverse_constrained == true)
- {
- // Query the correct face_index on face_1 respecting the given
- // orientation:
- const unsigned int j
- = cell_to_rotated_face_index[fe.face_to_cell_index(target,
- 0, /* It doesn't really matter, just assume
- * we're on the first face...
- */
- face_orientation, face_flip, face_rotation)];
-
- // if the two aren't already identity constrained (whichever way
- // around) or already identical (in case of rotated periodicity constraints),
- // then enter the constraint. otherwise there is nothing for us still to do
- bool enter_constraint = false;
- if (!constraint_matrix.is_constrained(dofs_1[j]))
- {
- if (dofs_2[i] != dofs_1[j])
- enter_constraint = true;
- }
- else //dofs_1[j] is constrained, is it identity or inverse constrained?
- {
- const std::vector<std::pair<types::global_dof_index, double > > *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<dofs_per_face; ++jj)
- {
- // Query the correct face_index on face_1 respecting the given
- // orientation:
- const unsigned int j =
- cell_to_rotated_face_index[fe.face_to_cell_index
- (jj, 0, face_orientation, face_flip, face_rotation)];
-
- // And finally constrain the two DoFs respecting component_mask:
- if (transformation(i,jj) != 0)
- constraint_matrix.add_entry(dofs_2[i], dofs_1[j],
- transformation(i,jj));
- }
- }
- }
- }
+ }
+ }
+ }
}
}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 2D case.
+
+
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <fstream>
+#include <iostream>
+
+#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<GridTools::PeriodicFacePair<typename DoFHandler<dim>::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<DoFHandler<dim> >(periodicity_vectorDof, constraints);
+
+ constraints.close();
+ /* std::map<types::global_dof_index, Point<dim> > support_points;
+ DoFTools::map_dofs_to_support_points (MappingQ<dim,dim>(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 <int dim>
+class PeriodicReference : public Function<dim>
+{
+public:
+ PeriodicReference () : Function<dim>() {}
+ virtual double value (const Point<dim> &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 <int dim>
+void get_point_value
+(const DoFHandler<dim> &dof_handler,
+ const Point<dim> &point,
+ const Vector<double> &solution,
+ Vector<double> &value)
+{
+ VectorTools::point_value (dof_handler, solution,
+ point, value);
+}
+
+
+void check_periodicity(const DoFHandler<2> &dof_handler, Vector<double> &solution, const unsigned int cycle)
+{
+ unsigned int n_points = 2;
+ for (unsigned int i = 0; i<cycle; i++)
+ n_points*=2;
+
+ //don't test exactly at the support points, since point_value is not stable there
+ const double eps = 1./(16.*n_points);
+
+ for (unsigned int i=1; i< n_points; i++)
+ {
+ Vector<double> value1(1);
+ Vector<double> 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<<point1<< "\t" << value1[0] << std::endl;
+ std::cout<<point2<< "\t" << value2[0] << std::endl;
+ AssertThrow(false, ExcInternalError());
+ }
+ else
+ {
+ std::cout << point1 << "\t" << "pass" << std::endl;
+ }
+ }
+ for (unsigned int i=1; i< n_points; i++)
+ {
+ Vector<double> value1(1);
+ Vector<double> 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<<point1<< "\t" << value1[0] << std::endl;
+ std::cout<<point2<< "\t" << value2[0] << std::endl;
+ Assert(false, ExcInternalError());
+ }
+ else
+ {
+ std::cout << point1 << "\t" << "pass" << std::endl;
+ }
+ }
+}
+
+
+int main (int argc, char *argv[])
+{
+ initlog();
+
+ constexpr int dim = 2;
+ const double L=numbers::PI;
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_cube (triangulation, -L, L, true);
+
+ triangulation.refine_global(1);
+ typename Triangulation<dim>::active_cell_iterator cellBegin = triangulation.begin_active();
+ cellBegin->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
+
+ FE_Q<dim> fe(1);
+ DoFHandler<dim> dof_handler (triangulation);
+ dof_handler.distribute_dofs(fe);
+
+ std::vector<ConstraintMatrix> constraints(8);
+
+ PeriodicReference<dim> periodic_function;
+
+ std::vector<Vector<double> > projection(8, Vector<double> (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<dim>(3), periodic_function, projection[i]);
+ check_periodicity(dof_handler, projection[i], i);
+ }
+}
+
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature.h>
+#include <deal.II/base/function.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_out.h>
+#include <fstream>
+#include <iostream>
+
+#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<GridTools::PeriodicFacePair<typename DoFHandler<dim>::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<DoFHandler<dim> >(periodicity_vectorDof, constraints);
+
+ constraints.close();
+ /*std::map<types::global_dof_index, Point<dim> > support_points;
+ DoFTools::map_dofs_to_support_points (MappingQ<dim,dim>(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 <int dim>
+class PeriodicReference : public Function<dim>
+{
+public:
+ PeriodicReference () : Function<dim>() {}
+ virtual double value (const Point<dim> &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 <int dim>
+void get_point_value
+(const DoFHandler<dim> &dof_handler,
+ const Point<dim> &point,
+ const Vector<double> &solution,
+ Vector<double> &value)
+{
+ VectorTools::point_value (dof_handler, solution,
+ point, value);
+}
+
+
+void check_periodicity(const DoFHandler<3> &dof_handler, Vector<double> &solution, const unsigned int cycle)
+{
+ unsigned int n_points = 2;
+ for (unsigned int i = 0; i<cycle; ++i)
+ n_points*=2;
+
+ //don't test exactly at the support points, since point_value is not stable there
+ const double eps = 1./(16.*n_points);
+
+ for (unsigned int i=1; i<n_points; ++i)
+ for (unsigned int j=1; j<n_points; ++j)
+ {
+ Vector<double> value1(1);
+ Vector<double> 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<<point1<< "\t" << value1[0] << std::endl;
+ std::cout<<point2<< "\t" << value2[0] << std::endl;
+ Assert(false, ExcInternalError());
+ }
+ else
+ {
+ std::cout << point1 << "\t" << "pass" << std::endl;
+ }
+ }
+
+ for (unsigned int i=1; i<n_points; ++i)
+ for (unsigned int j=1; j<n_points; ++j)
+ {
+ Vector<double> value1(1);
+ Vector<double> 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<<point1<< "\t" << value1[0] << std::endl;
+ std::cout<<point2<< "\t" << value2[0] << std::endl;
+ Assert(false, ExcInternalError());
+ }
+ else
+ {
+ std::cout << point1 << "\t" << "pass" << std::endl;
+ }
+ }
+
+ for (unsigned int i=1; i<n_points; ++i)
+ for (unsigned int j=1; j<n_points; ++j)
+ {
+ Vector<double> value1(1);
+ Vector<double> 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<<point1<< "\t" << value1[0] << std::endl;
+ std::cout<<point2<< "\t" << value2[0] << std::endl;
+ Assert(false, ExcInternalError());
+ }
+ else
+ {
+ std::cout << point1 << "\t" << "pass" << std::endl;
+ }
+ }
+}
+
+
+int main (int argc, char *argv[])
+{
+ initlog();
+
+ constexpr int dim = 3;
+ const double L=numbers::PI;
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_cube (triangulation, -L, L, true);
+
+ triangulation.refine_global(1);
+ typename Triangulation<dim>::active_cell_iterator cellBegin = triangulation.begin_active();
+ cellBegin->set_refine_flag();
+ triangulation.execute_coarsening_and_refinement();
+
+ FE_Q<dim> fe(1);
+ DoFHandler<dim> dof_handler (triangulation);
+ dof_handler.distribute_dofs(fe);
+
+ std::vector<ConstraintMatrix> constraints(4);
+
+ PeriodicReference<dim> periodic_function;
+
+ std::vector<Vector<double> > projection(4, Vector<double> (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<dim>(3), periodic_function, projection[i]);
+ check_periodicity(dof_handler, projection[i], i);
+ }
+}
+
--- /dev/null
+
+DEAL::Testing version 0
+DEAL::Testing version 1
+DEAL::Testing version 2
+DEAL::Testing version 3
projects / dealii.git / commitdiff