From 4d6efbc17820cb4d8b6c0ff47e6267be006b9226 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 1 Mar 2018 10:05:56 -0700 Subject: [PATCH] Update some code in step-12. Specifically: * Write the bilinear form with direction (-beta) obtained through integration by parts, rather than using beta and then doing -= on the matrix. * Use the names for FE*Values objects that we use everywhere else. * Some other minor cleanups. --- examples/step-12/step-12.cc | 112 ++++++++++++++++++------------------ 1 file changed, 57 insertions(+), 55 deletions(-) diff --git a/examples/step-12/step-12.cc b/examples/step-12/step-12.cc index ab59e5a21b..d91bb2d385 100644 --- a/examples/step-12/step-12.cc +++ b/examples/step-12/step-12.cc @@ -321,26 +321,27 @@ namespace Step12 // First, let us retrieve some of the objects used here from @p info. Note // that these objects can handle much more complex structures, thus the // access here looks more complicated than might seem necessary. - const FEValuesBase &fe_v = info.fe_values(); + const FEValuesBase &fe_values = info.fe_values(); FullMatrix &local_matrix = dinfo.matrix(0).matrix; - const std::vector &JxW = fe_v.get_JxW_values (); + const std::vector &JxW = fe_values.get_JxW_values (); // With these objects, we continue local integration like always. First, // we loop over the quadrature points and compute the advection vector in // the current point. - for (unsigned int point=0; point beta; - beta(0) = -fe_v.quadrature_point(point)(1); - beta(1) = fe_v.quadrature_point(point)(0); + beta(0) = -fe_values.quadrature_point(point)(1); + beta(1) = fe_values.quadrature_point(point)(0); beta /= beta.norm(); // We solve a homogeneous equation, thus no right hand side shows up // in the cell term. What's left is integrating the matrix entries. - for (unsigned int i=0; i::integrate_boundary_term (DoFInfo &dinfo, CellInfo &info) { - const FEValuesBase &fe_v = info.fe_values(); + const FEValuesBase &fe_face_values = info.fe_values(); FullMatrix &local_matrix = dinfo.matrix(0).matrix; Vector &local_vector = dinfo.vector(0).block(0); - const std::vector &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); + const std::vector &JxW = fe_face_values.get_JxW_values (); + const std::vector > &normals = fe_face_values.get_normal_vectors (); - std::vector g(fe_v.n_quadrature_points); + std::vector g(fe_face_values.n_quadrature_points); static BoundaryValues boundary_function; - boundary_function.value_list (fe_v.get_quadrature_points(), g); + boundary_function.value_list (fe_face_values.get_quadrature_points(), g); - for (unsigned int point=0; point beta; - beta(0) = -fe_v.quadrature_point(point)(1); - beta(1) = fe_v.quadrature_point(point)(0); + beta(0) = -fe_face_values.quadrature_point(point)(1); + beta(1) = fe_face_values.quadrature_point(point)(0); beta /= beta.norm(); const double beta_n=beta * normals[point]; if (beta_n>0) - for (unsigned int i=0; i &fe_v = info1.fe_values(); + const FEValuesBase &fe_face_values = info1.fe_values(); // For additional shape functions, we have to ask the neighbors // FEValuesBase. - const FEValuesBase &fe_v_neighbor = info2.fe_values(); + const FEValuesBase &fe_face_values_neighbor = info2.fe_values(); // Then we get references to the four local matrices. The letters u and v // refer to trial and test functions, respectively. The %numbers indicate @@ -420,53 +421,53 @@ namespace Step12 // hand side vectors. Fortunately, the interface terms only involve the // solution and the right hand side does not receive any contributions. - const std::vector &JxW = fe_v.get_JxW_values (); - const std::vector > &normals = fe_v.get_normal_vectors (); + const std::vector &JxW = fe_face_values.get_JxW_values (); + const std::vector > &normals = fe_face_values.get_normal_vectors (); - for (unsigned int point=0; point beta; - beta(0) = -fe_v.quadrature_point(point)(1); - beta(1) = fe_v.quadrature_point(point)(0); + beta(0) = -fe_face_values.quadrature_point(point)(1); + beta(1) = fe_face_values.quadrature_point(point)(0); beta /= beta.norm(); - const double beta_n=beta * normals[point]; - if (beta_n>0) + const double beta_dot_n = beta * normals[point]; + if (beta_dot_n>0) { // This term we've already seen: - for (unsigned int i=0; i void AdvectionProblem::refine_grid () { -- 2.39.5