From: Guido Kanschat Date: Fri, 12 Feb 2010 17:27:40 +0000 (+0000) Subject: updated doc, but not tested X-Git-Tag: v8.0.0~6485 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=24ee2a44b4b878272a0a2f5ac89427c9cdb62c21;p=dealii.git updated doc, but not tested git-svn-id: https://svn.dealii.org/trunk@20573 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-38/doc/intro.dox b/deal.II/examples/step-38/doc/intro.dox index 7eca8fa7bc..1f2de88b3e 100644 --- a/deal.II/examples/step-38/doc/intro.dox +++ b/deal.II/examples/step-38/doc/intro.dox @@ -4,14 +4,11 @@

Overview

This example is devoted to the MeshWorker framework and the discontinuous -Galerkin method, or in short: DG method. It solves the same problem as -@ref step_12 "step-12" (see there for a description of the problem and -discretization), but here we use the MeshWorker framework in order to save -programming the cell/face loops that are often rather. The aim of the -MeshWorker framework is to simplify this process, by putting the majority of -the boring setup into a framework class and leaving to user code only things -that are specific to the application. We have tried to strip this example of -peripheral information such that the structure becomes more clear. +Galerkin method, or in short: DG method. It includes the following topics. +
    +
  1. Discretization of the linear transport equation with the DG method. +
  2. Assembling of the system matrix using the MeshWorker::loop(). +
The particular concern of this program are the loops of DG methods. These turn out to be especially complex, primarily because for the face terms, we have to @@ -25,3 +22,95 @@ to write the local integrators for your problem. Second, you select classes from the MeshWorker namespace and combine them to achieve your goal. +

Problem

+ +The model problem solved in this example is the linear advection equation +@f[ + \nabla\cdot \left({\mathbf \beta} u\right)=f \qquad\mbox{in }\Omega, +\qquad\qquad\qquad\mathrm{[transport-equation]}@f] +subject to the boundary conditions +@f[ +u=g\quad\mbox{on }\Gamma_-, +@f] +on the inflow part $\Gamma_-$ of the boundary $\Gamma=\partial\Omega$ +of the domain. Here, ${\mathbf \beta}={\mathbf \beta}({\bf x})$ denotes a +vector field, $f$ a source function, $u$ the (scalar) solution +function, $g$ a boundary value function, +@f[ +\Gamma_-:=\{{\bf x}\in\Gamma, {\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})<0\} +@f] +the inflow part of the boundary of the domain and ${\bf n}$ denotes +the unit outward normal to the boundary $\Gamma$. Equation +[transport-equation] is the conservative version of the +transport equation already considered in step 9 of this tutorial. + +In particular, we consider problem [transport-equation] on +$\Omega=[0,1]^2$ with ${\mathbf \beta}=\frac{1}{|x|}(-x_2, x_1)$ +representing a circular counterclockwise flow field, $f=0$ and $g=1$ +on ${\bf x}\in\Gamma_-^1:=[0,0.5]\times\{0\}$ and $g=0$ on ${\bf x}\in +\Gamma_-\setminus \Gamma_-^1$. + + +

Discretization

+ +For deriving the DG +discretization we start with a variational, mesh-dependent +formulation of the problem, +@f[ + \sum_\kappa\left\{-\beta u,\nabla v)_\kappa+(u^+ \beta\cdot{\bf n}, v)_{\partial\kappa}\right\}=(f,v)_\Omega, +@f] + +that originates from [transport-equation] by +multiplication with a test function $v$ and integration by parts on +each cell $\kappa$ of the triangulation. Here $(\cdot, \cdot)_\kappa$ +and $(\cdot, \cdot)_{\partial\kappa}$ denote the +L2-inner products on the cell $\kappa$ and the +boundary $\partial\kappa$ of the cell, respectively. $u^+$ is the +value of u taken from the upwind cell with respect to $\beta$ +of the face, that is, the cell $\beta$ points away from. To discretize +the problem, the functions $u$ and $v$ are replaced by discrete +functions $u_h$ and $v_h$ that in the case of discontinuous Galerkin +methods belong to the space $V_h$ of discontinuous piecewise +polynomial functions of some degree $p$. + +Hence, the discontinuous Galerkin +scheme for the [transport-equation] is given +by: find $u_h\in V_h$ such that for all $v_h\in V_h$ following +equation holds: +@f[ + \sum_\kappa\left\{-(u_h,{\mathbf \beta}\cdot\nabla v_h)_\kappa + +({\mathbf \beta}\cdot{\bf n}\, u_h, v_h)_{\partial\kappa_+} + +({\mathbf \beta}\cdot{\bf n}\, u_h^-, v_h)_{\partial\kappa_-\setminus\Gamma}\right\} + =(f,v_h)_\Omega-({\mathbf \beta}\cdot{\bf n}\, g, v_h)_{\Gamma_-}, +\qquad\qquad\qquad\mathrm{[dg-transport1]}@f] +where $\partial\kappa_-:=\{{\bf x}\in\partial\kappa, +{\mathbf \beta}({\bf x})\cdot{\bf n}({\bf x})<0\}$ denotes the inflow boundary +and $\partial\kappa_+=\partial\kappa\setminus \partial \kappa_-$ the +outflow part of cell $\kappa$. Below, this equation will be referred +to as first version of the DG method. We note that after a +second integration by parts, we obtain: find $u_h\in V_h$ such that +@f[ + \sum_\kappa\left\{(\nabla\cdot\{{\mathbf \beta} u_h\},v_h)_\kappa + -({\mathbf \beta}\cdot{\bf n} [u_h], v_h)_{\partial\kappa_-}\right\} + =(f,v_h)_\Omega, \quad\forall v_h\in V_h, +@f] +where $[u_h]=u_h^+-u_h^-$ denotes the jump of the discrete function +between two neighboring cells and is defined to be $[u_h]=u_h^+-g$ on +the boundary of the domain. This is the discontinuous Galerkin scheme +for the transport equation given in its original notation. +Nevertheless, we will base the implementation of the scheme on the +form given by [dg-general1] and [upwind-flux], +or [dg-transport1], respectively. + +Finally, we rewrite [dg-general1] in terms of a summation over all +faces where each face $e=\partial \kappa\cap\partial \kappa'$ +between two neighboring cells $\kappa$ and $\kappa'$ occurs twice, obtaining + +@f[ + -\sum_\kappa(u_h,{\mathbf \beta}\cdot\nabla v_h)_\kappa + +\sum_{E\in\mathbb E_h^i} (u_h^-, \beta\cdot[v_h\mathbf n])_{E} + =(f,v_h)_\Omega-(g, \beta\cdot\mathbf n v_h)_{\Gamma_-}. +@f] + +In this form, we need to implement a diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc index 9d7a1ad814..21935f5fc6 100644 --- a/deal.II/examples/step-38/step-38.cc +++ b/deal.II/examples/step-38/step-38.cc @@ -11,12 +11,13 @@ /* further information on this license. */ // The first few files have already - // been covered in step-12 + // been covered in previous examples // and will thus not be further // commented on: #include #include #include +#include #include #include #include @@ -30,11 +31,32 @@ #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: the 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 -#include // Here come the new include files // for using the MeshWorker framework: @@ -42,19 +64,24 @@ #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 - // The last step is as in all - // previous programs: using namespace dealii; // @sect3{Equation data} // - // First, we need to describe the - // coefficients in the equation. Here, this - // concerns the boundary values, which we - // choose in the same way as for step-12: + // 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 { @@ -89,27 +116,29 @@ void BoundaryValues::value_list(const std::vector > &points, values[i]=0.; } } - - - // @sect3{Integrating cell and face matrices} - // @sect3{Class: DGMethod} + // @sect3{Class: Step12} // - // After these preparations, we - // proceed with the main part of this - // program. The main class, here - // called DGMethod is basically - // the main class of step-6. One of - // the differences is that there's no - // ConstraintMatrix object. This is, + // After this preparations, we + // proceed with the main class of + // this program, + // called Step12. It is basically + // the main class of step-6. We do + // not have a ConstraintMatrix, // 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 DGMethod +class Step12 { public: - DGMethod (); - ~DGMethod (); - + Step12 (); void run (); private: @@ -134,22 +163,19 @@ class DGMethod 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; - - // In step-12 we had two solution vectors - // that stored the solutions to the - // problems corresponding to the two - // different assembling routines - // assemble_system1 and - // assemble_system2. In this - // program, the goal is only to show the - // MeshWorker framework, so we only - // assemble the system in one of the two - // ways, and consequently we have only - // one solution vector along with the - // single assemble_system - // function declared above: + Vector solution; Vector right_hand_side; @@ -166,10 +192,11 @@ class DGMethod // then work on intermediate // objects for which first, we // here define typedefs to the - // two info objects handed to the + // info objects handed to the // local integration functions in // order to make our life easier // below. + typedef MeshWorker::DoFInfo DoFInfo; typedef typename MeshWorker::IntegrationWorker::CellInfo CellInfo; typedef typename MeshWorker::IntegrationWorker::FaceInfo FaceInfo; @@ -201,22 +228,18 @@ class DGMethod // types of arguments, but have // in fact other arguments // already bound. - static void integrate_cell_term (MeshWorker::DoFInfo& dinfo, CellInfo& info); - static void integrate_boundary_term (MeshWorker::DoFInfo& dinfo, FaceInfo& info); - static void integrate_face_term (MeshWorker::DoFInfo& dinfo1, - MeshWorker::DoFInfo& dinfo2, FaceInfo& info1, - FaceInfo& info2); + static void integrate_cell_term (DoFInfo& dinfo, CellInfo& info); + static void integrate_boundary_term (DoFInfo& dinfo, FaceInfo& info); + static void integrate_face_term (DoFInfo& dinfo1, DoFInfo& dinfo2, + FaceInfo& info1, FaceInfo& info2); }; - // We start with the - // constructor. This is the - // place to change the - // polynomial degree of the - // finite element shape - // functions. + // We start with the constructor. The + // 1 in the constructor call of #fe + // is the polynomial degree. template -DGMethod::DGMethod () +Step12::Step12 () : fe (1), dof_handler (triangulation) @@ -224,49 +247,37 @@ DGMethod::DGMethod () template -DGMethod::~DGMethod () +void Step12::setup_system () { - dof_handler.clear (); -} - - // In the function that sets up the usual // finite element data structures, we first // need to distribute the DoFs. -template -void DGMethod::setup_system () -{ dof_handler.distribute_dofs (fe); - // The DoFs of a cell are coupled with all - // DoFs of all neighboring cells, along - // with all of its siblings on the current - // cell. Therefore the maximum number of - // matrix entries per row is needed when - // all neighbors of a cell are once more - // refined than the cell under - // consideration. - sparsity_pattern.reinit (dof_handler.n_dofs(), - dof_handler.n_dofs(), - (GeometryInfo::faces_per_cell * - GeometryInfo::max_children_per_face - + - 1)*fe.dofs_per_cell); - + // We start by generating the + // sparsity pattern. To this end, + // we first fill an intermediate + // object of type + // CompressedSparsityPattern 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: - DoFTools::make_flux_sparsity_pattern (dof_handler, sparsity_pattern); - - // All following function calls are - // already known. - sparsity_pattern.compress(); + CompressedSparsityPattern c_sparsity(dof_handler.n_dofs()); + DoFTools::make_flux_sparsity_pattern (dof_handler, c_sparsity); + sparsity_pattern.copy_from(c_sparsity); + // 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()); } @@ -278,11 +289,11 @@ void DGMethod::setup_system () // 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 objects and use one of the + // integration functions and use one of the // classes in namespace MeshWorker::Assembler // to build the global system. template -void DGMethod::assemble_system () +void Step12::assemble_system () { // This is the magic object, which // knows everything about the data @@ -299,36 +310,6 @@ void DGMethod::assemble_system () // object distributes these into // the global sparse matrix and the // right hand side vector. - // - // MeshWorker::AssemblingIntegrator - // is not all that clever by - // itself, but its capabilities are - // provided the arguments provided - // to the constructor and by its - // second template argument. By - // exchanging - // MeshWorker::Assembler::SystemSimple, - // we could for instance assemble a - // BlockMatrix or just a Vector - // instead. - // - // 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 the - // result of std::bind if the local - // integrators were, for example, - // non-static member functions. MeshWorker::IntegrationWorker integration_worker; // First, we initialize the @@ -387,12 +368,30 @@ void DGMethod::assemble_system () // (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 the + // result of std::bind if the local + // integrators were, for example, + // non-static member functions. MeshWorker::integration_loop (dof_handler.begin_active(), dof_handler.end(), info_box, - &DGMethod::integrate_cell_term, - &DGMethod::integrate_boundary_term, - &DGMethod::integrate_face_term, + &Step12::integrate_cell_term, + &Step12::integrate_boundary_term, + &Step12::integrate_face_term, assembler, true); } @@ -415,7 +414,7 @@ void DGMethod::assemble_system () // added soon). template -void DGMethod::integrate_cell_term (MeshWorker::DoFInfo& dinfo, CellInfo& info) +void Step12::integrate_cell_term (DoFInfo& dinfo, CellInfo& info) { // First, let us retrieve some of // the objects used here from @@ -461,7 +460,7 @@ void DGMethod::integrate_cell_term (MeshWorker::DoFInfo& dinfo, CellIn // FESubfaceValues, in order to get access to // normal vectors. template -void DGMethod::integrate_boundary_term (MeshWorker::DoFInfo& dinfo, FaceInfo& info) +void Step12::integrate_boundary_term (DoFInfo& dinfo, FaceInfo& info) { const FEFaceValuesBase& fe_v = info.fe_values(); FullMatrix& local_matrix = dinfo.matrix(0).matrix; @@ -507,10 +506,8 @@ void DGMethod::integrate_boundary_term (MeshWorker::DoFInfo& dinfo, Fa // for each cell and two for coupling // back and forth. template -void DGMethod::integrate_face_term (MeshWorker::DoFInfo& dinfo1, - MeshWorker::DoFInfo& dinfo2, - FaceInfo& info1, - FaceInfo& info2) +void Step12::integrate_face_term (DoFInfo& dinfo1, DoFInfo& dinfo2, + FaceInfo& info1, FaceInfo& info2) { // For quadrature points, weights, // etc., we use the @@ -628,7 +625,7 @@ void DGMethod::integrate_face_term (MeshWorker::DoFInfo& dinfo1, // PreconditionBlockSOR class with // relaxation=1) does a much better job. template -void DGMethod::solve (Vector &solution) +void Step12::solve (Vector &solution) { SolverControl solver_control (1000, 1e-12, false, false); SolverRichardson<> solver (solver_control); @@ -688,7 +685,7 @@ void DGMethod::solve (Vector &solution) // (to be more precise, in $H^1_\beta$) // only. template -void DGMethod::refine_grid () +void Step12::refine_grid () { // The DerivativeApproximation // class computes the gradients to @@ -731,7 +728,7 @@ void DGMethod::refine_grid () // in previous examples and will not // be further commented on. template -void DGMethod::output_results (const unsigned int cycle) const +void Step12::output_results (const unsigned int cycle) const { // Write the grid in eps format. std::string filename = "grid-"; @@ -769,7 +766,7 @@ void DGMethod::output_results (const unsigned int cycle) const // The following run function is // similar to previous examples. template -void DGMethod::run () +void Step12::run () { for (unsigned int cycle=0; cycle<6; ++cycle) { @@ -809,7 +806,7 @@ int main () { try { - DGMethod<2> dgmethod; + Step12<2> dgmethod; dgmethod.run (); } catch (std::exception &exc)