--- /dev/null
+New: MGConstrainedDoFs now allows the user to define constraints
+on each multigrid level to be applied during prolongation in a
+matrix-free transfer.
+<br>
+(Conrad Clevenger, 2019/09/13)
const std::set<types::boundary_id> &boundary_ids,
const ComponentMask & component_mask = ComponentMask());
+ /**
+ * Add user defined constraints to be used on level @p level.
+ *
+ * The user can call this function multiple times and any new,
+ * conflicting constraints will overwrite the previous constraints
+ * for that DoF.
+ *
+ * Before the transfer, the user defined constraints will be distributed
+ * to the source vector, and then any DoF index set using
+ * make_zero_boundary_constraints() will be overwritten with
+ * value zero.
+ *
+ * @note This is currently only implemented for MGTransferMatrixFree.
+ */
+ void
+ add_user_constraints(const unsigned int level,
+ const AffineConstraints<double> &constraints_on_level);
+
/**
* Fill the internal data structures with information
* about no normal flux boundary dofs.
const AffineConstraints<double> &
get_level_constraint_matrix(const unsigned int level) const;
+ /**
+ * Return the user defined constraint matrix for a given level. These
+ * constraints are set using the function add_user_constraints() and
+ * should not contain constraints for DoF indices set in
+ * make_zero_boundary_constraints() as they will be overwritten during
+ * the transfer.
+ */
+ const AffineConstraints<double> &
+ get_user_constraint_matrix(const unsigned int level) const;
+
private:
/**
* The indices of boundary dofs for each level.
* periodic boundary conditions for each level .
*/
std::vector<AffineConstraints<double>> level_constraints;
+
+ /**
+ * Constraint matrices defined by user.
+ */
+ std::vector<AffineConstraints<double>> user_constraints;
};
boundary_indices.clear();
refinement_edge_indices.clear();
level_constraints.clear();
+ user_constraints.clear();
const unsigned int nlevels = dof.get_triangulation().n_global_levels();
// At this point level_constraint and refinement_edge_indices are empty.
- level_constraints.resize(nlevels);
refinement_edge_indices.resize(nlevels);
+ level_constraints.resize(nlevels);
+ user_constraints.resize(nlevels);
for (unsigned int l = 0; l < nlevels; ++l)
{
IndexSet relevant_dofs;
}
+inline void
+MGConstrainedDoFs::add_user_constraints(
+ const unsigned int level,
+ const AffineConstraints<double> &constraints_on_level)
+{
+ AssertIndexRange(level, user_constraints.size());
+
+ // Get the relevant DoFs from level_constraints if
+ // the user constraint matrix has not been initialized
+ if (user_constraints[level].get_local_lines().size() == 0)
+ user_constraints[level].reinit(level_constraints[level].get_local_lines());
+
+ user_constraints[level].merge(
+ constraints_on_level,
+ AffineConstraints<double>::MergeConflictBehavior::right_object_wins);
+ user_constraints[level].close();
+}
+
+
inline void
MGConstrainedDoFs::clear()
{
boundary_indices.clear();
refinement_edge_indices.clear();
+ user_constraints.clear();
}
+inline const AffineConstraints<double> &
+MGConstrainedDoFs::get_user_constraint_matrix(const unsigned int level) const
+{
+ AssertIndexRange(level, user_constraints.size());
+ return user_constraints[level];
+}
+
+
+
DEAL_II_NAMESPACE_CLOSE
#endif
T<float> &) const;
template void dealii::AffineConstraints<double>::distribute<T<float>>(
T<float> &) const;
+
+ template void dealii::AffineConstraints<double>::distribute<
+ LinearAlgebra::distributed::T<float>>(
+ LinearAlgebra::distributed::T<float> &) const;
}
Utilities::fixed_power<dim>(n_child_dofs_1d);
constexpr unsigned int three_to_dim = Utilities::pow(3, dim);
+ // If we have user defined MG constraints, we must create
+ // a non-const, ghosted version of the source vector to distribute
+ // constraints.
+ const LinearAlgebra::distributed::Vector<Number> *to_use = &src;
+ LinearAlgebra::distributed::Vector<Number> new_src;
+ if (this->mg_constrained_dofs != nullptr &&
+ this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+ .get_local_lines()
+ .size() > 0)
+ {
+ LinearAlgebra::distributed::Vector<Number> copy_src(src);
+
+ // Distribute any user defined constraints
+ this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+ .distribute(copy_src);
+
+ // Re-initialize new ghosted vector with correct constraints
+ new_src.reinit(copy_src);
+ new_src = copy_src;
+ new_src.update_ghost_values();
+
+ to_use = &new_src;
+ }
+
for (unsigned int cell = 0; cell < n_owned_level_cells[to_level - 1];
cell += vec_size)
{
for (unsigned int k = 0; k < (dim > 2 ? degree_size : 1); ++k)
for (unsigned int j = 0; j < (dim > 1 ? degree_size : 1); ++j)
for (unsigned int i = 0; i < degree_size; ++i, ++m)
- evaluation_data[m][v] = src.local_element(
+ evaluation_data[m][v] = to_use->local_element(
indices[c * n_scalar_cell_dofs +
k * n_child_dofs_1d * n_child_dofs_1d +
j * n_child_dofs_1d + i]);
Assert((to_level >= 1) && (to_level <= prolongation_matrices.size()),
ExcIndexRange(to_level, 1, prolongation_matrices.size() + 1));
+#ifdef DEBUG
+ if (this->mg_constrained_dofs != nullptr)
+ Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+ .get_local_lines()
+ .size() == 0,
+ ExcNotImplemented());
+#endif
+
prolongation_matrices[to_level - 1]->vmult(dst, src);
}
Assert(dst.n_blocks() == this->n_mg_blocks,
ExcDimensionMismatch(dst.n_blocks(), this->n_mg_blocks));
+#ifdef DEBUG
+ if (this->mg_constrained_dofs != nullptr)
+ Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+ .get_local_lines()
+ .size() == 0,
+ ExcNotImplemented());
+#endif
+
// Multiplicate with prolongation
// matrix, but only those blocks
// selected.
Assert((to_level >= 1) && (to_level <= prolongation_matrices.size()),
ExcIndexRange(to_level, 1, prolongation_matrices.size() + 1));
+#ifdef DEBUG
+ if (this->mg_constrained_dofs != nullptr)
+ Assert(this->mg_constrained_dofs->get_user_constraint_matrix(to_level - 1)
+ .get_local_lines()
+ .size() == 0,
+ ExcNotImplemented());
+#endif
+
prolongation_matrices[to_level - 1]
->block(selected_block, selected_block)
.vmult(dst, src);
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 - 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Check MGTransferMatrixFree with custom user constraints
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/multigrid/mg_transfer.h>
+#include <deal.II/multigrid/mg_transfer_matrix_free.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+check()
+{
+ Triangulation<dim> tr(Triangulation<dim>::limit_level_difference_at_vertices);
+ GridGenerator::hyper_cube(tr, 0, 1, false);
+
+ typename Triangulation<dim>::active_cell_iterator cell = tr.begin_active();
+ cell->face(0)->set_all_boundary_ids(1);
+
+ tr.refine_global(1);
+
+ FE_Q<dim> fe = FE_Q<dim>(1);
+
+ DoFHandler<dim> mgdof(tr);
+ mgdof.distribute_dofs(fe);
+ mgdof.distribute_mg_dofs();
+
+ MGConstrainedDoFs mg_constrained_dofs;
+ mg_constrained_dofs.initialize(mgdof);
+
+ IndexSet relevant_dofs;
+ DoFTools::extract_locally_relevant_level_dofs(mgdof, 0, relevant_dofs);
+ AffineConstraints<double> user_constraints;
+ user_constraints.reinit(relevant_dofs);
+
+ typename DoFHandler<dim>::level_face_iterator face0 = mgdof.begin(0)->face(0);
+ std::vector<types::global_dof_index> face_dofs(fe.dofs_per_face);
+ face0->get_mg_dof_indices(0, face_dofs);
+ for (unsigned int i = 0; i < face_dofs.size(); ++i)
+ {
+ if (user_constraints.can_store_line(face_dofs[i]))
+ {
+ user_constraints.add_line(face_dofs[i]);
+ user_constraints.set_inhomogeneity(face_dofs[i], 5.0);
+ }
+ }
+ user_constraints.close();
+ mg_constrained_dofs.add_user_constraints(0, user_constraints);
+
+ // build matrix-free transfer
+ MGTransferMatrixFree<dim, double> transfer(mg_constrained_dofs);
+ transfer.build(mgdof);
+
+ LinearAlgebra::distributed::Vector<double> src_level_0(mgdof.n_dofs(0));
+ for (unsigned int i = 0; i < mgdof.n_dofs(0); ++i)
+ deallog << src_level_0(i) << " ";
+ deallog << std::endl;
+
+ LinearAlgebra::distributed::Vector<double> dst_level_1(mgdof.n_dofs(1));
+ transfer.prolongate(1, dst_level_1, src_level_0);
+ for (unsigned int i = 0; i < mgdof.n_dofs(1); ++i)
+ deallog << dst_level_1(i) << " ";
+ deallog << std::endl;
+}
+
+
+int
+main()
+{
+ initlog();
+
+ check<2>();
+ deallog << std::endl;
+ check<3>();
+ deallog << std::endl;
+}
--- /dev/null
+
+DEAL::0.00000 0.00000 0.00000 0.00000
+DEAL::5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 5.00000 2.50000 0.00000
+DEAL::
+DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::5.00000 2.50000 5.00000 2.50000 5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 0.00000 0.00000 5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 5.00000 2.50000 5.00000 2.50000 0.00000 0.00000 5.00000 2.50000 0.00000
+DEAL::
projects / dealii.git / commitdiff