From 1001afeb494907d6826a99c1ea72b4c33000e72a Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 14 Aug 2008 16:25:46 +0000 Subject: [PATCH] Compute BDF-2 matrices. git-svn-id: https://svn.dealii.org/trunk@16549 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 104 +++++++++++++++++++--------- 1 file changed, 72 insertions(+), 32 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index a5fa1d70d1..308502274b 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -1765,7 +1765,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () { // TODO: right now, always do explicit // Euler - const bool use_bdf2_scheme = (timestep_number == 0 ? false : false); + const bool use_bdf2_scheme = (timestep_number == 0 ? true : false); system_matrix.block(2,2) = 0; @@ -1796,22 +1796,26 @@ void BoussinesqFlowProblem::assemble_rhs_T () // well-known extractor for // accessing the temperature // part of the FE system. - std::vector > old_solution_values(n_q_points, - Vector(dim+2)); std::vector > present_solution_values (n_q_points, Vector(dim+2)); - std::vector > > old_solution_grads( - n_q_points, - std::vector >(dim+2)); - - std::vector neighbor_temperature (n_face_q_points); + std::vector > old_solution_values(n_q_points, + Vector(dim+2)); + std::vector > old_old_solution_values(n_q_points, + Vector(dim+2)); + std::vector > > + old_solution_grads(n_q_points, + std::vector >(dim+2)); + std::vector > > + old_old_solution_grads(n_q_points, + std::vector >(dim+2)); TemperatureBoundaryValues temperature_boundary_values; + const FEValuesExtractors::Scalar temperature (dim+1); std::vector phi_T (dofs_per_cell); std::vector > grad_phi_T (dofs_per_cell); - + // Now, let's start the loop // over all cells in the // triangulation. The first @@ -1830,9 +1834,13 @@ void BoussinesqFlowProblem::assemble_rhs_T () local_rhs = 0; fe_values.reinit (cell); + fe_values.get_function_values (solution, present_solution_values); fe_values.get_function_values (old_solution, old_solution_values); + fe_values.get_function_values (old_old_solution, old_old_solution_values); + fe_values.get_function_gradients (old_solution, old_solution_grads); - fe_values.get_function_values (solution, present_solution_values); + fe_values.get_function_gradients (old_old_solution, old_old_solution_grads); + // build matrix contributions local_matrix = 0; @@ -1846,6 +1854,9 @@ void BoussinesqFlowProblem::assemble_rhs_T () const double kappa = std::max (5e-3 * cell->diameter(), 1e-5); + const double artificial_diffusion = 0; + + for (unsigned int q=0; q::assemble_rhs_T () const Point p = fe_values.quadrature_point(q); const double gamma = RightHandSide().value (p, dim+1); + + const double old_T = old_solution_values[q](dim+1); + const double old_old_T = old_old_solution_values[q](dim+1); + + const Tensor<1,dim> old_grad_T = old_solution_grads[q][dim+1]; + const Tensor<1,dim> old_old_grad_T = old_old_solution_grads[q][dim+1]; + + Tensor<1,dim> present_u; + for (unsigned int d=0; d old_grad_T = old_solution_grads[q][dim+1]; - - Tensor<1,dim> present_u; - for (unsigned int d=0; d