From 8a8de190e3f0a2b1514257dda5be083cd5453ba4 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 8 Jan 2010 09:12:39 +0000 Subject: [PATCH] Small edits. git-svn-id: https://svn.dealii.org/trunk@20331 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-38/doc/intro.dox | 28 +-- deal.II/examples/step-38/step-38.cc | 229 +++++++++++++------------ 2 files changed, 133 insertions(+), 124 deletions(-) diff --git a/deal.II/examples/step-38/doc/intro.dox b/deal.II/examples/step-38/doc/intro.dox index fdc96eed63..7eca8fa7bc 100644 --- a/deal.II/examples/step-38/doc/intro.dox +++ b/deal.II/examples/step-38/doc/intro.dox @@ -1,21 +1,23 @@ -

Example of the MeshWorker framework with an advection problem

+

An example of the MeshWorker framework with an advection problem

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 reprogramming the cell/face -loops. We have tried to strip this example of peripheral information -such that the structure becomes more clear. +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. -In particular the loops of DG methods turn out to be complex, because -for the face terms, we have to distinguish the cases of boundary, -regular interior faces and interior faces with hanging nodes, -respectively. The MeshWorker framework implements the standard loop -over all cells and faces in MeshWorker::loop() and takes care of +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 +distinguish the cases of boundary, regular interior faces and interior faces +with hanging nodes, respectively. The MeshWorker framework implements the +standard loop over all cells and faces in MeshWorker::loop() and takes care of distinguishing between all the different faces. There are two things left to do if you use MeshWorker: first, you need diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc index 090e4edb86..383a724020 100644 --- a/deal.II/examples/step-38/step-38.cc +++ b/deal.II/examples/step-38/step-38.cc @@ -3,7 +3,7 @@ /* $Id$ */ /* */ -/* Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010 by the deal.II authors */ +/* Copyright (C) 2010 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ @@ -13,7 +13,7 @@ // The first few files have already // been covered in step-12 // and will thus not be further - // commented on. + // commented on: #include #include #include @@ -37,7 +37,7 @@ #include // Here come the new include files - // for using the MeshWorker framework + // for using the MeshWorker framework: #include #include @@ -50,6 +50,10 @@ 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: template class BoundaryValues: public Function { @@ -60,15 +64,14 @@ class BoundaryValues: public Function const unsigned int component=0) const; }; - // 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. + // 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, @@ -93,7 +96,7 @@ void BoundaryValues::value_list(const std::vector > &points, // the MeshWorker framework. Since it // will be used by // MeshWorker::AssemblingIntegrator, - // it needs the functions for cell, + // it needs functions for cell, // boundary and interior face // integration specified exactly as // below. @@ -174,25 +177,24 @@ void DGIntegrator::cell(CellInfo& info) const beta(1) = fe_v.quadrature_point(point)(0); beta /= beta.norm(); + // 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 void DGIntegrator::bdry(FaceInfo& info) const @@ -257,7 +259,7 @@ void DGIntegrator::face(FaceInfo& info1, FaceInfo& info2) const // four local matrices. The letters // u and v refer to trial and test // functions, respectively. The - // numbers indicate the cells + // %numbers indicate the cells // provided by info1 and info2. By // convention, the two matrices in // each info object refer to the @@ -277,7 +279,7 @@ void DGIntegrator::face(FaceInfo& info1, FaceInfo& info2) const // vectors. Fortunately, the // interface terms only involve the // solution and the right hand side - // does not obtain a contribution. + // does not receive any contributions. const std::vector &JxW = fe_v.get_JxW_values (); const std::vector > &normals = fe_v.get_normal_vectors (); @@ -293,7 +295,7 @@ void DGIntegrator::face(FaceInfo& info1, FaceInfo& info2) const if (beta_n>0) { // This term we've already - // seen. + // seen: for (unsigned int i=0; i::face(FaceInfo& info1, FaceInfo& info2) const else { // This one we've already - // seen, too. + // seen, too: for (unsigned int i=0; i::face(FaceInfo& info1, FaceInfo& info2) const // And this is another new // one: $(\beta\cdot n \hat // u,\hat v)_{\partial - // \kappa_-}$. + // \kappa_-}$: for (unsigned int k=0; k fe; DoFHandler dof_handler; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; - // And there are two solution - // vectors, that store the - // solutions to the problems - // corresponding to the two + // 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; + // 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; }; + // We start with the + // constructor. This is the + // place to change the + // polynomial degree of the + // finite element shape + // functions. template DGMethod::DGMethod () : - // Change here for DG - // methods of - // different degrees. fe (1), dof_handler (triangulation) {} @@ -413,29 +424,35 @@ DGMethod::~DGMethod () } + // In the function that sets up the usual + // finite element data structures, we first + // need to distribute the DoFs. template void DGMethod::setup_system () { - // First we need to distribute the - // DoFs. dof_handler.distribute_dofs (fe); - // The DoFs of a cell are coupled - // with all DoFs of all neighboring - // cells. 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 + // 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); + (GeometryInfo::faces_per_cell * + GeometryInfo::max_children_per_face + + + 1)*fe.dofs_per_cell); - // For DG discretizations we call - // the function analogue to - // DoFTools::make_sparsity_pattern. + // 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 @@ -449,15 +466,15 @@ void DGMethod::setup_system () } // @sect4{Function: assemble_system} - // 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 objects and use - // one of the classes in Assembler to - // build the global system. + + // 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 objects and use one of the + // classes in namespace MeshWorker::Assembler + // to build the global system. template void DGMethod::assemble_system () { @@ -485,18 +502,21 @@ void DGMethod::assemble_system () // the global sparse matrix and the // right hand side vector. // - // It turns out, // MeshWorker::AssemblingIntegrator - // itself is not that clever at - // all, but all its capabilities - // are provided by the two later + // is not all that clever by itself, + // but its capabilities + // are provided by its two latter // template arguments. By // exchanging // MeshWorker::Assembler::SystemSimple, // we could for instance assemble a // BlockMatrix or just a Vector // instead. - MeshWorker::AssemblingIntegrator, Vector >, DGIntegrator > + MeshWorker::AssemblingIntegrator + , + Vector >, + DGIntegrator > integrator(dg); // First, we initialize the @@ -512,7 +532,9 @@ void DGMethod::assemble_system () // independently, we have to hand // over this value three times. const unsigned int n_gauss_points = dof_handler.get_fe().degree+1; - integrator.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points); + integrator.initialize_gauss_quadrature(n_gauss_points, + n_gauss_points, + n_gauss_points); // These are the types of values we // need for integrating our @@ -521,19 +543,21 @@ void DGMethod::assemble_system () // and interior faces, as well as // interior neighbor faces, which is // forced by the four @p true values. - UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; + UpdateFlags update_flags = update_quadrature_points | + update_values | + update_gradients; integrator.add_update_flags(update_flags, true, true, true, true); // Finally, we have to tell the - // assembler base class, where to + // assembler base class where to // put the local data. These will // be our system matrix and the // right hand side. integrator.initialize(system_matrix, right_hand_side); - // Finally, we get to the + // We are now ready to get to the // integration loop. @p info_box is - // an object, that generates the + // an object that generates the // extended iterators for cells and // faces of type // MeshWorker::IntegrationInfo. Since @@ -556,25 +580,21 @@ void DGMethod::assemble_system () // // 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 downstream the flow, - // then a block Gauss-Seidel + // 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) makes a much better - // job. + // relaxation=1) does a much better job. template void DGMethod::solve (Vector &solution) { @@ -585,8 +605,8 @@ void DGMethod::solve (Vector &solution) // preconditioner, PreconditionBlockSSOR > preconditioner; - // we assigned the matrix to it and - // set the right block size. + // then assign the matrix to it and + // set the right block size: preconditioner.initialize(system_matrix, fe.dofs_per_cell); // After these preparations we are @@ -715,18 +735,7 @@ void DGMethod::output_results (const unsigned int cycle) const // The following run function is - // similar to previous examples. The - // only difference is that the - // problem is assembled and solved - // twice on each refinement step; - // first by assemble_system1 that - // implements the first version and - // then by assemble_system2 that - // implements the second version of - // writing the DG - // discretization. Furthermore the - // time needed by each of the two - // assembling routines is measured. + // similar to previous examples. template void DGMethod::run () { @@ -757,15 +766,13 @@ void DGMethod::run () assemble_system (); solve (solution); - // Finally we perform the - // output. output_results (cycle); } } // The following main function is - // similar to previous examples and - // need not to be commented on. + // similar to previous examples as well, and + // need not be commented on. int main () { try -- 2.39.5