From f211638c30d3b0a8daed32cbbacedcc9d7c352cb Mon Sep 17 00:00:00 2001 From: hartmann Date: Tue, 20 Nov 2001 17:55:56 +0000 Subject: [PATCH] Introduce assemble_boundary_term function. New subsections. git-svn-id: https://svn.dealii.org/trunk@5221 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-12/step-12.cc | 349 +++++++++++++++------------- 1 file changed, 183 insertions(+), 166 deletions(-) diff --git a/deal.II/examples/step-12/step-12.cc b/deal.II/examples/step-12/step-12.cc index 5bea1272c8..c34c8ef336 100644 --- a/deal.II/examples/step-12/step-12.cc +++ b/deal.II/examples/step-12/step-12.cc @@ -185,6 +185,8 @@ void BoundaryValues::value_list(const std::vector > &points, } } + // @sect3{Class: DGTransportEquation} + // // Next we define the equation- // dependent and DG-method-dependent // class ``DGTransportEquation''. Its @@ -206,25 +208,32 @@ class DGTransportEquation FullMatrix &u_v_matrix, Vector &cell_vector); + void assemble_boundary_term(const FEFaceValues& 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); + FullMatrix &un_v_matrix); 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); + FullMatrix &un_vn_matrix); private: Beta beta_function; RHS rhs_function; BoundaryValues boundary_function; }; + // @sect4{Function: assemble_cell_term} + // + // The ``assemble_cell_term'' + // function assembles the cell terms + // of the discretization. // ``u_v_matrix'' is a cell matrix, // i.e. for a DG method of degree 1, // it is of size 4 times 4, and @@ -240,7 +249,7 @@ void DGTransportEquation::assemble_cell_term( Vector &cell_vector) { // First we ask ``fe_v'' for the - // shape grads, shape values and + // shape gradients, shape values and // quadrature weights, const vector > > &grad_v = fe_v.get_shape_grads (); const FullMatrix &v = fe_v.get_shape_values (); @@ -256,11 +265,9 @@ void DGTransportEquation::assemble_cell_term( 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. + // vector are assembled due to the + // terms $-(u,\beta\cdot\nabla + // v)_K$ and $(f,v)_K$. for (unsigned int point=0; point::assemble_cell_term( } } + // @sect4{Function: assemble_boundary_term} + // + // The ``assemble_boundary_term'' + // function assembles the face terms + // at boundary faces. When this + // function is invoked, ``fe_v'' is + // already reinited with the current + // cell and current face. Hence it + // provides the shape values on that + // boundary face. +template +void DGTransportEquation::assemble_boundary_term( + const FEFaceValues& fe_v, + FullMatrix &u_v_matrix, + Vector &cell_vector) +{ + // First we check whether the + // current face is really at the + // boundary. + Assert(fe_v.get_face()->at_boundary(), ExcInternalError()); + + // Again, as in the previous + // 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 (); + // but here also for the normals. + const 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); + + beta_function.value_list (fe_v.get_quadrature_points(), beta); + boundary_function.value_list (fe_v.get_quadrature_points(), g); + + // Then we assemble cell vector and + // cell matrix according to the DG + // method given in the + // introduction. + for (unsigned int point=0; point0) + // We assemble the term + // $(\beta\cdot n + // u,v)_{\partial K_+}$, + for (unsigned int i=0; i::assemble_cell_term( // on the face. // // In addition to the cell matrix - // ``u_v_matrix'' and the - // ``cell_vector'' this function has + // ``u_v_matrix'' this function has // got a new argument // ``un_v_matrix'', that stores // contributions to the system matrix @@ -309,41 +384,38 @@ 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) + FullMatrix &un_v_matrix) { - // Again, we ask the FEValues - // objects for the shape values and - // the quadrature weights + // First we check that the current + // face is not at the boundary by + // accident. + Assert(!fe_v.get_face()->at_boundary(), ExcInternalError()); + + // Again, as in the previous + // function, we ask the FEValues + // objects for the shape values, + // the quadrature weights and the + // 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 (); - // but also for the normals. const vector > &normals = fe_v.get_normal_vectors (); - // We also evaluate the flow field - // at the quadrature points + // and we evaluate the flow field + // at the quadrature points. vector > beta (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 + // Then we assemble the cell + // matrices according to the DG + // method given in the // introduction. for (unsigned int point=0; point0) - // The term $(\beta\cdot n + // We assemble the term + // $(\beta\cdot n // u,v)_{\partial K_+}$, for (unsigned int i=0; i::assemble_face_term1( v(i,point) * JxW[point]; else - { - // at the boundary the term - // $(\beta\cdot n - // g,v)_{\partial - // K_-\cap\partial\Omega}$, - if (face->at_boundary()) - 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, - Vector &cell_vector) + FullMatrix &un_vn_matrix) { // the first few lines are the same const FullMatrix &v = fe_v.get_shape_values (); @@ -411,17 +472,12 @@ void DGTransportEquation::assemble_face_term2( 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, + // This terms we've already seen. for (unsigned int i=0; i::assemble_face_term2( v(i,point) * JxW[point]; - // on inner faces we - // additionally have the - // term $(\beta\cdot n - // u,\hat v)_{\partial K_+}, - if (!face->at_boundary()) - for (unsigned int k=0; kat_boundary()) - for (unsigned int i=0; i::setup_system () }; + // @sect4{Function: assemble_system1} + // // 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_cell_term'', + // ``assemble_boundary_term'' and // ``assemble_face_term1'' functions // of the DGTransportEquation object. // The ``assemble_face_term1'' @@ -623,8 +673,7 @@ void DGMethod::setup_system () // introduction: // // 1. face is at boundary (current - // cell: FEFaceValues, neighboring - // cell does not exist); + // cell: FEFaceValues); // // 2. neighboring cell is finer // (current cell: FESubfaceValues, @@ -788,33 +837,17 @@ void DGMethod::assemble_system1 () // 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); + // terms. + dg.assemble_boundary_term(fe_v_face, + u_v_matrix, + cell_vector); } else { - // When we are not on the - // boundary of the - // domain then there - // must exist a + // Now we are not on + // the boundary of the + // domain, therefore + // there must exist a // neighboring cell. neighbor = cell->neighbor(face_no); @@ -911,8 +944,7 @@ void DGMethod::assemble_system1 () dg.assemble_face_term1(fe_v_subface, fe_v_face_neighbor, u_v_matrix, - un_v_matrix, - cell_vector); + un_v_matrix); // get dof // indices of @@ -985,8 +1017,7 @@ void DGMethod::assemble_system1 () dg.assemble_face_term1(fe_v_face, fe_v_face_neighbor, u_v_matrix, - un_v_matrix, - cell_vector); + un_v_matrix); // End of ``if // (neighbor->level() // == @@ -1041,8 +1072,7 @@ void DGMethod::assemble_system1 () dg.assemble_face_term1(fe_v_face, fe_v_subface_neighbor, u_v_matrix, - un_v_matrix, - cell_vector); + un_v_matrix); } // Get dof indices @@ -1085,7 +1115,8 @@ void DGMethod::assemble_system1 () }; - + // @sect4{Function: assemble_system2} + // // We proceed with the // ``assemble_system2'' function that // implements the DG discretization @@ -1169,14 +1200,6 @@ void DGMethod::assemble_system2 () 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. @@ -1205,13 +1228,9 @@ void DGMethod::assemble_system2 () { 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); + dg.assemble_boundary_term(fe_v_face, + u_v_matrix, + cell_vector); } else { @@ -1242,8 +1261,7 @@ void DGMethod::assemble_system2 () u_v_matrix, un_v_matrix, u_vn_matrix, - un_vn_matrix, - dummy_rhs); + un_vn_matrix); neighbor_child->get_dof_indices (dofs_neighbor); @@ -1281,8 +1299,7 @@ void DGMethod::assemble_system2 () u_v_matrix, un_v_matrix, u_vn_matrix, - un_vn_matrix, - dummy_rhs); + un_vn_matrix); neighbor->get_dof_indices (dofs_neighbor); -- 2.39.5