From 84483781743136274c7c075d10bcbe928ba3e38f Mon Sep 17 00:00:00 2001 From: Jiaqi Zhang Date: Mon, 9 Nov 2020 10:24:47 -0500 Subject: [PATCH] update --- doc/doxygen/references.bib | 2 +- examples/step-74/doc/intro.dox | 56 +++-- examples/step-74/doc/results.dox | 288 +++++++++++-------------- examples/step-74/doc/tooltip | 2 +- examples/step-74/step-74.cc | 351 +++++++++++++++---------------- 5 files changed, 334 insertions(+), 365 deletions(-) diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index 9421a8205c..d920803a5e 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -830,7 +830,7 @@ year = {2008}, @book{di2011mathematical, title={Mathematical aspects of discontinuous Galerkin methods}, - author={Di Pietro and Daniele Antonio and Ern, Alexandre}, + author={Di Pietro, Daniele Antonio and Ern, Alexandre}, volume={69}, year={2011}, publisher={Springer Science \& Business Media} diff --git a/examples/step-74/doc/intro.dox b/examples/step-74/doc/intro.dox index 54a00110b4..90965c30cf 100644 --- a/examples/step-74/doc/intro.dox +++ b/examples/step-74/doc/intro.dox @@ -28,7 +28,7 @@ This tutorial includes the following topics.

The equation

In this example, we consider Poisson's equation @f[ -- \nu \Delta u = f \qquad \mbox{in } \Omega, +- \nabla \cdot \left( \nu \nabla u\right) = f \qquad \mbox{in } \Omega, @f] subject to the boundary condition @f[ @@ -38,8 +38,8 @@ For simplicity, we assume that the diffusion coefficient $\nu$ is constant here. Note that if $\nu$ is discontinuous, we need to take this into account when computing jump terms on cell faces. -We denote the mesh by $\Gamma_h$, and $K\in\Gamma_h$ is a mesh cell. -The sets of interior and boundary faces are denoted by $F^i_h$ and $F^b_h$ +We denote the mesh by ${\mathbb T}_h$, and $K\in{\mathbb T}_h$ is a mesh cell. +The sets of interior and boundary faces are denoted by ${\mathbb F}^i_h$ and ${\mathbb F}^b_h$ respectively. Let $K^0$ and $K^1$ be the two cells sharing a face $f\in F_h^i$, and $\mathbf n$ be the outer normal vector of $K^0$. Then the jump and average operators are given by @@ -55,24 +55,24 @@ $\average{v}=v$. The discretization using the SIPG is given by the following weak formula (more details can be found in @cite di2011mathematical and the references therein) @f{multline*} - \sum_{K\in \Gamma_h} (\nabla v_h, \nu \nabla u_h)_K + \sum_{K\in {\mathbb T}_h} (\nabla v_h, \nu \nabla u_h)_K \\ - - \sum_{F \in F_h^i} \biggl\{ - \bigl< \jump{v_h}, \nu\average{ \nabla u_h} \cdot \mathbf n \bigr>_F - +\bigl<\average{ \nabla v_h }\cdot \mathbf n,\nu\jump{u_h}\bigr>_F - -\bigl<\jump{v_h},\nu \sigma \jump{u_h} \bigr>_F - \biggr\} + - \sum_{F \in F_h^i} \left\{ + \left< \jump{v_h}, \nu\average{ \nabla u_h} \cdot \mathbf n \right>_F + +\left<\average{ \nabla v_h }\cdot \mathbf n,\nu\jump{u_h}\right>_F + -\left<\jump{v_h},\nu \sigma \jump{u_h} \right>_F + \right\} \\ - - \sum_{F \in F_h^b} \biggl\{ - \bigl_F - + \bigl< \nabla v_h \cdot \mathbf n , \nu u_h\bigr>_F - - \bigl< v_h,\nu \sigma u_h\bigr>_F - \biggr\} + - \sum_{F \in F_h^b} \left\{ + \left_F + + \left< \nabla v_h \cdot \mathbf n , \nu u_h\right>_F + - \left< v_h,\nu \sigma u_h\right>_F + \right\} \\ = (v_h, f)_\Omega - - \sum_{F \in F_h^b} \biggl\{ - \bigl< \nabla v_h \cdot \mathbf n, \nu g_D\bigr>_F - \bigl_F - \biggr\}. + - \sum_{F \in F_h^b} \left\{ + \left< \nabla v_h \cdot \mathbf n, \nu g_D\right>_F - \left_F + \right\}. @f}

The penalty parameter

@@ -85,9 +85,9 @@ One can just pick a large constant, while other options could be the multiples o we follow step-39 and use $\gamma = p(p+1)$.

A posteriori error estimator

-In this example, we use the error estimator by Karakashian and Pascal @cite karakashian2003posteriori +In this example, we use the error estimator by Karakashian and Pascal @cite karakashian2003posteriori with a slight modification @f[ -\eta^2 = \sum_{K \in \Gamma_h} \eta^2_{K} + \sum_{f_i \in F^i_h} \eta^2_{f_i} + \sum_{f_b \in F^i_b}\eta^2_{f_b} +\eta^2 = \sum_{K \in {\mathbb T}_h} \eta^2_{K} + \sum_{f_i \in {\mathbb F}^i_h} \eta^2_{f_i} + \sum_{f_b \in F^i_b}\eta^2_{f_b} @f] where @f[ @@ -99,7 +99,7 @@ where @f[ \eta_{f_b}^2 = \sigma \left\| u_h-g_D \right\|_f^2 @f] -The only difference is that we use $\sigma = \gamma/h_f$ instead of $\gamma^2/h_f$ for the jump terms of $u$ (the first term in $\eta^2_{f_i}$ and $\eta_{f_b}^2$). +Here we use $\sigma = \gamma/h_f$ instead of $\gamma^2/h_f$ for the jump terms of $u_h$ (the first term in $\eta^2_{f_i}$ and $\eta_{f_b}^2$). In each cell $K$, we compute @f[ @@ -117,3 +117,19 @@ Then the error estimate square per cell is @f] Note that we compute $\eta_{local}^2$ instead of $\eta_{local}$ to simplify the implementation. The error estimate square per cell is stored in a global vector, whose $L_1$ norm is equal to $\eta^2$. + +

The test case

+In the first test problem, we run a convergence test using a smooth manufactured solution with $\nu =1$ in 2D +@f[ +u=\sin(2\pi x)\sin(2\pi y), (x,y)\in (0,1)\times (0,1), +@f] +correspondingly, +@f[ +u=0, \qquad \mbox{on } \partial \Omega. +@f] +and $f= 8\pi^2 u$. We compute errors against the manufactured solution and evaluate the convergence rate. + +In the second test, we choose Functions::LSingularityFunction on a L-shaped domain (GridGenerator::hyper_L) in 2D. +The solution is given in the polar coordinates by $u(r,\phi) = r^{\frac{2}{3}}\sin \left(\frac{2}{3}\phi \right)$, +which has a singularity at the origin. An error estimator is constructed to detect the region with large errors, +according to which the mesh is refined adaptively. diff --git a/examples/step-74/doc/results.dox b/examples/step-74/doc/results.dox index 9896f504e4..b14ec3acfd 100644 --- a/examples/step-74/doc/results.dox +++ b/examples/step-74/doc/results.dox @@ -3,184 +3,132 @@ 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 - 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 - 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 - 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 - 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 - 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 - 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: + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
cyclen_cellssn_dofsL2 rateH1rateEnergy
0162561.933e-03 1.061e-01 1.506e-01
16410249.605e-054.338.995e-033.561.133e-02
225640965.606e-064.109.018e-043.329.736e-04
31024163843.484e-074.011.071e-043.071.088e-04
44096655362.179e-084.001.327e-053.011.331e-05
5163842621441.363e-094.001.656e-063.001.657e-06
+ 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 - 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 @@ -193,3 +141,9 @@ and one order lower than $L_2$ error, which illustrates its ability to predict regions with large errors. + +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 diff --git a/examples/step-74/doc/tooltip b/examples/step-74/doc/tooltip index 321d3af5f4..fbb04bcfc1 100644 --- a/examples/step-74/doc/tooltip +++ b/examples/step-74/doc/tooltip @@ -1 +1 @@ -Symmetric interior penalty Galerkin for Poisson's equation. \ No newline at end of file +Symmetric interior penalty Galerkin for Poisson's equation. diff --git a/examples/step-74/step-74.cc b/examples/step-74/step-74.cc index 04080a4bc9..9d42c971ae 100644 --- a/examples/step-74/step-74.cc +++ b/examples/step-74/step-74.cc @@ -1,6 +1,6 @@ /* --------------------------------------------------------------------- * - * 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. * @@ -40,7 +40,6 @@ #include // Here the discontinuous finite elements and FEInterfaceValues are defined. #include -#include #include #include @@ -115,15 +114,21 @@ namespace Step74 {} virtual void value_list(const std::vector> &points, std::vector & 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 + void + SmoothRightHandSide::value_list(const std::vector> &points, + std::vector & 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 @@ -135,17 +140,23 @@ namespace Step74 {} virtual void value_list(const std::vector> &points, std::vector & 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 + void + SingularRightHandSide::value_list(const std::vector> &points, + std::vector & 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 @@ -155,8 +166,8 @@ namespace Step74 const Vector & solution, std::vector & jump) { - const unsigned n_q = fe_iv.n_quadrature_points; - std::vector face_values[2]; + const unsigned n_q = fe_iv.n_quadrature_points; + std::array, 2> face_values; jump.resize(n_q); for (unsigned i = 0; i < 2; ++i) { @@ -206,8 +217,8 @@ namespace Step74 { FullMatrix cell_matrix; std::vector joint_dof_indices; - double values[2]; - unsigned int cell_indices[2]; + std::array values; + std::array cell_indices; }; struct CopyData @@ -252,6 +263,9 @@ namespace Step74 double compute_energy_norm(); Triangulation triangulation; + const unsigned degree; + QGauss quadrature; + QGauss face_quadrature; const MappingQ1 mapping; using ScratchData = MeshWorker::ScratchData; @@ -283,12 +297,14 @@ namespace Step74 }; // 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 SIPGLaplace::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) { @@ -329,44 +345,41 @@ namespace Step74 template void SIPGLaplace::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 &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 &JxW = scratch_data.get_JxW_values(); - - std::vector 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 &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 &JxW = scratch_data.get_JxW_values(); + + std::vector 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 &fe_fv = scratch_data.reinit(cell, face_no); const auto & q_points = scratch_data.get_quadrature_points(); @@ -382,7 +395,7 @@ namespace Step74 const double extent1 = cell->extent_in_direction( GeometryInfo::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) { @@ -426,14 +439,14 @@ namespace Step74 // 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 &fe_iv = scratch_data.reinit(cell, f, sf, ncell, nf, nsf); @@ -453,7 +466,7 @@ namespace Step74 cell->extent_in_direction(GeometryInfo::unit_normal_direction[f]); const double extent2 = ncell->extent_in_direction( GeometryInfo::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) { @@ -483,10 +496,10 @@ namespace Step74 // 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 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, @@ -496,12 +509,9 @@ namespace Step74 // 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); } }; @@ -509,9 +519,6 @@ namespace Step74 // 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 quadrature(n_gauss_points); - QGauss face_quadrature(n_gauss_points); UpdateFlags cell_flags = update_values | update_gradients | update_quadrature_points | update_JxW_values; @@ -538,7 +545,6 @@ namespace Step74 template void SIPGLaplace::solve() { - std::cout << " Solving system..." << std::endl; SparseDirectUMFPACK A_direct; A_direct.initialize(system_matrix); A_direct.vmult(solution, system_rhs); @@ -547,10 +553,8 @@ namespace Step74 template void SIPGLaplace::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 data_out; @@ -565,45 +569,43 @@ namespace Step74 template void SIPGLaplace::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 &fe_v = scratch_data.reinit(cell); + const auto cell_worker = + [&](const auto &cell, auto &scratch_data, auto ©_data) { + const FEValues &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 &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 &JxW = fe_v.get_JxW_values(); - std::vector> hessians(n_q_points); - fe_v.get_function_hessians(solution, hessians); + std::vector> hessians(n_q_points); + fe_v.get_function_hessians(solution, hessians); - std::vector rhs(n_q_points); - rhs_function->value_list(q_points, rhs); + std::vector 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 &fe_fv = scratch_data.reinit(cell, face_no); const auto & q_points = fe_fv.get_quadrature_points(); @@ -619,7 +621,7 @@ namespace Step74 const double extent1 = cell->extent_in_direction( GeometryInfo::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) @@ -633,14 +635,14 @@ namespace Step74 // 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 &fe_iv = scratch_data.reinit(cell, f, sf, ncell, nf, nsf); @@ -668,7 +670,7 @@ namespace Step74 cell->extent_in_direction(GeometryInfo::unit_normal_direction[f]); const double extent2 = ncell->extent_in_direction( GeometryInfo::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; @@ -684,7 +686,7 @@ namespace Step74 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; @@ -693,10 +695,6 @@ namespace Step74 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 quadrature(n_gauss_points); - QGauss face_quadrature(n_gauss_points); - UpdateFlags cell_flags = update_hessians | update_quadrature_points | update_JxW_values; UpdateFlags face_flags = update_values | update_gradients | @@ -725,39 +723,37 @@ namespace Step74 template double SIPGLaplace::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 &fe_v = scratch_data.reinit(cell); + const auto cell_worker = + [&](const auto &cell, auto &scratch_data, auto ©_data) { + const FEValues &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 &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 &JxW = fe_v.get_JxW_values(); - std::vector> grad_u(n_q_points); - fe_v.get_function_gradients(solution, grad_u); + std::vector> grad_u(n_q_points); + fe_v.get_function_gradients(solution, grad_u); - std::vector> grad_exact(n_q_points); - exact_solution->gradient_list(q_points, grad_exact); + std::vector> 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 &fe_fv = scratch_data.reinit(cell, face_no); const auto & q_points = fe_fv.get_quadrature_points(); @@ -773,7 +769,7 @@ namespace Step74 const double extent1 = cell->extent_in_direction( GeometryInfo::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) @@ -784,14 +780,14 @@ namespace Step74 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 &fe_iv = scratch_data.reinit(cell, f, sf, ncell, nf, nsf); @@ -813,7 +809,7 @@ namespace Step74 cell->extent_in_direction(GeometryInfo::unit_normal_direction[f]); const double extent2 = ncell->extent_in_direction( GeometryInfo::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) @@ -824,7 +820,7 @@ namespace Step74 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) @@ -832,17 +828,17 @@ namespace Step74 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 quadrature(n_gauss_points); - QGauss 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(fe.degree + 2), + cell_flags, + QGauss(fe.degree + 2), + face_flags); CopyData cd; MeshWorker::mesh_loop(dof_handler.begin_active(), @@ -914,16 +910,16 @@ namespace Step74 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 void SIPGLaplace::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; @@ -949,7 +945,7 @@ namespace Step74 if (cycle == 0) { GridGenerator::hyper_L(triangulation); - triangulation.refine_global(2); + triangulation.refine_global(3); } else { @@ -962,11 +958,11 @@ namespace Step74 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(); @@ -982,6 +978,10 @@ namespace Step74 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())); @@ -1009,8 +1009,7 @@ namespace Step74 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); } -- 2.39.5