From ebb63b7332467b15f5cb45d8b0e88e4f3c6ff49f Mon Sep 17 00:00:00 2001 From: hartmann Date: Tue, 20 Nov 2001 16:23:23 +0000 Subject: [PATCH] Complete program, and major part of documentation. git-svn-id: https://svn.dealii.org/trunk@5220 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-12/step-12.cc | 1568 +++++++++++++++++++++------ 1 file changed, 1220 insertions(+), 348 deletions(-) diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc index 741eb6c9ab..5bea1272c8 100644 --- a/deal.II/examples/step-12/step-12.cc +++ b/deal.II/examples/step-12/step-12.cc @@ -21,28 +21,60 @@ #include #include #include -#include + // This is the first new file. It + // declares the MappingQ1 class that + // gives the standard bilinear + // mapping. For bilinear mappings use + // an object of this class rather + // than an object of the MappingQ(1) + // class, as the MappingQ1 class is + // optimized due to the + // pre-knowledge of the actual + // polynomial degree 1. #include + + // Here the discontinuous finite + // elements are defined. They are + // used as all other finite elements. #include -#include -#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 structur + // of system matrices arising from DG + // discretizations. +#include +#include -#include + // We are going to use gradients as + // refinement indicator. +#include -template -class Beta -{ - public: - Beta () {}; + // Finally we do some time comparison + // using the ``Timer'' class. +#include - void value_list (const std::vector > &points, - std::vector > &values) const; -}; + // And this again is C++: +#include + // First we define the class + // representing the equation-specific + // functions. Both classes, ``RHS'' + // and ``BoundaryValues'', are + // derived from the Function + // class. Only the ``value_list'' + // function are implemented because + // only lists of function values are + // computed rather than single + // values. template class RHS: public Function { @@ -56,90 +88,66 @@ class RHS: public Function template -class BoundaryFunction: public Function +class BoundaryValues: public Function { public: - BoundaryFunction() {}; + BoundaryValues() {}; virtual void value_list (const std::vector > &points, std::vector &values, const unsigned int component=0) const; }; + + // The class ``Beta'' that represents + // the vector valued flow field of + // the linear transport equation is + // not derived from the Function + // class as we prefer to get function + // values of type ``Point'' rather + // than of type + // ``Vector''. This, because + // there exist scalar products + // between ``Point'' and ``Point'' as + // well as between ``Point'' and + // ``Tensor'', simplifying terms like + // $\beta\cdot n$ and + // $\beta\cdot\nabla v$. template -class DGAssembler +class Beta { public: - DGAssembler() {}; + Beta () {}; - void assemble_cell_term(const FEValuesBase& fe_v, - FullMatrix &cell_matrix, - Vector &cell_vector); - - void assemble_face_term(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &cell_matrix, - FullMatrix &cell_inflow_matrix, - Vector &cell_vector); - - private: - Beta beta_function; - RHS rhs_function; - BoundaryFunction boundary_function; + void value_list (const std::vector > &points, + std::vector > &values) const; }; - // The main class is again almost - // unchanged. Two additions, however, - // are made: we have added the - // ``refine'' function, which is used - // to adaptively refine the grid - // (instead of the global refinement - // in the previous examples), and a - // variable which will hold the - // constraints associated to the - // hanging nodes. + + // The implementation of the + // ``value_list'' functions of these + // classes are rather simple. For + // simplicity the right hand side is + // set to be zero. template -class TransportProblem +void RHS::value_list(const std::vector > &, + std::vector &values, + const unsigned int) const { - public: - TransportProblem (); - ~TransportProblem (); - - void run (); - - private: - void setup_system (); - void assemble_system (); - void solve (); - void refine_grid (); - void output_results (const unsigned int cycle) const; - - Triangulation triangulation; - MappingQ1 mapping; - - // We need a finite element - // again. This time, we will want - // to use quadratic polynomials - // (but this is only specified in - // the constructor): - FE_DGQ fe; - DoFHandler dof_handler; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - Vector solution; - Vector right_hand_side; - - DGAssembler dg_assembler; -}; - + for (unsigned int i=0; i void Beta<2>::value_list(const std::vector > &points, std::vector > &values) const { - Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size())); + Assert(values.size()==points.size(), + ExcDimensionMismatch(values.size(),points.size())); + for (unsigned int i=0; i &p=points[i]; @@ -151,27 +159,23 @@ void Beta<2>::value_list(const std::vector > &points, } } - - - -template -void RHS::value_list(const std::vector > &, - std::vector &values, - const unsigned int) const -{ - for (unsigned int i=0; i -void BoundaryFunction::value_list(const std::vector > &points, +void BoundaryValues::value_list(const std::vector > &points, std::vector &values, const unsigned int) const { - Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size())); + Assert(values.size()==points.size(), + ExcDimensionMismatch(values.size(),points.size())); + for (unsigned int i=0; i::value_list(const std::vector > &points, } } + // Next we define the equation- + // dependent and DG-method-dependent + // class ``DGTransportEquation''. Its + // member functions were already + // mentioned in the Introduction and + // will be explained + // below. Furthermore it includes + // objects of the previously defined + // ``Beta'', ``RHS'' and + // ``BoundaryValues'' function + // classes. +template +class DGTransportEquation +{ + public: + DGTransportEquation() {}; + void assemble_cell_term(const FEValues& fe_v, + FullMatrix &u_v_matrix, + Vector &cell_vector); + + void assemble_face_term1(const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &u_v_matrix, + FullMatrix &un_v_matrix, + Vector &cell_vector); + + void assemble_face_term2(const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &u_v_matrix, + FullMatrix &un_v_matrix, + FullMatrix &u_vn_matrix, + FullMatrix &un_vn_matrix, + Vector &cell_vector); + private: + Beta beta_function; + RHS rhs_function; + BoundaryValues boundary_function; +}; - + // ``u_v_matrix'' is a cell matrix, + // i.e. for a DG method of degree 1, + // it is of size 4 times 4, and + // ``cell_vector'' is of size 4. + // When this function is invoked, + // ``fe_v'' was reinited with the + // current cell before and includes + // all shape values needed. template -void DGAssembler::assemble_cell_term(const FEValuesBase& fe_v, - FullMatrix &cell_matrix, - Vector &cell_vector) +void DGTransportEquation::assemble_cell_term( + const FEValues& fe_v, + FullMatrix &u_v_matrix, + Vector &cell_vector) { + // First we ask ``fe_v'' for the + // shape grads, shape values and + // quadrature weights, const vector > > &grad_v = fe_v.get_shape_grads (); const FullMatrix &v = fe_v.get_shape_values (); const vector &JxW = fe_v.get_JxW_values (); + // Then the flow field beta and the + // ``rhs_function'' are evaluated at + // the quadrature points, vector > beta (fe_v.n_quadrature_points); vector rhs (fe_v.n_quadrature_points); beta_function.value_list (fe_v.get_quadrature_points(), beta); rhs_function.value_list (fe_v.get_quadrature_points(), rhs); + // and the cell matrix and cell + // vector are assembled as in + // previous tutorial steps. Here, + // the terms $-(u,\beta\cdot\nabla + // v)_K$ and $(f,v)_K$ are + // assembled. for (unsigned int point=0; point::assemble_cell_term(const FEValuesBase& fe_v, } + // The ``assemble_face_term1'' + // function assembles the face terms + // corresponding to the first version + // of the DG method, cf. above. Then, + // the face terms are given as a sum + // of integrals over all cell + // boundaries. + // + // When this function is invoked, + // ``fe_v'' and ``fe_v_neighbor'' are + // already reinited with the current + // cell and the neighoring cell, + // respectively, as well as with the + // current face. Hence they provide + // the inner and outer shape values + // on the face. + // + // In addition to the cell matrix + // ``u_v_matrix'' and the + // ``cell_vector'' this function has + // got a new argument + // ``un_v_matrix'', that stores + // contributions to the system matrix + // that are based on outer values of + // u, see $\hat u_h$ in the + // Introduction, and inner values of + // v, see $v_h$. Here we note that + // ``un'' is the short notation for + // ``u_neighbor'' and represents + // $\hat u_h$. template -void DGAssembler::assemble_face_term(const FEFaceValuesBase& fe_v, - const FEFaceValuesBase& fe_v_neighbor, - FullMatrix &cell_matrix, - FullMatrix &cell_inflow_matrix, - Vector &cell_vector) +void DGTransportEquation::assemble_face_term1( + const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &u_v_matrix, + FullMatrix &un_v_matrix, + Vector &cell_vector) { - DoFHandler::face_iterator face=fe_v.get_face(); - + // Again, we ask the FEValues + // objects for the shape values and + // the quadrature weights const FullMatrix &v = fe_v.get_shape_values (); const FullMatrix &v_neighbor = fe_v_neighbor.get_shape_values (); const vector &JxW = fe_v.get_JxW_values (); + // but also for the normals. const vector > &normals = fe_v.get_normal_vectors (); + // We also evaluate the flow field + // at the quadrature points vector > beta (fe_v.n_quadrature_points); - vector g(fe_v.n_quadrature_points); beta_function.value_list (fe_v.get_quadrature_points(), beta); + // and the boundary values if the + // current face belongs to the + // boundary. + vector g(fe_v.n_quadrature_points); + DoFHandler::face_iterator face=fe_v.get_face(); if (face->at_boundary()) boundary_function.value_list (fe_v.get_quadrature_points(), g); + // Then we assemble the cell matrix + // and cell vector according to the + // DG method given in the + // introduction. for (unsigned int point=0; point0) + // The term $(\beta\cdot n + // u,v)_{\partial K_+}$, for (unsigned int i=0; iat_boundary()) for (unsigned int i=0; i::assemble_face_term(const FEFaceValuesBase& fe_v, v(i,point) * JxW[point]; else + // and on inner faces the + // term $(\beta\cdot n + // \hat u,v)_{\partial + // K_-}$ for (unsigned int i=0; i +void DGTransportEquation::assemble_face_term2( + const FEFaceValuesBase& fe_v, + const FEFaceValuesBase& fe_v_neighbor, + FullMatrix &u_v_matrix, + FullMatrix &un_v_matrix, + FullMatrix &u_vn_matrix, + FullMatrix &un_vn_matrix, + Vector &cell_vector) +{ + // the first few lines are the same + const FullMatrix &v = fe_v.get_shape_values (); + const FullMatrix &v_neighbor = fe_v_neighbor.get_shape_values (); + const vector &JxW = fe_v.get_JxW_values (); + const vector > &normals = fe_v.get_normal_vectors (); + + vector > beta (fe_v.n_quadrature_points); + beta_function.value_list (fe_v.get_quadrature_points(), beta); + + vector g(fe_v.n_quadrature_points); + DoFHandler::face_iterator face=fe_v.get_face(); + if (face->at_boundary()) + boundary_function.value_list (fe_v.get_quadrature_points(), g); + + for (unsigned int point=0; point0) + { + // This terms we've already seen, + for (unsigned int i=0; iat_boundary()) + for (unsigned int k=0; kat_boundary()) + for (unsigned int i=0; i +class DGMethod +{ + public: + DGMethod (); + ~DGMethod (); + + void run (); + + private: + void setup_system (); + void assemble_system1 (); + void assemble_system2 (); + void solve (Vector &solution); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + MappingQ1 mapping; + + // Furthermore we want to + // use DG elements of degree 1 + // (but this is only specified in + // the constructor): + FE_DGQ fe; + DoFHandler dof_handler; + + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + + // We define the quadrature + // formulae for the cell and the + // face terms of the + // discretization. + QGauss4 quadrature; + QGauss4 face_quadrature; + + // And there are two solution + // vectors, that store the + // solutions to the problems + // corresponding to the two + // different assembling routines + // ``assemble_system1'' and + // ``assemble_system2''; + Vector solution1; + Vector solution2; + Vector right_hand_side; + + // Finally this class includes an + // object of the + // DGTransportEquations class + // described above. + DGTransportEquation dg; +}; + + + + // Now for the implementation of the + // main class. Constructor and + // destructor follow the same + // pattern that was used previously, + // so we need not comment on these + // two functions: template -TransportProblem::TransportProblem () : +DGMethod::DGMethod () : fe (1), dof_handler (triangulation) {} template -TransportProblem::~TransportProblem () +DGMethod::~DGMethod () { dof_handler.clear (); }; @@ -280,172 +565,460 @@ TransportProblem::~TransportProblem () template -void TransportProblem::setup_system () +void DGMethod::setup_system () { - // To distribute degrees of - // freedom, the ``dof_handler'' - // variable takes only the finite - // element object. In this case, it - // will distribute four degrees of - // freedom per cell. + // 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 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(), - dof_handler.max_couplings_between_dofs()); + (GeometryInfo::faces_per_cell + *GeometryInfo::subfaces_per_face+1)*fe.dofs_per_cell); + + // For DG discretizations we call + // the function analogue to + // DoFTools::make_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()); + solution1.reinit (dof_handler.n_dofs()); + solution2.reinit (dof_handler.n_dofs()); right_hand_side.reinit (dof_handler.n_dofs()); }; - + // We proceed with the + // ``assemble_system1'' function that + // implements the DG discretization + // in its first version. This + // function repeatedly calls the + // ``assemble_cell_term'' and + // ``assemble_face_term1'' functions + // of the DGTransportEquation object. + // The ``assemble_face_term1'' + // function takes two + // FEFaceValuesBase objects; one for + // the shape functions on the current + // cell and the other for shape + // functions on the neighboring cell + // under consideration. Both objects + // are either of class FEFaceValues + // or of class FESubfaceValues (both + // derived from FEFaceValuesBase) + // according to following cases + // already mentioned in the + // introduction: + // + // 1. face is at boundary (current + // cell: FEFaceValues, neighboring + // cell does not exist); + // + // 2. neighboring cell is finer + // (current cell: FESubfaceValues, + // neighboring cell: FEFaceValues); + // + // 3. neighboring cell is of the same + // refinement level (both, current + // and neighboring cell: + // FEFaceValues); + // + // 4. neighboring cell is coarser + // (current cell: FEFaceValues, + // neighboring cell: + // FESubfaceValues). + // + // If we considered globally refined + // meshes then only cases 1 and 3 + // would occur. But as we consider + // also locally refined meshes we + // need to distinguish all four cases + // making the following assembling + // function a bit longish. template -void TransportProblem::assemble_system () +void DGMethod::assemble_system1 () { - // See Cockburn paper for the proper quadrature. - QGauss4 quadrature; - QGauss4 face_quadrature; - const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; vector dofs (dofs_per_cell); vector dofs_neighbor (dofs_per_cell); + // First we create the Update flags + // for the FEValues and the + // FEFaceValues objects. UpdateFlags update_flags = UpdateFlags(update_values | update_gradients | update_q_points | update_JxW_values); + // Note, that on faces we do not + // need gradients but we need + // normal vectors. UpdateFlags face_update_flags = UpdateFlags(update_values | update_q_points | update_JxW_values | update_normal_vectors); + + // On the neighboring cell we only + // need the shape values. Given a + // specific face, the quadrature + // points and `JxW values' are the + // same as for the current cells, + // the normal vectors are known to + // be the negative of the normal + // vectors of the current cell. + UpdateFlags neighbor_face_update_flags = UpdateFlags(update_values); - + // Then we create the FEValues + // object. Note, that since version + // 3.2.0 the constructor of this + // class takes a Mapping object as + // first argument. Although the + // constructor without Mapping + // argument is still supported it + // is recommended to use the new + // constructor. This reduces the + // effect of `hidden magic' (the + // old constructor implicitely + // assumes a MappingQ1 mapping) and + // makes it easier to change the + // Mapping object later. FEValues fe_v ( mapping, fe, quadrature, update_flags); + + // Similarly we create the + // FEFaceValues and FESubfaceValues + // objects for both, the current + // and the neighboring cell. Within + // the following nested loop over + // all cells and all faces of the + // cell they will be reinited to + // the current cell and the face + // (and subface) number. FEFaceValues fe_v_face ( mapping, fe, face_quadrature, face_update_flags); FESubfaceValues fe_v_subface ( mapping, fe, face_quadrature, face_update_flags); FEFaceValues fe_v_face_neighbor ( - mapping, fe, face_quadrature, UpdateFlags(update_values | update_default)); + mapping, fe, face_quadrature, neighbor_face_update_flags); FESubfaceValues fe_v_subface_neighbor ( - mapping, fe, face_quadrature, UpdateFlags(update_values | update_default)); - - // includes the u and v terms - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - // includes u_hat and v terms - FullMatrix cell_inflow_matrix (dofs_per_cell, dofs_per_cell); + mapping, fe, face_quadrature, neighbor_face_update_flags); + + // Now we create the cell matrices + // and vectors. Here we need two + // cell matrices, both for face + // terms that include test + // functions ``v'' (shape functions + // of the current cell). To be more + // precise, the first matrix will + // include the `u and v terms' and + // the second that will include the + // `un and v terms'. Here we recall + // our the convention that `un' is + // the short cut for `u_neighbor' + // and represents the $u_hat$, see + // introduction. + FullMatrix u_v_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix un_v_matrix (dofs_per_cell, dofs_per_cell); Vector cell_vector (dofs_per_cell); + // Furthermore we need some cell + // and face iterators DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); DoFHandler::face_iterator face; DoFHandler::cell_iterator neighbor; - DoFHandler::cell_iterator neighbor_child; + DoFHandler::active_cell_iterator neighbor_child; + // Now we start the loop over all + // active cells for (;cell!=endc; ++cell) { - // re-init fe values for this cell + // and reinit the FEValues + // object for the current cell, fe_v.reinit (cell); - cell_matrix.clear (); - cell_vector.clear (); - - dg_assembler.assemble_cell_term(fe_v, - cell_matrix, - cell_vector); - + // Call the function that + // assembles the cell + // terms. The first argument is + // the FEValues that was + // already reinited on current + // the cell. + dg.assemble_cell_term(fe_v, + u_v_matrix, + cell_vector); + + // As in previous example steps + // the vector `dofs' includes + // the dof_indices. cell->get_dof_indices (dofs); + // This is the start of the + // nested loop over all faces. for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { + // First we set the face + // iterator. face = cell->face(face_no); - - cell_inflow_matrix.clear(); + // Now we distinguish the + // four different cases in + // the ordering mentioned + // above. We start with + // faces belonging to the + // boundary of the domain. if (face->at_boundary()) { + // We reinit the + // FEFaceValues object + // to the current face fe_v_face.reinit (cell, face_no); - dg_assembler.assemble_face_term(fe_v_face, - fe_v_face, - cell_matrix, - cell_inflow_matrix, - cell_vector); + // and assemble the + // corresponding face + // terms. Here, the + // second and fourth + // arguments are only + // dummy arguments. On + // the boundary of the + // domain the + // ``assemble_face_term1'' + // function will not + // access to shape + // values on the + // non-existent + // neighboring + // cell. Also, + // ``un_v_matrix'' will + // be unchanged. + dg.assemble_face_term1(fe_v_face, + fe_v_face, + u_v_matrix, + un_v_matrix, + cell_vector); } - else // if (!face->at_boundary()) + else { - Assert (cell->neighbor(face_no).state() == valid, ExcInternalError()); + // When we are not on the + // boundary of the + // domain then there + // must exist a + // neighboring cell. neighbor = cell->neighbor(face_no); - - if (face->has_children()) // i.e. neighbor is one level more refined than cell + + // We proceed with the + // second and most + // complicated case: + // the neighboring cell + // is more refined than + // the current cell. As + // in deal.II + // neighboring cells + // are restricted to + // have a level + // difference of not + // more than one, the + // neighboring cell is + // known to be only + // ONCE more refined + // than the current + // cell. Furthermore + // also the face is + // once more refined, + // i.e. it has + // children. + if (face->has_children()) { - // store which number #cell# has in the - // list of neighbors of #neighbor# - const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); + // first we store + // which number the + // current cell has + // in the list of + // neighbors of the + // neighboring + // cell. Hence, + // neighbor->neighbor(neighbor2) + // equals the + // current cell + // ``cell''. + const unsigned int neighbor2= + cell->neighbor_of_neighbor(face_no); - // loop over all subfaces - for (unsigned int subface_no=0; subface_no::subfaces_per_face; + // We loop over + // subfaces + for (unsigned int subface_no=0; + subface_no::subfaces_per_face; ++subface_no) { - // get an iterator pointing to the - // cell behind the present subface + // and set the + // cell + // iterator + // ``neighbor_child'' + // to the cell + // placed + // `behind' the + // current + // subface. neighbor_child = neighbor->child(GeometryInfo:: child_cell_on_face(neighbor2,subface_no)); + + // As these are + // quite + // complicated + // indirections + // we check for + // the internal + // consistency. Assert (neighbor_child->face(neighbor2) == face->child(subface_no), ExcInternalError()); Assert (!neighbor_child->has_children(), ExcInternalError()); + // As already + // mentioned + // above for + // this case + // (case 2) we + // employ the + // FESubfaceValues + // of the + // current + // cell, here + // reinited for + // the current + // cell, face + // and subface, + // and we + // employ the + // FEFaceValues + // of the + // neighboring + // child cell. fe_v_subface.reinit (cell, face_no, subface_no); fe_v_face_neighbor.reinit (neighbor_child, neighbor2); - - cell_inflow_matrix.clear(); - - dg_assembler.assemble_face_term(fe_v_subface, - fe_v_face_neighbor, - cell_matrix, - cell_inflow_matrix, - cell_vector); - - // get indices of dofs of neighbor_child cell + + dg.assemble_face_term1(fe_v_subface, + fe_v_face_neighbor, + u_v_matrix, + un_v_matrix, + cell_vector); + + // get dof + // indices of + // the + // neighbor_child + // cell neighbor_child->get_dof_indices (dofs_neighbor); - // distribute cell matrix + // distribute + // cell matrix + // to the + // system_matrix for (unsigned int i=0; ihas_children()) + // End of ``if + // (face->has_children())'' + else { + // We proceed with + // case 3, + // i.e. neighboring + // cell is of the + // same refinement + // level as the + // current cell. if (neighbor->level() == cell->level()) { - // store which number #cell# has in the - // list of neighbors of #neighbor# + // Like before we + // store which + // number the + // current cell has + // in the list of + // neighbors of the + // neighboring + // cell. const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); + // We reinit + // the + // FEFaceValues + // of the + // current and + // neighboring + // cell to the + // current face + // and assemble + // the + // corresponding + // face terms. fe_v_face.reinit (cell, face_no); fe_v_face_neighbor.reinit (neighbor, neighbor2); - dg_assembler.assemble_face_term(fe_v_face, - fe_v_face_neighbor, - cell_matrix, - cell_inflow_matrix, - cell_vector); + dg.assemble_face_term1(fe_v_face, + fe_v_face_neighbor, + u_v_matrix, + un_v_matrix, + cell_vector); + // End of ``if + // (neighbor->level() + // == + // cell->level())'' } - else // if (neighbor->level() < cell->level()) i.e. neighbor is one level coarser than cell + else { + // Finally we + // consider + // case 4. When + // the + // neighboring + // cell is not + // finer and + // not of the + // same + // refinement + // level as the + // current cell + // it must be + // coarser. Assert(neighbor->level() < cell->level(), ExcInternalError()); + // Find out the + // how many'th + // face_no and + // subface_no + // the current + // face is + // w.r.t. the + // neighboring + // cell. const std::pair faceno_subfaceno= cell->neighbor_of_coarser_neighbor(face_no); const unsigned int neighbor_face_no=faceno_subfaceno.first, @@ -454,127 +1027,460 @@ void TransportProblem::assemble_system () Assert (neighbor->neighbor(neighbor_face_no) ->child(GeometryInfo::child_cell_on_face( face_no,neighbor_subface_no)) == cell, ExcInternalError()); - - // now 'neighbor_face_no' stores the number - // of a face in the list of faces of 'neighbor'. - // This face has got a subface that is - // between 'cell' and 'neighbor'. - // 'neighbor_subface_no' stores the number - // of this subface in the list of subfaces of this - // face 'neighbor->face(neighbor_face_no)' - // that is between 'cell' and 'neighbor' + + // Reinit the + // appropriate + // FEFaceValues + // and assemble + // the face + // terms. fe_v_face.reinit (cell, face_no); fe_v_subface_neighbor.reinit (neighbor, neighbor_face_no, neighbor_subface_no); - dg_assembler.assemble_face_term(fe_v_face, - fe_v_subface_neighbor, - cell_matrix, - cell_inflow_matrix, - cell_vector); - } // else // if (neighbor->level() < cell->level()) - - // get indices of dofs of neighbor_child cell + dg.assemble_face_term1(fe_v_face, + fe_v_subface_neighbor, + u_v_matrix, + un_v_matrix, + cell_vector); + } + + // Get dof indices + // of the + // neighbor_child + // cell, neighbor->get_dof_indices (dofs_neighbor); - // distribute cell_inflow_matrix + // distribute the + // un_v_matrix, for (unsigned int i=0; ihas_children()) - } // else // if (!face->at_boundary()) - } //for (face_no...) + un_v_matrix(i,k)); + + // and clear the + // ``un_v_matrix'' + // on each face. + un_v_matrix.clear(); + } + // End of ``face not at boundary'': + } + // End of loop over all faces: + } - // distribute cell matrix + // Finally we distribute the + // u_v_matrix, for (unsigned int i=0; i -void TransportProblem::solve () -{ - SolverControl solver_control (1000, 1e-12); +void DGMethod::assemble_system2 () +{ + const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell; + vector dofs (dofs_per_cell); + vector dofs_neighbor (dofs_per_cell); + + UpdateFlags update_flags = UpdateFlags(update_values + | update_gradients + | update_q_points + | update_JxW_values); + + UpdateFlags face_update_flags = UpdateFlags(update_values + | update_q_points + | update_JxW_values + | update_normal_vectors); + + UpdateFlags neighbor_face_update_flags = UpdateFlags(update_values); + + // Here we do not need + // ``fe_v_face_neighbor'' as case 4 + // does not occur. + FEValues fe_v ( + mapping, fe, quadrature, update_flags); + FEFaceValues fe_v_face ( + mapping, fe, face_quadrature, face_update_flags); + FESubfaceValues fe_v_subface ( + mapping, fe, face_quadrature, face_update_flags); + FEFaceValues fe_v_face_neighbor ( + mapping, fe, face_quadrature, neighbor_face_update_flags); + + + FullMatrix u_v_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix un_v_matrix (dofs_per_cell, dofs_per_cell); + + // Additionally we need following + // two cell matrices, both for face + // term that include test function + // ``vn'' (shape functions of the + // neighboring cell). To be more + // precise, the first matrix will + // include the `u and vn terms' and + // the second that will include the + // `un and vn terms'. + FullMatrix u_vn_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix un_vn_matrix (dofs_per_cell, dofs_per_cell); + + Vector cell_vector (dofs_per_cell); + + // Furthermore, here we define a + // dummy matrix and rhs to + // emphasize when arguments of the + // ``assemble_face_term2'' + // functions will not be access. + FullMatrix dummy_matrix; + Vector dummy_rhs; + + // The following lines are roughly + // the same as in the previous + // function. + DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + DoFHandler::face_iterator face; + DoFHandler::cell_iterator neighbor; + DoFHandler::cell_iterator neighbor_child; + + for (;cell!=endc; ++cell) + { + fe_v.reinit (cell); + + dg.assemble_cell_term(fe_v, + u_v_matrix, + cell_vector); + + cell->get_dof_indices (dofs); + + for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) + { + face = cell->face(face_no); + + // Case 1: + if (face->at_boundary()) + { + fe_v_face.reinit (cell, face_no); + + dg.assemble_face_term2(fe_v_face, + fe_v_face, + u_v_matrix, + dummy_matrix, + dummy_matrix, + dummy_matrix, + cell_vector); + } + else + { + Assert (cell->neighbor(face_no).state() == valid, ExcInternalError()); + neighbor = cell->neighbor(face_no); + + // Case 2: + if (face->has_children()) + { + const unsigned int neighbor2= + cell->neighbor_of_neighbor(face_no); + + for (unsigned int subface_no=0; + subface_no::subfaces_per_face; + ++subface_no) + { + neighbor_child = neighbor->child( + GeometryInfo::child_cell_on_face(neighbor2,subface_no)); + Assert (neighbor_child->face(neighbor2) == face->child(subface_no), + ExcInternalError()); + Assert (!neighbor_child->has_children(), ExcInternalError()); + + fe_v_subface.reinit (cell, face_no, subface_no); + fe_v_face_neighbor.reinit (neighbor_child, neighbor2); + + dg.assemble_face_term2(fe_v_subface, + fe_v_face_neighbor, + u_v_matrix, + un_v_matrix, + u_vn_matrix, + un_vn_matrix, + dummy_rhs); + + neighbor_child->get_dof_indices (dofs_neighbor); + + for (unsigned int i=0; ilevel() == cell->level() && + neighbor->index() > cell->index()) + { + const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); + + fe_v_face.reinit (cell, face_no); + fe_v_face_neighbor.reinit (neighbor, neighbor2); + + dg.assemble_face_term2(fe_v_face, + fe_v_face_neighbor, + u_v_matrix, + un_v_matrix, + u_vn_matrix, + un_vn_matrix, + dummy_rhs); + + neighbor->get_dof_indices (dofs_neighbor); + + for (unsigned int i=0; i +void DGMethod::solve (Vector &solution) +{ + SolverControl solver_control (1000, 1e-12, false, false); PrimitiveVectorMemory<> vector_memory; SolverRichardson<> solver (solver_control, vector_memory); + // Here we create the + // preconditioner, PreconditionBlockSSOR preconditioner; + + // we asigned the matrix to it and + // set the right block size. preconditioner.initialize(system_matrix, fe.dofs_per_cell); + + // As the inverses of the diagonal + // blocks are needed in each + // preconditioner step, it is wise + // to invert the diagonal blocks of + // the matrix before starting the + // solver. Otherwise, the diagonal + // blocks are inverted in each + // preconditioner step, + // significantly slowing down the + // linear solving process. preconditioner.invert_diagblocks(); - + + // 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 the gradients 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 + // simpliest 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. According + // to the argumentation in Step 9, + // here we consider + // $h^{1+d/2}|\nabla_h + // u_h|$. Futhermore 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 in H^1 (to be more precise, in + // H^1_\beta) only. template -void TransportProblem::refine_grid () +void DGMethod::refine_grid () { - Vector estimated_error_per_cell (triangulation.n_active_cells()); - - FunctionMap::type neumann_boundary; - - KellyErrorEstimator::estimate (dof_handler, - QGauss3(), - neumann_boundary, - solution, - estimated_error_per_cell); + // 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, + solution2, + gradient_indicator); + + // and they are cell-wise scaled by + // the factor $h^{1+d/2}$ + DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int cell_no=0; cell!=endc; ++cell, ++cell_no) + 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, - estimated_error_per_cell, - 0.3, 0.03); + 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. This was covered + // in previous examples and will not + // be further commented on. template -void TransportProblem::output_results (const unsigned int cycle) const +void DGMethod::output_results (const unsigned int cycle) const { - // We want to write the grid in - // each cycle. Here is another way - // to quickly produce a filename - // based on the cycle number. It - // assumes that the numbers `0' - // through `9' are represented - // consecutively in the character - // set (which is the case in all - // known character sets). However, - // this will only work if the cycle - // number is less than ten, which - // we check by an assertion. + // Write the grid in eps format. std::string filename = "grid-"; filename += ('0' + cycle); Assert (cycle < 10, ExcInternalError()); filename += ".eps"; + cout << "Writing grid to <" << filename << ">..." << endl; std::ofstream eps_output (filename.c_str()); - // Using this filename, we write - // each grid as a postscript file. GridOut grid_out; grid_out.write_eps (triangulation, eps_output); - - // output of the solution + + // Output of the solution in + // gnuplot format. filename = "sol-"; filename += ('0' + cycle); Assert (cycle < 10, ExcInternalError()); filename += ".gnuplot"; + cout << "Writing solution to <" << filename << ">..." << endl; std::ofstream gnuplot_output (filename.c_str()); DataOut data_out; data_out.attach_dof_handler (dof_handler); - data_out.add_data_vector (solution, "u"); + data_out.add_data_vector (solution2, "u"); data_out.build_patches (); @@ -582,11 +1488,23 @@ void TransportProblem::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. template -void TransportProblem::run () +void DGMethod::run () { - for (unsigned int cycle=0; cycle<3; ++cycle) + for (unsigned int cycle=0; cycle<6; ++cycle) { std::cout << "Cycle " << cycle << ':' << std::endl; @@ -597,18 +1515,7 @@ void TransportProblem::run () triangulation.refine_global (3); } else - // In case this is not the - // first cycle, we want to - // refine the grid. Unlike - // the global refinement - // employed in the last - // example, we now use the - // adaptive procedure - // described in the function - // which we now call: - { - refine_grid (); - }; + refine_grid (); std::cout << " Number of active cells: " @@ -620,97 +1527,62 @@ void TransportProblem::run () std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; - - assemble_system (); - solve (); + + // The constructor of the Timer + // class automatically starts + // the time measurement. + Timer assemble_timer; + // First assembling routine. + assemble_system1 (); + // The operator () accesses the + // current time without + // disturbing the time + // measurement. + cout << "Time of assemble_system1: " << assemble_timer() << endl; + solve (solution1); + + // As preparation for the + // second assembling routine we + // reinit the system matrix, the + // right hand side vector and + // the Timer object. + system_matrix.reinit(); + right_hand_side.clear(); + assemble_timer.reset(); + + // We start the Timer, + assemble_timer.start(); + // call the second assembling routine + assemble_system2 (); + // and access the current time. + cout << "Time of assemble_system2: " << assemble_timer() << endl; + solve (solution2); + + // To make sure that both + // versions of the DG method + // yield the same + // discretization and hence the + // same solution we check the + // two solutions for equality. + solution1-=solution2; + const double difference=solution1.linfty_norm(); + if (difference<1e-13) + cout << "solution1 and solution2 do not differ." << endl; + + // Finally we perform the + // output. output_results (cycle); } } + + int main () { + DGMethod<2> dgmethod_2d; + dgmethod_2d.run (); - // The general idea behind the - // layout of this function is as - // follows: let's try to run the - // program as we did before... - try - { - TransportProblem<2> Transport_problem_2d; - Transport_problem_2d.run (); - } - // ...and if this should fail, try - // to gather as much information as - // possible. Specifically, if the - // exception that was thrown is an - // object of a class that is - // derived from the C++ standard - // class ``exception'', then we can - // use the ``what'' member function - // to get a string which describes - // the reason why the exception was - // thrown. - // - // The deal.II exception classes - // are all derived from the - // standard class, and in - // particular, the ``exc.what()'' - // function will return - // approximately the same string as - // would be generated if the - // exception was thrown using the - // ``Assert'' macro. You have seen - // the output of such an exception - // in the previous example, and you - // then know that it contains the - // file and line number of where - // the exception occured, and some - // other information. This is also - // what would be printed in the - // following. - 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; - // We can't do much more than - // printing as much information - // as we can get to, so abort - // with error: - return 1; - } - // If the exception that was thrown - // somewhere was not an object of a - // class derived from the standard - // ``exception'' class, then we - // can't do anything at all. We - // then simply print an error - // message and exit. - catch (...) - { - std::cerr << std::endl << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - }; - - // If we got to this point, there - // was no exception which - // propagated up to the main - // function (maybe there were some, - // but they were caught somewhere - // in the program or the - // library). Therefore, the program - // performed as was expected and we - // can return without error. return 0; }; + + -- 2.39.5