From: bangerth Date: Wed, 20 Jan 2010 13:35:12 +0000 (+0000) Subject: Get rid again of the DGIntegrator class and re-merge everything into X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=96c0bab875eba1014a47f3243a8dff4aa6e06d9e;p=dealii-svn.git Get rid again of the DGIntegrator class and re-merge everything into the main class. git-svn-id: https://svn.dealii.org/trunk@20400 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc index 91a3816862..86af46cb78 100644 --- a/deal.II/examples/step-38/step-38.cc +++ b/deal.II/examples/step-38/step-38.cc @@ -91,44 +91,307 @@ void BoundaryValues::value_list(const std::vector > &points, // @sect3{Integrating cell and face matrices} + // @sect3{Class: DGMethod} // - // We define a class that fits into - // the MeshWorker framework. Since it - // will be used by - // MeshWorker::AssemblingIntegrator, - // it needs functions for cell, - // boundary and interior face - // integration specified exactly as - // below. - - // The base class Subscriptor is - // needed so that - // MeshWorker::AssemblingIntegrator - // can store a SmartPointer to an - // object of this class. + // 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, + // because there are no hanging node + // constraints in DG discretizations. template -class DGIntegrator : public Subscriptor +class DGMethod { public: - // First, we define the types of - // the two info objects handed to - // the local integration - // functions in order to make our - // life easier below. + DGMethod (); + ~DGMethod (); + + 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. + FE_DGQ fe; + DoFHandler dof_handler; + + 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; + + // 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 typedefs to the + // two info objects handed to the + // local integration functions in + // order to make our life easier + // below. typedef typename MeshWorker::IntegrationWorker::CellInfo CellInfo; typedef typename MeshWorker::IntegrationWorker::FaceInfo FaceInfo; // The following three functions - // are the ones that get called + // are then the ones that get called // inside the generic loop over all // cells and faces. They are the // ones doing the actual // integration. - static void cell(CellInfo& info); - static void bdry(FaceInfo& info); - static void face(FaceInfo& info1, FaceInfo& info2); + // + // 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 + // boost::bind (or std::bind, + // respectively) 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 (CellInfo& info); + static void integrate_boundary_term (FaceInfo& info); + static void integrate_face_term (FaceInfo& info1, + FaceInfo& info2); }; + + // We start with the + // constructor. This is the + // place to change the + // polynomial degree of the + // finite element shape + // functions. +template +DGMethod::DGMethod () + : + fe (1), + dof_handler (triangulation) +{} + + +template +DGMethod::~DGMethod () +{ + 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); + + // 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(); + + system_matrix.reinit (sparsity_pattern); + + solution.reinit (dof_handler.n_dofs()); + right_hand_side.reinit (dof_handler.n_dofs()); +} + + // @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 namespace MeshWorker::Assembler + // to build the global system. +template +void DGMethod::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::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::AssemblingIntegrator + , + Vector > > + integrator(&DGMethod::integrate_cell_term, + &DGMethod::integrate_boundary_term, + &DGMethod::integrate_face_term); + + // 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; + integrator.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. + 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 + // put the local data. These will + // be our system matrix and the + // right hand side. + integrator.initialize(system_matrix, right_hand_side); + + // We are now ready to get to the + // integration loop. @p info_box is + // an object that generates the + // extended iterators for cells and + // faces of type + // MeshWorker::IntegrationInfo. Since + // we need five different of them, + // this is a handy shortcut. It + // receives all the stuff we + // created so far. + MeshWorker::IntegrationInfoBox info_box(dof_handler); + info_box.initialize(integrator, fe, mapping); + + // Finally, the integration loop + // over all active cells + // (determined by the first + // argument, which is an active iterator). + MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator); +} + + // @sect4{The local integrators} // These functions are analogous to @@ -147,7 +410,7 @@ class DGIntegrator : public Subscriptor // added soon). template -void DGIntegrator::cell(CellInfo& info) +void DGMethod::integrate_cell_term (CellInfo& info) { // First, let us retrieve some of // the objects used here from @@ -193,7 +456,7 @@ void DGIntegrator::cell(CellInfo& info) // FESubfaceValues, in order to get access to // normal vectors. template -void DGIntegrator::bdry(FaceInfo& info) +void DGMethod::integrate_boundary_term (FaceInfo& info) { const FEFaceValuesBase& fe_v = info.fe(); FullMatrix& local_matrix = info.M1[0].matrix; @@ -239,7 +502,8 @@ void DGIntegrator::bdry(FaceInfo& info) // for each cell and two for coupling // back and forth. template -void DGIntegrator::face(FaceInfo& info1, FaceInfo& info2) +void DGMethod::integrate_face_term (FaceInfo& info1, + FaceInfo& info2) { // For quadrature points, weights, // etc., we use the @@ -337,243 +601,6 @@ void DGIntegrator::face(FaceInfo& info1, FaceInfo& info2) } - // @sect3{Class: DGMethod} - // - // 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, - // because there are no hanging node - // constraints in DG discretizations. -template -class DGMethod -{ - public: - DGMethod (); - ~DGMethod (); - - 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. - FE_DGQ fe; - DoFHandler dof_handler; - - 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; -}; - - - // We start with the - // constructor. This is the - // place to change the - // polynomial degree of the - // finite element shape - // functions. -template -DGMethod::DGMethod () - : - fe (1), - dof_handler (triangulation) -{} - - -template -DGMethod::~DGMethod () -{ - 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); - - // 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(); - - system_matrix.reinit (sparsity_pattern); - - solution.reinit (dof_handler.n_dofs()); - right_hand_side.reinit (dof_handler.n_dofs()); -} - - // @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 namespace MeshWorker::Assembler - // to build the global system. -template -void DGMethod::assemble_system () -{ - // Here we generate an object of - // our own integration class, which - // knows how to compute cell and - // face contributions for the - // matrix and the residual. - const DGIntegrator dg; - - // This is the magic object, which - // knows everything about the data - // structures and local integration - // (the latter through our object - // @p dg). This is the object doing - // the work in the function - // MeshWorker::loop(), which is - // implicitly called by - // MeshWorker::integration_loop() - // below. - // After @p dg did the local - // integration, the - // MeshWorker::Assembler::SystemSimple - // 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 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 > > - integrator(&DGIntegrator::cell, - &DGIntegrator::bdry, - &DGIntegrator::face); - - // 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; - integrator.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. - 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 - // put the local data. These will - // be our system matrix and the - // right hand side. - integrator.initialize(system_matrix, right_hand_side); - - // We are now ready to get to the - // integration loop. @p info_box is - // an object that generates the - // extended iterators for cells and - // faces of type - // MeshWorker::IntegrationInfo. Since - // we need five different of them, - // this is a handy shortcut. It - // receives all the stuff we - // created so far. - MeshWorker::IntegrationInfoBox info_box(dof_handler); - info_box.initialize(integrator, fe, mapping); - - // Finally, the integration loop - // over all active cells - // (determined by the first - // argument, which is an active iterator). - MeshWorker::integration_loop(dof_handler.begin_active(), dof_handler.end(), info_box, integrator); -} - - // @sect3{All the rest} // // For this simple problem we use the