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 MeshWorker::mesh_loop() with scratch and copy objects.
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
\]
\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
In our finite element program, we distinguish three different integrals, corresponding to the sums over cells, interior faces and boundary faces above. Since MeshWorker::mesh_loop() organizes these traversals for us, we only need to implement the integrals over each mesh element and store the local results in copy objects before they are transferred to the global data structures. The class MatrixIntegrator below has these three functions for the left hand side of the formula, the class RHSIntegrator for the right. This mirrors the WorkStream scratch/copy pattern already introduced in step-9, except that MeshWorker::mesh_loop() also manages the boundary and interior-face traversals needed by this discontinuous Galerkin formulation.
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.
Finally, we take our exact solution from the library as well as quadrature and additional tools.
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.
This is the function we use to set the boundary values and also the exact solution we compare to.
This organization is what makes parallel assembly practical: the worker stage must not touch shared global data directly, and it should also avoid repeatedly allocating expensive temporary objects. Consequently, each worker receives a scratch object that owns reusable FEValues-like data and other temporary arrays, together with a copy object that stores the local contributions produced on the current cell or face. If you would like to see the general WorkStream idea in a simpler setting, step-9 gives a more introductory discussion before we apply the same pattern here to cells, faces, and subfaces.
We start with a class that stores scratch data for assembling the global and multigrid matrices. This is the most elaborate scratch object because matrix assembly needs FEValues on cells, on boundary faces, on regular interior faces, and on subfaces at refinement edges.
Next comes the scratch object used for assembling the right-hand side. Here we only need access to boundary faces because the inhomogeneous terms come from the Nitsche boundary contributions.
The error estimator, in turn, needs the following scratch class. Besides FEValues-like objects on the relevant geometric entities, this object stores reusable buffers for solution values, gradients, Hessians, and exact boundary values so that the worker lambdas do not allocate these arrays repeatedly.
To compute the actual error norms, we use another scratch class. As above, we keep both FEValues-like objects and temporary storage for function values and gradients in one per-thread object that can be reused for many cells and faces.
The first copy-data structure represents the local data associated with one face contribution. These data are kept separately because each interior face contributes four blocks coupling the two adjacent cells.
Matrix assembly then uses the following copy-data class. Each worker first fills the cell matrix and, if necessary, appends additional face contributions; the copier then transfers all of this local data into the global sparse matrices.
For the right-hand side, a simpler copy-data class is sufficient. In this case the local result is just one cell vector together with the corresponding global DoF indices.
Finally, the estimator and error computations share the following copy-data class. Instead of matrices or vectors, the local results are a small fixed number of scalar indicators attached to the current cell and to the faces touching it. The template argument is the number of such indicators, which is known at compile time. Specifically, we will use this class below with n_values=1 for the error estimator, and with n_values=2 for the error computation (where we compute both the \(L_2\) and \(H^1\) errors).
The first namespace defining local integrators is responsible for assembling the global matrix as well as the level matrices. 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.
The second set of local integrators 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.
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). The cell contribution is the Laplacian of the discrete solution, since the right hand side is zero.
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.
Finally, on interior faces, the estimator consists of the jumps of the solution and its normal derivative, weighted appropriately.
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 MeshWorker::mesh_loop() 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 \]
Below, the first function is, as always, the integration on cells. The exact solution is evaluated directly in the quadrature points and then compared against the discrete solution.
This class does the main job, like in previous examples. For a description of the functions declared here, please refer to the implementation below.
The member objects related to the discretization are here.
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. These matrices are filled in the multigrid mesh loop below.
The flux matrix at the refinement edge, coupling fine level degrees of freedom to coarse level.
The transpose of the flux matrix at the refinement edge, coupling coarse level degrees of freedom to fine level.
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.
First, we use the finite element to distribute degrees of freedom over the mesh and number them.
Then, we already know the size of the vectors representing finite element functions.
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.
The global system is set up, now we attend to the level matrices. We resize all matrix objects to hold one matrix per level.
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.
These are roughly the same lines as above for the global matrix, now for each level.
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.
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.
None of these workers writes into the global sparse matrix directly. Instead, they only populate the copy object associated with the current cell. Once the local work is complete, the copier takes the data stored in that copy object and inserts it into the global matrix. This separation is what allows MeshWorker::mesh_loop() to parallelize the local integration safely.
Now, we do the same for the level matrices. Not too surprisingly, this function looks like a twin of the previous one. The mesh loop again traverses cells, boundary faces, and interior faces, and the worker functions compute local matrix contributions without touching global data. The scratch object provides the FEValues-like data needed for those local computations, and the copy object collects the local cell matrix together with the face blocks.
The main difference lies in the copier. Rather than assembling into a single global matrix, it dispatches the local data into the appropriate level matrix and, on refinement edges, into the up- and down-transfer matrices that are required by the multigrid algorithm.
Here we have another clone of the assembly function. The difference to assembling the system matrix consists in that we assemble a vector here.
The mesh loop still uses the same worker/copier split. The cell worker is only responsible for initializing the copy object for the current cell, whereas the actual local work happens on boundary faces: there, the worker evaluates the inhomogeneous boundary values and accumulates the associated Nitsche terms into the local right-hand-side vector. The copier then adds that local vector to the global right-hand side.
Now that we have coded all functions building the discrete linear system, it is about time that we actually solve it.
The solver of choice is conjugate gradient.
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.
Then, we need an exact solver for the matrix on the coarsest level.
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.
Do two smoothing steps on each level.
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.
The smoother class optionally implements the variable V-cycle, which we do not want here.
Finally, we must wrap our matrices in an object having the required multiplication functions.
Now, we are ready to set up the V-cycle operator and the multilevel preconditioner.
Let us not forget the edge matrices needed because of the adaptive refinement.
and use it to solve the system.
The next function estimates the error. The big difference to the previous mesh loop functions is that we now also read from the discrete solution vector. The results of the estimator are stored in a vector with one entry per cell.
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 close relative of the estimation function right above: the mesh loop again visits cells, boundary faces, and interior faces; the workers evaluate local quantities with the help of the scratch object; and the copier writes the resulting indicators into global data structures only after the local computation is finished.
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.
And finally the adaptive loop, more or less like in previous examples.
This log for instance shows that the number of conjugate gradient iteration steps is constant at approximately 15. This is the key qualitative result of the example: despite adaptive mesh refinement and the discontinuous Galerkin discretization, the multigrid preconditioner keeps the iteration count essentially mesh-independent.
#include <array>
#include <fstream>
#include <iostream>
namespace Step39
{
template <int dim>
struct MatrixScratchData
{
: fe_values(mapping, fe, cell_quadrature, cell_update_flags)
, boundary_fe_values(mapping, fe, face_quadrature, face_update_flags)
, face_fe_values(mapping, fe, face_quadrature, face_update_flags)
, subface_values(mapping, fe, face_quadrature, face_update_flags)
, neighbor_face_values(mapping, fe, face_quadrature, face_update_flags)
, neighbor_subface_values(mapping, fe, face_quadrature, face_update_flags)
{}
MatrixScratchData(const MatrixScratchData<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())
, boundary_fe_values(scratch_data.boundary_fe_values.get_mapping(),
scratch_data.boundary_fe_values.get_fe(),
scratch_data.boundary_fe_values.get_quadrature(),
scratch_data.boundary_fe_values.get_update_flags())
, face_fe_values(scratch_data.face_fe_values.get_mapping(),
scratch_data.face_fe_values.get_fe(),
scratch_data.face_fe_values.get_quadrature(),
scratch_data.face_fe_values.get_update_flags())
, subface_values(scratch_data.subface_values.get_mapping(),
scratch_data.subface_values.get_fe(),
scratch_data.subface_values.get_quadrature(),
scratch_data.subface_values.get_update_flags())
, neighbor_face_values(
scratch_data.neighbor_face_values.get_mapping(),
scratch_data.neighbor_face_values.get_fe(),
scratch_data.neighbor_face_values.get_quadrature(),
scratch_data.neighbor_face_values.get_update_flags())
, neighbor_subface_values(
scratch_data.neighbor_subface_values.get_mapping(),
scratch_data.neighbor_subface_values.get_fe(),
scratch_data.neighbor_subface_values.get_quadrature(),
scratch_data.neighbor_subface_values.get_update_flags())
{}
};
template <int dim>
struct RightHandSideScratchData
{
: boundary_fe_values(mapping, fe, face_quadrature, face_update_flags)
, boundary_values(face_quadrature.
size())
{}
RightHandSideScratchData(const RightHandSideScratchData<dim> &scratch_data)
: boundary_fe_values(scratch_data.boundary_fe_values.get_mapping(),
scratch_data.boundary_fe_values.get_fe(),
scratch_data.boundary_fe_values.get_quadrature(),
, boundary_values(scratch_data.boundary_values.
size())
{}
std::vector<double> boundary_values;
};
template <int dim>
struct EstimatorScratchData
{
: fe_values(mapping, fe, cell_quadrature, cell_update_flags)
, boundary_fe_values(mapping,
fe,
boundary_quadrature,
boundary_update_flags)
, face_fe_values(mapping, fe, face_quadrature, face_update_flags)
, subface_values(mapping, fe, face_quadrature, face_update_flags)
, neighbor_face_values(mapping, fe, face_quadrature, face_update_flags)
, neighbor_subface_values(mapping, fe, face_quadrature, face_update_flags)
, cell_hessians(cell_quadrature.
size())
, boundary_solution_values(boundary_quadrature.
size())
, boundary_exact_values(boundary_quadrature.
size())
, face_solution_values(face_quadrature.
size())
, neighbor_face_solution_values(face_quadrature.
size())
, face_solution_gradients(face_quadrature.
size())
, neighbor_face_solution_gradients(face_quadrature.
size())
{}
EstimatorScratchData(const EstimatorScratchData<dim> &scratch_data)
: fe_values(scratch_data.fe_values.get_mapping(),
scratch_data.fe_values.get_fe(),
scratch_data.fe_values.get_quadrature(),
, boundary_fe_values(scratch_data.boundary_fe_values.get_mapping(),
scratch_data.boundary_fe_values.get_fe(),
scratch_data.boundary_fe_values.get_quadrature(),
, face_fe_values(scratch_data.face_fe_values.get_mapping(),
scratch_data.face_fe_values.get_fe(),
scratch_data.face_fe_values.get_quadrature(),
, subface_values(scratch_data.subface_values.get_mapping(),
scratch_data.subface_values.get_fe(),
scratch_data.subface_values.get_quadrature(),
, neighbor_face_values(
scratch_data.neighbor_face_values.get_mapping(),
scratch_data.neighbor_face_values.get_fe(),
scratch_data.neighbor_face_values.get_quadrature(),
, neighbor_subface_values(
scratch_data.neighbor_subface_values.get_mapping(),
scratch_data.neighbor_subface_values.get_fe(),
scratch_data.neighbor_subface_values.get_quadrature(),
, cell_hessians(scratch_data.cell_hessians.
size())
, boundary_solution_values(scratch_data.boundary_solution_values.
size())
, boundary_exact_values(scratch_data.boundary_exact_values.
size())
, face_solution_values(scratch_data.face_solution_values.
size())
, neighbor_face_solution_values(
scratch_data.neighbor_face_solution_values.
size())
, face_solution_gradients(scratch_data.face_solution_gradients.
size())
, neighbor_face_solution_gradients(
scratch_data.neighbor_face_solution_gradients.
size())
{}
std::vector<Tensor<2, dim>> cell_hessians;
std::vector<double> boundary_solution_values;
std::vector<double> boundary_exact_values;
std::vector<double> face_solution_values;
std::vector<double> neighbor_face_solution_values;
std::vector<Tensor<1, dim>> face_solution_gradients;
std::vector<Tensor<1, dim>> neighbor_face_solution_gradients;
};
template <int dim>
struct ErrorScratchData
{
: fe_values(mapping, fe, cell_quadrature, cell_update_flags)
, boundary_fe_values(mapping,
fe,
boundary_quadrature,
boundary_update_flags)
, face_fe_values(mapping, fe, face_quadrature, face_update_flags)
, subface_values(mapping, fe, face_quadrature, face_update_flags)
, neighbor_face_values(mapping, fe, face_quadrature, face_update_flags)
, neighbor_subface_values(mapping, fe, face_quadrature, face_update_flags)
, cell_solution_values(cell_quadrature.
size())
, cell_solution_gradients(cell_quadrature.
size())
, cell_exact_values(cell_quadrature.
size())
, cell_exact_gradients(cell_quadrature.
size())
, boundary_solution_values(boundary_quadrature.
size())
, boundary_exact_values(boundary_quadrature.
size())
, face_solution_values(face_quadrature.
size())
, neighbor_face_solution_values(face_quadrature.
size())
{}
ErrorScratchData(const ErrorScratchData<dim> &scratch_data)
: fe_values(scratch_data.fe_values.get_mapping(),
scratch_data.fe_values.get_fe(),
scratch_data.fe_values.get_quadrature(),
, boundary_fe_values(scratch_data.boundary_fe_values.get_mapping(),
scratch_data.boundary_fe_values.get_fe(),
scratch_data.boundary_fe_values.get_quadrature(),
, face_fe_values(scratch_data.face_fe_values.get_mapping(),
scratch_data.face_fe_values.get_fe(),
scratch_data.face_fe_values.get_quadrature(),
, subface_values(scratch_data.subface_values.get_mapping(),
scratch_data.subface_values.get_fe(),
scratch_data.subface_values.get_quadrature(),
, neighbor_face_values(
scratch_data.neighbor_face_values.get_mapping(),
scratch_data.neighbor_face_values.get_fe(),
scratch_data.neighbor_face_values.get_quadrature(),
, neighbor_subface_values(
scratch_data.neighbor_subface_values.get_mapping(),
scratch_data.neighbor_subface_values.get_fe(),
scratch_data.neighbor_subface_values.get_quadrature(),
, cell_solution_values(scratch_data.cell_solution_values.
size())
, cell_solution_gradients(scratch_data.cell_solution_gradients.
size())
, cell_exact_values(scratch_data.cell_exact_values.
size())
, cell_exact_gradients(scratch_data.cell_exact_gradients.
size())
, boundary_solution_values(scratch_data.boundary_solution_values.
size())
, boundary_exact_values(scratch_data.boundary_exact_values.
size())
, face_solution_values(scratch_data.face_solution_values.
size())
, neighbor_face_solution_values(
scratch_data.neighbor_face_solution_values.
size())
{}
std::vector<double> cell_solution_values;
std::vector<Tensor<1, dim>> cell_solution_gradients;
std::vector<double> cell_exact_values;
std::vector<Tensor<1, dim>> cell_exact_gradients;
std::vector<double> boundary_solution_values;
std::vector<double> boundary_exact_values;
std::vector<double> face_solution_values;
std::vector<double> neighbor_face_solution_values;
};
struct FaceCopyData
{
FaceCopyData()
{}
unsigned int level_1;
unsigned int level_2;
std::vector<types::global_dof_index> dof_indices_1;
std::vector<types::global_dof_index> dof_indices_2;
};
template <int dim>
struct MatrixCopyData
{
MatrixCopyData(const unsigned int dofs_per_cell = 0)
, local_dof_indices(dofs_per_cell)
{}
template <typename CellIterator>
void reinit(
const CellIterator &cell,
const bool use_level_dofs =
false)
{
const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
local_dof_indices.resize(dofs_per_cell);
if (use_level_dofs)
cell->get_mg_dof_indices(local_dof_indices);
else
cell->get_dof_indices(local_dof_indices);
face_data.clear();
}
template <typename CellIterator>
FaceCopyData &emplace_face_data(const CellIterator &cell_1,
const CellIterator &cell_2,
const bool use_level_dofs = false)
{
const unsigned int dofs_per_cell_1 = cell_1->get_fe().n_dofs_per_cell();
const unsigned int dofs_per_cell_2 = cell_2->get_fe().n_dofs_per_cell();
face_data.emplace_back();
FaceCopyData &face_copy = face_data.back();
face_copy.level_1 = cell_1->level();
face_copy.level_2 = cell_2->level();
face_copy.matrix_11.reinit(dofs_per_cell_1, dofs_per_cell_1);
face_copy.matrix_12.reinit(dofs_per_cell_1, dofs_per_cell_2);
face_copy.matrix_21.reinit(dofs_per_cell_2, dofs_per_cell_1);
face_copy.matrix_22.reinit(dofs_per_cell_2, dofs_per_cell_2);
face_copy.dof_indices_1.resize(dofs_per_cell_1);
face_copy.dof_indices_2.resize(dofs_per_cell_2);
if (use_level_dofs)
{
cell_1->get_mg_dof_indices(face_copy.dof_indices_1);
cell_2->get_mg_dof_indices(face_copy.dof_indices_2);
}
else
{
cell_1->get_dof_indices(face_copy.dof_indices_1);
cell_2->get_dof_indices(face_copy.dof_indices_2);
}
return face_copy;
}
std::vector<types::global_dof_index> local_dof_indices;
std::vector<FaceCopyData> face_data;
};
struct RightHandSideCopyData
{
RightHandSideCopyData(const unsigned int dofs_per_cell = 0)
: cell_rhs(dofs_per_cell)
, local_dof_indices(dofs_per_cell)
{}
template <typename CellIterator>
void reinit(
const CellIterator &cell)
{
const unsigned int dofs_per_cell = cell->get_fe().n_dofs_per_cell();
cell_rhs.reinit(dofs_per_cell);
local_dof_indices.resize(dofs_per_cell);
cell->get_dof_indices(local_dof_indices);
}
std::vector<types::global_dof_index> local_dof_indices;
};
template <unsigned int n_values>
struct ErrorCopyData
{
struct FaceContribution
{
FaceContribution()
{
}
unsigned int cell_index_1;
unsigned int cell_index_2;
std::array<double, n_values>
values;
};
ErrorCopyData()
{
cell_values.fill(0.);
}
template <typename CellIterator>
void reinit(
const CellIterator &cell)
{
cell_values.fill(0.);
face_data.clear();
}
template <typename CellIterator>
FaceContribution &emplace_face_data(const CellIterator &cell_1,
const CellIterator &cell_2)
{
face_data.emplace_back();
FaceContribution &face_contribution = face_data.back();
face_contribution.cell_index_1 = cell_1->active_cell_index();
face_contribution.cell_index_2 = cell_2->active_cell_index();
face_contribution.values.fill(0.);
return face_contribution;
}
std::array<double, n_values> cell_values;
std::vector<FaceContribution> face_data;
};
namespace MatrixIntegrator
{
template <int dim, typename CellIterator>
double ip_penalty_factor(const CellIterator &cell1,
const unsigned int face1,
const unsigned int deg1,
const CellIterator &cell2,
const unsigned int face2,
const 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 / cell1->extent_in_direction(normal1);
double penalty2 = deg2sq / cell2->extent_in_direction(normal2);
if (cell1->has_children() && !cell2->has_children())
penalty1 *= 2;
else if (!cell1->has_children() && cell2->has_children())
penalty2 *= 2;
const double penalty = 0.5 * (penalty1 + penalty2);
return penalty;
}
template <int dim>
{
{
const double dx = fe_values.
JxW(k);
{
const double Mii =
M(i, i) += Mii;
{
M(i, j) += Mij;
M(j, i) += Mij;
}
}
}
}
template <int dim, typename CellIterator>
void boundary(const CellIterator &cell,
const unsigned int face_no,
{
const unsigned int polynomial_degree =
const double ip_penalty = ip_penalty_factor<dim>(
cell, face_no, polynomial_degree, cell, face_no, polynomial_degree);
{
const double dx = fe_face_values.
JxW(k);
M(i, j) += (2. * fe_face_values.
shape_value(i, k) * ip_penalty *
dx;
}
}
template <int dim, typename CellIterator>
void face(const CellIterator &cell_1,
const unsigned int face_no_1,
&fe_face_values_1,
const CellIterator &cell_2,
const unsigned int face_no_2,
{
const unsigned int polynomial_degree =
const double ip_penalty = ip_penalty_factor<dim>(cell_1,
face_no_1,
polynomial_degree,
cell_2,
face_no_2,
polynomial_degree);
const double nui = 1.;
const double nue = 1.;
const double nu = .5 * (nui + nue);
{
const double dx = fe_face_values_1.
JxW(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 dnvi = n * fe_face_values_1.
shape_grad(i, k);
const double dnve = n * fe_face_values_2.
shape_grad(i, k);
const double dnui = n * fe_face_values_1.
shape_grad(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;
}
}
}
}
}
namespace RHSIntegrator
{
template <int dim, typename CellIterator>
void boundary(const CellIterator &cell,
const unsigned int face_no,
const std::vector<double> &boundary_values,
{
const double penalty = 2. * degree * (degree + 1) *
cell->face(face_no)->measure() / cell->measure();
local_vector(i) +=
+
* boundary_values[k] * fe.
JxW(k);
}
}
namespace Estimator
{
template <int dim, typename CellIterator>
double cell(const CellIterator &cell,
{
{
const double t = cell->diameter() *
trace(DDuh[k]);
}
}
template <int dim, typename CellIterator>
double boundary(const CellIterator &cell,
const unsigned int face_no,
const std::vector<double> &uh,
const std::vector<double> &boundary_values)
{
const double penalty = 2. * degree * (degree + 1) *
cell->face(face_no)->measure() / cell->measure();
{
const double diff = boundary_values[k] - uh[k];
value += penalty * diff * diff * fe.
JxW(k);
}
}
template <int dim, typename CellIterator>
double face(const CellIterator &cell_1,
const unsigned int face_no_1,
&fe,
const std::vector<double> &uh1,
const CellIterator &cell_2,
const unsigned int face_no_2,
const std::vector<double> &uh2,
{
const double penalty1 = degree * (degree + 1) *
cell_1->face(face_no_1)->measure() /
cell_1->measure();
const double penalty2 = degree * (degree + 1) *
cell_2->face(face_no_2)->measure() /
cell_2->measure();
const double penalty = penalty1 + penalty2;
const double h = cell_1->face(face_no_1)->measure();
{
const double diff1 = uh1[k] - uh2[k];
const double diff2 =
value += (penalty * diff1 * diff1 + h * diff2 * diff2) * fe.
JxW(k);
}
}
}
namespace ErrorIntegrator
{
template <int dim>
std::array<double, 2>
const std::vector<double> &uh,
const std::vector<double> &exact_values,
{
std::array<double, 2>
values = {};
{
for (
unsigned int d = 0;
d < dim; ++
d)
{
const double diff = exact_gradients[k][
d] - Duh[k][
d];
}
const double diff = exact_values[k] - uh[k];
}
}
template <int dim, typename CellIterator>
double boundary(const CellIterator &cell,
const unsigned int face_no,
const std::vector<double> &uh,
const std::vector<double> &exact_values)
{
const double penalty = 2. * degree * (degree + 1) *
cell->face(face_no)->measure() / cell->measure();
{
const double diff = exact_values[k] - uh[k];
value += penalty * diff * diff * fe.
JxW(k);
}
}
template <int dim, typename CellIterator>
double face(const CellIterator &cell_1,
const unsigned int face_no_1,
&fe,
const std::vector<double> &uh1,
const CellIterator &cell_2,
const unsigned int face_no_2,
const std::vector<double> &uh2)
{
const double penalty1 = degree * (degree + 1) *
cell_1->face(face_no_1)->measure() /
cell_1->measure();
const double penalty2 = degree * (degree + 1) *
cell_2->face(face_no_2)->measure() /
cell_2->measure();
const double penalty = penalty1 + penalty2;
{
const double diff = uh1[k] - uh2[k];
value += penalty * diff * diff * fe.
JxW(k);
}
}
}
template <int dim>
class InteriorPenaltyProblem
{
public:
InteriorPenaltyProblem();
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;
};
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);
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();
{
mg_sparsity[
level].copy_from(dsp);
{
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()
{
const MatrixScratchData<dim> scratch(mapping,
fe,
const MatrixCopyData<dim> copy_data(fe.n_dofs_per_cell());
dof_handler.begin_active(),
dof_handler.end(),
[&](const CellIterator &cell,
MatrixScratchData<dim> &scratch_data,
MatrixCopyData<dim> ©) {
copy.reinit(cell);
scratch_data.fe_values.reinit(cell);
MatrixIntegrator::cell<dim>(scratch_data.fe_values, copy.cell_matrix);
},
[&](const MatrixCopyData<dim> ©) {
matrix.add(copy.local_dof_indices, copy.cell_matrix);
for (const auto &face : copy.face_data)
{
matrix.add(face.dof_indices_1, face.dof_indices_1, face.matrix_11);
matrix.add(face.dof_indices_1, face.dof_indices_2, face.matrix_12);
matrix.add(face.dof_indices_2, face.dof_indices_1, face.matrix_21);
matrix.add(face.dof_indices_2, face.dof_indices_2, face.matrix_22);
}
},
scratch,
copy_data,
[&](const CellIterator &cell,
const unsigned int face_no,
MatrixScratchData<dim> &scratch_data,
MatrixCopyData<dim> &
copy) {
scratch_data.boundary_fe_values.reinit(cell, face_no);
MatrixIntegrator::boundary<dim>(cell,
face_no,
scratch_data.boundary_fe_values,
},
[&](const CellIterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const CellIterator &neighbor,
const unsigned int neighbor_face_no,
const unsigned int neighbor_subface_no,
MatrixScratchData<dim> &scratch_data,
MatrixCopyData<dim> &
copy) {
FaceCopyData &face_copy =
copy.emplace_face_data(cell, neighbor);
{
scratch_data.face_fe_values.reinit(cell, face_no);
{
scratch_data.neighbor_face_values.reinit(neighbor,
neighbor_face_no);
MatrixIntegrator::face<dim>(cell,
face_no,
scratch_data.face_fe_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_values,
face_copy.matrix_11,
face_copy.matrix_12,
face_copy.matrix_21,
face_copy.matrix_22);
}
else
{
scratch_data.neighbor_subface_values.reinit(
neighbor, neighbor_face_no, neighbor_subface_no);
MatrixIntegrator::face<dim>(
cell,
face_no,
scratch_data.face_fe_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_subface_values,
face_copy.matrix_11,
face_copy.matrix_12,
face_copy.matrix_21,
face_copy.matrix_22);
}
}
else
{
scratch_data.subface_values.reinit(cell, face_no, subface_no);
{
scratch_data.neighbor_face_values.reinit(neighbor,
neighbor_face_no);
MatrixIntegrator::face<dim>(cell,
face_no,
scratch_data.subface_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_values,
face_copy.matrix_11,
face_copy.matrix_12,
face_copy.matrix_21,
face_copy.matrix_22);
}
else
{
scratch_data.neighbor_subface_values.reinit(
neighbor, neighbor_face_no, neighbor_subface_no);
MatrixIntegrator::face<dim>(
cell,
face_no,
scratch_data.subface_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_subface_values,
face_copy.matrix_11,
face_copy.matrix_12,
face_copy.matrix_21,
face_copy.matrix_22);
}
}
});
}
template <int dim>
void InteriorPenaltyProblem<dim>::assemble_mg_matrix()
{
const MatrixScratchData<dim> scratch(mapping,
fe,
const MatrixCopyData<dim> copy_data(fe.n_dofs_per_cell());
dof_handler.begin_mg(),
dof_handler.end_mg(),
[&](const CellIterator &cell,
MatrixScratchData<dim> &scratch_data,
MatrixCopyData<dim> ©) {
copy.reinit(cell, true);
scratch_data.fe_values.reinit(cell);
MatrixIntegrator::cell<dim>(scratch_data.fe_values, copy.cell_matrix);
},
[&](const MatrixCopyData<dim> ©) {
mg_matrix[copy.level].add(copy.local_dof_indices, copy.cell_matrix);
for (const auto &face : copy.face_data)
if (face.level_1 == face.level_2)
{
mg_matrix[face.level_1].add(face.dof_indices_1,
face.dof_indices_1,
face.matrix_11);
mg_matrix[face.level_1].add(face.dof_indices_1,
face.dof_indices_2,
face.matrix_12);
mg_matrix[face.level_1].add(face.dof_indices_2,
face.dof_indices_1,
face.matrix_21);
mg_matrix[face.level_1].add(face.dof_indices_2,
face.dof_indices_2,
face.matrix_22);
}
else
{
Assert(face.level_1 > face.level_2, ExcInternalError());
mg_matrix[face.level_1].add(face.dof_indices_1,
face.dof_indices_1,
face.matrix_11);
for (unsigned int j = 0; j < face.dof_indices_2.size(); ++j)
for (unsigned int k = 0; k < face.dof_indices_1.size(); ++k)
{
mg_matrix_dg_up[face.level_1].add(face.dof_indices_2[j],
face.dof_indices_1[k],
face.matrix_12(k, j));
mg_matrix_dg_down[face.level_1].add(face.dof_indices_2[j],
face.dof_indices_1[k],
face.matrix_21(j, k));
}
}
},
scratch,
copy_data,
[&](const CellIterator &cell,
const unsigned int face_no,
MatrixScratchData<dim> &scratch_data,
MatrixCopyData<dim> &
copy) {
scratch_data.boundary_fe_values.reinit(cell, face_no);
MatrixIntegrator::boundary<dim>(cell,
face_no,
scratch_data.boundary_fe_values,
},
[&](const CellIterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const CellIterator &neighbor,
const unsigned int neighbor_face_no,
const unsigned int neighbor_subface_no,
MatrixScratchData<dim> &scratch_data,
MatrixCopyData<dim> &
copy) {
FaceCopyData &face_copy =
copy.emplace_face_data(cell, neighbor,
true);
{
scratch_data.face_fe_values.reinit(cell, face_no);
{
scratch_data.neighbor_face_values.reinit(neighbor,
neighbor_face_no);
MatrixIntegrator::face<dim>(cell,
face_no,
scratch_data.face_fe_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_values,
face_copy.matrix_11,
face_copy.matrix_12,
face_copy.matrix_21,
face_copy.matrix_22);
}
else
{
scratch_data.neighbor_subface_values.reinit(
neighbor, neighbor_face_no, neighbor_subface_no);
MatrixIntegrator::face<dim>(
cell,
face_no,
scratch_data.face_fe_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_subface_values,
face_copy.matrix_11,
face_copy.matrix_12,
face_copy.matrix_21,
face_copy.matrix_22);
}
}
else
{
scratch_data.subface_values.reinit(cell, face_no, subface_no);
{
scratch_data.neighbor_face_values.reinit(neighbor,
neighbor_face_no);
MatrixIntegrator::face<dim>(cell,
face_no,
scratch_data.subface_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_values,
face_copy.matrix_11,
face_copy.matrix_12,
face_copy.matrix_21,
face_copy.matrix_22);
}
else
{
scratch_data.neighbor_subface_values.reinit(
neighbor, neighbor_face_no, neighbor_subface_no);
MatrixIntegrator::face<dim>(
cell,
face_no,
scratch_data.subface_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_subface_values,
face_copy.matrix_11,
face_copy.matrix_12,
face_copy.matrix_21,
face_copy.matrix_22);
}
}
});
}
template <int dim>
void InteriorPenaltyProblem<dim>::assemble_right_hand_side()
{
const RightHandSideScratchData<dim> scratch(mapping,
fe,
const RightHandSideCopyData copy_data(fe.n_dofs_per_cell());
dof_handler.begin_active(),
dof_handler.end(),
[&](const CellIterator &cell,
RightHandSideScratchData<dim> &,
RightHandSideCopyData ©) { copy.reinit(cell); },
[&](const RightHandSideCopyData ©) {
right_hand_side.add(copy.local_dof_indices, copy.cell_rhs);
},
scratch,
copy_data,
[&](const CellIterator &cell,
const unsigned int face_no,
RightHandSideScratchData<dim> &scratch_data,
RightHandSideCopyData ©) {
scratch_data.boundary_fe_values.reinit(cell, face_no);
exact_solution.value_list(
scratch_data.boundary_fe_values.get_quadrature_points(),
scratch_data.boundary_values);
RHSIntegrator::boundary<dim>(cell,
face_no,
scratch_data.boundary_fe_values,
scratch_data.boundary_values,
copy.cell_rhs);
});
right_hand_side *= -1.;
}
template <int dim>
void InteriorPenaltyProblem<dim>::solve()
{
mg_transfer.
build(dof_handler);
RELAXATION::AdditionalData smoother_data(1.);
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);
std::cout << "Converged in " << control.last_step() << " iterations"
<< std::endl;
}
template <int dim>
double InteriorPenaltyProblem<dim>::estimate()
{
estimates.block(0).reinit(triangulation.n_active_cells());
const unsigned int n_gauss_points =
dof_handler.get_fe().tensor_degree() + 1;
const EstimatorScratchData<dim> scratch(mapping,
fe,
const ErrorCopyData<1> copy_data;
dof_handler.begin_active(),
dof_handler.end(),
[&](const CellIterator &cell,
EstimatorScratchData<dim> &scratch_data,
ErrorCopyData<1> ©) {
copy.reinit(cell);
scratch_data.fe_values.reinit(cell);
const FEValues<dim> &fe_values = scratch_data.fe_values;
fe_values.get_function_hessians(solution, scratch_data.cell_hessians);
copy.cell_values[0] =
Estimator::cell<dim>(cell, fe_values, scratch_data.cell_hessians);
},
[&](const ErrorCopyData<1> ©) {
estimates.block(0)(copy.cell_index) += copy.cell_values[0];
for (const auto &face : copy.face_data)
{
estimates.block(0)(face.cell_index_1) += 0.5 * face.values[0];
estimates.block(0)(face.cell_index_2) += 0.5 * face.values[0];
}
},
scratch,
copy_data,
[&](const CellIterator &cell,
const unsigned int face_no,
EstimatorScratchData<dim> &scratch_data,
ErrorCopyData<1> &
copy) {
scratch_data.boundary_fe_values.reinit(cell, face_no);
scratch_data.boundary_fe_values;
solution, scratch_data.boundary_solution_values);
scratch_data.boundary_exact_values);
Estimator::boundary<dim>(cell,
face_no,
fe_face_values,
scratch_data.boundary_solution_values,
scratch_data.boundary_exact_values);
},
[&](const CellIterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const CellIterator &neighbor,
const unsigned int neighbor_face_no,
const unsigned int neighbor_subface_no,
EstimatorScratchData<dim> &scratch_data,
ErrorCopyData<1> &
copy) {
auto &face_data =
copy.emplace_face_data(cell, neighbor);
{
scratch_data.face_fe_values.reinit(cell, face_no);
scratch_data.face_fe_values;
{
scratch_data.neighbor_face_values.
reinit(neighbor,
neighbor_face_no);
scratch_data.neighbor_face_values;
solution, scratch_data.face_solution_values);
solution, scratch_data.neighbor_face_solution_values);
solution, scratch_data.face_solution_gradients);
solution, scratch_data.neighbor_face_solution_gradients);
face_data.values[0] = Estimator::face<dim>(
cell,
face_no,
fe_face_values,
scratch_data.face_solution_values,
scratch_data.face_solution_gradients,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_solution_values,
scratch_data.neighbor_face_solution_gradients);
}
else
{
scratch_data.neighbor_subface_values.reinit(
neighbor, neighbor_face_no, neighbor_subface_no);
scratch_data.neighbor_subface_values;
solution, scratch_data.face_solution_values);
solution, scratch_data.neighbor_face_solution_values);
solution, scratch_data.face_solution_gradients);
solution, scratch_data.neighbor_face_solution_gradients);
face_data.values[0] = Estimator::face<dim>(
cell,
face_no,
fe_face_values,
scratch_data.face_solution_values,
scratch_data.face_solution_gradients,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_solution_values,
scratch_data.neighbor_face_solution_gradients);
}
}
else
{
scratch_data.subface_values.reinit(cell, face_no, subface_no);
scratch_data.subface_values;
{
scratch_data.neighbor_face_values.reinit(neighbor,
neighbor_face_no);
scratch_data.neighbor_face_values;
solution, scratch_data.face_solution_values);
solution, scratch_data.neighbor_face_solution_values);
solution, scratch_data.face_solution_gradients);
solution, scratch_data.neighbor_face_solution_gradients);
face_data.values[0] = Estimator::face<dim>(
cell,
face_no,
fe_face_values,
scratch_data.face_solution_values,
scratch_data.face_solution_gradients,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_solution_values,
scratch_data.neighbor_face_solution_gradients);
}
else
{
scratch_data.neighbor_subface_values.reinit(
neighbor, neighbor_face_no, neighbor_subface_no);
scratch_data.neighbor_subface_values;
solution, scratch_data.face_solution_values);
solution, scratch_data.neighbor_face_solution_values);
solution, scratch_data.face_solution_gradients);
solution, scratch_data.neighbor_face_solution_gradients);
face_data.values[0] = Estimator::face<dim>(
cell,
face_no,
fe_face_values,
scratch_data.face_solution_values,
scratch_data.face_solution_gradients,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_solution_values,
scratch_data.neighbor_face_solution_gradients);
}
}
});
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());
const unsigned int n_gauss_points =
dof_handler.get_fe().tensor_degree() + 1;
const ErrorScratchData<dim> scratch(mapping,
fe,
const ErrorCopyData<2> copy_data;
dof_handler.begin_active(),
dof_handler.end(),
[&](const CellIterator &cell,
ErrorScratchData<dim> &scratch_data,
ErrorCopyData<2> ©) {
copy.reinit(cell);
scratch_data.fe_values.reinit(cell);
const FEValues<dim> &fe_values = scratch_data.fe_values;
fe_values.get_function_values(solution,
scratch_data.cell_solution_values);
fe_values.get_function_gradients(solution,
scratch_data.cell_solution_gradients);
exact_solution.value_list(fe_values.get_quadrature_points(),
scratch_data.cell_exact_values);
exact_solution.gradient_list(fe_values.get_quadrature_points(),
scratch_data.cell_exact_gradients);
copy.cell_values =
ErrorIntegrator::cell<dim>(fe_values,
scratch_data.cell_solution_values,
scratch_data.cell_solution_gradients,
scratch_data.cell_exact_values,
scratch_data.cell_exact_gradients);
},
[&](const ErrorCopyData<2> ©) {
errors.block(0)(copy.cell_index) += copy.cell_values[0];
errors.block(1)(copy.cell_index) += copy.cell_values[1];
for (const auto &face : copy.face_data)
{
errors.block(0)(face.cell_index_1) += 0.5 * face.values[0];
errors.block(0)(face.cell_index_2) += 0.5 * face.values[0];
}
},
scratch,
copy_data,
[&](const CellIterator &cell,
const unsigned int face_no,
ErrorScratchData<dim> &scratch_data,
ErrorCopyData<2> &
copy) {
scratch_data.boundary_fe_values.reinit(cell, face_no);
scratch_data.boundary_fe_values;
solution, scratch_data.boundary_solution_values);
scratch_data.boundary_exact_values);
ErrorIntegrator::boundary<dim>(cell,
face_no,
fe_face_values,
scratch_data.boundary_solution_values,
scratch_data.boundary_exact_values);
},
[&](const CellIterator &cell,
const unsigned int face_no,
const unsigned int subface_no,
const CellIterator &neighbor,
const unsigned int neighbor_face_no,
const unsigned int neighbor_subface_no,
ErrorScratchData<dim> &scratch_data,
ErrorCopyData<2> &
copy) {
auto &face_data =
copy.emplace_face_data(cell, neighbor);
{
scratch_data.face_fe_values.reinit(cell, face_no);
scratch_data.face_fe_values;
{
scratch_data.neighbor_face_values.
reinit(neighbor,
neighbor_face_no);
scratch_data.neighbor_face_values;
solution, scratch_data.face_solution_values);
solution, scratch_data.neighbor_face_solution_values);
face_data.values[0] = ErrorIntegrator::face<dim>(
cell,
face_no,
fe_face_values,
scratch_data.face_solution_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_solution_values);
}
else
{
scratch_data.neighbor_subface_values.reinit(
neighbor, neighbor_face_no, neighbor_subface_no);
scratch_data.neighbor_subface_values;
solution, scratch_data.face_solution_values);
solution, scratch_data.neighbor_face_solution_values);
face_data.values[0] = ErrorIntegrator::face<dim>(
cell,
face_no,
fe_face_values,
scratch_data.face_solution_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_solution_values);
}
}
else
{
scratch_data.subface_values.reinit(cell, face_no, subface_no);
scratch_data.subface_values;
{
scratch_data.neighbor_face_values.reinit(neighbor,
neighbor_face_no);
scratch_data.neighbor_face_values;
solution, scratch_data.face_solution_values);
solution, scratch_data.neighbor_face_solution_values);
face_data.values[0] = ErrorIntegrator::face<dim>(
cell,
face_no,
fe_face_values,
scratch_data.face_solution_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_solution_values);
}
else
{
scratch_data.neighbor_subface_values.reinit(
neighbor, neighbor_face_no, neighbor_subface_no);
scratch_data.neighbor_subface_values;
solution, scratch_data.face_solution_values);
solution, scratch_data.neighbor_face_solution_values);
face_data.values[0] = ErrorIntegrator::face<dim>(
cell,
face_no,
fe_face_values,
scratch_data.face_solution_values,
neighbor,
neighbor_face_no,
scratch_data.neighbor_face_solution_values);
}
}
});
std::cout << "energy-error: " << errors.block(0).l2_norm() << std::endl;
std::cout << "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 =
std::cout << "Writing solution to <" << filename << ">..." << std::endl
<< std::endl;
std::ofstream gnuplot_output(filename);
}
template <int dim>
void InteriorPenaltyProblem<dim>::run(
unsigned int n_steps)
{
std::cout << "Element: " << fe.get_name() << std::endl;
for (
unsigned int s = 0; s <
n_steps; ++s)
{
std::cout << "Step " << s << std::endl;
if (estimates.block(0).empty())
triangulation.refine_global(1);
else
{
triangulation, estimates.block(0), 0.5, 0.0);
triangulation.execute_coarsening_and_refinement();
}
std::cout << "Triangulation " << triangulation.n_active_cells()
<< " cells, " << triangulation.n_levels() << " levels"
<< std::endl;
setup_system();
std::cout << "DoFHandler " << dof_handler.n_dofs()
<< " dofs, level dofs";
for (
unsigned int l = 0;
l < triangulation.n_levels(); ++
l)
std::cout << ' ' << dof_handler.n_dofs(l);
std::cout << std::endl;
std::cout << "Assemble matrix" << std::endl;
assemble_matrix();
std::cout << "Assemble multilevel matrix" << std::endl;
assemble_mg_matrix();
std::cout << "Assemble right hand side" << std::endl;
assemble_right_hand_side();
std::cout << "Solve" << std::endl;
solve();
error();
std::cout << "Estimate " << estimate() << std::endl;
output_results(s);
}
}
}
{
try
{
using namespace Step39;
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)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell, const unsigned int face_no)
const std::vector< Point< spacedim > > & get_quadrature_points() const
const Tensor< 1, spacedim > & normal_vector(const unsigned int q_point) const
void get_function_gradients(const ReadVector< Number > &fe_function, std::vector< Tensor< 1, spacedim, Number > > &gradients) const
const FiniteElement< dim, spacedim > & get_fe() const
const double & shape_value(const unsigned int i, const unsigned int q_point) const
unsigned int tensor_degree() const
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(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
@ matrix
Contents is actually a matrix.
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
* * * * * * TimeRateUpdateFlags TimeRateRequest< ValueType, dim, Number > get_update_flags() const
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Using the perl script 
1.9.8