From 9a624030a477adcb9aca21050bce5a869195ec30 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Wed, 13 Aug 2008 20:03:08 +0000
Subject: [PATCH] Don't integrate by parts the terms on the temperature rhs.

git-svn-id: https://svn.dealii.org/trunk@16536 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/examples/step-31/step-31.cc | 29 ++++++++++-------------------
 1 file changed, 10 insertions(+), 19 deletions(-)

diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc
index cd1118a02b..2e7a91ea4f 100644
--- a/deal.II/examples/step-31/step-31.cc
+++ b/deal.II/examples/step-31/step-31.cc
@@ -1801,7 +1801,7 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
 							n_face_q_points, 
 							Vector<double>(dim+2));
 
-  std::vector<std::vector<Tensor<1,dim> > >  present_solution_grads(
+  std::vector<std::vector<Tensor<1,dim> > >  old_solution_grads(
 				  n_q_points,
 				  std::vector<Tensor<1,dim> >(dim+2));
 
@@ -1829,35 +1829,26 @@ void BoussinesqFlowProblem<dim>::assemble_rhs_T ()
       fe_values.reinit (cell);
 
       fe_values.get_function_values (old_solution, 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 (solution, present_solution_grads);
 
       for (unsigned int q=0; q<n_q_points; ++q)
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           {
-            const double old_T = old_solution_values[q](dim+1);
+            const double        old_T      = old_solution_values[q](dim+1);
+	    const Tensor<1,dim> old_grad_T = old_solution_grads[q][dim+1];
+	    
             Tensor<1,dim> present_u;
             for (unsigned int d=0; d<dim; ++d)
               present_u[d] = present_solution_values[q](d);
 
-	    double present_div_u = 0;
-            for (unsigned int d=0; d<dim; ++d)
-              present_div_u += present_solution_grads[q][d][d];
-
-            const double        phi_i_T      = fe_values[temperature].value (i, q);
-            const Tensor<1,dim> grad_phi_i_T = fe_values[temperature].gradient (i, q);
-
+            const double     phi_i_T      = fe_values[temperature].value (i, q);
 	    const Point<dim> p = fe_values.quadrature_point(q);
 
-            local_rhs(i) += (time_step *
-                             old_T *
-                             (present_u *
-			      grad_phi_i_T
-			      +
-			      present_div_u *
-			      phi_i_T)
-                             +
-                             old_T * phi_i_T
+            local_rhs(i) += (old_T * phi_i_T
+			     -
+			     time_step *
+                             present_u * old_grad_T * phi_i_T
 			     +
 			     time_step *
 			     RightHandSide<dim>().value (p, dim+1)
-- 
2.39.5