From 6b96ba02fcc551bd967b62cc110be5705bc576ca Mon Sep 17 00:00:00 2001 From: hartmann Date: Thu, 22 Nov 2001 15:37:47 +0000 Subject: [PATCH] A couple of changes and corrections. git-svn-id: https://svn.dealii.org/trunk@5241 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-12/step-12.cc | 699 ++++++++++++++++------------ 1 file changed, 389 insertions(+), 310 deletions(-) diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc index f3dddad75d..8d568d396a 100644 --- a/deal.II/examples/step-12/step-12.cc +++ b/deal.II/examples/step-12/step-12.cc @@ -1,5 +1,5 @@ /* $Id$ */ -/* Author: Ralf Hartmann, University of Heidelberg, 2000 */ +/* Author: Ralf Hartmann, University of Heidelberg, 2001 */ // The first few files have already // been covered in previous examples @@ -38,7 +38,7 @@ // used as all other finite elements. #include // We are going to use the simplest - // possible solver, called richardson + // possible solver, called Richardson // iteration, that represents a simple // defect correction. This, in // combination with a block SSOR @@ -60,11 +60,13 @@ #include - // First we define the class + // @sect3{Equation data} + // + // First we define the classes // representing the equation-specific // functions. Both classes, ``RHS'' // and ``BoundaryValues'', are - // derived from the Function + // derived from the ``Function'' // class. Only the ``value_list'' // function are implemented because // only lists of function values are @@ -74,8 +76,6 @@ template class RHS: public Function { public: - RHS() {}; - virtual void value_list (const std::vector > &points, std::vector &values, const unsigned int component=0) const; @@ -86,18 +86,16 @@ template class BoundaryValues: public Function { public: - 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 + // The class ``Beta'' represents the + // vector valued flow field of the + // linear transport equation and is + // not derived from the ``Function'' // class as we prefer to get function // values of type ``Point'' rather // than of type @@ -112,8 +110,6 @@ template class Beta { public: - Beta () {}; - void value_list (const std::vector > &points, std::vector > &values) const; }; @@ -123,35 +119,41 @@ class Beta // ``value_list'' functions of these // classes are rather simple. For // simplicity the right hand side is - // set to be zero. + // set to be zero but will be + // assembled anyway. template -void RHS::value_list(const std::vector > &, +void RHS::value_list(const std::vector > &points, std::vector &values, const unsigned int) const { + // Usually we check whether input + // parameter have the right sizes. + Assert(values.size()==points.size(), + ExcDimensionMismatch(values.size(),points.size())); + for (unsigned int i=0; i -void Beta<2>::value_list(const std::vector > &points, - std::vector > &values) const +template +void Beta::value_list(const std::vector > &points, + std::vector > &values) const { Assert(values.size()==points.size(), ExcDimensionMismatch(values.size(),points.size())); for (unsigned int i=0; i &p=points[i]; - Point<2> &beta=values[i]; + const Point &p=points[i]; + Point &beta=values[i]; beta(0)=-p(1); beta(1)=p(0); - beta/=sqrt(beta*beta); + beta/=sqrt(beta.square()); } } @@ -185,9 +187,10 @@ void BoundaryValues::value_list(const std::vector > &points, // @sect3{Class: DGTransportEquation} // - // Next we define the equation- - // dependent and DG-method-dependent - // class ``DGTransportEquation''. Its + // 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 @@ -204,27 +207,27 @@ class DGTransportEquation void assemble_cell_term(const FEValues& fe_v, FullMatrix &u_v_matrix, - Vector &cell_vector); + Vector &cell_vector) const; void assemble_boundary_term(const FEFaceValues& fe_v, FullMatrix &u_v_matrix, - Vector &cell_vector); + Vector &cell_vector) const; void assemble_face_term1(const FEFaceValuesBase& fe_v, const FEFaceValuesBase& fe_v_neighbor, FullMatrix &u_v_matrix, - FullMatrix &un_v_matrix); + FullMatrix &un_v_matrix) const; 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); + FullMatrix &un_vn_matrix) const; private: - Beta beta_function; - RHS rhs_function; - BoundaryValues boundary_function; + const Beta beta_function; + const RHS rhs_function; + const BoundaryValues boundary_function; }; @@ -238,27 +241,27 @@ class DGTransportEquation // 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 + // ``fe_v'' is already reinit'ed with the // current cell before and includes // all shape values needed. template void DGTransportEquation::assemble_cell_term( - const FEValues& fe_v, + const FEValues &fe_v, FullMatrix &u_v_matrix, - Vector &cell_vector) + Vector &cell_vector) const { // First we ask ``fe_v'' for the // shape gradients, shape values and // quadrature weights, - const vector > > &grad_v = fe_v.get_shape_grads (); + const std::vector > > &grad_v = fe_v.get_shape_grads (); const FullMatrix &v = fe_v.get_shape_values (); - const vector &JxW = fe_v.get_JxW_values (); + const std::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); + std::vector > beta (fe_v.n_quadrature_points); + std::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); @@ -286,7 +289,7 @@ void DGTransportEquation::assemble_cell_term( // function assembles the face terms // at boundary faces. When this // function is invoked, ``fe_v'' is - // already reinited with the current + // already reinit'ed with the current // cell and current face. Hence it // provides the shape values on that // boundary face. @@ -294,7 +297,7 @@ template void DGTransportEquation::assemble_boundary_term( const FEFaceValues& fe_v, FullMatrix &u_v_matrix, - Vector &cell_vector) + Vector &cell_vector) const { // First we check whether the // current face is really at the @@ -302,19 +305,19 @@ void DGTransportEquation::assemble_boundary_term( Assert(fe_v.get_face()->at_boundary(), ExcInternalError()); // Again, as in the previous - // function, we ask the FEValues + // function, we ask the ``FEValues'' // object for the shape values and // the quadrature weights const FullMatrix &v = fe_v.get_shape_values (); - const vector &JxW = fe_v.get_JxW_values (); + const std::vector &JxW = fe_v.get_JxW_values (); // but here also for the normals. - const vector > &normals = fe_v.get_normal_vectors (); + const std::vector > &normals = fe_v.get_normal_vectors (); // We evaluate the flow field // and the boundary values at the // quadrature points. - vector > beta (fe_v.n_quadrature_points); - vector g(fe_v.n_quadrature_points); + std::vector > beta (fe_v.n_quadrature_points); + std::vector g(fe_v.n_quadrature_points); beta_function.value_list (fe_v.get_quadrature_points(), beta); boundary_function.value_list (fe_v.get_quadrature_points(), g); @@ -325,11 +328,11 @@ void DGTransportEquation::assemble_boundary_term( // introduction. for (unsigned int point=0; point0) + const double beta_n=beta[point] * normals[point]; // We assemble the term // $(\beta\cdot n // u,v)_{\partial K_+}$, + if (beta_n>0) for (unsigned int i=0; i::assemble_boundary_term( // // When this function is invoked, // ``fe_v'' and ``fe_v_neighbor'' are - // already reinited with the current + // already reinit'ed with the current // cell and the neighoring cell, // respectively, as well as with the // current face. Hence they provide @@ -375,7 +378,7 @@ void DGTransportEquation::assemble_boundary_term( // contributions to the system matrix // that are based on outer values of // u, see $\hat u_h$ in the - // Introduction, and inner values of + // introduction, and inner values of // v, see $v_h$. Here we note that // ``un'' is the short notation for // ``u_neighbor'' and represents @@ -385,7 +388,7 @@ void DGTransportEquation::assemble_face_term1( const FEFaceValuesBase& fe_v, const FEFaceValuesBase& fe_v_neighbor, FullMatrix &u_v_matrix, - FullMatrix &un_v_matrix) + FullMatrix &un_v_matrix) const { // First we check that the current // face is not at the boundary by @@ -399,12 +402,12 @@ void DGTransportEquation::assemble_face_term1( // normals 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 (); + const std::vector &JxW = fe_v.get_JxW_values (); + const std::vector > &normals = fe_v.get_normal_vectors (); // and we evaluate the flow field // at the quadrature points. - vector > beta (fe_v.n_quadrature_points); + std::vector > beta (fe_v.n_quadrature_points); beta_function.value_list (fe_v.get_quadrature_points(), beta); // Then we assemble the cell @@ -413,11 +416,11 @@ void DGTransportEquation::assemble_face_term1( // introduction. for (unsigned int point=0; point0) + const double beta_n=beta[point] * normals[point]; // We assemble the term // $(\beta\cdot n // u,v)_{\partial K_+}$, + if (beta_n>0) for (unsigned int i=0; i::assemble_face_term2( FullMatrix &u_v_matrix, FullMatrix &un_v_matrix, FullMatrix &u_vn_matrix, - FullMatrix &un_vn_matrix) + FullMatrix &un_vn_matrix) const { // 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 (); + const std::vector &JxW = fe_v.get_JxW_values (); + const std::vector > &normals = fe_v.get_normal_vectors (); - vector > beta (fe_v.n_quadrature_points); + std::vector > beta (fe_v.n_quadrature_points); beta_function.value_list (fe_v.get_quadrature_points(), beta); for (unsigned int point=0; point0) { // This terms we've already seen. @@ -533,8 +536,8 @@ void DGTransportEquation::assemble_face_term2( // the main class of step 6. One of // the differences is that there's no // ConstraintMatrix object. This is, - // because there are no hanging nodes - // in DG discretizations. + // because there are no hanging node + // constraints in DG discretizations. template class DGMethod { @@ -553,12 +556,17 @@ class DGMethod void output_results (const unsigned int cycle) const; Triangulation triangulation; - MappingQ1 mapping; + const MappingQ1 mapping; - // Furthermore we want to - // use DG elements of degree 1 - // (but this is only specified in - // the constructor): + // 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 wanted + // degree. FE_DGQ fe; DoFHandler dof_handler; @@ -569,8 +577,8 @@ class DGMethod // formulae for the cell and the // face terms of the // discretization. - QGauss4 quadrature; - QGauss4 face_quadrature; + const QGauss4 quadrature; + const QGauss4 face_quadrature; // And there are two solution // vectors, that store the @@ -587,18 +595,15 @@ class DGMethod // object of the // DGTransportEquations class // described above. - DGTransportEquation dg; + const 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 DGMethod::DGMethod () : + // Change here for DG + // methods of + // different degrees. fe (1), dof_handler (triangulation) {} @@ -621,10 +626,11 @@ void DGMethod::setup_system () // 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. + // 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 @@ -657,67 +663,75 @@ void DGMethod::setup_system () // ``assemble_cell_term'', // ``assemble_boundary_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: + // of the ``DGTransportEquation'' + // object. The + // ``assemble_boundary_term'' covers + // the first case mentioned in the + // introduction. + // + // 1. face is at boundary // - // 1. face is at boundary (current - // cell: FEFaceValues); + // This function takes a + // ``FEFaceValues'' object as + // argument. In contrast to that + // ``assemble_face_term1'' + // 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 + // the remaining cases mentioned + // in the introduction: // // 2. neighboring cell is finer - // (current cell: FESubfaceValues, - // neighboring cell: FEFaceValues); + // (current cell: ``FESubfaceValues'', + // neighboring cell: ``FEFaceValues''); // // 3. neighboring cell is of the same // refinement level (both, current // and neighboring cell: - // FEFaceValues); + // ``FEFaceValues''); // // 4. neighboring cell is coarser - // (current cell: FEFaceValues, + // (current cell: ``FEFaceValues'', // neighboring cell: - // FESubfaceValues). + // ``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. + // meshes then only case 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 DGMethod::assemble_system1 () { 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); - + std::vector dofs (dofs_per_cell); + std::vector dofs_neighbor (dofs_per_cell); + + // First we create the + // ``UpdateFlags'' for the + // ``FEValues'' and the + // ``FEFaceValues'' objects. + UpdateFlags update_flags = 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); + UpdateFlags face_update_flags = update_values + | update_q_points + | update_JxW_values + | update_normal_vectors; // On the neighboring cell we only // need the shape values. Given a @@ -727,34 +741,36 @@ void DGMethod::assemble_system1 () // 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); + UpdateFlags neighbor_face_update_flags = update_values; - // Then we create the FEValues + // 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 + // 3.2.0 of deal.II 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. + // 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'' 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 ( @@ -774,8 +790,8 @@ void DGMethod::assemble_system1 () // 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' + // the convention that `un' is + // the shortcut for `u_neighbor' // and represents the $u_hat$, see // introduction. FullMatrix u_v_matrix (dofs_per_cell, dofs_per_cell); @@ -784,34 +800,45 @@ void DGMethod::assemble_system1 () 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::active_cell_iterator neighbor_child; + // iterators. + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); // Now we start the loop over all - // active cells + // active cells. for (;cell!=endc; ++cell) { - // and reinit the FEValues - // object for the current cell, + // In the + // ``assemble_face_term1'' + // function contributions to + // the cell matrices and the + // cell vector are only + // ADDED. Therefore on each + // cell we need to reset the + // ``u_v_matrix'' and + // ``cell_vector'' to zero, + // before assembling the cell terms. + u_v_matrix.clear (); + cell_vector.clear (); + + // Now we reinit the ``FEValues'' + // object for the current cell fe_v.reinit (cell); - // Call the function that - // assembles the cell + // and call the function + // that assembles the cell // terms. The first argument is - // the FEValues that was - // already reinited on current - // the cell. + // the ``FEValues'' that was + // previously reinit'ed on the + // current cell. dg.assemble_cell_term(fe_v, u_v_matrix, cell_vector); - // As in previous example steps - // the vector `dofs' includes - // the dof_indices. + // As in previous examples the + // vector `dofs' includes the + // dof_indices. cell->get_dof_indices (dofs); // This is the start of the @@ -819,9 +846,14 @@ void DGMethod::assemble_system1 () for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { // First we set the face - // iterator. - face = cell->face(face_no); - + // iterator + typename DoFHandler::face_iterator face=cell->face(face_no); + + // and clear the + // ``un_v_matrix'' on each + // face. + un_v_matrix.clear(); + // Now we distinguish the // four different cases in // the ordering mentioned @@ -831,8 +863,9 @@ void DGMethod::assemble_system1 () if (face->at_boundary()) { // We reinit the - // FEFaceValues object - // to the current face + // ``FEFaceValues'' + // object to the + // current face fe_v_face.reinit (cell, face_no); // and assemble the @@ -849,7 +882,8 @@ void DGMethod::assemble_system1 () // domain, therefore // there must exist a // neighboring cell. - neighbor = cell->neighbor(face_no); + typename DoFHandler::cell_iterator neighbor= + cell->neighbor(face_no);; // We proceed with the // second and most @@ -864,17 +898,21 @@ void DGMethod::assemble_system1 () // difference of not // more than one, the // neighboring cell is - // known to be only + // known to be at most // ONCE more refined // than the current // cell. Furthermore // also the face is - // once more refined, + // more refined, // i.e. it has - // children. + // children. Here we + // note that the + // following part of + // code will not work + // for ``dim==1''. if (face->has_children()) { - // first we store + // First we store // which number the // current cell has // in the list of @@ -904,34 +942,56 @@ void DGMethod::assemble_system1 () // `behind' the // current // subface. - neighbor_child = neighbor->child(GeometryInfo:: - child_cell_on_face(neighbor2,subface_no)); - + typename DoFHandler::active_cell_iterator neighbor_child= + neighbor->child(GeometryInfo:: + child_cell_on_face(neighbor2,subface_no)); + // As these are // quite // complicated // indirections - // we check for + // which one + // does not + // usually get + // right at + // first + // attempt we + // check for // the internal // consistency. Assert (neighbor_child->face(neighbor2) == face->child(subface_no), ExcInternalError()); Assert (!neighbor_child->has_children(), ExcInternalError()); + // We need to + // reset the + // ``un_v_matrix'' + // on each + // subface + // because on + // each subface + // the ``un'' + // belong to + // different + // neighboring + // cells. + un_v_matrix.clear(); + // As already // mentioned // above for - // this case - // (case 2) we - // employ the - // FESubfaceValues + // the current + // case (case + // 2) we employ + // the + // ``FESubfaceValues'' // of the // current - // cell, here + // cell (here // reinited for // the current // cell, face - // and subface, + // and subface) // and we // employ the // FEFaceValues @@ -946,38 +1006,27 @@ void DGMethod::assemble_system1 () u_v_matrix, un_v_matrix); - // get dof + // Then we get + // the dof // indices of // the // neighbor_child // cell neighbor_child->get_dof_indices (dofs_neighbor); + // and // distribute - // cell matrix + // ``un_v_matrix'' // to the // system_matrix for (unsigned int i=0; ihas_children())'' } - // End of ``if - // (face->has_children())'' else { // We proceed with @@ -989,19 +1038,21 @@ void DGMethod::assemble_system1 () // current cell. if (neighbor->level() == cell->level()) { - // Like before we - // store which - // number the - // current cell has - // in the list of - // neighbors of the - // neighboring - // cell. + // 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 + // ``FEFaceValues'' // of the // current and // neighboring @@ -1061,7 +1112,7 @@ void DGMethod::assemble_system1 () // Reinit the // appropriate - // FEFaceValues + // ``FEFaceValues'' // and assemble // the face // terms. @@ -1075,23 +1126,19 @@ void DGMethod::assemble_system1 () un_v_matrix); } - // Get dof indices - // of the - // neighbor_child + // Now we get the + // dof indices of + // the + // ``neighbor_child'' // cell, neighbor->get_dof_indices (dofs_neighbor); - // distribute the - // un_v_matrix, + // and distribute the + // ``un_v_matrix''. for (unsigned int i=0; i::assemble_system1 () } // Finally we distribute the - // u_v_matrix, + // ``u_v_matrix'' for (unsigned int i=0; i 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); + std::vector dofs (dofs_per_cell); + std::vector dofs_neighbor (dofs_per_cell); - UpdateFlags update_flags = UpdateFlags(update_values - | update_gradients - | update_q_points - | update_JxW_values); + UpdateFlags update_flags = 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); + UpdateFlags face_update_flags = update_values + | update_q_points + | update_JxW_values + | update_normal_vectors; + + UpdateFlags neighbor_face_update_flags = update_values; // Here we do not need // ``fe_v_face_neighbor'' as case 4 @@ -1186,15 +1229,16 @@ void DGMethod::assemble_system2 () 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'. + // Additionally we need the + // 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); @@ -1203,14 +1247,14 @@ void DGMethod::assemble_system2 () // 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; - + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); for (;cell!=endc; ++cell) { + u_v_matrix.clear (); + cell_vector.clear (); + fe_v.reinit (cell); dg.assemble_cell_term(fe_v, @@ -1221,7 +1265,8 @@ void DGMethod::assemble_system2 () for (unsigned int face_no=0; face_no::faces_per_cell; ++face_no) { - face = cell->face(face_no); + typename DoFHandler::face_iterator face= + cell->face(face_no); // Case 1: if (face->at_boundary()) @@ -1235,8 +1280,8 @@ void DGMethod::assemble_system2 () else { Assert (cell->neighbor(face_no).state() == valid, ExcInternalError()); - neighbor = cell->neighbor(face_no); - + typename DoFHandler::cell_iterator neighbor= + cell->neighbor(face_no); // Case 2: if (face->has_children()) { @@ -1247,12 +1292,17 @@ void DGMethod::assemble_system2 () subface_no::subfaces_per_face; ++subface_no) { - neighbor_child = neighbor->child( - GeometryInfo::child_cell_on_face(neighbor2,subface_no)); + typename DoFHandler::cell_iterator 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()); - + + un_v_matrix.clear(); + u_vn_matrix.clear(); + un_vn_matrix.clear(); + fe_v_subface.reinit (cell, face_no, subface_no); fe_v_face_neighbor.reinit (neighbor_child, neighbor2); @@ -1275,10 +1325,6 @@ void DGMethod::assemble_system2 () system_matrix.add(dofs_neighbor[i], dofs_neighbor[j], un_vn_matrix(i,j)); } - - un_v_matrix.clear(); - u_vn_matrix.clear(); - un_vn_matrix.clear(); } } else @@ -1290,7 +1336,11 @@ void DGMethod::assemble_system2 () neighbor->index() > cell->index()) { const unsigned int neighbor2=cell->neighbor_of_neighbor(face_no); - + + un_v_matrix.clear(); + u_vn_matrix.clear(); + un_vn_matrix.clear(); + fe_v_face.reinit (cell, face_no); fe_v_face_neighbor.reinit (neighbor, neighbor2); @@ -1313,10 +1363,6 @@ void DGMethod::assemble_system2 () system_matrix.add(dofs_neighbor[i], dofs_neighbor[j], un_vn_matrix(i,j)); } - - un_v_matrix.clear(); - u_vn_matrix.clear(); - un_vn_matrix.clear(); } // Due to rule b) @@ -1333,15 +1379,15 @@ void DGMethod::assemble_system2 () for (unsigned int i=0; i::assemble_system2 () // structur of system matrices // arising from DG // discretizations. The size of these - // blocks are the number of DoFs - // per cell. Here, we use a SSOR + // 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 @@ -1380,8 +1426,9 @@ void DGMethod::solve (Vector &solution) // 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 + // solver. Otherwise, it takes less + // memory, but the diagonal blocks + // are inverted in each // preconditioner step, // significantly slowing down the // linear solving process. @@ -1396,42 +1443,47 @@ void DGMethod::solve (Vector &solution) // 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 + // namely an approximation to the + // gradient of the solution. As here + // we consider the DG(1) method + // (i.e. we use piecewise bilinear + // shape functions) we could simply + // compute the gradients on each + // cell. But we do not want to base + // our refinement indicator on the + // gradients on each cell only, but + // want to base them also on jumps of + // the discontinuous solution + // function over faces between + // neighboring cells. The 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 + // ``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 + // ``GradientEstimation'' described + // in Step 9 of this tutorial. In + // fact, the + // ``DerivativeApproximation'' class + // was developed following the + // ``GradientEstimation'' class of + // Step 9. Relating to the + // discussion in Step 9, here we + // consider $h^{1+d/2}|\nabla_h // u_h|$. 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. + // 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 DGMethod::refine_grid () { - // The DerivativeApproximation + // The ``DerivativeApproximation'' // class computes the gradients to // float precision. This is // sufficient as they are @@ -1492,7 +1544,7 @@ void DGMethod::output_results (const unsigned int cycle) const Assert (cycle < 10, ExcInternalError()); filename += ".gnuplot"; - cout << "Writing solution to <" << filename << ">..." << endl; + cout << "Writing solution to <" << filename << ">..." << endl << endl; std::ofstream gnuplot_output (filename.c_str()); DataOut data_out; @@ -1583,8 +1635,8 @@ void DGMethod::run () // two solutions for equality. solution1-=solution2; const double difference=solution1.linfty_norm(); - if (difference<1e-13) - cout << "solution1 and solution2 do not differ." << endl; + if (difference>1e-13) + cout << "solution1 and solution2 differ!!" << endl; // Finally we perform the // output. @@ -1592,13 +1644,40 @@ void DGMethod::run () } } - - + // The following ``main'' function is + // similar to previous examples and + // need not to be commented on. int main () { - DGMethod<2> dgmethod_2d; - dgmethod_2d.run (); - + try + { + DGMethod<2> dgmethod; + dgmethod.run (); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; + return 0; }; -- 2.39.5