--- /dev/null
+Fixed: There was an inconsistency between the way
+DoFTools::make_periodicity_constraints() and MGConstrainedDoFs
+set up the constraints for periodic constraints. This is fixed
+now.
+<br>
+(Peter Munch, Laura Prieto Saavedra, 2023/10/19)
const ComponentMask &component_mask = {},
const number periodicity_factor = 1.);
+ namespace internal
+ {
+ /**
+ * Internally used in make_periodicity_constraints.
+ *
+ * Enter constraints for periodicity into the given AffineConstraints
+ * object. this function is called when at least one of the two face
+ * iterators corresponds to an active object without further children
+ *
+ * The matrix @p transformation maps degrees of freedom from one face
+ * to another. If the DoFs on the two faces are supposed to match exactly,
+ * then the matrix so provided will be the identity matrix. if face 2 is
+ * once refined from face 1, then the matrix needs to be the interpolation
+ * matrix from a face to this particular child
+ *
+ * @note: @p face_1 is supposed to be active.
+ *
+ * @note We have to be careful not to accidentally create constraint
+ * cycles when adding periodic constraints: For example, as the
+ * corresponding testcase bits/periodicity_05 demonstrates, we can
+ * occasionally get into trouble if we already have the constraint
+ * x1 == x2 and want to insert x2 == x1. We avoid this by skipping
+ * such "identity constraints" if the opposite constraint already
+ * exists.
+ */
+ template <typename FaceIterator, typename number>
+ void
+ set_periodicity_constraints(
+ const FaceIterator &face_1,
+ const std_cxx20::type_identity_t<FaceIterator> &face_2,
+ const FullMatrix<double> &transformation,
+ AffineConstraints<number> &affine_constraints,
+ const ComponentMask &component_mask,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation,
+ const number periodicity_factor,
+ const unsigned int level = numbers::invalid_unsigned_int);
+ } // namespace internal
+
/**
* @}
*/
#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/subscriptor.h>
+#include <deal.II/dofs/dof_tools.h>
+
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/multigrid/mg_tools.h>
level_constraints[l].reinit(dof.locally_owned_mg_dofs(l),
relevant_dofs);
}
+ }
+
+ // TODO: currently we only consider very basic periodic constraints
+ const IdentityMatrix transformation(dof.get_fe().n_dofs_per_face());
+ const ComponentMask component_mask;
+ const double periodicity_factor = 1.0;
+
+ for (const auto &[first_cell, second_cell] :
+ dof.get_triangulation().get_periodic_face_map())
+ {
+ // only consider non-artificial cells
+ if (first_cell.first->is_artificial_on_level())
+ continue;
+ if (second_cell.first.first->is_artificial_on_level())
+ continue;
+
+ // consider cell pairs with the same level
+ if (first_cell.first->level() != second_cell.first.first->level())
+ continue;
+
+ DoFTools::internal::set_periodicity_constraints(
+ first_cell.first->as_dof_handler_level_iterator(dof)->face(
+ first_cell.second),
+ second_cell.first.first->as_dof_handler_level_iterator(dof)->face(
+ second_cell.first.second),
+ transformation,
+ level_constraints[first_cell.first->level()],
+ component_mask,
+ second_cell.second[0],
+ second_cell.second[1],
+ second_cell.second[2],
+ periodicity_factor,
+ first_cell.first->level());
+ }
- // Loop through relevant cells and faces finding those which are periodic
- // neighbors.
- for (const auto &cell : dof.cell_iterators_on_level(l))
- if (cell->is_artificial_on_level() == false)
- for (auto f : cell->face_indices())
- if (cell->has_periodic_neighbor(f) &&
- cell->periodic_neighbor(f)->level() == cell->level())
- {
- if (cell->is_locally_owned_on_level())
- {
- Assert(
- cell->periodic_neighbor(f)->level_subdomain_id() !=
- numbers::artificial_subdomain_id,
- ExcMessage(
- "Periodic neighbor of a locally owned cell must either be owned or ghost."));
- }
- // Cell is a level-ghost and its neighbor is a
- // level-artificial cell nothing to do here
- else if (cell->periodic_neighbor(f)->level_subdomain_id() ==
- numbers::artificial_subdomain_id)
- {
- Assert(cell->is_locally_owned_on_level() == false,
- ExcInternalError());
- continue;
- }
-
- const unsigned int dofs_per_face =
- dof.get_fe(0).n_dofs_per_face(f);
- std::vector<types::global_dof_index> dofs_1(dofs_per_face);
- std::vector<types::global_dof_index> dofs_2(dofs_per_face);
-
- cell->periodic_neighbor(f)
- ->face(cell->periodic_neighbor_face_no(f))
- ->get_mg_dof_indices(l, dofs_1, 0);
- cell->face(f)->get_mg_dof_indices(l, dofs_2, 0);
- // Store periodicity information in the level
- // AffineConstraints object. Skip DoFs for which we've
- // previously entered periodicity constraints already; this
- // can happen, for example, for a vertex dof at a periodic
- // boundary that we visit from more than one cell
- for (unsigned int i = 0; i < dofs_per_face; ++i)
- if (level_constraints[l].can_store_line(dofs_2[i]) &&
- level_constraints[l].can_store_line(dofs_1[i]) &&
- !level_constraints[l].is_constrained(dofs_2[i]) &&
- !level_constraints[l].is_constrained(dofs_1[i]))
- {
- level_constraints[l].add_constraint(dofs_2[i],
- {{dofs_1[i], 1.}},
- 0.);
- }
- }
+ for (unsigned int l = min_level; l <= max_level; ++l)
+ {
level_constraints[l].close();
// Initialize with empty IndexSet of correct size
- namespace
+ namespace internal
{
- /**
- * @internal
- *
- * Internally used in make_periodicity_constraints.
- *
- * enter constraints for periodicity into the given AffineConstraints
- * object. this function is called when at least one of the two face
- * iterators corresponds to an active object without further children
- *
- * @param transformation A matrix that maps degrees of freedom from one face
- * to another. If the DoFs on the two faces are supposed to match exactly,
- * then the matrix so provided will be the identity matrix. if face 2 is
- * once refined from face 1, then the matrix needs to be the interpolation
- * matrix from a face to this particular child
- *
- * @precondition: face_1 is supposed to be active
- *
- * @note We have to be careful not to accidentally create constraint
- * cycles when adding periodic constraints: For example, as the
- * corresponding testcase bits/periodicity_05 demonstrates, we can
- * occasionally get into trouble if we already have the constraint
- * x1 == x2 and want to insert x2 == x1. We avoid this by skipping
- * such "identity constraints" if the opposite constraint already
- * exists.
- */
template <typename FaceIterator, typename number>
void
set_periodicity_constraints(
const bool face_orientation,
const bool face_flip,
const bool face_rotation,
- const number periodicity_factor)
+ const number periodicity_factor,
+ const unsigned int level)
{
static const int dim = FaceIterator::AccessorType::dimension;
static const int spacedim = FaceIterator::AccessorType::space_dimension;
+ const bool use_mg = level != numbers::invalid_unsigned_int;
+
// we should be in the case where face_1 is active, i.e. has no children:
- Assert(!face_1->has_children(), ExcInternalError());
+ Assert(use_mg || (!face_1->has_children()), ExcInternalError());
Assert(face_1->n_active_fe_indices() == 1, ExcInternalError());
// If face_2 does have children, then we need to iterate over these
// children and set periodic constraints in the inverse direction:
- if (face_2->has_children())
+ if ((!use_mg) && face_2->has_children())
{
Assert(face_2->n_children() ==
GeometryInfo<dim>::max_children_per_face,
std::vector<types::global_dof_index> dofs_1(dofs_per_face);
std::vector<types::global_dof_index> dofs_2(dofs_per_face);
- face_1->get_dof_indices(dofs_1, face_1_index);
- face_2->get_dof_indices(dofs_2, face_2_index);
+ if (use_mg)
+ face_1->get_mg_dof_indices(level, dofs_1, face_1_index);
+ else
+ face_1->get_dof_indices(dofs_1, face_1_index);
+
+ if (use_mg)
+ face_2->get_mg_dof_indices(level, dofs_2, face_2_index);
+ else
+ face_2->get_dof_indices(dofs_2, face_2_index);
// If either of the two faces has an invalid dof index, stop. This is
// so that there is no attempt to match artificial cells of parallel
}
} /* for dofs_per_face */
}
+ } // namespace internal
-
+ namespace
+ {
// Internally used in make_periodicity_constraints.
//
// Build up a (possibly rotated) interpolation matrix that is used in
// matrix is the identity matrix.
if (first_vector_components.empty() && matrix.m() == 0)
{
- set_periodicity_constraints(face_2,
- face_1,
- transformation,
- affine_constraints,
- component_mask,
- face_orientation,
- face_flip,
- face_rotation,
- periodicity_factor);
+ internal::set_periodicity_constraints(face_2,
+ face_1,
+ transformation,
+ affine_constraints,
+ component_mask,
+ face_orientation,
+ face_flip,
+ face_rotation,
+ periodicity_factor);
}
else
{
FullMatrix<double> inverse(transformation.m());
inverse.invert(transformation);
- set_periodicity_constraints(face_2,
- face_1,
- inverse,
- affine_constraints,
- component_mask,
- face_orientation,
- face_flip,
- face_rotation,
- periodicity_factor);
+ internal::set_periodicity_constraints(face_2,
+ face_1,
+ inverse,
+ affine_constraints,
+ component_mask,
+ face_orientation,
+ face_flip,
+ face_rotation,
+ periodicity_factor);
}
}
else
// the rotation when constraining face_2 to face_1. Therefore
// face_flip has to be toggled if face_rotation is true: In case of
// inverted orientation, nothing has to be done.
- set_periodicity_constraints(face_1,
- face_2,
- transformation,
- affine_constraints,
- component_mask,
- face_orientation,
- face_orientation ?
- face_rotation ^ face_flip :
- face_flip,
- face_rotation,
- periodicity_factor);
+ internal::set_periodicity_constraints(face_1,
+ face_2,
+ transformation,
+ affine_constraints,
+ component_mask,
+ face_orientation,
+ face_orientation ?
+ face_rotation ^ face_flip :
+ face_flip,
+ face_rotation,
+ periodicity_factor);
}
}
}
const FullMatrix<double> &,
const std::vector<unsigned int> &,
const S);
+
+ template void DoFTools::internal::set_periodicity_constraints(
+ const DoFHandler<deal_II_dimension,
+ deal_II_space_dimension>::face_iterator &face_1,
+ const DoFHandler<deal_II_dimension,
+ deal_II_space_dimension>::face_iterator &face_2,
+ const FullMatrix<double> &transformation,
+ AffineConstraints<S> &affine_constraints,
+ const ComponentMask &component_mask,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation,
+ const S periodicity_factor,
+ const unsigned int level);
+
+ template void DoFTools::internal::set_periodicity_constraints(
+ const DoFHandler<deal_II_dimension,
+ deal_II_space_dimension>::level_face_iterator &face_1,
+ const DoFHandler<deal_II_dimension,
+ deal_II_space_dimension>::level_face_iterator &face_2,
+ const FullMatrix<double> &transformation,
+ AffineConstraints<S> &affine_constraints,
+ const ComponentMask &component_mask,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation,
+ const S periodicity_factor,
+ const unsigned int level);
#endif
}
DEAL::no. cells: 16
DEAL::Transfer matrices:
Level 0
-(3,3) 0.250000
-(5,3) 0.500000
-(7,3) 0.500000
-(8,3) 1.00000
+(0,0) 1.00000
+(1,0) 0.500000
+(2,0) 0.500000
+(3,0) 0.250000
Level 1
+(0,0) 1.00000
+(1,0) 0.500000
+(1,1) 0.500000
+(2,0) 0.500000
+(2,2) 0.500000
+(3,0) 0.250000
+(3,1) 0.250000
+(3,2) 0.250000
(3,3) 0.250000
-(3,5) 0.250000
-(3,7) 0.250000
-(3,8) 0.250000
+(4,1) 1.00000
+(5,1) 0.500000
(5,3) 0.500000
-(5,7) 0.500000
+(6,2) 1.00000
+(7,2) 0.500000
(7,3) 0.500000
-(7,5) 0.500000
(8,3) 1.00000
+(9,0) 0.500000
+(9,1) 0.500000
+(10,0) 0.250000
+(10,1) 0.250000
+(10,2) 0.250000
(10,3) 0.250000
-(10,5) 0.250000
-(10,7) 0.250000
-(10,8) 0.250000
-(12,5) 0.500000
-(12,8) 0.500000
+(13,2) 0.500000
(13,3) 0.500000
-(13,5) 0.500000
-(14,5) 1.00000
+(15,0) 0.500000
+(15,2) 0.500000
+(16,0) 0.250000
+(16,1) 0.250000
+(16,2) 0.250000
(16,3) 0.250000
-(16,5) 0.250000
-(16,7) 0.250000
-(16,8) 0.250000
+(17,1) 0.500000
(17,3) 0.500000
-(17,7) 0.500000
-(19,7) 0.500000
-(19,8) 0.500000
-(20,7) 1.00000
+(21,0) 0.250000
+(21,1) 0.250000
+(21,2) 0.250000
(21,3) 0.250000
-(21,5) 0.250000
-(21,7) 0.250000
-(21,8) 0.250000
-(22,5) 0.500000
-(22,8) 0.500000
-(23,7) 0.500000
-(23,8) 0.500000
-(24,8) 1.00000
DEAL::
-DEAL::Diff prolongate l1: 0.823013
+DEAL::Diff prolongate l1: 0.866046
Process #0
Local range: [0, 9), global size: 9
Vector data:
-0.000e+00 0.000e+00 0.000e+00 7.984e-01 0.000e+00 7.984e-01 0.000e+00 7.984e-01 7.984e-01
+8.402e-01 8.402e-01 8.402e-01 8.402e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
Process #0
Local range: [0, 9), global size: 9
Vector data:
-0.000e+00 0.000e+00 0.000e+00 -5.988e-01 0.000e+00 -3.992e-01 0.000e+00 -3.992e-01 0.000e+00
-DEAL::Diff prolongate l2: 0
-DEAL::Diff restrict l1: 0.555735
-DEAL::Diff restrict add l1: 0.555735
-DEAL::Diff restrict l2: 0
-DEAL::Diff restrict add l2: 0
+0.000e+00 -4.201e-01 -4.201e-01 -6.301e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff restrict l1: 1.41100
+DEAL::Diff restrict add l1: 1.41100
+DEAL::Diff restrict l2: 5.43896e-16
+DEAL::Diff restrict add l2: 4.44089e-16
DEAL::
DEAL:: 0, 2 <-> 1, 3 via periodic boundary
DEAL:: 1, 3 <-> 0, 2 via periodic boundary
DEAL:: Constraint matrix:
- 0 1: 1.00000
- 2 3: 1.00000
+ 1 0: 1.00000
+ 3 2: 1.00000
DEAL::Level 1:
DEAL:: Faces neighboring coarser cells:
DEAL:: Refinement edge DoFs: {}
DEAL:: 0, 2 <-> 4, 5 via periodic boundary
DEAL:: 2, 6 <-> 5, 8 via periodic boundary
DEAL:: Constraint matrix:
- 0 4: 1.00000
- 2 5: 1.00000
- 6 8: 1.00000
+ 4 0: 1.00000
+ 5 2: 1.00000
+ 8 6: 1.00000
DEAL::Level 2:
DEAL:: Faces neighboring coarser cells:
DEAL:: 0, 2 via periodic boundary
DEAL:: 0, 2, 4, 6 <-> 1, 3, 5, 7 via periodic boundary
DEAL:: 1, 3, 5, 7 <-> 0, 2, 4, 6 via periodic boundary
DEAL:: Constraint matrix:
- 0 1: 1.00000
- 2 3: 1.00000
- 4 5: 1.00000
- 6 7: 1.00000
+ 1 0: 1.00000
+ 3 2: 1.00000
+ 5 4: 1.00000
+ 7 6: 1.00000
DEAL::Level 1:
DEAL:: Faces neighboring coarser cells:
DEAL:: Refinement edge DoFs: {}
DEAL:: 4, 6, 18, 20 <-> 10, 11, 22, 23 via periodic boundary
DEAL:: 6, 14, 20, 24 <-> 11, 17, 23, 26 via periodic boundary
DEAL:: Constraint matrix:
- 0 8: 1.00000
- 2 9: 1.00000
- 4 10: 1.00000
- 6 11: 1.00000
- 12 16: 1.00000
- 14 17: 1.00000
- 18 22: 1.00000
- 20 23: 1.00000
- 24 26: 1.00000
+ 8 0: 1.00000
+ 9 2: 1.00000
+ 10 4: 1.00000
+ 11 6: 1.00000
+ 16 12: 1.00000
+ 17 14: 1.00000
+ 22 18: 1.00000
+ 23 20: 1.00000
+ 26 24: 1.00000
DEAL::Level 2:
DEAL:: Faces neighboring coarser cells:
DEAL:: 0, 2, 4, 6 via periodic boundary
--- /dev/null
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2023 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.md at
+ * the top level directory of deal.II.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * Test that DoFTools::make_periodicity_constraints() and MGConstrainedDoFs
+ * creates the same constrains for periodic constraints in the case of
+ * globally refined meshes.
+ */
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+int
+main()
+{
+ initlog();
+
+ const unsigned int dim = 2;
+ const unsigned int n_refinements = 2;
+
+ Triangulation<dim> tria(
+ Triangulation<dim>::MeshSmoothing::limit_level_difference_at_vertices);
+ GridGenerator::hyper_cube(tria, 0, 1, true);
+
+ std::vector<
+ GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator>>
+ periodic_faces;
+
+ GridTools::collect_periodic_faces(tria, 0, 1, 0, periodic_faces);
+
+ tria.add_periodicity(periodic_faces);
+ tria.refine_global(n_refinements);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(FE_Q<dim>(1));
+ dof_handler.distribute_mg_dofs();
+
+ MGConstrainedDoFs mg_constraints;
+ mg_constraints.initialize(dof_handler);
+
+ const auto &af_level = mg_constraints.get_level_constraints(n_refinements);
+
+ AffineConstraints<double> af;
+ DoFTools::make_periodicity_constraints<dim, dim, double>(
+ dof_handler, 0, 1, 0, af);
+
+ for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
+ {
+ Assert(af_level.is_constrained(i) == af.is_constrained(i),
+ ExcInternalError());
+
+ if (af_level.is_constrained(i) == false)
+ continue;
+
+ Assert(*af_level.get_constraint_entries(i) ==
+ *af.get_constraint_entries(i),
+ ExcInternalError());
+ }
+
+#if 0
+ af_level.print(std::cout);
+ std::cout << std::endl;
+ af.print(std::cout);
+ std::cout << std::endl;
+#endif
+
+ deallog << "OK!" << std::endl;
+}
--- /dev/null
+
+DEAL::OK!
DEAL::no. cells: 16
DEAL::Transfer matrices:
Level 0
-(3,3) 0.250000
-(5,3) 0.500000
-(7,3) 0.500000
-(8,3) 1.00000
+(0,0) 1.00000
+(1,0) 0.500000
+(2,0) 0.500000
+(3,0) 0.250000
Level 1
+(0,0) 1.00000
+(1,0) 0.500000
+(1,1) 0.500000
+(2,0) 0.500000
+(2,2) 0.500000
+(3,0) 0.250000
+(3,1) 0.250000
+(3,2) 0.250000
(3,3) 0.250000
-(3,5) 0.250000
-(3,7) 0.250000
-(3,8) 0.250000
+(4,1) 1.00000
+(5,1) 0.500000
(5,3) 0.500000
-(5,7) 0.500000
+(6,2) 1.00000
+(7,2) 0.500000
(7,3) 0.500000
-(7,5) 0.500000
(8,3) 1.00000
+(9,0) 0.500000
+(9,1) 0.500000
+(10,0) 0.250000
+(10,1) 0.250000
+(10,2) 0.250000
(10,3) 0.250000
-(10,5) 0.250000
-(10,7) 0.250000
-(10,8) 0.250000
-(12,5) 0.500000
-(12,8) 0.500000
+(13,2) 0.500000
(13,3) 0.500000
-(13,5) 0.500000
-(14,5) 1.00000
+(15,0) 0.500000
+(15,2) 0.500000
+(16,0) 0.250000
+(16,1) 0.250000
+(16,2) 0.250000
(16,3) 0.250000
-(16,5) 0.250000
-(16,7) 0.250000
-(16,8) 0.250000
+(17,1) 0.500000
(17,3) 0.500000
-(17,7) 0.500000
-(19,7) 0.500000
-(19,8) 0.500000
-(20,7) 1.00000
+(21,0) 0.250000
+(21,1) 0.250000
+(21,2) 0.250000
(21,3) 0.250000
-(21,5) 0.250000
-(21,7) 0.250000
-(21,8) 0.250000
-(22,5) 0.500000
-(22,8) 0.500000
-(23,7) 0.500000
-(23,8) 0.500000
-(24,8) 1.00000
DEAL::
-DEAL::Diff prolongate l1: 0.823013
+DEAL::Diff prolongate l1: 0.866046
Process #0
Local range: [0, 9), global size: 9
Vector data:
-0.000e+00 0.000e+00 0.000e+00 7.984e-01 0.000e+00 7.984e-01 0.000e+00 7.984e-01 7.984e-01
+8.402e-01 8.402e-01 8.402e-01 8.402e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
Process #0
Local range: [0, 9), global size: 9
Vector data:
-0.000e+00 0.000e+00 0.000e+00 -5.988e-01 0.000e+00 -3.992e-01 0.000e+00 -3.992e-01 0.000e+00
-DEAL::Diff prolongate l2: 0
-DEAL::Diff restrict l1: 0.555735
-DEAL::Diff restrict add l1: 0.555735
-DEAL::Diff restrict l2: 0
-DEAL::Diff restrict add l2: 0
+0.000e+00 -4.201e-01 -4.201e-01 -6.301e-01 0.000e+00 0.000e+00 0.000e+00 0.000e+00 0.000e+00
+DEAL::Diff prolongate l2: 0.00000
+DEAL::Diff restrict l1: 1.41100
+DEAL::Diff restrict add l1: 1.41100
+DEAL::Diff restrict l2: 6.66134e-16
+DEAL::Diff restrict add l2: 6.28037e-16
DEAL::
projects / dealii.git / commitdiff