From f8cc05a619ee86ff86fda0e53d61f374359f38a6 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Mon, 21 Apr 2008 15:42:19 +0000
Subject: [PATCH] Implemented a faster matrix assembly routine, following the
 one in step-22.

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

diff --git a/deal.II/examples/step-31/step-31.cc b/deal.II/examples/step-31/step-31.cc
index e4f3b1368c..c7b2bbb6cf 100644
--- a/deal.II/examples/step-31/step-31.cc
+++ b/deal.II/examples/step-31/step-31.cc
@@ -501,6 +501,13 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 
   const double Raleigh_number = 10;
 
+  std::vector<Tensor<1,dim> >          phi_u       (dofs_per_cell);
+  std::vector<SymmetricTensor<2,dim> > phi_grads_u (dofs_per_cell);
+  std::vector<double>                  div_phi_u   (dofs_per_cell);
+  std::vector<double>                  phi_p       (dofs_per_cell);
+  std::vector<double>                  phi_T       (dofs_per_cell);
+  std::vector<Tensor<1,dim> >          grad_phi_T  (dofs_per_cell);
+
   const FEValuesExtractors::Vector velocities (0);
   const FEValuesExtractors::Scalar pressure (dim);
   const FEValuesExtractors::Scalar temperature (dim+1);
@@ -520,47 +527,54 @@ void BoussinesqFlowProblem<dim>::assemble_system ()
 	{
 	  const double old_temperature = old_solution_values[q](dim+1);
 
+				// Extract the basis relevant
+				// terms in the inner products
+				// once in advance as shown
+				// in step-22. This accelerates
+				// the assembly process,
+	  for (unsigned int k=0; k<dofs_per_cell; ++k)
+	    {
+	      phi_u[k] = fe_values[velocities].value (k,q);
+	      if (rebuild_matrices)
+	        {
+		  phi_grads_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);
+		}
+	    }
+
 	  for (unsigned int i=0; i<dofs_per_cell; ++i)
 	    {
 
 	      const Tensor<1,dim> phi_i_u = fe_values[velocities].value (i, q);
 
 	      if (rebuild_matrices)
-		{
-		  const SymmetricTensor<2,dim>
-		    phi_i_grads_u = fe_values[velocities].symmetric_gradient (i, q);
-		  const Tensor<1,dim> phi_i_u       = fe_values[velocities].value (i, q);
-		  const double        div_phi_i_u   = fe_values[velocities].divergence (i, q);
-		  const double        phi_i_p       = fe_values[pressure].value (i, q);
-		  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);
-
-		  for (unsigned int j=0; j<dofs_per_cell; ++j)
-		    {
-		      const SymmetricTensor<2,dim>
-			phi_j_grads_u = fe_values[velocities].symmetric_gradient (j, q);
-		      const double        div_phi_j_u = fe_values[velocities].divergence (j, q);
-		      const double        phi_j_p     = fe_values[pressure].value (j, q);
-		      const double        phi_j_T     = fe_values[temperature].value (j, q);
-
-		      local_matrix(i,j) += (phi_i_grads_u * phi_j_grads_u
-					    - div_phi_i_u * phi_j_p
-					    - phi_i_p * div_phi_j_u
-					    + phi_i_p * phi_j_p
-					    + phi_i_T * phi_j_T)
-					   * fe_values.JxW(q);
-		    }
-		}
+		for (unsigned int j=0; j<dofs_per_cell; ++j)
+		  local_matrix(i,j) += (phi_grads_u[i] * phi_grads_u[j]
+					- div_phi_u[i] * phi_p[j]
+					- phi_p[i] * div_phi_u[j]
+					+ phi_p[i] * phi_p[j]
+					+ phi_T[i] * phi_T[j])
+				       * fe_values.JxW(q);
 
 	      const Point<dim> gravity (0,1);
 
 	      local_rhs(i) += (Raleigh_number *
-			       gravity * phi_i_u * old_temperature)*
+			       gravity * phi_u[i] * old_temperature)*
 			      fe_values.JxW(q);
           }
 	}
 
 
+				// The assembly of the face
+				// cells which enters the
+				// right hand sides cannot
+				// be accelerated with the
+				// above technique, since
+				// all the basis functions are
+				// only evaluated once.
       for (unsigned int face_no=0;
            face_no<GeometryInfo<dim>::faces_per_cell;
            ++face_no)
-- 
2.39.5