From e27565f1c4e7e728350a7a44bee957b4d4595dd0 Mon Sep 17 00:00:00 2001 From: bangerth Date: Thu, 14 Aug 2008 15:59:02 +0000 Subject: [PATCH] Move building the temperature matrix to the function that also builds the right hand side. This ensures that we use the correct time step, and also that we can next build the BDF-2 matrix. git-svn-id: https://svn.dealii.org/trunk@16546 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-31/step-31.cc | 76 ++++++++++++++++++++--------- 1 file changed, 54 insertions(+), 22 deletions(-) diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc index 3a452c2fd4..ec62724e9c 100644 --- a/deal.II/examples/step-31/step-31.cc +++ b/deal.II/examples/step-31/step-31.cc @@ -1450,8 +1450,6 @@ void BoussinesqFlowProblem::assemble_system () std::vector > grads_phi_u (dofs_per_cell); std::vector div_phi_u (dofs_per_cell); std::vector phi_p (dofs_per_cell); - std::vector phi_T (dofs_per_cell); - std::vector > grad_phi_T (dofs_per_cell); const FEValuesExtractors::Vector velocities (0); const FEValuesExtractors::Scalar pressure (dim); @@ -1518,8 +1516,6 @@ void BoussinesqFlowProblem::assemble_system () grads_phi_u[k] = fe_values[velocities].symmetric_gradient(k,q); div_phi_u[k] = fe_values[velocities].divergence (k, q); phi_p[k] = fe_values[pressure].value (k, q); - phi_T[k] = fe_values[temperature].value (k, q); - grad_phi_T[k] = fe_values[temperature].gradient (k, q); } } @@ -1527,27 +1523,22 @@ void BoussinesqFlowProblem::assemble_system () { const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q); - // define viscosity and - // diffusion. for the - // latter, take the - // maximum of what we - // really want and the - // minimal amount of - // diffusion - // (determined - // impirically) to keep - // the scheme stable - const double eta = 1, - kappa = std::max (5e-4 * cell->diameter(), - 1e-6); - + // define viscosity + const double eta = 1; + + // build Stokes part of + // the matrix. we have + // to rebuild the + // temperature part of + // it in each time step + // due to the fact that + // it depends on the + // time step size if (rebuild_matrices) for (unsigned int j=0; j gravity = ( (dim == 2) ? (Point (0,1)) : @@ -1775,7 +1766,8 @@ void BoussinesqFlowProblem::assemble_rhs_T () // TODO: right now, always do explicit // Euler const bool is_first_timestep = (timestep_number == 0 ? true : true); - + + system_matrix.block(2,2) = 0; QGauss quadrature_formula(degree+2); QGauss face_quadrature_formula(degree+2); @@ -1788,6 +1780,7 @@ void BoussinesqFlowProblem::assemble_rhs_T () const unsigned int n_face_q_points = face_quadrature_formula.size(); Vector local_rhs (dofs_per_cell); + FullMatrix local_matrix (dofs_per_cell, dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); @@ -1816,6 +1809,9 @@ void BoussinesqFlowProblem::assemble_rhs_T () 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 @@ -1838,6 +1834,35 @@ void BoussinesqFlowProblem::assemble_rhs_T () fe_values.get_function_gradients (old_solution, old_solution_grads); fe_values.get_function_values (solution, present_solution_values); + // build matrix contributions + local_matrix = 0; + + // define diffusion. take the + // maximum of what we really + // want and the minimal amount + // of diffusion (determined + // impirically) to keep the + // scheme stable + const double kappa = std::max (5e-4 * cell->diameter(), + 1e-6); + + for (unsigned int q=0; q::assemble_rhs_T () } cell->get_dof_indices (local_dof_indices); + + for (unsigned int i=0; i