* Authors: Guido Kanschat, University of Heidelberg, 2003
* Baerbel Janssen, University of Heidelberg, 2010
* Wolfgang Bangerth, Texas A&M University, 2010
+ * Timo Heister, Clemson University, 2018
*/
#include <deal.II/multigrid/mg_smoother.h>
#include <deal.II/multigrid/mg_matrix.h>
-// Finally we include the MeshWorker framework. This framework through its
-// function loop() and integration_loop(), automates loops over cells and
-// assembling of data into vectors, matrices, etc. It obeys constraints
-// automatically. Since we have to build several matrices and have to be aware
-// of several sets of constraints, this will save us a lot of headache.
-#include <deal.II/meshworker/dof_info.h>
-#include <deal.II/meshworker/integration_info.h>
-#include <deal.II/meshworker/simple.h>
-#include <deal.II/meshworker/output.h>
-#include <deal.II/meshworker/loop.h>
-
-// In order to save effort, we use the pre-implemented Laplacian found in
-#include <deal.II/integrators/laplace.h>
-#include <deal.II/integrators/l2.h>
+// We will be using MeshWorker::mesh_loop to loop over the cells, so include it
+// here:
+#include <deal.II/meshworker/mesh_loop.h>
+
// This is C++:
#include <iostream>
namespace Step16
{
- // @sect3{The integrator on each cell}
-
- // The MeshWorker::integration_loop() expects a class that provides functions
- // for integration on cells and boundary and interior faces. This is done by
- // the following class. In the constructor, we tell the loop that cell
- // integrals should be computed (the 'true'), but integrals should not be
- // computed on boundary and interior faces (the two 'false'). Accordingly, we
- // only need a cell function, but none for the faces.
+ // @sect3{The Scratch and Copy objects}
+ //
+ // We use MeshWorker::mesh_loop() to assemble our matrices. For this, we need
+ // a ScratchData object to store temporary data on each cell (this is just the
+ // FEValues object) and a CopyData object that will contain the output of each
+ // cell assembly.
template <int dim>
- class LaplaceIntegrator : public MeshWorker::LocalIntegrator<dim>
+ struct ScratchData
{
- public:
- LaplaceIntegrator();
- virtual void cell(MeshWorker::DoFInfo<dim> & dinfo,
- MeshWorker::IntegrationInfo<dim> &info) const override;
+ ScratchData(const Mapping<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const unsigned int quadrature_degree,
+ const UpdateFlags update_flags)
+ : fe_values(mapping, fe, QGauss<dim>(quadrature_degree), update_flags)
+ {}
+
+ ScratchData(const ScratchData<dim> &scratch_data)
+ : fe_values(scratch_data.fe_values.get_mapping(),
+ scratch_data.fe_values.get_fe(),
+ scratch_data.fe_values.get_quadrature(),
+ scratch_data.fe_values.get_update_flags())
+ {}
+
+ FEValues<dim> fe_values;
};
-
- template <int dim>
- LaplaceIntegrator<dim>::LaplaceIntegrator()
- : MeshWorker::LocalIntegrator<dim>(true, false, false)
- {}
-
-
- // Next the actual integrator on each cell. We solve a Poisson problem with a
- // coefficient one in the right half plane and one tenth in the left
- // half plane.
-
- // The MeshWorker::LocalResults base class of MeshWorker::DoFInfo contains
- // objects that can be filled in this local integrator. How many objects are
- // created is determined inside the MeshWorker framework by the assembler
- // class. Here, we test for instance that one matrix is required
- // (MeshWorker::LocalResults::n_matrices()). The matrices are accessed through
- // MeshWorker::LocalResults::matrix(), which takes the number of the matrix as
- // its first argument. The second argument is only used for integrals over
- // faces when there are two matrices for each test function used. Then, a
- // second matrix with indicator 'true' would exist with the same index.
-
- // MeshWorker::IntegrationInfo provides one or several FEValues objects, which
- // below are used by LocalIntegrators::Laplace::cell_matrix() or
- // LocalIntegrators::L2::L2(). Since we are assembling only a single PDE,
- // there is also only one of these objects with index zero.
-
- // In addition, we note that this integrator serves to compute the matrices
- // for the multilevel preconditioner as well as the matrix and the right hand
- // side for the global system. Since the assembler for a system requires an
- // additional vector, MeshWorker::LocalResults::n_vectors() is returning a
- // nonzero value. Accordingly, we fill a right hand side vector at the end of
- // this function. Since LocalResults can deal with several BlockVector
- // objects, but we are again in the simplest case here, we enter the
- // information into block zero of vector zero.
- template <int dim>
- void
- LaplaceIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> & dinfo,
- MeshWorker::IntegrationInfo<dim> &info) const
+ struct CopyData
{
- AssertDimension(dinfo.n_matrices(), 1);
- const double coefficient = (dinfo.cell->center()(0) > 0.) ? .1 : 1.;
+ unsigned int level;
+ FullMatrix<double> cell_matrix;
+ Vector<double> cell_rhs;
+ std::vector<types::global_dof_index> local_dof_indices;
- LocalIntegrators::Laplace::cell_matrix(dinfo.matrix(0, false).matrix,
- info.fe_values(0),
- coefficient);
-
- if (dinfo.n_vectors() > 0)
- {
- std::vector<double> rhs(info.fe_values(0).n_quadrature_points, 1.);
- LocalIntegrators::L2::L2(dinfo.vector(0).block(0),
- info.fe_values(0),
- rhs);
- }
- }
+ template <class Iterator>
+ void reinit(const Iterator &cell, unsigned int dofs_per_cell)
+ {
+ cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+ cell_rhs.reinit(dofs_per_cell);
+ local_dof_indices.resize(dofs_per_cell);
+ cell->get_active_or_mg_dof_indices(local_dof_indices);
+ level = cell->level();
+ }
+ };
// @sect3{The <code>LaplaceProblem</code> class template}
- // This main class is basically the same class as in step-6. As far as
- // member functions is concerned, the only addition is the
- // <code>assemble_multigrid</code> function that assembles the matrices that
- // correspond to the discrete operators on intermediate levels:
+ // This main class is similar to the same class in step-6. As far as
+ // member functions is concerned, the only additions are:
+ // - The <code>assemble_multigrid</code> function that assembles the matrices
+ // that correspond to the discrete operators on intermediate levels.
+ // - The <code>cell_worker</code> function that assembles our PDE on a single
+ // cell.
template <int dim>
class LaplaceProblem
{
void run();
private:
+ template <class Iterator>
+ void cell_worker(const Iterator & cell,
+ ScratchData<dim> &scratch_data,
+ CopyData & copy_data);
+
void setup_system();
void assemble_system();
void assemble_multigrid();
const unsigned int degree;
// The following members are the essential data structures for the multigrid
- // method. The first two represent the sparsity patterns and the matrices on
- // individual levels of the multilevel hierarchy, very much like the objects
- // for the global mesh above.
+ // method. The first four represent the sparsity patterns and the matrices
+ // on individual levels of the multilevel hierarchy, very much like the
+ // objects for the global mesh above.
//
// Then we have two new matrices only needed for multigrid methods with
// local smoothing on adaptive meshes. They convey data between the interior
// level and information about indices lying on a refinement edge between
// two different refinement levels. It thus serves a similar purpose as
// AffineConstraints, but on each level.
- MGLevelObject<SparsityPattern> mg_sparsity_patterns;
+ MGLevelObject<SparsityPattern> mg_sparsity_patterns;
+ MGLevelObject<SparsityPattern> mg_interface_sparsity_patterns;
+
MGLevelObject<SparseMatrix<double>> mg_matrices;
- MGLevelObject<SparseMatrix<double>> mg_interface_in;
- MGLevelObject<SparseMatrix<double>> mg_interface_out;
+ MGLevelObject<SparseMatrix<double>> mg_interface_matrices;
MGConstrainedDoFs mg_constrained_dofs;
};
dof_handler.distribute_dofs(fe);
dof_handler.distribute_mg_dofs();
- deallog << " Number of degrees of freedom: " << dof_handler.n_dofs()
- << " (by level: ";
+ std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
+ << " (by level: ";
for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
- deallog << dof_handler.n_dofs(level)
- << (level == triangulation.n_levels() - 1 ? ")" : ", ");
- deallog << std::endl;
+ std::cout << dof_handler.n_dofs(level)
+ << (level == triangulation.n_levels() - 1 ? ")" : ", ");
+ std::cout << std::endl;
- DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
- DoFTools::make_sparsity_pattern(dof_handler, dsp);
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
const std::map<types::boundary_id, const Function<dim> *>
dirichlet_boundary_functions = {
{types::boundary_id(0), &homogeneous_dirichlet_bc}};
- VectorTools::interpolate_boundary_values(
- static_cast<const DoFHandler<dim> &>(dof_handler),
- dirichlet_boundary_functions,
- constraints);
+ VectorTools::interpolate_boundary_values(dof_handler,
+ dirichlet_boundary_functions,
+ constraints);
constraints.close();
- constraints.condense(dsp);
- sparsity_pattern.copy_from(dsp);
+
+ {
+ DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+ DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints);
+ sparsity_pattern.copy_from(dsp);
+ }
system_matrix.reinit(sparsity_pattern);
// The multigrid constraints have to be initialized. They need to know
- // about the boundary values as well, so we pass the
- // <code>dirichlet_boundary</code> here as well.
+ // where Dirichlet boundary conditions are prescribed.
mg_constrained_dofs.clear();
mg_constrained_dofs.initialize(dof_handler);
mg_constrained_dofs.make_zero_boundary_constraints(dof_handler,
// their SparsityPattern before the can be destroyed upon resizing.
const unsigned int n_levels = triangulation.n_levels();
- mg_interface_in.resize(0, n_levels - 1);
- mg_interface_in.clear_elements();
- mg_interface_out.resize(0, n_levels - 1);
- mg_interface_out.clear_elements();
+ mg_interface_matrices.resize(0, n_levels - 1);
mg_matrices.resize(0, n_levels - 1);
- mg_matrices.clear_elements();
mg_sparsity_patterns.resize(0, n_levels - 1);
+ mg_interface_sparsity_patterns.resize(0, n_levels - 1);
// Now, we have to provide a matrix on each level. To this end, we first use
// the MGTools::make_sparsity_pattern function to generate a preliminary
// compressed sparsity pattern on each level (see the @ref Sparsity module
// for more information on this topic) and then copy it over to the one we
- // really want. The next step is to initialize both kinds of level matrices
- // with these sparsity patterns.
+ // really want. The next step is to initialize the interface matrices with
+ // the fitting sparsity pattern.
//
// It may be worth pointing out that the interface matrices only have
// entries for degrees of freedom that sit at or next to the interface
// between coarser and finer levels of the mesh. They are therefore even
// sparser than the matrices on the individual levels of our multigrid
- // hierarchy. If we were more concerned about memory usage (and possibly the
- // speed with which we can multiply with these matrices), we should use
- // separate and different sparsity patterns for these two kinds of matrices.
+ // hierarchy. Therefore, we use a function specifically build for this
+ // purpose to generate it.
for (unsigned int level = 0; level < n_levels; ++level)
{
- DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
- dof_handler.n_dofs(level));
- MGTools::make_sparsity_pattern(dof_handler, dsp, level);
+ {
+ DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
+ dof_handler.n_dofs(level));
+ MGTools::make_sparsity_pattern(dof_handler, dsp, level);
+
+ mg_sparsity_patterns[level].copy_from(dsp);
+ mg_matrices[level].reinit(mg_sparsity_patterns[level]);
+ }
+ {
+ DynamicSparsityPattern dsp(dof_handler.n_dofs(level),
+ dof_handler.n_dofs(level));
+ MGTools::make_interface_sparsity_pattern(dof_handler,
+ mg_constrained_dofs,
+ dsp,
+ level);
+ mg_interface_sparsity_patterns[level].copy_from(dsp);
+ mg_interface_matrices[level].reinit(
+ mg_interface_sparsity_patterns[level]);
+ }
+ }
+ }
+
+
+ // @sect4{LaplaceProblem::cell_worker}
- mg_sparsity_patterns[level].copy_from(dsp);
+ // The cell_worker function is used to assemble the matrix and right-hand side
+ // on the given cell. This function is used for the active cells to generate
+ // the system_matrix and on each level to build the level matrices.
+ //
+ // Note that we also assemble a right-hand side when called from
+ // assemble_multigrid() even though it is not used.
+ template <int dim>
+ template <class Iterator>
+ void LaplaceProblem<dim>::cell_worker(const Iterator & cell,
+ ScratchData<dim> &scratch_data,
+ CopyData & copy_data)
+ {
+ FEValues<dim> &fe_values = scratch_data.fe_values;
+ fe_values.reinit(cell);
+
+ const unsigned int dofs_per_cell = fe_values.get_fe().dofs_per_cell;
+ const unsigned int n_q_points = fe_values.get_quadrature().size();
+
+ copy_data.reinit(cell, dofs_per_cell);
+
+ const std::vector<double> &JxW = fe_values.get_JxW_values();
+
+ for (unsigned int q = 0; q < n_q_points; ++q)
+ {
+ const double coefficient =
+ (fe_values.get_quadrature_points()[q][0] < 0.0) ? 1.0 : 0.1;
+ //(cell->center().square() < 0.5 * 0.5) ? 10.0:1.0;
- mg_matrices[level].reinit(mg_sparsity_patterns[level]);
- mg_interface_in[level].reinit(mg_sparsity_patterns[level]);
- mg_interface_out[level].reinit(mg_sparsity_patterns[level]);
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ {
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ {
+ copy_data.cell_matrix(i, j) +=
+ coefficient *
+ (fe_values.shape_grad(i, q) * fe_values.shape_grad(j, q)) *
+ JxW[q];
+ }
+ copy_data.cell_rhs(i) += 1.0 * fe_values.shape_value(i, q) * JxW[q];
+ }
}
}
+
// @sect4{LaplaceProblem::assemble_system}
- // The following function assembles the linear system on the finest level of
- // the mesh. Since we want to reuse the code here for the level assembling
- // below, we use the local integrator class LaplaceIntegrator and leave the
- // loops to the MeshWorker framework. Thus, this function first sets up the
- // objects necessary for this framework, namely
- // - a MeshWorker::IntegrationInfoBox object, which will provide all the
- // required data in quadrature points on the cell. This object can be seen
- // as an extension of FEValues, providing a lot more useful information,
- // - a MeshWorker::DoFInfo object, which on the one hand side extends the
- // functionality of cell iterators, but also provides space for return
- // values in its base class LocalResults,
- // - an assembler, in this case for the whole system. The term 'simple' here
- // refers to the fact that the global system does not have a block
- // structure,
- // - the local integrator, which implements the actual forms.
- //
- // After the loop has combined all of these into a matrix and a right hand
- // side, there is one thing left to do: the assemblers leave matrix rows and
- // columns of constrained degrees of freedom untouched. Therefore, we put a
- // one on the diagonal to make the whole system well posed. The value one, or
- // any fixed value has the advantage, that its effect on the spectrum of the
- // matrix is easily understood. Since the corresponding eigenvectors form an
- // invariant subspace, the value chosen does not affect the convergence of
- // Krylov space solvers.
+ // The following function assembles the linear system on the active cells of
+ // the mesh. For this, we pass two lambda functions to the mesh_loop()
+ // function. The cell_worker function redirects to the class member function
+ // of the same name, while the copier is specific to this function and copies
+ // local matrix and vector to the corresponding global ones using the
+ // constraints.
template <int dim>
void LaplaceProblem<dim>::assemble_system()
{
- MappingQ1<dim> mapping;
- MeshWorker::IntegrationInfoBox<dim> info_box;
- UpdateFlags update_flags =
- update_values | update_gradients | update_hessians;
- info_box.add_update_flags_all(update_flags);
- info_box.initialize(fe, mapping);
-
- MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
- MeshWorker::Assembler::SystemSimple<SparseMatrix<double>, Vector<double>>
- assembler;
- assembler.initialize(constraints);
- assembler.initialize(system_matrix, system_rhs);
-
- LaplaceIntegrator<dim> matrix_integrator;
- MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
- dof_handler.end(),
- dof_info,
- info_box,
- matrix_integrator,
- assembler);
-
- for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
- if (constraints.is_constrained(i))
- system_matrix.set(i, i, 1.);
+ MappingQ1<dim> mapping;
+
+ typedef decltype(dof_handler.begin_active()) Iterator;
+
+ auto cell_worker = [&](const Iterator & cell,
+ ScratchData<dim> &scratch_data,
+ CopyData & copy_data) {
+ this->cell_worker(cell, scratch_data, copy_data);
+ };
+
+ auto copier = [&](const CopyData &cd) {
+ this->constraints.distribute_local_to_global(cd.cell_matrix,
+ cd.cell_rhs,
+ cd.local_dof_indices,
+ system_matrix,
+ system_rhs);
+ };
+
+ const unsigned int n_gauss_points = degree + 1;
+
+ ScratchData<dim> scratch_data(mapping,
+ fe,
+ n_gauss_points,
+ update_values | update_gradients |
+ update_JxW_values |
+ update_quadrature_points);
+
+ MeshWorker::mesh_loop(dof_handler.begin_active(),
+ dof_handler.end(),
+ cell_worker,
+ copier,
+ scratch_data,
+ CopyData(),
+ MeshWorker::assemble_own_cells);
}
// @sect4{LaplaceProblem::assemble_multigrid}
- // The next function is the one that builds the linear operators (matrices)
+ // The next function is the one that builds the matrices
// that define the multigrid method on each level of the mesh. The integration
// core is the same as above, but the loop below will go over all existing
// cells instead of just the active ones, and the results must be entered into
// the correct level matrices. Fortunately, MeshWorker hides most of that from
// us, and thus the difference between this function and the previous lies
// only in the setup of the assembler and the different iterators in the loop.
- // Also, fixing up the matrices in the end is a little more complicated.
+ //
+ // We generate an AffineConstraints<> object
template <int dim>
void LaplaceProblem<dim>::assemble_multigrid()
{
- MappingQ1<dim> mapping;
- MeshWorker::IntegrationInfoBox<dim> info_box;
- UpdateFlags update_flags =
- update_values | update_gradients | update_hessians;
- info_box.add_update_flags_all(update_flags);
- info_box.initialize(fe, mapping);
-
- MeshWorker::DoFInfo<dim> dof_info(dof_handler);
-
- MeshWorker::Assembler::MGMatrixSimple<SparseMatrix<double>> assembler;
- assembler.initialize(mg_constrained_dofs);
- assembler.initialize(mg_matrices);
- assembler.initialize_interfaces(mg_interface_in, mg_interface_out);
-
- LaplaceIntegrator<dim> matrix_integrator;
- MeshWorker::integration_loop<dim, dim>(dof_handler.begin_mg(),
- dof_handler.end_mg(),
- dof_info,
- info_box,
- matrix_integrator,
- assembler);
-
- const unsigned int nlevels = triangulation.n_levels();
- for (unsigned int level = 0; level < nlevels; ++level)
+ MappingQ1<dim> mapping;
+ const unsigned int n_levels = triangulation.n_levels();
+
+ std::vector<AffineConstraints<>> boundary_constraints(n_levels);
+ for (unsigned int level = 0; level < n_levels; ++level)
{
- for (unsigned int i = 0; i < dof_handler.n_dofs(level); ++i)
- if (mg_constrained_dofs.is_boundary_index(level, i) ||
- mg_constrained_dofs.at_refinement_edge(level, i))
- mg_matrices[level].set(i, i, 1.);
+ IndexSet dofset;
+ DoFTools::extract_locally_relevant_level_dofs(dof_handler,
+ level,
+ dofset);
+ boundary_constraints[level].reinit(dofset);
+ boundary_constraints[level].add_lines(
+ mg_constrained_dofs.get_refinement_edge_indices(level));
+ boundary_constraints[level].add_lines(
+ mg_constrained_dofs.get_boundary_indices(level));
+ boundary_constraints[level].close();
}
+
+ typedef decltype(dof_handler.begin_mg()) Iterator;
+
+ auto cell_worker = [&](const Iterator & cell,
+ ScratchData<dim> &scratch_data,
+ CopyData & copy_data) {
+ this->cell_worker(cell, scratch_data, copy_data);
+ };
+
+ auto copier = [&](const CopyData &cd) {
+ boundary_constraints[cd.level].distribute_local_to_global(
+ cd.cell_matrix, cd.local_dof_indices, mg_matrices[cd.level]);
+
+ const unsigned int dofs_per_cell = cd.local_dof_indices.size();
+
+ // TODO EXPLAIN:
+
+ for (unsigned int i = 0; i < dofs_per_cell; ++i)
+ for (unsigned int j = 0; j < dofs_per_cell; ++j)
+ if (mg_constrained_dofs.is_interface_matrix_entry(
+ cd.level, cd.local_dof_indices[i], cd.local_dof_indices[j]))
+ {
+ mg_interface_matrices[cd.level].add(cd.local_dof_indices[i],
+ cd.local_dof_indices[j],
+ cd.cell_matrix(i, j));
+ }
+ };
+
+ const unsigned int n_gauss_points = degree + 1;
+
+ ScratchData<dim> scratch_data(mapping,
+ fe,
+ n_gauss_points,
+ update_values | update_gradients |
+ update_JxW_values |
+ update_quadrature_points);
+
+ MeshWorker::mesh_loop(dof_handler.begin_mg(),
+ dof_handler.end_mg(),
+ cell_worker,
+ copier,
+ scratch_data,
+ CopyData(),
+ MeshWorker::assemble_own_cells);
}
// initialize both up and down versions of the operator with the matrices
// we already built:
mg::Matrix<Vector<double>> mg_matrix(mg_matrices);
- mg::Matrix<Vector<double>> mg_interface_up(mg_interface_in);
- mg::Matrix<Vector<double>> mg_interface_down(mg_interface_out);
+ mg::Matrix<Vector<double>> mg_interface_up(mg_interface_matrices);
+ mg::Matrix<Vector<double>> mg_interface_down(mg_interface_matrices);
// Now, we are ready to set up the V-cycle operator and the multilevel
// preconditioner.
solution = 0;
solver.solve(system_matrix, solution, system_rhs, preconditioner);
+ std::cout << " Number of CG iterations: " << solver_control.last_step()
+ << "\n"
+ << std::endl;
constraints.distribute(solution);
}
// The following two functions postprocess a solution once it is
// computed. In particular, the first one refines the mesh at the beginning
// of each cycle while the second one outputs results at the end of each
- // such cycle. The functions are almost unchanged from those in step-6, with
- // the exception of one minor difference: we generate output in VTK
- // format, to use the more modern visualization programs available today
- // compared to those that were available when step-6 was written.
+ // such cycle. The functions are almost unchanged from those in step-6.
template <int dim>
void LaplaceProblem<dim>::refine_grid()
{
KellyErrorEstimator<dim>::estimate(
dof_handler,
- QGauss<dim - 1>(3),
+ QGauss<dim - 1>(degree + 2),
std::map<types::boundary_id, const Function<dim> *>(),
solution,
estimated_error_per_cell);
{
for (unsigned int cycle = 0; cycle < 8; ++cycle)
{
- deallog << "Cycle " << cycle << std::endl;
+ std::cout << "Cycle " << cycle << std::endl;
if (cycle == 0)
{
GridGenerator::hyper_ball(triangulation);
- triangulation.refine_global(1);
+ triangulation.refine_global(2);
}
else
refine_grid();
- deallog << " Number of active cells: "
- << triangulation.n_active_cells() << std::endl;
+ std::cout << " Number of active cells: "
+ << triangulation.n_active_cells() << std::endl;
setup_system();
{
using namespace Step16;
- deallog.depth_console(2);
-
LaplaceProblem<2> laplace_problem(1);
laplace_problem.run();
}