The output of this program consist of the console output and
solutions in vtu format.
-Convergence rate for the smooth case with polynomial degree 3:
+In the first test case, when you run the program, the screen output should look like the following:
@code
-DEAL::Cycle 0
-DEAL::Number of active cells: 16
-DEAL::Number of degrees of freedom: 256
- Solving system...
-DEAL::Writing solution to <sol_Q3-00.vtu>
- Error in the L2 norm : 0.00193285
- Error in the H1 seminorm : 0.106087
- Error in the enery norm : 0.150625
-
-DEAL::Cycle 1
-DEAL::Number of active cells: 64
-DEAL::Number of degrees of freedom: 1024
- Solving system...
-DEAL::Writing solution to <sol_Q3-01.vtu>
- Error in the L2 norm : 9.60497e-05
- Error in the H1 seminorm : 0.0089954
- Error in the enery norm : 0.0113265
-
-DEAL::Cycle 2
-DEAL::Number of active cells: 256
-DEAL::Number of degrees of freedom: 4096
- Solving system...
-DEAL::Writing solution to <sol_Q3-02.vtu>
- Error in the L2 norm : 5.60638e-06
- Error in the H1 seminorm : 0.000901791
- Error in the enery norm : 0.000973568
-
-DEAL::Cycle 3
-DEAL::Number of active cells: 1024
-DEAL::Number of degrees of freedom: 16384
- Solving system...
-DEAL::Writing solution to <sol_Q3-03.vtu>
- Error in the L2 norm : 3.48399e-07
- Error in the H1 seminorm : 0.000107055
- Error in the enery norm : 0.000108757
-
-DEAL::Cycle 4
-DEAL::Number of active cells: 4096
-DEAL::Number of degrees of freedom: 65536
- Solving system...
-DEAL::Writing solution to <sol_Q3-04.vtu>
- Error in the L2 norm : 2.17928e-08
- Error in the H1 seminorm : 1.32657e-05
- Error in the enery norm : 1.33065e-05
-
-DEAL::Cycle 5
-DEAL::Number of active cells: 16384
-DEAL::Number of degrees of freedom: 262144
- Solving system...
-DEAL::Writing solution to <sol_Q3-05.vtu>
- Error in the L2 norm : 1.36274e-09
- Error in the H1 seminorm : 1.65576e-06
- Error in the enery norm : 1.65683e-06
-
-degree = 3
-| cycle | cells | dofs | L2 | L2...red.rate.log2 | H1 | H1...red.rate.log2 | Energy |
-| 0 | 16 | 256 | 1.933e-03 | - | 1.061e-01 | - | 1.506e-01 |
-| 1 | 64 | 1024 | 9.605e-05 | 4.33 | 8.995e-03 | 3.56 | 1.133e-02 |
-| 2 | 256 | 4096 | 5.606e-06 | 4.10 | 9.018e-04 | 3.32 | 9.736e-04 |
-| 3 | 1024 | 16384 | 3.484e-07 | 4.01 | 1.071e-04 | 3.07 | 1.088e-04 |
-| 4 | 4096 | 65536 | 2.179e-08 | 4.00 | 1.327e-05 | 3.01 | 1.331e-05 |
-| 5 | 16384 | 262144 | 1.363e-09 | 4.00 | 1.656e-06 | 3.00 | 1.657e-06 |
+Cycle 0
+ Number of active cells : 16
+ Number of degrees of freedom : 256
+ Error in the L2 norm : 0.00193285
+ Error in the H1 seminorm : 0.106087
+ Error in the energy norm : 0.150625
+
+Cycle 1
+ Number of active cells : 64
+ Number of degrees of freedom : 1024
+ Error in the L2 norm : 9.60497e-05
+ Error in the H1 seminorm : 0.0089954
+ Error in the energy norm : 0.0113265
+
+Cycle 2
+.
+.
+.
@endcode
+Convergence rate for the smooth case with polynomial degree 3:
+<table align="center" class="doxtable">
+ <tr>
+ <th>cycle</th>
+ <th>n_cellss</th>
+ <th>n_dofs</th>
+ <th>L2 </th>
+ <th>rate</th>
+ <th>H1</th>
+ <th>rate</th>
+ <th>Energy</th>
+ </tr>
+ <tr>
+ <td align="center">0</td>
+ <td align="right">16</td>
+ <td align="right">256</td>
+ <td align="center">1.933e-03</td>
+ <td> </td>
+ <td align="center">1.061e-01</td>
+ <td> </td>
+ <td align="center">1.506e-01</td>
+ </tr>
+ <tr>
+ <td align="center">1</td>
+ <td align="right">64</td>
+ <td align="right">1024</td>
+ <td align="center">9.605e-05</td>
+ <td align="center">4.33</td>
+ <td align="center">8.995e-03</td>
+ <td align="center">3.56</td>
+ <td align="center">1.133e-02</td>
+ </tr>
+ <tr>
+ <td align="center">2</td>
+ <td align="right">256</td>
+ <td align="right">4096</td>
+ <td align="center">5.606e-06</td>
+ <td align="center">4.10</td>
+ <td align="center">9.018e-04</td>
+ <td align="center">3.32</td>
+ <td align="center">9.736e-04</td>
+ </tr>
+ <tr>
+ <td align="center">3</td>
+ <td align="right">1024</td>
+ <td align="right">16384</td>
+ <td align="center">3.484e-07</td>
+ <td align="center">4.01</td>
+ <td align="center">1.071e-04</td>
+ <td align="center">3.07</td>
+ <td align="center">1.088e-04</td>
+ </tr>
+ <tr>
+ <td align="center">4</td>
+ <td align="right">4096</td>
+ <td align="right">65536</td>
+ <td align="center">2.179e-08</td>
+ <td align="center">4.00</td>
+ <td align="center">1.327e-05</td>
+ <td align="center">3.01</td>
+ <td align="center">1.331e-05</td>
+ </tr>
+ <tr>
+ <td align="center">5</td>
+ <td align="right">16384</td>
+ <td align="right">262144</td>
+ <td align="center">1.363e-09</td>
+ <td align="center">4.00</td>
+ <td align="center">1.656e-06</td>
+ <td align="center">3.00</td>
+ <td align="center">1.657e-06</td>
+ </tr>
+</table>
+
Theoretically, for polynomial degree $p$, the order of convergence in $L_2$
norm and $H_1$ seminorm should be $p+1$ and $p$, respectively. Our numerical
results are in good agreement with theory.
-Errors of the singular case with polynomial degree 3.
+In the second test case, when you run the program, the screen output should look like the following:
@code
-DEAL::Cycle 0
-DEAL::Number of active cells: 12
-DEAL::Number of degrees of freedom: 192
- Solving system...
-DEAL::Writing solution to <sol_Q3-00.vtu>
- Error in the L2 norm : 0.00278279
- Error in the H1 seminorm : 0.0748987
- Error in the enery norm : 0.106218
-
-DEAL::Cycle 1
-DEAL::Number of active cells: 15
-DEAL::Number of degrees of freedom: 240
- Solving system...
-DEAL::Writing solution to <sol_Q3-01.vtu>
- Error in the L2 norm : 0.00179741
- Error in the H1 seminorm : 0.0568531
- Error in the enery norm : 0.0815378
-
-DEAL::Cycle 2
-DEAL::Number of active cells: 18
-DEAL::Number of degrees of freedom: 288
- Solving system...
-DEAL::Writing solution to <sol_Q3-02.vtu>
- Error in the L2 norm : 0.00171775
- Error in the H1 seminorm : 0.0598664
- Error in the enery norm : 0.0850871
-
-DEAL::Cycle 3
-DEAL::Number of active cells: 21
-DEAL::Number of degrees of freedom: 336
- Solving system...
-DEAL::Writing solution to <sol_Q3-03.vtu>
- Error in the L2 norm : 0.000939329
- Error in the H1 seminorm : 0.0470787
- Error in the enery norm : 0.0668118
-
-DEAL::Cycle 4
-DEAL::Number of active cells: 27
-DEAL::Number of degrees of freedom: 432
- Solving system...
-DEAL::Writing solution to <sol_Q3-04.vtu>
- Error in the L2 norm : 0.000339722
- Error in the H1 seminorm : 0.0273542
- Error in the enery norm : 0.0394557
-
-DEAL::Cycle 5
-DEAL::Number of active cells: 33
-DEAL::Number of degrees of freedom: 528
- Solving system...
-DEAL::Writing solution to <sol_Q3-05.vtu>
- Error in the L2 norm : 0.000217819
- Error in the H1 seminorm : 0.0226131
- Error in the enery norm : 0.0324468
-
-DEAL::Cycle 6
-DEAL::Number of active cells: 42
-DEAL::Number of degrees of freedom: 672
- Solving system...
-DEAL::Writing solution to <sol_Q3-06.vtu>
- Error in the L2 norm : 9.95869e-05
- Error in the H1 seminorm : 0.0143484
- Error in the enery norm : 0.0205855
-
-DEAL::Cycle 7
-DEAL::Number of active cells: 54
-DEAL::Number of degrees of freedom: 864
- Solving system...
-DEAL::Writing solution to <sol_Q3-07.vtu>
- Error in the L2 norm : 6.45972e-05
- Error in the H1 seminorm : 0.00914956
- Error in the enery norm : 0.0131191
-
-DEAL::Cycle 8
-DEAL::Number of active cells: 69
-DEAL::Number of degrees of freedom: 1104
- Solving system...
-DEAL::Writing solution to <sol_Q3-08.vtu>
- Error in the L2 norm : 4.16999e-05
- Error in the H1 seminorm : 0.00584172
- Error in the enery norm : 0.0083634
-
-DEAL::Cycle 9
-DEAL::Number of active cells: 87
-DEAL::Number of degrees of freedom: 1392
- Solving system...
-DEAL::Writing solution to <sol_Q3-09.vtu>
- Error in the L2 norm : 2.7367e-05
- Error in the H1 seminorm : 0.00372588
- Error in the enery norm : 0.00531894
-
-degree = 3
-| cycle | cells | dofs | L2 | H1 | Energy | Estimator |
-| 0 | 12 | 192 | 2.783e-03 | 7.490e-02 | 1.062e-01 | 3.455e-01 |
-| 1 | 15 | 240 | 1.797e-03 | 5.685e-02 | 8.154e-02 | 3.155e-01 |
-| 2 | 18 | 288 | 1.718e-03 | 5.987e-02 | 8.509e-02 | 2.582e-01 |
-| 3 | 21 | 336 | 9.393e-04 | 4.708e-02 | 6.681e-02 | 2.174e-01 |
-| 4 | 27 | 432 | 3.397e-04 | 2.735e-02 | 3.946e-02 | 2.014e-01 |
-| 5 | 33 | 528 | 2.178e-04 | 2.261e-02 | 3.245e-02 | 1.277e-01 |
-| 6 | 42 | 672 | 9.959e-05 | 1.435e-02 | 2.059e-02 | 8.352e-02 |
-| 7 | 54 | 864 | 6.460e-05 | 9.150e-03 | 1.312e-02 | 5.609e-02 |
-| 8 | 69 | 1104 | 4.170e-05 | 5.842e-03 | 8.363e-03 | 3.807e-02 |
-| 9 | 87 | 1392 | 2.737e-05 | 3.726e-03 | 5.319e-03 | 2.590e-02 |
+Cycle 0
+ Number of active cells : 192
+ Number of degrees of freedom : 3072
+ Error in the L2 norm : 0.000323585
+ Error in the H1 seminorm : 0.0296202
+ Error in the energy norm : 0.0420478
+ Estimated error : 0.136067
+
+Cycle 1
+ Number of active cells : 249
+ Number of degrees of freedom : 3984
+ Error in the L2 norm : 0.000114739
+ Error in the H1 seminorm : 0.0186571
+ Error in the energy norm : 0.0264879
+ Estimated error : 0.0857186
+
+Cycle 2
+.
+.
+.
@endcode
The following figure provides a log-log plot of the errors versus
the number of degrees of freedom. Let $n$ be the number of degrees of
-freedom, then $h$ is of order $1/\sqrt{n}$ in 2d. Combining the theoretical
+freedom, then $h$ is of order $1/\sqrt{n}$ in 2D. Combining the theoretical
results in the previous case,
we see that the error in $L_2$ norm is of order $O(n^{-\frac{p+1}{2}})$
and in $H_1$ seminorm is $O(n^{-\frac{p}{2}})$. From the figure, we see
its ability to predict regions with large errors.
<img width="600px" src="https://www.dealii.org/images/steps/developer/step-74-log-log-plot.png" alt="">
+
+While this tutorial is focused on the implementation, the step-59 tutorial program achieves an efficient
+large-scale solver in terms of computing time with matrix-free solution techniques.
+Note that the step-59 tutorial does not work with meshes containing hanging nodes at this moment,
+because the multigrid interface matrices are not as easily determined,
+but that is merely the lack of some interfaces in deal.II, nothing fundamental.
\ No newline at end of file
/* ---------------------------------------------------------------------
*
- * Copyright (C) 2009 - 2015 by the deal.II authors
+ * Copyright (C) 2020 by the deal.II authors
*
* This file is part of the deal.II library.
*
#include <deal.II/fe/mapping_q1.h>
// Here the discontinuous finite elements and FEInterfaceValues are defined.
#include <deal.II/fe/fe_dgq.h>
-#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_interface_values.h>
#include <deal.II/numerics/derivative_approximation.h>
{}
virtual void value_list(const std::vector<Point<dim>> &points,
std::vector<double> & values,
- const unsigned int /*component*/) const override
- {
- using numbers::PI;
- for (unsigned int i = 0; i < values.size(); ++i)
- values[i] = 8. * PI * PI * std::sin(2. * PI * points[i][0]) *
- std::sin(2. * PI * points[i][1]);
- }
+ const unsigned int /*component*/) const override;
};
+ template <int dim>
+ void
+ SmoothRightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
+ std::vector<double> & values,
+ const unsigned int /*component*/) const
+ {
+ using numbers::PI;
+ for (unsigned int i = 0; i < values.size(); ++i)
+ values[i] = 8. * PI * PI * std::sin(2. * PI * points[i][0]) *
+ std::sin(2. * PI * points[i][1]);
+ }
+
// The right-hand side corresponds to the function
// Functions::LSingularityFunction.
template <int dim>
{}
virtual void value_list(const std::vector<Point<dim>> &points,
std::vector<double> & values,
- const unsigned int /*component*/) const override
- {
- for (unsigned int i = 0; i < values.size(); ++i)
- // We assume that the diffusion coefficient $\nu$ = 1.
- values[i] = -ref.laplacian(points[i]);
- }
+ const unsigned int /*component*/) const override;
private:
Functions::LSingularityFunction ref;
};
+ template <int dim>
+ void
+ SingularRightHandSide<dim>::value_list(const std::vector<Point<dim>> &points,
+ std::vector<double> & values,
+ const unsigned int /*component*/) const
+ {
+ for (unsigned int i = 0; i < values.size(); ++i)
+ // We assume that the diffusion coefficient $\nu$ = 1.
+ values[i] = -ref.laplacian(points[i]);
+ }
+
// @sect3{Auxiliary functions}
// The following two auxiliary functions are used to compute
// jump terms for $u_h$ and $\nabla u_h$ on the
const Vector<double> & solution,
std::vector<double> & jump)
{
- const unsigned n_q = fe_iv.n_quadrature_points;
- std::vector<double> face_values[2];
+ const unsigned n_q = fe_iv.n_quadrature_points;
+ std::array<std::vector<double>, 2> face_values;
jump.resize(n_q);
for (unsigned i = 0; i < 2; ++i)
{
{
FullMatrix<double> cell_matrix;
std::vector<types::global_dof_index> joint_dof_indices;
- double values[2];
- unsigned int cell_indices[2];
+ std::array<double, 2> values;
+ std::array<unsigned int, 2> cell_indices;
};
struct CopyData
double compute_energy_norm();
Triangulation<dim> triangulation;
+ const unsigned degree;
+ QGauss<dim> quadrature;
+ QGauss<dim - 1> face_quadrature;
const MappingQ1<dim> mapping;
using ScratchData = MeshWorker::ScratchData<dim>;
};
// The constructor here reads the test case as an input and then determines
- // the correct solution and right-hand side classes. The 3 in the constructor
- // call of fe is the polynomial degree.
+ // the correct solution and right-hand side classes.
template <int dim>
SIPGLaplace<dim>::SIPGLaplace(const Test_Case &test_case)
- : mapping()
- , fe(3)
+ : degree(3)
+ , quadrature(degree + 1)
+ , face_quadrature(degree + 1)
+ , mapping()
+ , fe(degree)
, dof_handler(triangulation)
, test_case(test_case)
{
template <int dim>
void SIPGLaplace<dim>::assemble_system()
{
- typedef decltype(dof_handler.begin_active()) Iterator;
-
// This function assembles the cell integrals.
- auto cell_worker = [&](const Iterator &cell,
- ScratchData & scratch_data,
- CopyData & copy_data) {
- const FEValues<dim> &fe_v = scratch_data.reinit(cell);
- const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
- copy_data.reinit(cell, dofs_per_cell);
-
- const auto & q_points = scratch_data.get_quadrature_points();
- const unsigned int n_q_points = q_points.size();
- const std::vector<double> &JxW = scratch_data.get_JxW_values();
-
- std::vector<double> rhs(n_q_points);
- rhs_function->value_list(q_points, rhs);
-
- for (unsigned int point = 0; point < n_q_points; ++point)
- for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
- {
- for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
- copy_data.cell_matrix(i, j) +=
- diffusion_coefficient * // nu
- fe_v.shape_grad(i, point) * // grad v_h
- fe_v.shape_grad(j, point) * // grad u_h
- JxW[point]; // dx
-
- copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h
- rhs[point] * // f
- JxW[point]; // dx
- }
- };
+ const auto cell_worker =
+ [&](const auto &cell, auto &scratch_data, auto ©_data) {
+ const FEValues<dim> &fe_v = scratch_data.reinit(cell);
+ const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
+ copy_data.reinit(cell, dofs_per_cell);
+
+ const auto & q_points = scratch_data.get_quadrature_points();
+ const unsigned int n_q_points = q_points.size();
+ const std::vector<double> &JxW = scratch_data.get_JxW_values();
+
+ std::vector<double> rhs(n_q_points);
+ rhs_function->value_list(q_points, rhs);
+
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
+ {
+ for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
+ copy_data.cell_matrix(i, j) +=
+ diffusion_coefficient * // nu
+ fe_v.shape_grad(i, point) * // grad v_h
+ fe_v.shape_grad(j, point) * // grad u_h
+ JxW[point]; // dx
+
+ copy_data.cell_rhs(i) += fe_v.shape_value(i, point) * // v_h
+ rhs[point] * // f
+ JxW[point]; // dx
+ }
+ };
// This function assembles face integrals on the boundary.
- auto boundary_worker = [&](const Iterator & cell,
- const unsigned int &face_no,
- ScratchData & scratch_data,
- CopyData & copy_data) {
+ const auto boundary_worker = [&](const auto & cell,
+ const unsigned int &face_no,
+ auto & scratch_data,
+ auto & copy_data) {
const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
const auto & q_points = scratch_data.get_quadrature_points();
const double extent1 = cell->extent_in_direction(
GeometryInfo<dim>::unit_normal_direction[face_no]);
- const double penalty = compute_penalty(fe.get_degree(), extent1, extent1);
+ const double penalty = compute_penalty(degree, extent1, extent1);
for (unsigned int point = 0; point < n_q_points; ++point)
{
// To reinitialize FEInterfaceValues, we need to pass cells,
// face and subface indices (for adaptive refinement)
// to the reinit() function of FEInterfaceValues.
- auto face_worker = [&](const Iterator & cell,
- const unsigned int &f,
- const unsigned int &sf,
- const Iterator & ncell,
- const unsigned int &nf,
- const unsigned int &nsf,
- ScratchData & scratch_data,
- CopyData & copy_data) {
+ const auto face_worker = [&](const auto & cell,
+ const unsigned int &f,
+ const unsigned int &sf,
+ const auto & ncell,
+ const unsigned int &nf,
+ const unsigned int &nsf,
+ auto & scratch_data,
+ auto & copy_data) {
const FEInterfaceValues<dim> &fe_iv =
scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
cell->extent_in_direction(GeometryInfo<dim>::unit_normal_direction[f]);
const double extent2 = ncell->extent_in_direction(
GeometryInfo<dim>::unit_normal_direction[nf]);
- const double penalty = compute_penalty(fe.get_degree(), extent1, extent2);
+ const double penalty = compute_penalty(degree, extent1, extent2);
for (unsigned int point = 0; point < n_q_points; ++point)
{
// the global matrix and right-hand side.
// Though there are no hanging node constraints in DG discretization,
// we define an empty AffineConstraints oject that
- // allows us to use copy_local_to_global functionality.
+ // allows us to use distribute_local_to_global functionality.
AffineConstraints<double> constraints;
constraints.close();
- auto copier = [&](const CopyData &c) {
+ const auto copier = [&](const auto &c) {
constraints.distribute_local_to_global(c.cell_matrix,
c.cell_rhs,
c.local_dof_indices,
// Copy data from interior face assembly to the global matrix.
for (auto &cdf : c.face_data)
{
- const unsigned int joint_dofs_per_face = cdf.joint_dof_indices.size();
- for (unsigned int i = 0; i < joint_dofs_per_face; ++i)
- for (unsigned int k = 0; k < joint_dofs_per_face; ++k)
- system_matrix.add(cdf.joint_dof_indices[i],
- cdf.joint_dof_indices[k],
- cdf.cell_matrix(i, k));
+ constraints.distribute_local_to_global(cdf.cell_matrix,
+ cdf.joint_dof_indices,
+ system_matrix);
}
};
// and pass them together with the lambda functions
// above to MeshWorker::mesh_loop. In addition, we
// need to specify that we want to assemble interior faces once.
- const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
- QGauss<dim> quadrature(n_gauss_points);
- QGauss<dim - 1> face_quadrature(n_gauss_points);
UpdateFlags cell_flags = update_values | update_gradients |
update_quadrature_points | update_JxW_values;
template <int dim>
void SIPGLaplace<dim>::solve()
{
- std::cout << " Solving system..." << std::endl;
SparseDirectUMFPACK A_direct;
A_direct.initialize(system_matrix);
A_direct.vmult(solution, system_rhs);
template <int dim>
void SIPGLaplace<dim>::output_results(const unsigned int cycle) const
{
- std::string filename = "sol_Q" +
- Utilities::int_to_string(fe.get_degree(), 1) + "-" +
+ std::string filename = "sol_Q" + Utilities::int_to_string(degree, 1) + "-" +
Utilities::int_to_string(cycle, 2) + ".vtu";
- std::cout << "Writing solution to <" << filename << ">" << std::endl;
std::ofstream output(filename);
DataOut<dim> data_out;
template <int dim>
void SIPGLaplace<dim>::compute_error_estimate()
{
- typedef decltype(dof_handler.begin_active()) Iterator;
estimated_error_square_per_cell.reinit(triangulation.n_active_cells());
// Assemble cell residual $h_K^2 \left\| f + \nu \Delta u_h \right\|_K^2$.
- auto cell_worker = [&](const Iterator &cell,
- ScratchData & scratch_data,
- CopyData & copy_data) {
- const FEValues<dim> &fe_v = scratch_data.reinit(cell);
+ const auto cell_worker =
+ [&](const auto &cell, auto &scratch_data, auto ©_data) {
+ const FEValues<dim> &fe_v = scratch_data.reinit(cell);
- copy_data.cell_index = cell->active_cell_index();
+ copy_data.cell_index = cell->active_cell_index();
- const auto & q_points = fe_v.get_quadrature_points();
- const unsigned int n_q_points = q_points.size();
- const std::vector<double> &JxW = fe_v.get_JxW_values();
+ const auto & q_points = fe_v.get_quadrature_points();
+ const unsigned int n_q_points = q_points.size();
+ const std::vector<double> &JxW = fe_v.get_JxW_values();
- std::vector<Tensor<2, dim>> hessians(n_q_points);
- fe_v.get_function_hessians(solution, hessians);
+ std::vector<Tensor<2, dim>> hessians(n_q_points);
+ fe_v.get_function_hessians(solution, hessians);
- std::vector<double> rhs(n_q_points);
- rhs_function->value_list(q_points, rhs);
+ std::vector<double> rhs(n_q_points);
+ rhs_function->value_list(q_points, rhs);
- const double hk = cell->diameter();
- double residual_norm_square = 0;
+ const double hk = cell->diameter();
+ double residual_norm_square = 0;
- for (unsigned int point = 0; point < n_q_points; ++point)
- {
- const double residual =
- rhs[point] + diffusion_coefficient * trace(hessians[point]);
- residual_norm_square += residual * residual * JxW[point];
- }
- copy_data.value = hk * hk * residual_norm_square;
- };
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ {
+ const double residual =
+ rhs[point] + diffusion_coefficient * trace(hessians[point]);
+ residual_norm_square += residual * residual * JxW[point];
+ }
+ copy_data.value = hk * hk * residual_norm_square;
+ };
// Assemble boundary terms $\sum_{f\in \partial K \cap \partial \Omega}
// \sigma \left\| [ u_h-g_D ] \right\|_f^2 $.
- auto boundary_worker = [&](const Iterator & cell,
- const unsigned int &face_no,
- ScratchData & scratch_data,
- CopyData & copy_data) {
+ const auto boundary_worker = [&](const auto & cell,
+ const unsigned int &face_no,
+ auto & scratch_data,
+ auto & copy_data) {
const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
const auto & q_points = fe_fv.get_quadrature_points();
const double extent1 = cell->extent_in_direction(
GeometryInfo<dim>::unit_normal_direction[face_no]);
- const double penalty = compute_penalty(fe.get_degree(), extent1, extent1);
+ const double penalty = compute_penalty(degree, extent1, extent1);
double difference_norm_square = 0.;
for (unsigned int point = 0; point < q_points.size(); ++point)
// Assemble interior face terms $\sum_{f\in \partial K}\lbrace \sigma
// \left\| [u_h] \right\|_f^2 + h_f \left\| [\nu \nabla u_h \cdot
// \mathbf n ] \right\|_f^2 \rbrace$.
- auto face_worker = [&](const Iterator & cell,
- const unsigned int &f,
- const unsigned int &sf,
- const Iterator & ncell,
- const unsigned int &nf,
- const unsigned int &nsf,
- ScratchData & scratch_data,
- CopyData & copy_data) {
+ const auto face_worker = [&](const auto & cell,
+ const unsigned int &f,
+ const unsigned int &sf,
+ const auto & ncell,
+ const unsigned int &nf,
+ const unsigned int &nsf,
+ auto & scratch_data,
+ auto & copy_data) {
const FEInterfaceValues<dim> &fe_iv =
scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
cell->extent_in_direction(GeometryInfo<dim>::unit_normal_direction[f]);
const double extent2 = ncell->extent_in_direction(
GeometryInfo<dim>::unit_normal_direction[nf]);
- const double penalty = compute_penalty(fe.get_degree(), extent1, extent2);
+ const double penalty = compute_penalty(degree, extent1, extent2);
double flux_jump_square = 0;
double u_jump_square = 0;
copy_data_face.values[1] = copy_data_face.values[0];
};
- auto copier = [&](const CopyData ©_data) {
+ const auto copier = [&](const auto ©_data) {
if (copy_data.cell_index != numbers::invalid_unsigned_int)
estimated_error_square_per_cell[copy_data.cell_index] +=
copy_data.value;
estimated_error_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
};
- const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
- QGauss<dim> quadrature(n_gauss_points);
- QGauss<dim - 1> face_quadrature(n_gauss_points);
-
UpdateFlags cell_flags =
update_hessians | update_quadrature_points | update_JxW_values;
UpdateFlags face_flags = update_values | update_gradients |
template <int dim>
double SIPGLaplace<dim>::compute_energy_norm()
{
- typedef decltype(dof_handler.begin_active()) Iterator;
energy_norm_square_per_cell.reinit(triangulation.n_active_cells());
- auto cell_worker = [&](const Iterator &cell,
- ScratchData & scratch_data,
- CopyData & copy_data) {
- const FEValues<dim> &fe_v = scratch_data.reinit(cell);
+ const auto cell_worker =
+ [&](const auto &cell, auto &scratch_data, auto ©_data) {
+ const FEValues<dim> &fe_v = scratch_data.reinit(cell);
- copy_data.cell_index = cell->active_cell_index();
+ copy_data.cell_index = cell->active_cell_index();
- const auto & q_points = fe_v.get_quadrature_points();
- const unsigned int n_q_points = q_points.size();
- const std::vector<double> &JxW = fe_v.get_JxW_values();
+ const auto & q_points = fe_v.get_quadrature_points();
+ const unsigned int n_q_points = q_points.size();
+ const std::vector<double> &JxW = fe_v.get_JxW_values();
- std::vector<Tensor<1, dim>> grad_u(n_q_points);
- fe_v.get_function_gradients(solution, grad_u);
+ std::vector<Tensor<1, dim>> grad_u(n_q_points);
+ fe_v.get_function_gradients(solution, grad_u);
- std::vector<Tensor<1, dim>> grad_exact(n_q_points);
- exact_solution->gradient_list(q_points, grad_exact);
+ std::vector<Tensor<1, dim>> grad_exact(n_q_points);
+ exact_solution->gradient_list(q_points, grad_exact);
- double norm_square = 0;
- for (unsigned int point = 0; point < n_q_points; ++point)
- {
- norm_square +=
- (grad_u[point] - grad_exact[point]).norm_square() * JxW[point];
- }
- copy_data.value = norm_square;
- };
+ double norm_square = 0;
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ {
+ norm_square +=
+ (grad_u[point] - grad_exact[point]).norm_square() * JxW[point];
+ }
+ copy_data.value = norm_square;
+ };
- auto boundary_worker = [&](const Iterator & cell,
- const unsigned int &face_no,
- ScratchData & scratch_data,
- CopyData & copy_data) {
+ const auto boundary_worker = [&](const auto & cell,
+ const unsigned int &face_no,
+ auto & scratch_data,
+ auto & copy_data) {
const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
const auto & q_points = fe_fv.get_quadrature_points();
const double extent1 = cell->extent_in_direction(
GeometryInfo<dim>::unit_normal_direction[face_no]);
- const double penalty = compute_penalty(fe.get_degree(), extent1, extent1);
+ const double penalty = compute_penalty(degree, extent1, extent1);
double difference_norm_square = 0.;
for (unsigned int point = 0; point < q_points.size(); ++point)
copy_data.value += penalty * difference_norm_square;
};
- auto face_worker = [&](const Iterator & cell,
- const unsigned int &f,
- const unsigned int &sf,
- const Iterator & ncell,
- const unsigned int &nf,
- const unsigned int &nsf,
- ScratchData & scratch_data,
- CopyData & copy_data) {
+ const auto face_worker = [&](const auto & cell,
+ const unsigned int &f,
+ const unsigned int &sf,
+ const auto & ncell,
+ const unsigned int &nf,
+ const unsigned int &nsf,
+ auto & scratch_data,
+ auto & copy_data) {
const FEInterfaceValues<dim> &fe_iv =
scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
cell->extent_in_direction(GeometryInfo<dim>::unit_normal_direction[f]);
const double extent2 = ncell->extent_in_direction(
GeometryInfo<dim>::unit_normal_direction[nf]);
- const double penalty = compute_penalty(fe.get_degree(), extent1, extent2);
+ const double penalty = compute_penalty(degree, extent1, extent2);
double u_jump_square = 0;
for (unsigned int point = 0; point < n_q_points; ++point)
copy_data_face.values[1] = copy_data_face.values[0];
};
- auto copier = [&](const CopyData ©_data) {
+ const auto copier = [&](const auto ©_data) {
if (copy_data.cell_index != numbers::invalid_unsigned_int)
energy_norm_square_per_cell[copy_data.cell_index] += copy_data.value;
for (auto &cdf : copy_data.face_data)
energy_norm_square_per_cell[cdf.cell_indices[j]] += cdf.values[j];
};
- const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1;
- QGauss<dim> quadrature(n_gauss_points);
- QGauss<dim - 1> face_quadrature(n_gauss_points);
-
UpdateFlags cell_flags =
update_gradients | update_quadrature_points | update_JxW_values;
UpdateFlags face_flags =
update_values | update_quadrature_points | update_JxW_values;
- ScratchData scratch_data(
- mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
+ ScratchData scratch_data(mapping,
+ fe,
+ QGauss<dim>(fe.degree + 2),
+ cell_flags,
+ QGauss<dim - 1>(fe.degree + 2),
+ face_flags);
CopyData cd;
MeshWorker::mesh_loop(dof_handler.begin_active(),
const double energy_error = compute_energy_norm();
convergence_table.add_value("Energy", energy_error);
- std::cout << " Error in the L2 norm : " << L2_error << std::endl
- << " Error in the H1 seminorm : " << H1_error << std::endl
- << " Error in the energy norm : " << energy_error
+ std::cout << " Error in the L2 norm : " << L2_error << std::endl
+ << " Error in the H1 seminorm : " << H1_error << std::endl
+ << " Error in the energy norm : " << energy_error
<< std::endl;
}
template <int dim>
void SIPGLaplace<dim>::run()
{
- unsigned int max_cycle = test_case == Test_Case::convergence_rate ? 6 : 10;
+ unsigned int max_cycle = test_case == Test_Case::convergence_rate ? 6 : 20;
for (unsigned int cycle = 0; cycle < max_cycle; ++cycle)
{
std::cout << "Cycle " << cycle << std::endl;
if (cycle == 0)
{
GridGenerator::hyper_L(triangulation);
- triangulation.refine_global(2);
+ triangulation.refine_global(3);
}
else
{
Assert(false, ExcNotImplemented());
}
}
- std::cout << "Number of active cells: "
+ std::cout << " Number of active cells : "
<< triangulation.n_active_cells() << std::endl;
setup_system();
- std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs()
+ std::cout << " Number of degrees of freedom : " << dof_handler.n_dofs()
<< std::endl;
assemble_system();
if (test_case == Test_Case::l_singularity)
{
compute_error_estimate();
+ std::cout << " Estimated error : "
+ << std::sqrt(estimated_error_square_per_cell.l1_norm())
+ << std::endl;
+
convergence_table.add_value(
"Estimator",
std::sqrt(estimated_error_square_per_cell.l1_norm()));
convergence_table.evaluate_convergence_rates(
"H1", ConvergenceTable::reduction_rate_log2);
}
-
- std::cout << "degree = " << fe.get_degree() << std::endl;
+ std::cout << "degree = " << degree << std::endl;
convergence_table.write_text(
std::cout, TableHandler::TextOutputFormat::org_mode_table);
}