From: Wolfgang Bangerth Date: Thu, 9 May 2024 16:28:59 +0000 (+0530) Subject: Zap step-12b. X-Git-Tag: v9.6.0-rc1~269^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c3e046301479b76302de46bc5086e8cbbcf454c6;p=dealii.git Zap step-12b. --- diff --git a/doc/doxygen/headers/io.h b/doc/doxygen/headers/io.h index f60fe51b8e..45a64e78e1 100644 --- a/doc/doxygen/headers/io.h +++ b/doc/doxygen/headers/io.h @@ -108,7 +108,7 @@ * documentation of that class lists the supported formats. * * Several of the tutorial programs, notably step-1, step-6, step-9, step-10, - * step-12b, and step-14 demonstrate the use of the GridOut class. + * and step-14 demonstrate the use of the GridOut class. * * *

Matrix output

diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index 89351c0cd5..8ea384b5a9 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -201,12 +201,6 @@ * * * - * step-12b - * A variant of step-12 (for linear advection problems), but - * using the MeshWorker framework for the assembly of linear systems. - * - * - * * step-13 * Software design questions and * how to write a modular, extensible finite element program. @@ -939,7 +933,6 @@ * * * step-12, - * step-12b, * step-21, * step-39, * step-46, diff --git a/doc/news/9.1.1-vs-9.2.0.h b/doc/news/9.1.1-vs-9.2.0.h index b4075c1ab4..867121b5fa 100644 --- a/doc/news/9.1.1-vs-9.2.0.h +++ b/doc/news/9.1.1-vs-9.2.0.h @@ -734,7 +734,8 @@ inconvenience this causes.
  • New: The step-12 tutorial program has been changed to use - FEInterfaceValues. The old version of step-12 is still available as step-12b. + FEInterfaceValues. The old version of step-12 was still available + until deal.II 9.5 as step 12b.
    (Timo Heister, 2019/08/27)
    • diff --git a/examples/step-12/doc/intro.dox b/examples/step-12/doc/intro.dox index 5e55e023a3..57de566848 100644 --- a/examples/step-12/doc/intro.dox +++ b/examples/step-12/doc/intro.dox @@ -1,10 +1,5 @@
    - Note: A variant called step-12b of this tutorial exists, using -MeshWorker and LocalIntegrators instead of assembling matrices using -FEInterfaceValues as is done in this tutorial. - -

    An example of an advection problem using the Discountinuous Galerkin method

    diff --git a/examples/step-12b/CMakeLists.txt b/examples/step-12b/CMakeLists.txt deleted file mode 100644 index f0be3ec66a..0000000000 --- a/examples/step-12b/CMakeLists.txt +++ /dev/null @@ -1,39 +0,0 @@ -## -# CMake script for the step-12b tutorial program: -## - -# Set the name of the project and target: -set(TARGET "step-12b") - -# Declare all source files the target consists of. Here, this is only -# the one step-X.cc file, but as you expand your project you may wish -# to add other source files as well. If your project becomes much larger, -# you may want to either replace the following statement by something like -# file(GLOB_RECURSE TARGET_SRC "source/*.cc") -# file(GLOB_RECURSE TARGET_INC "include/*.h") -# set(TARGET_SRC ${TARGET_SRC} ${TARGET_INC}) -# or switch altogether to the large project CMakeLists.txt file discussed -# in the "CMake in user projects" page accessible from the "User info" -# page of the documentation. -set(TARGET_SRC - ${TARGET}.cc - ) - -# Usually, you will not need to modify anything beyond this point... - -cmake_minimum_required(VERSION 3.13.4) - -find_package(deal.II 9.6.0 - HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} - ) -if(NOT ${deal.II_FOUND}) - message(FATAL_ERROR "\n" - "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" - "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" - "or set an environment variable \"DEAL_II_DIR\" that contains this path." - ) -endif() - -deal_ii_initialize_cached_variables() -project(${TARGET}) -deal_ii_invoke_autopilot() diff --git a/examples/step-12b/doc/builds-on b/examples/step-12b/doc/builds-on deleted file mode 100644 index 943473db87..0000000000 --- a/examples/step-12b/doc/builds-on +++ /dev/null @@ -1 +0,0 @@ -step-7 diff --git a/examples/step-12b/doc/intro.dox b/examples/step-12b/doc/intro.dox deleted file mode 100644 index 06effea3a1..0000000000 --- a/examples/step-12b/doc/intro.dox +++ /dev/null @@ -1,13 +0,0 @@ - -

    Introduction

    - -This is a variant of step-12 with the only change that we are using the -MeshWorker framework with the pre-made LocalIntegrator helper classes instead -of assembling the face terms using FEInterfaceValues. - -The details of this framework on how it is used in practice will be explained -as part of this tutorial program. - -

    The testcase

    - -The problem we solve here is the same as the one in step-12. diff --git a/examples/step-12b/doc/kind b/examples/step-12b/doc/kind deleted file mode 100644 index c1d9154931..0000000000 --- a/examples/step-12b/doc/kind +++ /dev/null @@ -1 +0,0 @@ -techniques diff --git a/examples/step-12b/doc/results.dox b/examples/step-12b/doc/results.dox deleted file mode 100644 index ee39c3ffe5..0000000000 --- a/examples/step-12b/doc/results.dox +++ /dev/null @@ -1,32 +0,0 @@ -

    Results

    - - -The output of this program is very similar to step-16 and we are not repeating the output here. - -We show the solutions on the initial mesh, the mesh after two -and after five adaptive refinement steps. - - - - - - -Then we show the final grid (after 5 refinement steps) and the solution again, -this time with a nicer 3d rendering (obtained using the DataOutBase::write_vtk -function and the VTK-based VisIt visualization program) that better shows the -sharpness of the jump on the refined mesh and the over- and undershoots of the -solution along the interface: - - - - - -And finally we show a plot of a 3d computation. - - - - - -

    Possibilities for extensions

    - -For ideas for further extensions, please see see step-12. diff --git a/examples/step-12b/doc/tooltip b/examples/step-12b/doc/tooltip deleted file mode 100644 index 06e84e3b17..0000000000 --- a/examples/step-12b/doc/tooltip +++ /dev/null @@ -1 +0,0 @@ -Discontinuous Galerkin for linear advection, MeshWorker version. diff --git a/examples/step-12b/step-12b.cc b/examples/step-12b/step-12b.cc deleted file mode 100644 index 2bfb222e79..0000000000 --- a/examples/step-12b/step-12b.cc +++ /dev/null @@ -1,674 +0,0 @@ -/* ------------------------------------------------------------------------ - * - * SPDX-License-Identifier: LGPL-2.1-or-later - * Copyright (C) 2009 - 2024 by the deal.II authors - * - * This file is part of the deal.II library. - * - * Part of the source code is dual licensed under Apache-2.0 WITH - * LLVM-exception OR LGPL-2.1-or-later. Detailed license information - * governing the source code and code contributions can be found in - * LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. - * - * ------------------------------------------------------------------------ - * - * Author: Guido Kanschat, Texas A&M University, 2009 - */ - - -// The first few files have already been covered in previous examples and will -// thus not be further commented on: -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -// Here the discontinuous finite elements are defined. They are used in the same -// way as all other finite elements, though -- as you have seen in previous -// tutorial programs -- there isn't much user interaction with finite element -// classes at all: they are passed to DoFHandler and -// FEValues objects, and that is about it. -#include -// We are going to use the simplest possible solver, called Richardson -// iteration, that represents a simple defect correction. This, in combination -// with a block SSOR preconditioner (defined in precondition_block.h), that -// uses the special block matrix structure of system matrices arising from DG -// discretizations. -#include -#include -// We are going to use gradients as refinement indicator. -#include - -// Here come the new include files for using the MeshWorker framework. The first -// contains the class MeshWorker::DoFInfo, which provides local integrators with -// a mapping between local and global degrees of freedom. It stores the results -// of local integrals as well in its base class MeshWorker::LocalResults. -// In the second of these files, we find an object of type -// MeshWorker::IntegrationInfo, which is mostly a wrapper around a group of -// FEValues objects. The file meshworker/simple.h contains classes -// assembling locally integrated data into a global system containing only a -// single matrix. Finally, we will need the file that runs the loop over all -// mesh cells and faces. -#include -#include -#include -#include - -// Like in all programs, we finish this section by including the needed C++ -// headers and declaring we want to use objects in the dealii namespace without -// prefix. -#include -#include - - -namespace Step12 -{ - using namespace dealii; - - // @sect3{Equation data} - // - // First, we define a class describing the inhomogeneous boundary data. Since - // only its values are used, we implement value_list(), but leave all other - // functions of Function undefined. - template - class BoundaryValues : public Function - { - public: - BoundaryValues() = default; - virtual void value_list(const std::vector> &points, - std::vector &values, - const unsigned int component = 0) const override; - }; - - // Given the flow direction, the inflow boundary of the unit square $[0,1]^2$ - // are the right and the lower boundaries. We prescribe discontinuous boundary - // values 1 and 0 on the x-axis and value 0 on the right boundary. The values - // of this function on the outflow boundaries will not be used within the DG - // scheme. - template - void BoundaryValues::value_list(const std::vector> &points, - std::vector &values, - const unsigned int component) const - { - (void)component; - AssertIndexRange(component, 1); - AssertDimension(values.size(), points.size()); - - for (unsigned int i = 0; i < values.size(); ++i) - { - if (points[i][0] < 0.5) - values[i] = 1.; - else - values[i] = 0.; - } - } - - - // Finally, a function that computes and returns the wind field - // $\beta=\beta(\mathbf x)$. As explained in the introduction, we will use a - // rotational field around the origin in 2d. In 3d, we simply leave the - // $z$-component unset (i.e., at zero), whereas the function can not be used - // in 1d in its current implementation: - template - Tensor<1, dim> beta(const Point &p) - { - Assert(dim >= 2, ExcNotImplemented()); - - Tensor<1, dim> wind_field; - wind_field[0] = -p[1]; - wind_field[1] = p[0]; - wind_field /= wind_field.norm(); - - return wind_field; - } - - - // @sect3{The AdvectionProblem class} - // - // After this preparations, we proceed with the main class of this program, - // called AdvectionProblem. It is basically the main class of step-6. We do - // not have an AffineConstraints object, because there are no hanging node - // constraints in DG discretizations. - - // Major differences will only come up in the implementation of the assemble - // functions, since here, we not only need to cover the flux integrals over - // faces, we also use the MeshWorker interface to simplify the loops - // involved. - template - class AdvectionProblem - { - public: - AdvectionProblem(); - void run(); - - private: - void setup_system(); - void assemble_system(); - void solve(Vector &solution); - void refine_grid(); - void output_results(const unsigned int cycle) const; - - Triangulation triangulation; - const MappingQ1 mapping; - - // Furthermore we want to use DG elements of degree 1 (but this is only - // specified in the constructor). If you want to use a DG method of a - // different degree the whole program stays the same, only replace 1 in - // the constructor by the desired polynomial degree. - const FE_DGQ fe; - DoFHandler dof_handler; - - // The next four members represent the linear system to be solved. - // system_matrix and right_hand_side are generated - // by assemble_system(), the solution is computed - // in solve(). The sparsity_pattern is used to - // determine the location of nonzero elements in system_matrix. - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector right_hand_side; - - // Finally, we have to provide functions that assemble the cell, boundary, - // and inner face terms. Within the MeshWorker framework, the loop over all - // cells and much of the setup of operations will be done outside this - // class, so all we have to provide are these three operations. They will - // then work on intermediate objects for which first, we here define - // alias to the info objects handed to the local integration functions - // in order to make our life easier below. - using DoFInfo = MeshWorker::DoFInfo; - using CellInfo = MeshWorker::IntegrationInfo; - - // The following three functions are then the ones that get called inside - // the generic loop over all cells and faces. They are the ones doing the - // actual integration. - // - // In our code below, these functions do not access member variables of the - // current class, so we can mark them as static and simply pass - // pointers to these functions to the MeshWorker framework. If, however, - // these functions would want to access member variables (or needed - // additional arguments beyond the ones specified below), we could use the - // facilities of lambda functions to provide the - // MeshWorker framework with objects that act as if they had the required - // number and types of arguments, but have in fact other arguments already - // bound. - static void integrate_cell_term(DoFInfo &dinfo, CellInfo &info); - static void integrate_boundary_term(DoFInfo &dinfo, CellInfo &info); - static void integrate_face_term(DoFInfo &dinfo1, - DoFInfo &dinfo2, - CellInfo &info1, - CellInfo &info2); - }; - - - // We start with the constructor. The 1 in the constructor call of - // fe is the polynomial degree. - template - AdvectionProblem::AdvectionProblem() - : mapping() - , fe(1) - , dof_handler(triangulation) - {} - - - template - void AdvectionProblem::setup_system() - { - // In the function that sets up the usual finite element data structures, we - // first need to distribute the DoFs. - dof_handler.distribute_dofs(fe); - - // We start by generating the sparsity pattern. To this end, we first fill - // an intermediate object of type DynamicSparsityPattern with the couplings - // appearing in the system. After building the pattern, this object is - // copied to sparsity_pattern and can be discarded. - - // To build the sparsity pattern for DG discretizations, we can call the - // function analogue to DoFTools::make_sparsity_pattern, which is called - // DoFTools::make_flux_sparsity_pattern: - DynamicSparsityPattern dsp(dof_handler.n_dofs()); - DoFTools::make_flux_sparsity_pattern(dof_handler, dsp); - sparsity_pattern.copy_from(dsp); - - // Finally, we set up the structure of all components of the linear system. - system_matrix.reinit(sparsity_pattern); - solution.reinit(dof_handler.n_dofs()); - right_hand_side.reinit(dof_handler.n_dofs()); - } - - // @sect4{The assemble_system function} - - // Here we see the major difference to assembling by hand. Instead of writing - // loops over cells and faces, we leave all this to the MeshWorker framework. - // In order to do so, we just have to define local integration functions and - // use one of the classes in namespace MeshWorker::Assembler to build the - // global system. - template - void AdvectionProblem::assemble_system() - { - // This is the magic object, which knows everything about the data - // structures and local integration. This is the object doing the work in - // the function MeshWorker::loop(), which is implicitly called by - // MeshWorker::integration_loop() below. After the functions to which we - // provide pointers did the local integration, the - // MeshWorker::Assembler::SystemSimple object distributes these into the - // global sparse matrix and the right hand side vector. - MeshWorker::IntegrationInfoBox info_box; - - // First, we initialize the quadrature formulae and the update flags in the - // worker base class. For quadrature, we play safe and use a QGauss formula - // with number of points one higher than the polynomial degree used. Since - // the quadratures for cells, boundary and interior faces can be selected - // independently, we have to hand over this value three times. - const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1; - info_box.initialize_gauss_quadrature(n_gauss_points, - n_gauss_points, - n_gauss_points); - - // These are the types of values we need for integrating our system. They - // are added to the flags used on cells, boundary and interior faces, as - // well as interior neighbor faces, which is forced by the four @p true - // values. - info_box.initialize_update_flags(); - UpdateFlags update_flags = - update_quadrature_points | update_values | update_gradients; - info_box.add_update_flags(update_flags, true, true, true, true); - - // After preparing all data in info_box, we initialize the FEValues - // objects in there. - info_box.initialize(fe, mapping); - - // The object created so far helps us do the local integration on each cell - // and face. Now, we need an object which receives the integrated (local) - // data and forwards them to the assembler. - MeshWorker::DoFInfo dof_info(dof_handler); - - // Now, we have to create the assembler object and tell it, where to put the - // local data. These will be our system matrix and the right hand side. - MeshWorker::Assembler::SystemSimple, Vector> - assembler; - assembler.initialize(system_matrix, right_hand_side); - - // Finally, the integration loop over all active cells (determined by the - // first argument, which is an active iterator). - // - // As noted in the discussion when declaring the local integration functions - // in the class declaration, the arguments expected by the assembling - // integrator class are not actually function pointers. Rather, they are - // objects that can be called like functions with a certain number of - // arguments. Consequently, we could also pass objects with appropriate - // operator() implementations here, or lambda functions if the local - // integrators were, for example, non-static member functions. - MeshWorker::loop, - MeshWorker::IntegrationInfoBox>( - dof_handler.begin_active(), - dof_handler.end(), - dof_info, - info_box, - &AdvectionProblem::integrate_cell_term, - &AdvectionProblem::integrate_boundary_term, - &AdvectionProblem::integrate_face_term, - assembler); - } - - - // @sect4{The local integrators} - - // These are the functions given to the MeshWorker::integration_loop() called - // just above. They compute the local contributions to the system matrix and - // right hand side on cells and faces. - template - void AdvectionProblem::integrate_cell_term(DoFInfo &dinfo, - CellInfo &info) - { - // First, let us retrieve some of the objects used here from @p info. Note - // that these objects can handle much more complex structures, thus the - // access here looks more complicated than might seem necessary. - const FEValuesBase &fe_values = info.fe_values(); - FullMatrix &local_matrix = dinfo.matrix(0).matrix; - const std::vector &JxW = fe_values.get_JxW_values(); - - // With these objects, we continue local integration like always. First, we - // loop over the quadrature points and compute the advection vector in the - // current point. - for (unsigned int point = 0; point < fe_values.n_quadrature_points; ++point) - { - const Tensor<1, dim> beta_at_q_point = - beta(fe_values.quadrature_point(point)); - - // We solve a homogeneous equation, thus no right hand side shows up in - // the cell term. What's left is integrating the matrix entries. - for (unsigned int i = 0; i < fe_values.dofs_per_cell; ++i) - for (unsigned int j = 0; j < fe_values.dofs_per_cell; ++j) - local_matrix(i, j) += -beta_at_q_point * // - fe_values.shape_grad(i, point) * // - fe_values.shape_value(j, point) * // - JxW[point]; - } - } - - // Now the same for the boundary terms. Note that now we use FEValuesBase, the - // base class for both FEFaceValues and FESubfaceValues, in order to get - // access to normal vectors. - template - void AdvectionProblem::integrate_boundary_term(DoFInfo &dinfo, - CellInfo &info) - { - const FEValuesBase &fe_face_values = info.fe_values(); - FullMatrix &local_matrix = dinfo.matrix(0).matrix; - Vector &local_vector = dinfo.vector(0).block(0); - - const std::vector &JxW = fe_face_values.get_JxW_values(); - const std::vector> &normals = - fe_face_values.get_normal_vectors(); - - std::vector g(fe_face_values.n_quadrature_points); - - static BoundaryValues boundary_function; - boundary_function.value_list(fe_face_values.get_quadrature_points(), g); - - for (unsigned int point = 0; point < fe_face_values.n_quadrature_points; - ++point) - { - const double beta_dot_n = - beta(fe_face_values.quadrature_point(point)) * normals[point]; - if (beta_dot_n > 0) - for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i) - for (unsigned int j = 0; j < fe_face_values.dofs_per_cell; ++j) - local_matrix(i, j) += beta_dot_n * // - fe_face_values.shape_value(j, point) * // - fe_face_values.shape_value(i, point) * // - JxW[point]; - else - for (unsigned int i = 0; i < fe_face_values.dofs_per_cell; ++i) - local_vector(i) += -beta_dot_n * // - g[point] * // - fe_face_values.shape_value(i, point) * // - JxW[point]; - } - } - - // Finally, the interior face terms. The difference here is that we receive - // two info objects, one for each cell adjacent to the face and we assemble - // four matrices, one for each cell and two for coupling back and forth. - template - void AdvectionProblem::integrate_face_term(DoFInfo &dinfo1, - DoFInfo &dinfo2, - CellInfo &info1, - CellInfo &info2) - { - // For quadrature points, weights, etc., we use the FEValuesBase object of - // the first argument. - const FEValuesBase &fe_face_values = info1.fe_values(); - const unsigned int dofs_per_cell = fe_face_values.dofs_per_cell; - - // For additional shape functions, we have to ask the neighbors - // FEValuesBase. - const FEValuesBase &fe_face_values_neighbor = info2.fe_values(); - const unsigned int neighbor_dofs_per_cell = - fe_face_values_neighbor.dofs_per_cell; - - // Then we get references to the four local matrices. The letters u and v - // refer to trial and test functions, respectively. The %numbers indicate - // the cells provided by info1 and info2. By convention, the two matrices - // in each info object refer to the test functions on the respective cell. - // The first matrix contains the interior couplings of that cell, while the - // second contains the couplings between cells. - FullMatrix &u1_v1_matrix = dinfo1.matrix(0, false).matrix; - FullMatrix &u2_v1_matrix = dinfo1.matrix(0, true).matrix; - FullMatrix &u1_v2_matrix = dinfo2.matrix(0, true).matrix; - FullMatrix &u2_v2_matrix = dinfo2.matrix(0, false).matrix; - - // Here, following the previous functions, we would have the local right - // hand side vectors. Fortunately, the interface terms only involve the - // solution and the right hand side does not receive any contributions. - - const std::vector &JxW = fe_face_values.get_JxW_values(); - const std::vector> &normals = - fe_face_values.get_normal_vectors(); - - for (unsigned int point = 0; point < fe_face_values.n_quadrature_points; - ++point) - { - const double beta_dot_n = - beta(fe_face_values.quadrature_point(point)) * normals[point]; - if (beta_dot_n > 0) - { - // This term we've already seen: - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int j = 0; j < dofs_per_cell; ++j) - u1_v1_matrix(i, j) += beta_dot_n * // - fe_face_values.shape_value(j, point) * // - fe_face_values.shape_value(i, point) * // - JxW[point]; - - // We additionally assemble the term $(\beta\cdot n u,\hat - // v)_{\partial \kappa_+}$, - for (unsigned int k = 0; k < neighbor_dofs_per_cell; ++k) - for (unsigned int j = 0; j < dofs_per_cell; ++j) - u1_v2_matrix(k, j) += - -beta_dot_n * // - fe_face_values.shape_value(j, point) * // - fe_face_values_neighbor.shape_value(k, point) * // - JxW[point]; - } - else - { - // This one we've already seen, too: - for (unsigned int i = 0; i < dofs_per_cell; ++i) - for (unsigned int l = 0; l < neighbor_dofs_per_cell; ++l) - u2_v1_matrix(i, l) += - beta_dot_n * // - fe_face_values_neighbor.shape_value(l, point) * // - fe_face_values.shape_value(i, point) * // - JxW[point]; - - // And this is another new one: $(\beta\cdot n \hat u,\hat - // v)_{\partial \kappa_-}$: - for (unsigned int k = 0; k < neighbor_dofs_per_cell; ++k) - for (unsigned int l = 0; l < neighbor_dofs_per_cell; ++l) - u2_v2_matrix(k, l) += - -beta_dot_n * // - fe_face_values_neighbor.shape_value(l, point) * // - fe_face_values_neighbor.shape_value(k, point) * // - JxW[point]; - } - } - } - - - // @sect3{All the rest} - // - // For this simple problem we use the simplest possible solver, called - // Richardson iteration, that represents a simple defect correction. This, in - // combination with a block SSOR preconditioner, that uses the special block - // matrix structure of system matrices arising from DG discretizations. The - // size of these blocks are the number of DoFs per cell. Here, we use a SSOR - // preconditioning as we have not renumbered the DoFs according to the flow - // field. If the DoFs are renumbered in the downstream direction of the flow, - // then a block Gauss-Seidel preconditioner (see the PreconditionBlockSOR - // class with relaxation=1) does a much better job. - template - void AdvectionProblem::solve(Vector &solution) - { - SolverControl solver_control(1000, 1e-12); - SolverRichardson> solver(solver_control); - - // Here we create the preconditioner, - PreconditionBlockSSOR> preconditioner; - - // then assign the matrix to it and set the right block size: - preconditioner.initialize(system_matrix, fe.n_dofs_per_cell()); - - // After these preparations we are ready to start the linear solver. - solver.solve(system_matrix, solution, right_hand_side, preconditioner); - } - - - // We refine the grid according to a very simple refinement criterion, namely - // an approximation to the gradient of the solution. As here we consider the - // DG(1) method (i.e. we use piecewise bilinear shape functions) we could - // simply compute the gradients on each cell. But we do not want to base our - // refinement indicator on the gradients on each cell only, but want to base - // them also on jumps of the discontinuous solution function over faces - // between neighboring cells. The simplest way of doing that is to compute - // approximative gradients by difference quotients including the cell under - // consideration and its neighbors. This is done by the - // DerivativeApproximation class that computes the approximate - // gradients in a way similar to the GradientEstimation described - // in step-9 of this tutorial. In fact, the - // DerivativeApproximation class was developed following the - // GradientEstimation class of step-9. Relating to the discussion - // in step-9, here we consider $h^{1+d/2}|\nabla_h u_h|$. Furthermore we note - // that we do not consider approximate second derivatives because solutions to - // the linear advection equation are in general not in $H^2$ but only in $H^1$ - // (or, to be more precise: in $H^1_\beta$, i.e., the space of functions whose - // derivatives in direction $\beta$ are square integrable). - template - void AdvectionProblem::refine_grid() - { - // The DerivativeApproximation class computes the gradients to - // float precision. This is sufficient as they are approximate and serve as - // refinement indicators only. - Vector gradient_indicator(triangulation.n_active_cells()); - - // Now the approximate gradients are computed - DerivativeApproximation::approximate_gradient(mapping, - dof_handler, - solution, - gradient_indicator); - - // and they are cell-wise scaled by the factor $h^{1+d/2}$ - unsigned int cell_no = 0; - for (const auto &cell : dof_handler.active_cell_iterators()) - gradient_indicator(cell_no++) *= - std::pow(cell->diameter(), 1 + 1.0 * dim / 2); - - // Finally they serve as refinement indicator. - GridRefinement::refine_and_coarsen_fixed_number(triangulation, - gradient_indicator, - 0.3, - 0.1); - - triangulation.execute_coarsening_and_refinement(); - } - - - // The output of this program consists of eps-files of the adaptively refined - // grids and the numerical solutions given in gnuplot format. - template - void AdvectionProblem::output_results(const unsigned int cycle) const - { - // First write the grid in eps format. - { - const std::string filename = "grid-" + std::to_string(cycle) + ".eps"; - deallog << "Writing grid to <" << filename << '>' << std::endl; - std::ofstream eps_output(filename); - - GridOut grid_out; - grid_out.write_eps(triangulation, eps_output); - } - - // Then output the solution in gnuplot format. - { - const std::string filename = "sol-" + std::to_string(cycle) + ".gnuplot"; - deallog << "Writing solution to <" << filename << '>' << std::endl; - std::ofstream gnuplot_output(filename); - - DataOut data_out; - data_out.attach_dof_handler(dof_handler); - data_out.add_data_vector(solution, "u"); - - data_out.build_patches(); - - data_out.write_gnuplot(gnuplot_output); - } - } - - - // The following run function is similar to previous examples. - template - void AdvectionProblem::run() - { - for (unsigned int cycle = 0; cycle < 6; ++cycle) - { - deallog << "Cycle " << cycle << std::endl; - - if (cycle == 0) - { - GridGenerator::hyper_cube(triangulation); - - triangulation.refine_global(3); - } - else - refine_grid(); - - - deallog << "Number of active cells: " - << triangulation.n_active_cells() << std::endl; - - setup_system(); - - deallog << "Number of degrees of freedom: " << dof_handler.n_dofs() - << std::endl; - - assemble_system(); - solve(solution); - - output_results(cycle); - } - } -} // namespace Step12 - - -// The following main function is similar to previous examples as -// well, and need not be commented on. -int main() -{ - try - { - using namespace dealii; - deallog.depth_console(5); - - Step12::AdvectionProblem<2> dgmethod; - dgmethod.run(); - } - 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; -} diff --git a/examples/step-21/doc/intro.dox b/examples/step-21/doc/intro.dox index 79a80650f3..283bcea995 100644 --- a/examples/step-21/doc/intro.dox +++ b/examples/step-21/doc/intro.dox @@ -287,7 +287,7 @@ The quantities $S_+,\mathbf{u}_+$ then correspond to the values of these variables on the present cell, whereas $S_-,\mathbf{u}_-$ (needed on the inflow part of the boundary of $K$) are quantities taken from the neighboring cell. Some more context on discontinuous element techniques and evaluation of -fluxes can also be found in step-12 and step-12b. +fluxes can also be found in step-12.

    Linear solvers

    diff --git a/examples/step-39/doc/builds-on b/examples/step-39/doc/builds-on index f1775c1cbd..20c3064395 100644 --- a/examples/step-39/doc/builds-on +++ b/examples/step-39/doc/builds-on @@ -1 +1 @@ -step-12b +step-12 diff --git a/examples/step-39/doc/intro.dox b/examples/step-39/doc/intro.dox index 0b4debc58a..319e3efa40 100644 --- a/examples/step-39/doc/intro.dox +++ b/examples/step-39/doc/intro.dox @@ -79,9 +79,9 @@ Thus, the functions for assembling matrices, right hand side and error estimates below exhibit that these loops are all generic and can be programmed in the same way. -This program is related to step-12b, in that it uses MeshWorker and -discontinuous Galerkin methods. While there, we solved an advection -problem, here it is a diffusion problem. Here, we also use multigrid +This program is related to step-12, in that it uses MeshWorker and +discontinuous Galerkin methods. There we solved an advection +problem, while here it is a diffusion problem. Here, we also use multigrid preconditioning and a theoretically justified error estimator, see Karakashian and Pascal (2003). The multilevel scheme was discussed in detail in Kanschat (2004). The adaptive iteration and its convergence