--- /dev/null
+New: MGConstrainedDoFs object now stores information for periodicity
+constraints which is used during the multigrid transfer by
+MGTransferMatrixFree and MGTransferPrebuilt.
+<br>
+(Conrad Clevenger, Martin Kronbichler, 2017/07/08)
* edges of a refinement level are treated corrected but leaves degrees of
* freedom at the boundary of the domain untouched assuming that no
* Dirichlet boundary conditions for them exist.
+ *
+ * Furthermore, this call sets up a ConstraintMatrix on each level that
+ * contains possible periodicity constraints in case those have been added to
+ * the underlying triangulation. The constraint matrix can be queried by
+ * get_level_constraint_matrix(level). Note that the current implementation of
+ * periodicity constraints in this class does not support rotation matrices in the
+ * periodicity definition, i.e., the respective argument in the
+ * GridTools::collect_periodic_faces() may not be different from the identity
+ * matrix.
*/
template <int dim, int spacedim>
void initialize (const DoFHandler<dim,spacedim> &dof);
*/
bool have_boundary_indices () const;
+ /**
+ * Return the level constraint matrix for a given level, containing
+ * periodicity constraints (if enabled on the triangulation).
+ */
+ const ConstraintMatrix &
+ get_level_constraint_matrix (const unsigned int level) const;
+
private:
/**
* between the level and cells on a coarser level.
*/
std::vector<IndexSet> refinement_edge_indices;
+
+ /**
+ * Constraint matrices containing information regarding potential
+ * periodic boundary conditions for each level .
+ */
+ std::vector<ConstraintMatrix> level_constraints;
+
};
{
boundary_indices.clear();
refinement_edge_indices.clear();
+ level_constraints.clear();
const unsigned int nlevels = dof.get_triangulation().n_global_levels();
- //At this point refinement_edge_indices is empty.
+ //At this point level_constraint and refinement_edge_indices are empty.
+ level_constraints.resize(nlevels);
refinement_edge_indices.resize(nlevels);
for (unsigned int l=0; l<nlevels; ++l)
- refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
+ {
+ IndexSet relevant_dofs;
+ DoFTools::extract_locally_relevant_level_dofs(dof,l,relevant_dofs);
+ level_constraints[l].reinit(relevant_dofs);
+
+ // Loop through relevant cells and faces finding those which are periodic neighbors.
+ typename DoFHandler<dim,spacedim>::cell_iterator
+ cell = dof.begin(l),
+ endc = dof.end(l);
+ for (; cell!=endc; ++cell)
+ if (cell->level_subdomain_id() != numbers::artificial_subdomain_id)
+ {
+ for (unsigned int f = 0;
+ f < GeometryInfo<dim>::faces_per_cell;
+ ++f)
+ if (cell->has_periodic_neighbor(f))
+ {
+ 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 = cell->face(f)->get_fe(0).dofs_per_face;
+ 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 constraint matrix
+ // 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_line(dofs_2[i]);
+ level_constraints[l].add_entry(dofs_2[i], dofs_1[i], 1.);
+ }
+ }
+ }
+ level_constraints[l].close();
+
+ // Initialize with empty IndexSet of correct size
+ refinement_edge_indices[l] = IndexSet(dof.n_dofs(l));
+ }
MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices);
}
}
+
+
+inline
+const ConstraintMatrix &
+MGConstrainedDoFs::get_level_constraint_matrix (const unsigned int level) const
+{
+ AssertIndexRange(level, level_constraints.size());
+ return level_constraints[level];
+}
+
+
+
DEAL_II_NAMESPACE_CLOSE
#endif
fe.get_prolongation_matrix(c)(renumbering[j],renumbering[i]);
}
+ namespace
+ {
+ /**
+ * Helper function for setup_transfer. Checks for identity constrained dofs
+ * and replace with the indices of the dofs to which they are constrained
+ */
+ void replace (const MGConstrainedDoFs *mg_constrained_dofs,
+ const unsigned int level,
+ std::vector<types::global_dof_index> &dof_indices)
+ {
+ if (mg_constrained_dofs != nullptr &&
+ mg_constrained_dofs->get_level_constraint_matrix(level).n_constraints() > 0)
+ for (auto &ind : dof_indices)
+ if (mg_constrained_dofs->get_level_constraint_matrix(level).is_identity_constrained(ind))
+ {
+ Assert(mg_constrained_dofs->get_level_constraint_matrix(level).get_constraint_entries(ind)->size() == 1,
+ ExcInternalError());
+ ind = mg_constrained_dofs->get_level_constraint_matrix(level).get_constraint_entries(ind)->front().first;
+ }
+ }
+ }
// Sets up most of the internal data structures of the MGTransferMatrixFree
// class
continue;
cell->child(c)->get_mg_dof_indices(local_dof_indices);
+ replace(mg_constrained_dofs, level, local_dof_indices);
+
const IndexSet &owned_level_dofs = mg_dof.locally_owned_mg_dofs(level);
for (unsigned int i=0; i<local_dof_indices.size(); ++i)
if (!owned_level_dofs.is_element(local_dof_indices[i]))
{
cell->get_mg_dof_indices(local_dof_indices);
+ replace(mg_constrained_dofs, level-1, local_dof_indices);
+
const IndexSet &owned_level_dofs_l0 = mg_dof.locally_owned_mg_dofs(0);
for (unsigned int i=0; i<local_dof_indices.size(); ++i)
if (!owned_level_dofs_l0.is_element(local_dof_indices[i]))
}
+namespace
+{
+ /**
+ * Helper function for build_matrices. Checks for identity constrained dofs
+ * and replace with the indices of the dofs to which they are constrained
+ */
+ void replace (const MGConstrainedDoFs *mg_constrained_dofs,
+ const unsigned int level,
+ std::vector<types::global_dof_index> &dof_indices)
+ {
+ if (mg_constrained_dofs != nullptr &&
+ mg_constrained_dofs->get_level_constraint_matrix(level).n_constraints() > 0)
+ for (auto &ind : dof_indices)
+ if (mg_constrained_dofs->get_level_constraint_matrix(level).is_identity_constrained(ind))
+ {
+ Assert(mg_constrained_dofs->get_level_constraint_matrix(level).get_constraint_entries(ind)->size() == 1,
+ ExcInternalError());
+ ind = mg_constrained_dofs->get_level_constraint_matrix(level).get_constraint_entries(ind)->front().first;
+ }
+ }
+}
template <typename VectorType>
template <int dim, int spacedim>
{
cell->get_mg_dof_indices (dof_indices_parent);
+ replace(this->mg_constrained_dofs, level, dof_indices_parent);
+
Assert(cell->n_children()==GeometryInfo<dim>::max_children_per_cell,
ExcNotImplemented());
for (unsigned int child=0; child<cell->n_children(); ++child)
cell->child(child)->get_mg_dof_indices (dof_indices_child);
+ replace(this->mg_constrained_dofs, level+1, dof_indices_child);
+
// now tag the entries in the
// matrix which will be used
// for this pair of parent/child
{
cell->get_mg_dof_indices (dof_indices_parent);
+ replace(this->mg_constrained_dofs, level, dof_indices_parent);
+
Assert(cell->n_children()==GeometryInfo<dim>::max_children_per_cell,
ExcNotImplemented());
for (unsigned int child=0; child<cell->n_children(); ++child)
cell->child(child)->get_mg_dof_indices (dof_indices_child);
+ replace(this->mg_constrained_dofs, level+1, dof_indices_child);
+
// now set the entries in the matrix
for (unsigned int i=0; i<dofs_per_cell; ++i)
prolongation_matrices[level]->set (dof_indices_child[i],
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Check MGTransferMatrixFree and MGTransferPrebuilt for periodic boundary
+// conditions
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+
+
+template <int dim, typename Number>
+void check(const unsigned int fe_degree)
+{
+ deallog.threshold_double(std::max(5e2*(double)std::numeric_limits<Number>::epsilon(),
+ 1e-11));
+ FE_Q<dim> fe(fe_degree);
+ deallog << "FE: " << fe.get_name() << std::endl;
+
+ parallel::distributed::Triangulation<dim>
+ tr(MPI_COMM_WORLD,
+ Triangulation<dim>::limit_level_difference_at_vertices,
+ parallel::distributed::Triangulation<dim>::construct_multigrid_hierarchy);
+ GridGenerator::hyper_cube(tr, -1, 1, true);
+ std::vector<GridTools::PeriodicFacePair<typename Triangulation<dim>::cell_iterator> >
+ periodicity_vector;
+ GridTools::collect_periodic_faces(tr, 0, 1, 0,
+ periodicity_vector);
+ GridTools::collect_periodic_faces(tr, 2, 3, 1,
+ periodicity_vector);
+ tr.add_periodicity(periodicity_vector);
+ tr.refine_global(2);
+ deallog << "no. cells: " << tr.n_global_active_cells() << std::endl;
+
+ DoFHandler<dim> mgdof(tr);
+ mgdof.distribute_dofs(fe);
+ mgdof.distribute_mg_dofs(fe);
+
+ MGConstrainedDoFs mg_constrained_dofs;
+ mg_constrained_dofs.initialize(mgdof);
+
+ // build reference
+ MGTransferPrebuilt<LinearAlgebra::distributed::Vector<double> >
+ transfer_ref(mg_constrained_dofs);
+ transfer_ref.build_matrices(mgdof);
+ deallog << "Transfer matrices: " << std::endl;
+ transfer_ref.print_matrices(deallog.get_file_stream());
+ deallog << std::endl;
+
+ // build matrix-free transfer
+ MGTransferMatrixFree<dim, Number> transfer(mg_constrained_dofs);
+ transfer.build(mgdof);
+
+ // check prolongation for all levels using random vector
+ for (unsigned int level=1; level<mgdof.get_triangulation().n_global_levels(); ++level)
+ {
+ LinearAlgebra::distributed::Vector<Number> v1, v2;
+ LinearAlgebra::distributed::Vector<double> v1_cpy, v2_cpy, v3;
+ v1.reinit(mgdof.locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+ v2.reinit(mgdof.locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+ v3.reinit(mgdof.locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+ for (unsigned int i=0; i<v1.local_size(); ++i)
+ v1.local_element(i) = (double)Testing::rand()/RAND_MAX;
+ v1_cpy = v1;
+ transfer.prolongate(level, v2, v1);
+ transfer_ref.prolongate(level, v3, v1_cpy);
+ v2_cpy = v2;
+ v3 -= v2_cpy;
+ deallog << "Diff prolongate l" << level << ": " << v3.l2_norm() << std::endl;
+ if (v3.l2_norm() > 1e-12)
+ {
+ // On level 0, we expect the matrix-based constraints to be wrong
+ // because it cannot capture the periodicity connections with a
+ // single cell
+ v2_cpy.print(deallog.get_file_stream());
+ v3.print(deallog.get_file_stream());
+ }
+ }
+
+ // check restriction for all levels using random vector
+ for (unsigned int level=1; level<mgdof.get_triangulation().n_global_levels(); ++level)
+ {
+ LinearAlgebra::distributed::Vector<Number> v1, v2;
+ LinearAlgebra::distributed::Vector<double> v1_cpy, v2_cpy, v3;
+ v1.reinit(mgdof.locally_owned_mg_dofs(level), MPI_COMM_WORLD);
+ v2.reinit(mgdof.locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+ v3.reinit(mgdof.locally_owned_mg_dofs(level-1), MPI_COMM_WORLD);
+ for (unsigned int i=0; i<v1.local_size(); ++i)
+ v1.local_element(i) = (double)Testing::rand()/RAND_MAX;
+ v1_cpy = v1;
+ transfer.restrict_and_add(level, v2, v1);
+ transfer_ref.restrict_and_add(level, v3, v1_cpy);
+ v2_cpy = v2;
+ v3 -= v2_cpy;
+ deallog << "Diff restrict l" << level << ": " << v3.l2_norm() << std::endl;
+
+ v2 = 1.;
+ v3 = 1.;
+ transfer.restrict_and_add(level, v2, v1);
+ transfer_ref.restrict_and_add(level, v3, v1_cpy);
+ v2_cpy = v2;
+ v3 -= v2_cpy;
+ deallog << "Diff restrict add l" << level << ": " << v3.l2_norm() << std::endl;
+ }
+ deallog << std::endl;
+}
+
+
+int main(int argc, char **argv)
+{
+ // no threading in this test...
+ Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+ mpi_initlog();
+
+ check<2,double>(1);
+}
--- /dev/null
+
+DEAL::FE: FE_Q<2>(1)
+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
+
+Level 1
+(3,3) 0.250000
+(3,5) 0.250000
+(3,7) 0.250000
+(3,8) 0.250000
+(5,3) 0.500000
+(5,7) 0.500000
+(7,3) 0.500000
+(7,5) 0.500000
+(8,3) 1.00000
+(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,3) 0.500000
+(13,5) 0.500000
+(14,5) 1.00000
+(16,3) 0.250000
+(16,5) 0.250000
+(16,7) 0.250000
+(16,8) 0.250000
+(17,3) 0.500000
+(17,7) 0.500000
+(19,7) 0.500000
+(19,8) 0.500000
+(20,7) 1.00000
+(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
+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
+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
+DEAL::
projects / dealii.git / commitdiff