Reference documentation for deal.II version GIT relicensing-669-gbd74b205fb 2024-05-18 16:30:02+00:00
Searching...
No Matches
The step-39 tutorial program

This tutorial depends on step-12.

1. Introduction
2. The commented program
1. Results
2. The plain program

In this program, we use the interior penalty method and Nitsche's weak boundary conditions to solve Poisson's equation. We use multigrid methods on locally refined meshes, which are generated using a bulk criterion and a standard error estimator based on cell and face residuals. All operators are implemented using the MeshWorker interface.

Like in step-12, the discretization relies on finite element spaces, which are polynomial inside the mesh cells $$K\in \mathbb T_h$$, but have no continuity between cells. Since such functions have two values on each interior face $$F\in \mathbb F_h^i$$, one from each side, we define mean value and jump operators as follows: let K1 and K2 be the two cells sharing a face, and let the traces of functions ui and the outer normal vectors ni be labeled accordingly. Then, on the face, we let

$\average{ u } = \frac{u_1 + u_2}2$

Note, that if such an expression contains a normal vector, the averaging operator turns into a jump. The interior penalty method for the problem

$-\Delta u = f \text{ in }\Omega \qquad u = u^D \text{ on } \partial\Omega$

becomes

\begin{multline*} \sum_{K\in \mathbb T_h} (\nabla u, \nabla v)_K \\ + \sum_{F \in F_h^i} \biggl\{4\sigma_F (\average{ u \mathbf n}, \average{ v \mathbf n })_F - 2 (\average{ \nabla u },\average{ v\mathbf n })_F - 2 (\average{ \nabla v },\average{ u\mathbf n })_F \biggr\} \\ + \sum_{F \in F_h^b} \biggl\{2\sigma_F (u, v)_F - (\partial_n u,v)_F - (\partial_n v,u)_F \biggr\} \\ = (f, v)_\Omega + \sum_{F \in F_h^b} \biggl\{ 2\sigma_F (u^D, v)_F - (\partial_n v,u^D)_F \biggr\}. \end{multline*}

Here, $$\sigma_F$$ is the penalty parameter, which is chosen as follows: for a face F of a cell K, compute the value

$\sigma_{F,K} = p(p+1) \frac{|F|_{d-1}}{|K|_d},$

where p is the polynomial degree of the finite element functions and $$|\cdot|_d$$ and $$|\cdot|_{d-1}$$ denote the $$d$$ and $$d-1$$ dimensional Hausdorff measure of the corresponding object. If the face is at the boundary, choose $$\sigma_F = \sigma_{F,K}$$. For an interior face, we take the average of the two values at this face.

In our finite element program, we distinguish three different integrals, corresponding to the sums over cells, interior faces and boundary faces above. Since the MeshWorker::loop organizes the sums for us, we only need to implement the integrals over each mesh element. The class MatrixIntegrator below has these three functions for the left hand side of the formula, the class RHSIntegrator for the right.

As we will see below, even the error estimate is of the same structure, since it can be written as

\begin{align*} \eta^2 &= \eta_K^2 + \eta_F^2 + \eta_B^2 \\ \eta_K^2 &= \sum_{K\in \mathbb T_h} h^2 \|f + \Delta u_h\|^2 \\ \eta_F^2 &= \sum_{F \in F_h^i} \biggl\{ 4 \sigma_F \| \average{u_h\mathbf n} \|^2 + h \|\average{\partial_n u_h}\|^2 \biggr\} \\ \eta_B^2 &= \sum_{F \in F_h^b} 2\sigma_F \| u_h-u^D \|^2. \end{align*}

Thus, the functions for assembling matrices, right hand side and error estimates below exhibit that these loops are all generic and can be programmed in the same way.

This program is related to step-12, in that it uses MeshWorker and discontinuous Galerkin methods. There we solved an advection problem, while here it is a diffusion problem. Here, we also use multigrid preconditioning and a theoretically justified error estimator, see Karakashian and Pascal (2003). The multilevel scheme was discussed in detail in Kanschat (2004). The adaptive iteration and its convergence have been discussed (for triangular meshes) in Hoppe, Kanschat, and Warburton (2009).

# The commented program

The include files for the linear algebra: A regular SparseMatrix, which in turn will include the necessary files for SparsityPattern and Vector classes.

#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/precondition_block.h>
#include <deal.II/lac/block_vector.h>

Include files for setting up the mesh

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>

Include files for FiniteElement classes and DoFHandler.

#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/mapping_q1.h>

The include files for using the MeshWorker framework

#include <deal.II/meshworker/dof_info.h>
#include <deal.II/meshworker/integration_info.h>
#include <deal.II/meshworker/assembler.h>
#include <deal.II/meshworker/loop.h>

The include file for local integrators associated with the Laplacian

#include <deal.II/integrators/laplace.h>

Support for multigrid methods

#include <deal.II/multigrid/mg_tools.h>
#include <deal.II/multigrid/multigrid.h>
#include <deal.II/multigrid/mg_matrix.h>
#include <deal.II/multigrid/mg_transfer.h>
#include <deal.II/multigrid/mg_coarse.h>
#include <deal.II/multigrid/mg_smoother.h>

Finally, we take our exact solution from the library as well as quadrature and additional tools.

#include <deal.II/base/function_lib.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/data_out.h>

#include <iostream>
#include <fstream>

All classes of the deal.II library are in the namespace dealii. In order to save typing, we tell the compiler to search names in there as well.

namespace Step39
{
using namespace dealii;

This is the function we use to set the boundary values and also the exact solution we compare to.

### The local integrators

MeshWorker separates local integration from the loops over cells and faces. Thus, we have to write local integration classes for generating matrices, the right hand side and the error estimator.

All these classes have the same three functions for integrating over cells, boundary faces and interior faces, respectively. All the information needed for the local integration is provided by MeshWorker::IntegrationInfo<dim>. Note that the signature of the functions cannot be changed, because it is expected by MeshWorker::integration_loop().

The first class defining local integrators is responsible for computing cell and face matrices. It is used to assemble the global matrix as well as the level matrices.

template <int dim>
class MatrixIntegrator : public MeshWorker::LocalIntegrator<dim>
{
public:
MeshWorker::IntegrationInfo<dim> &info) const override;
MeshWorker::IntegrationInfo<dim> &info) const override;
MeshWorker::IntegrationInfo<dim> &info2) const override;
};

virtual void cell(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const
virtual void face(DoFInfo< dim, spacedim, number > &dinfo1, DoFInfo< dim, spacedim, number > &dinfo2, IntegrationInfo< dim, spacedim > &info1, IntegrationInfo< dim, spacedim > &info2) const
virtual void boundary(DoFInfo< dim, spacedim, number > &dinfo, IntegrationInfo< dim, spacedim > &info) const

On each cell, we integrate the Dirichlet form as well as the Nitsche boundary conditions and the interior penalty fluxes between cells.

The boundary and flux terms need a penalty parameter, which should be adjusted to the cell size and the polynomial degree. We compute it in two steps: First, we compute on each cell $$K_i$$ the value $$P_i = p_i(p_i+1)/h_i$$, where $$p_i$$ is the polynomial degree on cell $$K_i$$ and $$h_i$$ is the length of $$K_i$$ orthogonal to the current face. Second, if exactly one of the two cells adjacent to the face has children, its penalty is multiplied by two (to account for the fact that the mesh size $$h_i$$ there is only half that previously computed); it is possible that both adjacent cells are refined, in which case we are integrating over a non-active face and no adjustment is necessary. Finally, we return the average of the two penalty values.

template <int dim>
double ip_penalty_factor(const MeshWorker::DoFInfo<dim> &dinfo1,
const MeshWorker::DoFInfo<dim> &dinfo2,
unsigned int deg1,
unsigned int deg2)
{
const unsigned int normal1 =
const unsigned int normal2 =
const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1 + 1);
const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2 + 1);

double penalty1 = deg1sq / dinfo1.cell->extent_in_direction(normal1);
double penalty2 = deg2sq / dinfo2.cell->extent_in_direction(normal2);
if (dinfo1.cell->has_children() && !dinfo2.cell->has_children())
penalty1 *= 2;
else if (!dinfo1.cell->has_children() && dinfo2.cell->has_children())
penalty2 *= 2;

const double penalty = 0.5 * (penalty1 + penalty2);
return penalty;
}

template <int dim>
void MatrixIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo,
{
FullMatrix<double> &M = dinfo.matrix(0, false).matrix;

for (unsigned int k = 0; k < info.fe_values().n_quadrature_points; ++k)
{
const double dx = info.fe_values().JxW(k);

for (unsigned int i = 0; i < info.fe_values().dofs_per_cell; ++i)
{
const double Mii = (info.fe_values().shape_grad(i, k) *

M(i, i) += Mii;

for (unsigned int j = i + 1; j < info.fe_values().dofs_per_cell;
++j)
{
const double Mij = info.fe_values().shape_grad(j, k) *

M(i, j) += Mij;
M(j, i) += Mij;
}
}
}
}

unsigned int face_number
Definition dof_info.h:89
Triangulation< dim, spacedim >::cell_iterator cell
The current cell.
Definition dof_info.h:78

Boundary faces use the Nitsche method to impose boundary values:

template <int dim>
void
MatrixIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe_face_values = info.fe_values(0);

FullMatrix<double> &M = dinfo.matrix(0, false).matrix;
AssertDimension(M.n(), fe_face_values.dofs_per_cell);
AssertDimension(M.m(), fe_face_values.dofs_per_cell);

const unsigned int polynomial_degree =
info.fe_values(0).get_fe().tensor_degree();

const double ip_penalty =
ip_penalty_factor(dinfo, dinfo, polynomial_degree, polynomial_degree);

for (unsigned int k = 0; k < fe_face_values.n_quadrature_points; ++k)
{
const double dx = fe_face_values.JxW(k);
const Tensor<1, dim> &n = fe_face_values.normal_vector(k);

for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i)
for (unsigned int j = 0; j < fe_face_values.dofs_per_cell; ++j)
M(i, j) += (2. * fe_face_values.shape_value(i, k) * ip_penalty *
fe_face_values.shape_value(j, k) -
fe_face_values.shape_value(j, k) -
fe_face_values.shape_value(i, k)) *
dx;
}
}

#define AssertDimension(dim1, dim2)

Interior faces use the interior penalty method:

template <int dim>
void
MatrixIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,
{
const FEValuesBase<dim> &fe_face_values_1 = info1.fe_values(0);
const FEValuesBase<dim> &fe_face_values_2 = info2.fe_values(0);

FullMatrix<double> &M11 = dinfo1.matrix(0, false).matrix;
FullMatrix<double> &M12 = dinfo1.matrix(0, true).matrix;
FullMatrix<double> &M21 = dinfo2.matrix(0, true).matrix;
FullMatrix<double> &M22 = dinfo2.matrix(0, false).matrix;

AssertDimension(M11.n(), fe_face_values_1.dofs_per_cell);
AssertDimension(M11.m(), fe_face_values_1.dofs_per_cell);
AssertDimension(M12.n(), fe_face_values_1.dofs_per_cell);
AssertDimension(M12.m(), fe_face_values_1.dofs_per_cell);
AssertDimension(M21.n(), fe_face_values_1.dofs_per_cell);
AssertDimension(M21.m(), fe_face_values_1.dofs_per_cell);
AssertDimension(M22.n(), fe_face_values_1.dofs_per_cell);
AssertDimension(M22.m(), fe_face_values_1.dofs_per_cell);

const unsigned int polynomial_degree =
info1.fe_values(0).get_fe().tensor_degree();
const double ip_penalty =
ip_penalty_factor(dinfo1, dinfo2, polynomial_degree, polynomial_degree);

const double nui = 1.;
const double nue = 1.;
const double nu = .5 * (nui + nue);

for (unsigned int k = 0; k < fe_face_values_1.n_quadrature_points; ++k)
{
const double dx = fe_face_values_1.JxW(k);
const Tensor<1, dim> &n = fe_face_values_1.normal_vector(k);

for (unsigned int i = 0; i < fe_face_values_1.dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < fe_face_values_1.dofs_per_cell; ++j)
{
const double vi = fe_face_values_1.shape_value(i, k);
const double dnvi = n * fe_face_values_1.shape_grad(i, k);
const double ve = fe_face_values_2.shape_value(i, k);
const double dnve = n * fe_face_values_2.shape_grad(i, k);
const double ui = fe_face_values_1.shape_value(j, k);
const double dnui = n * fe_face_values_1.shape_grad(j, k);
const double ue = fe_face_values_2.shape_value(j, k);
const double dnue = n * fe_face_values_2.shape_grad(j, k);

M11(i, j) += (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
nu * ip_penalty * ui * vi) *
dx;
M12(i, j) += (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
nu * ip_penalty * vi * ue) *
dx;
M21(i, j) += (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
nu * ip_penalty * ui * ve) *
dx;
M22(i, j) += (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
nu * ip_penalty * ue * ve) *
dx;
}
}
}
}

MatrixType matrix
MatrixBlock< FullMatrix< number > > & matrix(const unsigned int i, const bool external=false)

The second local integrator builds the right hand side. In our example, the right hand side function is zero, such that only the boundary condition is set here in weak form.

template <int dim>
class RHSIntegrator : public MeshWorker::LocalIntegrator<dim>
{
public:
MeshWorker::IntegrationInfo<dim> &info) const override;
MeshWorker::IntegrationInfo<dim> &info) const override;
MeshWorker::IntegrationInfo<dim> &info2) const override;
};

template <int dim>
void RHSIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &,
{}

template <int dim>
void
RHSIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();
Vector<double> &local_vector = dinfo.vector(0).block(0);

const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() /
dinfo.cell->measure();

for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
local_vector(i) +=
(-penalty * fe.shape_value(i, k) // (-sigma * v_i(x_k)
* boundary_values[k] * fe.JxW(k); // u^D(x_k) * dx
}

template <int dim>
void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &,
{}

The third local integrator is responsible for the contributions to the error estimate. This is the standard energy estimator due to Karakashian and Pascal (2003).

template <int dim>
class Estimator : public MeshWorker::LocalIntegrator<dim>
{
public:
MeshWorker::IntegrationInfo<dim> &info) const override;
MeshWorker::IntegrationInfo<dim> &info) const override;
MeshWorker::IntegrationInfo<dim> &info2) const override;
};

The cell contribution is the Laplacian of the discrete solution, since the right hand side is zero.

template <int dim>
void Estimator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();

const std::vector<Tensor<2, dim>> &DDuh = info.hessians[0][0];
for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double t = dinfo.cell->diameter() * trace(DDuh[k]);
dinfo.value(0) += t * t * fe.JxW(k);
}
dinfo.value(0) = std::sqrt(dinfo.value(0));
}

::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
DEAL_II_HOST constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)

At the boundary, we use simply a weighted form of the boundary residual, namely the norm of the difference between the finite element solution and the correct boundary condition.

template <int dim>
void Estimator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();

const std::vector<double> &uh = info.values[0][0];

const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() /
dinfo.cell->measure();

for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double diff = boundary_values[k] - uh[k];
dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
}
dinfo.value(0) = std::sqrt(dinfo.value(0));
}

Finally, on interior faces, the estimator consists of the jumps of the solution and its normal derivative, weighted appropriately.

template <int dim>
void Estimator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,
{
const FEValuesBase<dim> &fe = info1.fe_values();
const std::vector<double> &uh1 = info1.values[0][0];
const std::vector<double> &uh2 = info2.values[0][0];
const std::vector<Tensor<1, dim>> &Duh1 = info1.gradients[0][0];
const std::vector<Tensor<1, dim>> &Duh2 = info2.gradients[0][0];

const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty1 =
degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure();
const double penalty2 =
degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure();
const double penalty = penalty1 + penalty2;
const double h = dinfo1.face->measure();

for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double diff1 = uh1[k] - uh2[k];
const double diff2 =
fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];
dinfo1.value(0) +=
(penalty * diff1 * diff1 + h * diff2 * diff2) * fe.JxW(k);
}
dinfo1.value(0) = std::sqrt(dinfo1.value(0));
dinfo2.value(0) = dinfo1.value(0);
}

const FiniteElement< dim, spacedim > & get_fe() const
unsigned int tensor_degree() const
Triangulation< dim, spacedim >::face_iterator face
The current face.
Definition dof_info.h:81
number & value(const unsigned int i)

Finally we have an integrator for the error. Since the energy norm for discontinuous Galerkin problems not only involves the difference of the gradient inside the cells, but also the jump terms across faces and at the boundary, we cannot just use VectorTools::integrate_difference(). Instead, we use the MeshWorker interface to compute the error ourselves.

There are several different ways to define this energy norm, but all of them are equivalent to each other uniformly with mesh size (some not uniformly with polynomial degree). Here, we choose

$\|u\|_{1,h} = \sum_{K\in \mathbb T_h} \|\nabla u\|_K^2 + \sum_{F \in F_h^i} 4\sigma_F\|\average{ u \mathbf n}\|^2_F + \sum_{F \in F_h^b} 2\sigma_F\|u\|^2_F$

template <int dim>
class ErrorIntegrator : public MeshWorker::LocalIntegrator<dim>
{
public:
MeshWorker::IntegrationInfo<dim> &info) const override;
MeshWorker::IntegrationInfo<dim> &info) const override;
MeshWorker::IntegrationInfo<dim> &info2) const override;
};

Here we have the integration on cells. There is currently no good interface in MeshWorker that would allow us to access values of regular functions in the quadrature points. Thus, we have to create the vectors for the exact function's values and gradients inside the cell integrator. After that, everything is as before and we just add up the squares of the differences.

Additionally to computing the error in the energy norm, we use the capability of the mesh worker to compute two functionals at the same time and compute the L2-error in the same loop. Obviously, this one does not have any jump terms and only appears in the integration on cells.

template <int dim>
void ErrorIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();

const std::vector<Tensor<1, dim>> &Duh = info.gradients[0][0];
const std::vector<double> &uh = info.values[0][0];

for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
double sum = 0;
for (unsigned int d = 0; d < dim; ++d)
{
const double diff = exact_gradients[k][d] - Duh[k][d];
sum += diff * diff;
}
const double diff = exact_values[k] - uh[k];
dinfo.value(0) += sum * fe.JxW(k);
dinfo.value(1) += diff * diff * fe.JxW(k);
}
dinfo.value(0) = std::sqrt(dinfo.value(0));
dinfo.value(1) = std::sqrt(dinfo.value(1));
}

template <int dim>
void
ErrorIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();

const std::vector<double> &uh = info.values[0][0];

const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() /
dinfo.cell->measure();

for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double diff = exact_values[k] - uh[k];
dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
}
dinfo.value(0) = std::sqrt(dinfo.value(0));
}

template <int dim>
void ErrorIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,
{
const FEValuesBase<dim> &fe = info1.fe_values();
const std::vector<double> &uh1 = info1.values[0][0];
const std::vector<double> &uh2 = info2.values[0][0];

const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty1 =
degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure();
const double penalty2 =
degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure();
const double penalty = penalty1 + penalty2;

for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double diff = uh1[k] - uh2[k];
dinfo1.value(0) += (penalty * diff * diff) * fe.JxW(k);
}
dinfo1.value(0) = std::sqrt(dinfo1.value(0));
dinfo2.value(0) = dinfo1.value(0);
}

### The main class

This class does the main job, like in previous examples. For a description of the functions declared here, please refer to the implementation below.

template <int dim>
class InteriorPenaltyProblem
{
public:

InteriorPenaltyProblem();

void run(unsigned int n_steps);

private:
void setup_system();
void assemble_matrix();
void assemble_mg_matrix();
void assemble_right_hand_side();
void error();
double estimate();
void solve();
void output_results(const unsigned int cycle) const;

The member objects related to the discretization are here.

const MappingQ1<dim> mapping;
const FE_DGQ<2> fe;
DoFHandler<dim> dof_handler;

const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation

Then, we have the matrices and vectors related to the global discrete system.

Finally, we have a group of sparsity patterns and sparse matrices related to the multilevel preconditioner. First, we have a level matrix and its sparsity pattern.

When we perform multigrid with local smoothing on locally refined meshes, additional matrices are required; see Kanschat (2004). Here is the sparsity pattern for these edge matrices. We only need one, because the pattern of the up matrix is the transpose of that of the down matrix. Actually, we do not care too much about these details, since the MeshWorker is filling these matrices.

MGLevelObject<SparsityPattern> mg_sparsity_dg_interface;

The flux matrix at the refinement edge, coupling fine level degrees of freedom to coarse level.

MGLevelObject<SparseMatrix<double>> mg_matrix_dg_down;

The transpose of the flux matrix at the refinement edge, coupling coarse level degrees of freedom to fine level.

};

The constructor simply sets up the coarse grid and the DoFHandler.

template <int dim>
InteriorPenaltyProblem<dim>::InteriorPenaltyProblem()
: triangulation(Triangulation<dim>::limit_level_difference_at_vertices)
, mapping()
, fe(3)
, dof_handler(triangulation)
, estimates(1)
{
}

void hyper_cube_slit(Triangulation< dim > &tria, const double left=0., const double right=1., const bool colorize=false)

In this function, we set up the dimension of the linear system and the sparsity patterns for the global matrix as well as the level matrices.

template <int dim>
void InteriorPenaltyProblem<dim>::setup_system()
{

First, we use the finite element to distribute degrees of freedom over the mesh and number them.

dof_handler.distribute_dofs(fe);
dof_handler.distribute_mg_dofs();
unsigned int n_dofs = dof_handler.n_dofs();

Then, we already know the size of the vectors representing finite element functions.

solution.reinit(n_dofs);
right_hand_side.reinit(n_dofs);

Next, we set up the sparsity pattern for the global matrix. Since we do not know the row sizes in advance, we first fill a temporary DynamicSparsityPattern object and copy it to the regular SparsityPattern once it is complete.

DynamicSparsityPattern dsp(n_dofs);
sparsity.copy_from(dsp);
matrix.reinit(sparsity);

const unsigned int n_levels = triangulation.n_levels();
void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern)

The global system is set up, now we attend to the level matrices. We resize all matrix objects to hold one matrix per level.

mg_matrix.resize(0, n_levels - 1);
mg_matrix.clear_elements();
mg_matrix_dg_up.resize(0, n_levels - 1);
mg_matrix_dg_up.clear_elements();
mg_matrix_dg_down.resize(0, n_levels - 1);
mg_matrix_dg_down.clear_elements();

It is important to update the sparsity patterns after clear() was called for the level matrices, since the matrices lock the sparsity pattern through the SmartPointer and Subscriptor mechanism.

mg_sparsity.resize(0, n_levels - 1);
mg_sparsity_dg_interface.resize(0, n_levels - 1);

Now all objects are prepared to hold one sparsity pattern or matrix per level. What's left is setting up the sparsity patterns on each level.

for (unsigned int level = mg_sparsity.min_level();
level <= mg_sparsity.max_level();
++level)
{
unsigned int level
Definition grid_out.cc:4626

These are roughly the same lines as above for the global matrix, now for each level.

DynamicSparsityPattern dsp(dof_handler.n_dofs(level));
mg_sparsity[level].copy_from(dsp);
mg_matrix[level].reinit(mg_sparsity[level]);

void make_flux_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity, const unsigned int level, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true)
Definition mg_tools.cc:605

Additionally, we need to initialize the transfer matrices at the refinement edge between levels. They are stored at the index referring to the finer of the two indices, thus there is no such object on level 0.

if (level > 0)
{
dsp.reinit(dof_handler.n_dofs(level - 1),
dof_handler.n_dofs(level));
mg_sparsity_dg_interface[level].copy_from(dsp);
mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]);
mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]);
}
}
}

void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
void make_flux_sparsity_pattern_edge(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity, const unsigned int level)
Definition mg_tools.cc:679

In this function, we assemble the global system matrix, where by global we indicate that this is the matrix of the discrete system we solve and it is covering the whole mesh.

template <int dim>
void InteriorPenaltyProblem<dim>::assemble_matrix()
{

First, we need t set up the object providing the values we integrate. This object contains all FEValues and FEFaceValues objects needed and also maintains them automatically such that they always point to the current cell. To this end, we need to tell it first, where and what to compute. Since we are not doing anything fancy, we can rely on their standard choice for quadrature rules.

Since their default update flags are minimal, we add what we need additionally, namely the values and gradients of shape functions on all objects (cells, boundary and interior faces). Afterwards, we are ready to initialize the container, which will create all necessary FEValuesBase objects for integration.

info_box.initialize(fe, mapping);

UpdateFlags
@ update_values
Shape function values.

This is the object into which we integrate local data. It is filled by the local integration routines in MatrixIntegrator and then used by the assembler to distribute the information into the global matrix.

MeshWorker::DoFInfo<dim> dof_info(dof_handler);

Furthermore, we need an object that assembles the local matrix into the global matrix. These assembler objects have all the knowledge of the structures of the target object, in this case a SparseMatrix, possible constraints and the mesh structure.

assembler.initialize(matrix);

void initialize(MatrixType &m)
Definition simple.h:584

Now comes the part we coded ourselves, the local integrator. This is the only part which is problem dependent.

MatrixIntegrator<dim> integrator;

Now, we throw everything into a MeshWorker::loop(), which here traverses all active cells of the mesh, computes cell and face matrices and assembles them into the global matrix. We use the variable dof_handler here in order to use the global numbering of degrees of freedom.

MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
dof_handler.end(),
dof_info,
info_box,
integrator,
assembler);
}

Now, we do the same for the level matrices. Not too surprisingly, this function looks like a twin of the previous one. Indeed, there are only two minor differences.

template <int dim>
void InteriorPenaltyProblem<dim>::assemble_mg_matrix()
{
info_box.initialize(fe, mapping);

MeshWorker::DoFInfo<dim> dof_info(dof_handler);

Obviously, the assembler needs to be replaced by one filling level matrices. Note that it automatically fills the edge matrices as well.

assembler.initialize(mg_matrix);
assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);

MatrixIntegrator<dim> integrator;
void initialize(MGLevelObject< MatrixType > &m)
Definition simple.h:767

Here is the other difference to the previous function: we run over all cells, not only the active ones. And we use functions ending on _mg since we need the degrees of freedom on each level, not the global numbering.

MeshWorker::integration_loop<dim, dim>(dof_handler.begin_mg(),
dof_handler.end_mg(),
dof_info,
info_box,
integrator,
assembler);
}

Here we have another clone of the assemble function. The difference to assembling the system matrix consists in that we assemble a vector here.

template <int dim>
void InteriorPenaltyProblem<dim>::assemble_right_hand_side()
{
UpdateFlags update_flags =
info_box.initialize(fe, mapping);

MeshWorker::DoFInfo<dim> dof_info(dof_handler);

Since this assembler allows us to fill several vectors, the interface is a little more complicated as above. The pointers to the vectors have to be stored in an AnyData object. While this seems to cause two extra lines of code here, it actually comes handy in more complex applications.

AnyData data;
assembler.initialize(data);

RHSIntegrator<dim> integrator;
MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
dof_handler.end(),
dof_info,
info_box,
integrator,
assembler);

right_hand_side *= -1.;
}

void add(type entry, const std::string &name)
Definition any_data.h:430

Now that we have coded all functions building the discrete linear system, it is about time that we actually solve it.

template <int dim>
void InteriorPenaltyProblem<dim>::solve()
{

The solver of choice is conjugate gradient.

SolverControl control(1000, 1.e-12);
SolverCG<Vector<double>> solver(control);

Now we are setting up the components of the multilevel preconditioner. First, we need transfer between grid levels. The object we are using here generates sparse matrices for these transfers.

mg_transfer.build(dof_handler);

void build(const DoFHandler< dim, spacedim > &dof_handler)

Then, we need an exact solver for the matrix on the coarsest level.

FullMatrix<double> coarse_matrix;
coarse_matrix.copy_from(mg_matrix[0]);
mg_coarse.initialize(coarse_matrix);

void copy_from(const MatrixType &)
void initialize(const FullMatrix< number > &A)

While transfer and coarse grid solver are pretty much generic, more flexibility is offered for the smoother. First, we choose Gauss-Seidel as our smoothing method.

using RELAXATION = PreconditionSOR<SparseMatrix<double>>;
mg_smoother.initialize(mg_matrix, smoother_data);

Do two smoothing steps on each level.

mg_smoother.set_steps(2);

Since the SOR method is not symmetric, but we use conjugate gradient iteration below, here is a trick to make the multilevel preconditioner a symmetric operator even for nonsymmetric smoothers.

mg_smoother.set_symmetric(true);

The smoother class optionally implements the variable V-cycle, which we do not want here.

mg_smoother.set_variable(false);

Finally, we must wrap our matrices in an object having the required multiplication functions.

mg::Matrix<Vector<double>> mgmatrix(mg_matrix);
mg::Matrix<Vector<double>> mgdown(mg_matrix_dg_down);
mg::Matrix<Vector<double>> mgup(mg_matrix_dg_up);

Now, we are ready to set up the V-cycle operator and the multilevel preconditioner.

mgmatrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
Definition mg.h:81

Let us not forget the edge matrices needed because of the adaptive refinement.

mg.set_edge_flux_matrices(mgdown, mgup);

After all preparations, wrap the Multigrid object into another object, which can be used as a regular preconditioner,

preconditioner(dof_handler, mg, mg_transfer);

and use it to solve the system.

solver.solve(matrix, solution, right_hand_side, preconditioner);
}

Another clone of the assemble function. The big difference to the previous ones is here that we also have an input vector.

template <int dim>
double InteriorPenaltyProblem<dim>::estimate()
{

The results of the estimator are stored in a vector with one entry per cell. Since cells in deal.II are not numbered, we have to create our own numbering in order to use this vector. For the assembler used below the information in which component of a vector the result is stored is transmitted by the user_index variable for each cell. We need to set this numbering up here.

On the other hand, somebody might have used the user indices already. So, let's be good citizens and save them before tampering with them.

std::vector<unsigned int> old_user_indices;
triangulation.save_user_indices(old_user_indices);

estimates.block(0).reinit(triangulation.n_active_cells());
unsigned int i = 0;
for (const auto &cell : triangulation.active_cell_iterators())
cell->set_user_index(i++);

This starts like before,

const unsigned int n_gauss_points =
dof_handler.get_fe().tensor_degree() + 1;
n_gauss_points + 1,
n_gauss_points);

void initialize_gauss_quadrature(unsigned int n_cell_points, unsigned int n_boundary_points, unsigned int n_face_points, const bool force=true)

but now we need to notify the info box of the finite element function we want to evaluate in the quadrature points. First, we create an AnyData object with this vector, which is the solution we just computed.

AnyData solution_data;

Then, we tell the Meshworker::VectorSelector for cells, that we need the second derivatives of this solution (to compute the Laplacian). Therefore, the Boolean arguments selecting function values and first derivatives a false, only the last one selecting second derivatives is true.

On interior and boundary faces, we need the function values and the first derivatives, but not second derivatives.

And we continue as before, with the exception that the default update flags are already adjusted to the values and derivatives we requested above.

info_box.initialize(fe, mapping, solution_data, solution);

MeshWorker::DoFInfo<dim> dof_info(dof_handler);

The assembler stores one number per cell, but else this is the same as in the computation of the right hand side.

AnyData out_data;
assembler.initialize(out_data, false);

Estimator<dim> integrator;
MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
dof_handler.end(),
dof_info,
info_box,
integrator,
assembler);

Right before we return the result of the error estimate, we restore the old user indices.

return estimates.block(0).l2_norm();
}

BlockType & block(const unsigned int i)

Here we compare our finite element solution with the (known) exact solution and compute the mean quadratic error of the gradient and the function itself. This function is a clone of the estimation function right above.

Since we compute the error in the energy and the L2-norm, respectively, our block vector needs two blocks here.

template <int dim>
void InteriorPenaltyProblem<dim>::error()
{
BlockVector<double> errors(2);
errors.block(0).reinit(triangulation.n_active_cells());
errors.block(1).reinit(triangulation.n_active_cells());

std::vector<unsigned int> old_user_indices;
triangulation.save_user_indices(old_user_indices);
unsigned int i = 0;
for (const auto &cell : triangulation.active_cell_iterators())
cell->set_user_index(i++);

const unsigned int n_gauss_points =
dof_handler.get_fe().tensor_degree() + 1;
n_gauss_points + 1,
n_gauss_points);

AnyData solution_data;

info_box.initialize(fe, mapping, solution_data, solution);

MeshWorker::DoFInfo<dim> dof_info(dof_handler);

AnyData out_data;
assembler.initialize(out_data, false);

ErrorIntegrator<dim> integrator;
MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
dof_handler.end(),
dof_info,
info_box,
integrator,
assembler);

deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl;
deallog << "L2-error: " << errors.block(1).l2_norm() << std::endl;
}

void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
LogStream deallog
Definition logstream.cc:36

Create graphical output. We produce the filename by collating the name from its various components, including the refinement cycle that we output with two digits.

template <int dim>
void
InteriorPenaltyProblem<dim>::output_results(const unsigned int cycle) const
{
const std::string filename =
"sol-" + Utilities::int_to_string(cycle, 2) + ".gnuplot";

deallog << "Writing solution to <" << filename << ">..." << std::endl
<< std::endl;
std::ofstream gnuplot_output(filename);

DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);

data_out.build_patches();

data_out.write_gnuplot(gnuplot_output);
}

void attach_dof_handler(const DoFHandler< dim, spacedim > &)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470

And finally the adaptive loop, more or less like in previous examples.

template <int dim>
void InteriorPenaltyProblem<dim>::run(unsigned int n_steps)
{
deallog << "Element: " << fe.get_name() << std::endl;
for (unsigned int s = 0; s < n_steps; ++s)
{
deallog << "Step " << s << std::endl;
if (estimates.block(0).size() == 0)
triangulation.refine_global(1);
else
{
triangulation, estimates.block(0), 0.5, 0.0);
triangulation.execute_coarsening_and_refinement();
}

deallog << "Triangulation " << triangulation.n_active_cells()
<< " cells, " << triangulation.n_levels() << " levels"
<< std::endl;

setup_system();
deallog << "DoFHandler " << dof_handler.n_dofs() << " dofs, level dofs";
for (unsigned int l = 0; l < triangulation.n_levels(); ++l)
deallog << ' ' << dof_handler.n_dofs(l);
deallog << std::endl;

deallog << "Assemble matrix" << std::endl;
assemble_matrix();
deallog << "Assemble multilevel matrix" << std::endl;
assemble_mg_matrix();
deallog << "Assemble right hand side" << std::endl;
assemble_right_hand_side();
deallog << "Solve" << std::endl;
solve();
error();
deallog << "Estimate " << estimate() << std::endl;
output_results(s);
}
}
} // namespace Step39

int main()
{
try
{
using namespace dealii;
using namespace Step39;

std::ofstream logfile("deallog");
deallog.attach(logfile);

InteriorPenaltyProblem<2> test1;
test1.run(12);
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}

return 0;
}
void attach(std::ostream &o, const bool print_job_id=true, const std::ios_base::fmtflags flags=std::ios::showpoint|std::ios::left)
Definition logstream.cc:218
unsigned int depth_console(const unsigned int n)
Definition logstream.cc:349
void refine_and_coarsen_fixed_fraction(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double top_fraction, const double bottom_fraction, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max(), const VectorTools::NormType norm_type=VectorTools::L1_norm)

# Results

### Logfile output

First, the program produces the usual logfile here stored in deallog. It reads (with omission of intermediate steps)

DEAL::Element: FE_DGQ<2>(3)
DEAL::Step 0
DEAL::Triangulation 16 cells, 2 levels
DEAL::DoFHandler 256 dofs, level dofs 64 256
DEAL::Assemble matrix
DEAL::Assemble multilevel matrix
DEAL::Assemble right hand side
DEAL::Solve
DEAL:cg::Starting value 37.4071
DEAL:cg::Convergence step 13 value 1.64974e-13
DEAL::energy-error: 0.297419
DEAL::L2-error: 0.00452447
DEAL::Estimate 0.990460
DEAL::Writing solution to <sol-00.gnuplot>...
DEAL::
DEAL::Step 1
DEAL::Triangulation 25 cells, 3 levels
DEAL::DoFHandler 400 dofs, level dofs 64 256 192
DEAL::Assemble matrix
DEAL::Assemble multilevel matrix
DEAL::Assemble right hand side
DEAL::Solve
DEAL:cg::Starting value 37.4071
DEAL:cg::Convergence step 14 value 3.72262e-13
DEAL::energy-error: 0.258559
DEAL::L2-error: 0.00288510
DEAL::Estimate 0.738624
DEAL::Writing solution to <sol-01.gnuplot>...
DEAL::
DEAL::Step 2
DEAL::Triangulation 34 cells, 4 levels
DEAL::DoFHandler 544 dofs, level dofs 64 256 256 128
DEAL::Assemble matrix
DEAL::Assemble multilevel matrix
DEAL::Assemble right hand side
DEAL::Solve
DEAL:cg::Starting value 37.4071
DEAL:cg::Convergence step 15 value 1.91610e-13
DEAL::energy-error: 0.189234
DEAL::L2-error: 0.00147954
DEAL::Estimate 0.657507
DEAL::Writing solution to <sol-02.gnuplot>...
...
DEAL::Step 10
DEAL::Triangulation 232 cells, 11 levels
DEAL::DoFHandler 3712 dofs, level dofs 64 256 896 768 768 640 512 256 256 256 256
DEAL::Assemble matrix
DEAL::Assemble multilevel matrix
DEAL::Assemble right hand side
DEAL::Solve
DEAL:cg::Starting value 51.1571
DEAL:cg::Convergence step 15 value 7.19599e-13
DEAL::energy-error: 0.0132475
DEAL::L2-error: 1.00423e-05
DEAL::Estimate 0.0470724
DEAL::Writing solution to <sol-10.gnuplot>...
DEAL::
DEAL::Step 11
DEAL::Triangulation 322 cells, 12 levels
DEAL::DoFHandler 5152 dofs, level dofs 64 256 1024 1024 896 768 768 640 448 320 320 320
DEAL::Assemble matrix
DEAL::Assemble multilevel matrix
DEAL::Assemble right hand side
DEAL::Solve
DEAL:cg::Starting value 52.2226
DEAL:cg::Convergence step 15 value 8.15195e-13
DEAL::energy-error: 0.00934891
DEAL::L2-error: 5.41095e-06
DEAL::Estimate 0.0329102
DEAL::Writing solution to <sol-11.gnuplot>...
DEAL::

This log for instance shows that the number of conjugate gradient iteration steps is constant at approximately 15.

### Postprocessing of the logfile

Using the perl script postprocess.pl, we extract relevant data into output.dat, which can be used to plot graphs with gnuplot. The graph above for instance was produced using the gnuplot script plot_errors.gpl via

./step-39 | perl postprocess.pl >output.dat
gnuplot plot_errors.gpl

Reference data can be found in output.reference.dat.

# The plain program

/* ------------------------------------------------------------------------
*
* Copyright (C) 2010 - 2023 by the deal.II authors
*
* This file is part of the deal.II library.
*
* Part of the source code is dual licensed under Apache-2.0 WITH
* LLVM-exception OR LGPL-2.1-or-later. Detailed license information
* governing the source code and code contributions can be found in
* LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
*
* ------------------------------------------------------------------------
*
* Author: Guido Kanschat, Texas A&M University, 2009
*/
#include <iostream>
#include <fstream>
namespace Step39
{
using namespace dealii;
template <int dim>
class MatrixIntegrator : public MeshWorker::LocalIntegrator<dim>
{
public:
void cell(MeshWorker::DoFInfo<dim> &dinfo,
MeshWorker::IntegrationInfo<dim> &info) const override;
void boundary(MeshWorker::DoFInfo<dim> &dinfo,
MeshWorker::IntegrationInfo<dim> &info) const override;
void face(MeshWorker::DoFInfo<dim> &dinfo1,
MeshWorker::IntegrationInfo<dim> &info2) const override;
};
template <int dim>
double ip_penalty_factor(const MeshWorker::DoFInfo<dim> &dinfo1,
const MeshWorker::DoFInfo<dim> &dinfo2,
unsigned int deg1,
unsigned int deg2)
{
const unsigned int normal1 =
const unsigned int normal2 =
const unsigned int deg1sq = (deg1 == 0) ? 1 : deg1 * (deg1 + 1);
const unsigned int deg2sq = (deg2 == 0) ? 1 : deg2 * (deg2 + 1);
double penalty1 = deg1sq / dinfo1.cell->extent_in_direction(normal1);
double penalty2 = deg2sq / dinfo2.cell->extent_in_direction(normal2);
if (dinfo1.cell->has_children() && !dinfo2.cell->has_children())
penalty1 *= 2;
else if (!dinfo1.cell->has_children() && dinfo2.cell->has_children())
penalty2 *= 2;
const double penalty = 0.5 * (penalty1 + penalty2);
return penalty;
}
template <int dim>
void MatrixIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo,
{
FullMatrix<double> &M = dinfo.matrix(0, false).matrix;
for (unsigned int k = 0; k < info.fe_values().n_quadrature_points; ++k)
{
const double dx = info.fe_values().JxW(k);
for (unsigned int i = 0; i < info.fe_values().dofs_per_cell; ++i)
{
const double Mii = (info.fe_values().shape_grad(i, k) *
M(i, i) += Mii;
for (unsigned int j = i + 1; j < info.fe_values().dofs_per_cell;
++j)
{
const double Mij = info.fe_values().shape_grad(j, k) *
M(i, j) += Mij;
M(j, i) += Mij;
}
}
}
}
template <int dim>
void
MatrixIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe_face_values = info.fe_values(0);
FullMatrix<double> &M = dinfo.matrix(0, false).matrix;
AssertDimension(M.n(), fe_face_values.dofs_per_cell);
AssertDimension(M.m(), fe_face_values.dofs_per_cell);
const unsigned int polynomial_degree =
info.fe_values(0).get_fe().tensor_degree();
const double ip_penalty =
ip_penalty_factor(dinfo, dinfo, polynomial_degree, polynomial_degree);
for (unsigned int k = 0; k < fe_face_values.n_quadrature_points; ++k)
{
const double dx = fe_face_values.JxW(k);
const Tensor<1, dim> &n = fe_face_values.normal_vector(k);
for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i)
for (unsigned int j = 0; j < fe_face_values.dofs_per_cell; ++j)
M(i, j) += (2. * fe_face_values.shape_value(i, k) * ip_penalty *
fe_face_values.shape_value(j, k) -
fe_face_values.shape_value(j, k) -
fe_face_values.shape_value(i, k)) *
dx;
}
}
template <int dim>
void
MatrixIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,
{
const FEValuesBase<dim> &fe_face_values_1 = info1.fe_values(0);
const FEValuesBase<dim> &fe_face_values_2 = info2.fe_values(0);
FullMatrix<double> &M11 = dinfo1.matrix(0, false).matrix;
FullMatrix<double> &M12 = dinfo1.matrix(0, true).matrix;
FullMatrix<double> &M21 = dinfo2.matrix(0, true).matrix;
FullMatrix<double> &M22 = dinfo2.matrix(0, false).matrix;
AssertDimension(M11.n(), fe_face_values_1.dofs_per_cell);
AssertDimension(M11.m(), fe_face_values_1.dofs_per_cell);
AssertDimension(M12.n(), fe_face_values_1.dofs_per_cell);
AssertDimension(M12.m(), fe_face_values_1.dofs_per_cell);
AssertDimension(M21.n(), fe_face_values_1.dofs_per_cell);
AssertDimension(M21.m(), fe_face_values_1.dofs_per_cell);
AssertDimension(M22.n(), fe_face_values_1.dofs_per_cell);
AssertDimension(M22.m(), fe_face_values_1.dofs_per_cell);
const unsigned int polynomial_degree =
info1.fe_values(0).get_fe().tensor_degree();
const double ip_penalty =
ip_penalty_factor(dinfo1, dinfo2, polynomial_degree, polynomial_degree);
const double nui = 1.;
const double nue = 1.;
const double nu = .5 * (nui + nue);
for (unsigned int k = 0; k < fe_face_values_1.n_quadrature_points; ++k)
{
const double dx = fe_face_values_1.JxW(k);
const Tensor<1, dim> &n = fe_face_values_1.normal_vector(k);
for (unsigned int i = 0; i < fe_face_values_1.dofs_per_cell; ++i)
{
for (unsigned int j = 0; j < fe_face_values_1.dofs_per_cell; ++j)
{
const double vi = fe_face_values_1.shape_value(i, k);
const double dnvi = n * fe_face_values_1.shape_grad(i, k);
const double ve = fe_face_values_2.shape_value(i, k);
const double dnve = n * fe_face_values_2.shape_grad(i, k);
const double ui = fe_face_values_1.shape_value(j, k);
const double dnui = n * fe_face_values_1.shape_grad(j, k);
const double ue = fe_face_values_2.shape_value(j, k);
const double dnue = n * fe_face_values_2.shape_grad(j, k);
M11(i, j) += (-.5 * nui * dnvi * ui - .5 * nui * dnui * vi +
nu * ip_penalty * ui * vi) *
dx;
M12(i, j) += (.5 * nui * dnvi * ue - .5 * nue * dnue * vi -
nu * ip_penalty * vi * ue) *
dx;
M21(i, j) += (-.5 * nue * dnve * ui + .5 * nui * dnui * ve -
nu * ip_penalty * ui * ve) *
dx;
M22(i, j) += (.5 * nue * dnve * ue + .5 * nue * dnue * ve +
nu * ip_penalty * ue * ve) *
dx;
}
}
}
}
template <int dim>
class RHSIntegrator : public MeshWorker::LocalIntegrator<dim>
{
public:
void cell(MeshWorker::DoFInfo<dim> &dinfo,
MeshWorker::IntegrationInfo<dim> &info) const override;
void boundary(MeshWorker::DoFInfo<dim> &dinfo,
MeshWorker::IntegrationInfo<dim> &info) const override;
void face(MeshWorker::DoFInfo<dim> &dinfo1,
MeshWorker::IntegrationInfo<dim> &info2) const override;
};
template <int dim>
void RHSIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &,
{}
template <int dim>
void
RHSIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();
Vector<double> &local_vector = dinfo.vector(0).block(0);
const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() /
dinfo.cell->measure();
for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
for (unsigned int i = 0; i < fe.dofs_per_cell; ++i)
local_vector(i) +=
(-penalty * fe.shape_value(i, k) // (-sigma * v_i(x_k)
* boundary_values[k] * fe.JxW(k); // u^D(x_k) * dx
}
template <int dim>
void RHSIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &,
{}
template <int dim>
class Estimator : public MeshWorker::LocalIntegrator<dim>
{
public:
void cell(MeshWorker::DoFInfo<dim> &dinfo,
MeshWorker::IntegrationInfo<dim> &info) const override;
void boundary(MeshWorker::DoFInfo<dim> &dinfo,
MeshWorker::IntegrationInfo<dim> &info) const override;
void face(MeshWorker::DoFInfo<dim> &dinfo1,
MeshWorker::IntegrationInfo<dim> &info2) const override;
};
template <int dim>
void Estimator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();
const std::vector<Tensor<2, dim>> &DDuh = info.hessians[0][0];
for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double t = dinfo.cell->diameter() * trace(DDuh[k]);
dinfo.value(0) += t * t * fe.JxW(k);
}
dinfo.value(0) = std::sqrt(dinfo.value(0));
}
template <int dim>
void Estimator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();
const std::vector<double> &uh = info.values[0][0];
const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() /
dinfo.cell->measure();
for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double diff = boundary_values[k] - uh[k];
dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
}
dinfo.value(0) = std::sqrt(dinfo.value(0));
}
template <int dim>
void Estimator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,
{
const FEValuesBase<dim> &fe = info1.fe_values();
const std::vector<double> &uh1 = info1.values[0][0];
const std::vector<double> &uh2 = info2.values[0][0];
const std::vector<Tensor<1, dim>> &Duh1 = info1.gradients[0][0];
const std::vector<Tensor<1, dim>> &Duh2 = info2.gradients[0][0];
const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty1 =
degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure();
const double penalty2 =
degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure();
const double penalty = penalty1 + penalty2;
const double h = dinfo1.face->measure();
for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double diff1 = uh1[k] - uh2[k];
const double diff2 =
fe.normal_vector(k) * Duh1[k] - fe.normal_vector(k) * Duh2[k];
dinfo1.value(0) +=
(penalty * diff1 * diff1 + h * diff2 * diff2) * fe.JxW(k);
}
dinfo1.value(0) = std::sqrt(dinfo1.value(0));
dinfo2.value(0) = dinfo1.value(0);
}
template <int dim>
class ErrorIntegrator : public MeshWorker::LocalIntegrator<dim>
{
public:
void cell(MeshWorker::DoFInfo<dim> &dinfo,
MeshWorker::IntegrationInfo<dim> &info) const override;
void boundary(MeshWorker::DoFInfo<dim> &dinfo,
MeshWorker::IntegrationInfo<dim> &info) const override;
void face(MeshWorker::DoFInfo<dim> &dinfo1,
MeshWorker::IntegrationInfo<dim> &info2) const override;
};
template <int dim>
void ErrorIntegrator<dim>::cell(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();
const std::vector<Tensor<1, dim>> &Duh = info.gradients[0][0];
const std::vector<double> &uh = info.values[0][0];
for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
double sum = 0;
for (unsigned int d = 0; d < dim; ++d)
{
const double diff = exact_gradients[k][d] - Duh[k][d];
sum += diff * diff;
}
const double diff = exact_values[k] - uh[k];
dinfo.value(0) += sum * fe.JxW(k);
dinfo.value(1) += diff * diff * fe.JxW(k);
}
dinfo.value(0) = std::sqrt(dinfo.value(0));
dinfo.value(1) = std::sqrt(dinfo.value(1));
}
template <int dim>
void
ErrorIntegrator<dim>::boundary(MeshWorker::DoFInfo<dim> &dinfo,
{
const FEValuesBase<dim> &fe = info.fe_values();
const std::vector<double> &uh = info.values[0][0];
const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty = 2. * degree * (degree + 1) * dinfo.face->measure() /
dinfo.cell->measure();
for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double diff = exact_values[k] - uh[k];
dinfo.value(0) += penalty * diff * diff * fe.JxW(k);
}
dinfo.value(0) = std::sqrt(dinfo.value(0));
}
template <int dim>
void ErrorIntegrator<dim>::face(MeshWorker::DoFInfo<dim> &dinfo1,
{
const FEValuesBase<dim> &fe = info1.fe_values();
const std::vector<double> &uh1 = info1.values[0][0];
const std::vector<double> &uh2 = info2.values[0][0];
const unsigned int degree = fe.get_fe().tensor_degree();
const double penalty1 =
degree * (degree + 1) * dinfo1.face->measure() / dinfo1.cell->measure();
const double penalty2 =
degree * (degree + 1) * dinfo2.face->measure() / dinfo2.cell->measure();
const double penalty = penalty1 + penalty2;
for (unsigned k = 0; k < fe.n_quadrature_points; ++k)
{
const double diff = uh1[k] - uh2[k];
dinfo1.value(0) += (penalty * diff * diff) * fe.JxW(k);
}
dinfo1.value(0) = std::sqrt(dinfo1.value(0));
dinfo2.value(0) = dinfo1.value(0);
}
template <int dim>
class InteriorPenaltyProblem
{
public:
InteriorPenaltyProblem();
void run(unsigned int n_steps);
private:
void setup_system();
void assemble_matrix();
void assemble_mg_matrix();
void assemble_right_hand_side();
void error();
double estimate();
void solve();
void output_results(const unsigned int cycle) const;
const FE_DGQ<2> fe;
DoFHandler<dim> dof_handler;
SparsityPattern sparsity;
Vector<double> solution;
Vector<double> right_hand_side;
MGLevelObject<SparsityPattern> mg_sparsity_dg_interface;
};
template <int dim>
InteriorPenaltyProblem<dim>::InteriorPenaltyProblem()
: triangulation(Triangulation<dim>::limit_level_difference_at_vertices)
, mapping()
, fe(3)
, dof_handler(triangulation)
, estimates(1)
{
}
template <int dim>
void InteriorPenaltyProblem<dim>::setup_system()
{
dof_handler.distribute_dofs(fe);
dof_handler.distribute_mg_dofs();
unsigned int n_dofs = dof_handler.n_dofs();
solution.reinit(n_dofs);
right_hand_side.reinit(n_dofs);
sparsity.copy_from(dsp);
matrix.reinit(sparsity);
const unsigned int n_levels = triangulation.n_levels();
mg_matrix.resize(0, n_levels - 1);
mg_matrix.clear_elements();
mg_matrix_dg_up.resize(0, n_levels - 1);
mg_matrix_dg_up.clear_elements();
mg_matrix_dg_down.resize(0, n_levels - 1);
mg_matrix_dg_down.clear_elements();
mg_sparsity.resize(0, n_levels - 1);
mg_sparsity_dg_interface.resize(0, n_levels - 1);
for (unsigned int level = mg_sparsity.min_level();
level <= mg_sparsity.max_level();
++level)
{
DynamicSparsityPattern dsp(dof_handler.n_dofs(level));
mg_sparsity[level].copy_from(dsp);
mg_matrix[level].reinit(mg_sparsity[level]);
if (level > 0)
{
dsp.reinit(dof_handler.n_dofs(level - 1),
dof_handler.n_dofs(level));
mg_sparsity_dg_interface[level].copy_from(dsp);
mg_matrix_dg_up[level].reinit(mg_sparsity_dg_interface[level]);
mg_matrix_dg_down[level].reinit(mg_sparsity_dg_interface[level]);
}
}
}
template <int dim>
void InteriorPenaltyProblem<dim>::assemble_matrix()
{
info_box.initialize(fe, mapping);
MeshWorker::DoFInfo<dim> dof_info(dof_handler);
assembler.initialize(matrix);
MatrixIntegrator<dim> integrator;
MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
dof_handler.end(),
dof_info,
info_box,
integrator,
assembler);
}
template <int dim>
void InteriorPenaltyProblem<dim>::assemble_mg_matrix()
{
info_box.initialize(fe, mapping);
MeshWorker::DoFInfo<dim> dof_info(dof_handler);
assembler.initialize(mg_matrix);
assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down);
MatrixIntegrator<dim> integrator;
MeshWorker::integration_loop<dim, dim>(dof_handler.begin_mg(),
dof_handler.end_mg(),
dof_info,
info_box,
integrator,
assembler);
}
template <int dim>
void InteriorPenaltyProblem<dim>::assemble_right_hand_side()
{
UpdateFlags update_flags =
info_box.initialize(fe, mapping);
MeshWorker::DoFInfo<dim> dof_info(dof_handler);
AnyData data;
assembler.initialize(data);
RHSIntegrator<dim> integrator;
MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
dof_handler.end(),
dof_info,
info_box,
integrator,
assembler);
right_hand_side *= -1.;
}
template <int dim>
void InteriorPenaltyProblem<dim>::solve()
{
SolverControl control(1000, 1.e-12);
SolverCG<Vector<double>> solver(control);
mg_transfer.build(dof_handler);
FullMatrix<double> coarse_matrix;
coarse_matrix.copy_from(mg_matrix[0]);
mg_coarse.initialize(coarse_matrix);
mg_smoother.initialize(mg_matrix, smoother_data);
mg_smoother.set_steps(2);
mg_smoother.set_symmetric(true);
mg_smoother.set_variable(false);
mg::Matrix<Vector<double>> mgmatrix(mg_matrix);
mg::Matrix<Vector<double>> mgdown(mg_matrix_dg_down);
mg::Matrix<Vector<double>> mgup(mg_matrix_dg_up);
mgmatrix, mg_coarse, mg_transfer, mg_smoother, mg_smoother);
mg.set_edge_flux_matrices(mgdown, mgup);
preconditioner(dof_handler, mg, mg_transfer);
solver.solve(matrix, solution, right_hand_side, preconditioner);
}
template <int dim>
double InteriorPenaltyProblem<dim>::estimate()
{
std::vector<unsigned int> old_user_indices;
triangulation.save_user_indices(old_user_indices);
estimates.block(0).reinit(triangulation.n_active_cells());
unsigned int i = 0;
for (const auto &cell : triangulation.active_cell_iterators())
cell->set_user_index(i++);
const unsigned int n_gauss_points =
dof_handler.get_fe().tensor_degree() + 1;
n_gauss_points + 1,
n_gauss_points);
AnyData solution_data;
info_box.initialize(fe, mapping, solution_data, solution);
MeshWorker::DoFInfo<dim> dof_info(dof_handler);
AnyData out_data;
assembler.initialize(out_data, false);
Estimator<dim> integrator;
MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
dof_handler.end(),
dof_info,
info_box,
integrator,
assembler);
return estimates.block(0).l2_norm();
}
template <int dim>
void InteriorPenaltyProblem<dim>::error()
{
errors.block(0).reinit(triangulation.n_active_cells());
errors.block(1).reinit(triangulation.n_active_cells());
std::vector<unsigned int> old_user_indices;
triangulation.save_user_indices(old_user_indices);
unsigned int i = 0;
for (const auto &cell : triangulation.active_cell_iterators())
cell->set_user_index(i++);
const unsigned int n_gauss_points =
dof_handler.get_fe().tensor_degree() + 1;
n_gauss_points + 1,
n_gauss_points);
AnyData solution_data;
info_box.initialize(fe, mapping, solution_data, solution);
MeshWorker::DoFInfo<dim> dof_info(dof_handler);
AnyData out_data;
assembler.initialize(out_data, false);
ErrorIntegrator<dim> integrator;
MeshWorker::integration_loop<dim, dim>(dof_handler.begin_active(),
dof_handler.end(),
dof_info,
info_box,
integrator,
assembler);
deallog << "energy-error: " << errors.block(0).l2_norm() << std::endl;
deallog << "L2-error: " << errors.block(1).l2_norm() << std::endl;
}
template <int dim>
void
InteriorPenaltyProblem<dim>::output_results(const unsigned int cycle) const
{
const std::string filename =
"sol-" + Utilities::int_to_string(cycle, 2) + ".gnuplot";
deallog << "Writing solution to <" << filename << ">..." << std::endl
<< std::endl;
std::ofstream gnuplot_output(filename);
DataOut<dim> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.build_patches();
data_out.write_gnuplot(gnuplot_output);
}
template <int dim>
void InteriorPenaltyProblem<dim>::run(unsigned int n_steps)
{
deallog << "Element: " << fe.get_name() << std::endl;
for (unsigned int s = 0; s < n_steps; ++s)
{
deallog << "Step " << s << std::endl;
if (estimates.block(0).size() == 0)
triangulation.refine_global(1);
else
{
triangulation, estimates.block(0), 0.5, 0.0);
triangulation.execute_coarsening_and_refinement();
}
deallog << "Triangulation " << triangulation.n_active_cells()
<< " cells, " << triangulation.n_levels() << " levels"
<< std::endl;
setup_system();
deallog << "DoFHandler " << dof_handler.n_dofs() << " dofs, level dofs";
for (unsigned int l = 0; l < triangulation.n_levels(); ++l)
deallog << ' ' << dof_handler.n_dofs(l);
deallog << std::endl;
deallog << "Assemble matrix" << std::endl;
assemble_matrix();
deallog << "Assemble multilevel matrix" << std::endl;
assemble_mg_matrix();
deallog << "Assemble right hand side" << std::endl;
assemble_right_hand_side();
deallog << "Solve" << std::endl;
solve();
error();
deallog << "Estimate " << estimate() << std::endl;
output_results(s);
}
}
} // namespace Step39
int main()
{
try
{
using namespace dealii;
using namespace Step39;
std::ofstream logfile("deallog");
deallog.attach(logfile);
InteriorPenaltyProblem<2> test1;
test1.run(12);
}
catch (std::exception &exc)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl
<< std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
void write_gnuplot(std::ostream &out) const
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition data_out.cc:1062
const std::vector< Point< spacedim > > & get_quadrature_points() const
const unsigned int dofs_per_cell
const Tensor< 1, spacedim > & normal_vector(const unsigned int q_point) const
const Tensor< 1, spacedim > & shape_grad(const unsigned int i, const unsigned int q_point) const
double JxW(const unsigned int q_point) const
const double & shape_value(const unsigned int i, const unsigned int q_point) const
size_type n() const
size_type m() const
virtual void gradient_list(const std::vector< Point< dim > > &points, std::vector< Tensor< 1, dim > > &gradients, const unsigned int component=0) const override
virtual void value_list(const std::vector< Point< dim > > &points, std::vector< double > &values, const unsigned int component=0) const override
void set_steps(const unsigned int)
void set_symmetric(const bool)
void set_variable(const bool)
void initialize(AnyData &results, bool separate_faces=true)
Definition functional.h:230
void initialize_fluxes(MGLevelObject< MatrixType > &flux_up, MGLevelObject< MatrixType > &flux_down)
Definition simple.h:782
void initialize(AnyData &results)
Definition simple.h:512
void initialize(const FiniteElement< dim, spacedim > &el, const Mapping< dim, spacedim > &mapping, const BlockInfo *block_info=nullptr)
MeshWorker::VectorSelector boundary_selector
MeshWorker::VectorSelector face_selector
MeshWorker::VectorSelector cell_selector
std::vector< std::vector< std::vector< Tensor< 2, spacedim > > > > hessians
const FEValuesBase< dim, spacedim > & fe_values() const