From d5fd2325fe159f3e03489674a2e6bf8c950654c7 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler.martin@gmail.com>
Date: Mon, 1 Sep 2014 22:10:58 +0200
Subject: [PATCH] Avoid setting the templated 'Number' to '0' (double) in
 trace(Tensor<2,dim,Number>).

Enables using trace also with e.g. Tensor<2,dim,VectorizedArray<double> >.
---
 include/deal.II/base/symmetric_tensor.h       |  4 ++--
 include/deal.II/base/tensor.h                 |  4 ++--
 tests/matrix_free/matrix_vector_stokes.cc     | 23 +++++++++----------
 .../matrix_vector_stokes_noflux.cc            |  2 +-
 4 files changed, 16 insertions(+), 17 deletions(-)

diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h
index 67710ad64e..aa913b69f8 100644
--- a/include/deal.II/base/symmetric_tensor.h
+++ b/include/deal.II/base/symmetric_tensor.h
@@ -2250,8 +2250,8 @@ double third_invariant (const SymmetricTensor<2,dim,Number> &t)
 template <int dim, typename Number>
 Number trace (const SymmetricTensor<2,dim,Number> &d)
 {
-  Number t=0;
-  for (unsigned int i=0; i<dim; ++i)
+  Number t = d.data[0];
+  for (unsigned int i=1; i<dim; ++i)
     t += d.data[i];
   return t;
 }
diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index 05e71b4982..cab217ad2b 100644
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -1758,8 +1758,8 @@ Number determinant (const Tensor<2,dim,Number> &t)
 template <int dim, typename Number>
 Number trace (const Tensor<2,dim,Number> &d)
 {
-  Number t=0;
-  for (unsigned int i=0; i<dim; ++i)
+  Number t=d[0][0];
+  for (unsigned int i=1; i<dim; ++i)
     t += d[i][i];
   return t;
 }
diff --git a/tests/matrix_free/matrix_vector_stokes.cc b/tests/matrix_free/matrix_vector_stokes.cc
index 46e1463c3a..e2cc02804c 100644
--- a/tests/matrix_free/matrix_vector_stokes.cc
+++ b/tests/matrix_free/matrix_vector_stokes.cc
@@ -85,17 +85,16 @@ public:
 
         for (unsigned int q=0; q<velocity.n_q_points; ++q)
           {
-            SymmetricTensor<2,dim,vector_t> sym_grad_u =
-              velocity.get_symmetric_gradient (q);
+            Tensor<2,dim,vector_t> grad_u = velocity.get_gradient (q);
             vector_t pres = pressure.get_value(q);
-            vector_t div = -velocity.get_divergence(q);
+            vector_t div = -trace(grad_u);
             pressure.submit_value   (div, q);
 
             // subtract p * I
             for (unsigned int d=0; d<dim; ++d)
-              sym_grad_u[d][d] -= pres;
+              grad_u[d][d] -= pres;
 
-            velocity.submit_symmetric_gradient(sym_grad_u, q);
+            velocity.submit_gradient(grad_u, q);
           }
 
         velocity.integrate (false,true);
@@ -221,9 +220,9 @@ void test ()
     const FEValuesExtractors::Vector velocities (0);
     const FEValuesExtractors::Scalar pressure (dim);
 
-    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<Tensor<2,dim> > phi_grad_u (dofs_per_cell);
+    std::vector<double>         div_phi_u  (dofs_per_cell);
+    std::vector<double>         phi_p      (dofs_per_cell);
 
     typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
@@ -237,16 +236,16 @@ void test ()
           {
             for (unsigned int k=0; k<dofs_per_cell; ++k)
               {
-                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_grad_u[k] = fe_values[velocities].gradient (k, q);
+                div_phi_u[k]  = fe_values[velocities].divergence (k, q);
+                phi_p[k]      = fe_values[pressure].value (k, q);
               }
 
             for (unsigned int i=0; i<dofs_per_cell; ++i)
               {
                 for (unsigned int j=0; j<=i; ++j)
                   {
-                    local_matrix(i,j) += (phi_grads_u[i] * phi_grads_u[j]
+                    local_matrix(i,j) += (scalar_product(phi_grad_u[i], phi_grad_u[j])
                                           - div_phi_u[i] * phi_p[j]
                                           - phi_p[i] * div_phi_u[j])
                                          * fe_values.JxW(q);
diff --git a/tests/matrix_free/matrix_vector_stokes_noflux.cc b/tests/matrix_free/matrix_vector_stokes_noflux.cc
index 1860f6797e..6da9d1f55a 100644
--- a/tests/matrix_free/matrix_vector_stokes_noflux.cc
+++ b/tests/matrix_free/matrix_vector_stokes_noflux.cc
@@ -87,7 +87,7 @@ public:
             SymmetricTensor<2,dim,vector_t> sym_grad_u =
               velocity.get_symmetric_gradient (q);
             vector_t pres = pressure.get_value(q);
-            vector_t div = -velocity.get_divergence(q);
+            vector_t div = -trace(sym_grad_u);
             pressure.submit_value   (div, q);
 
             // subtract p * I
-- 
2.39.5