From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 8 Jul 2013 19:43:27 +0000 (+0000)
Subject: Patch by Denis Davydov: Add the missing functions for the Tensor extractors.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0c4171a908097a4a36c9d139a3846bf81069ff82;p=dealii-svn.git

Patch by Denis Davydov: Add the missing functions for the Tensor extractors.

git-svn-id: https://svn.dealii.org/trunk@29955 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/include/deal.II/fe/fe_values.h b/deal.II/include/deal.II/fe/fe_values.h
index 4da0529a08..b9da03e4df 100644
--- a/deal.II/include/deal.II/fe/fe_values.h
+++ b/deal.II/include/deal.II/fe/fe_values.h
@@ -4619,6 +4619,114 @@ namespace FEValuesViews
       }
   }
 
+  template <int dim, int spacedim>
+  inline
+  typename Tensor<2, dim, spacedim>::value_type
+  Tensor<2, dim, spacedim>::value (const unsigned int shape_function,
+                                   const unsigned int q_point) const
+  {
+    typedef FEValuesBase<dim,spacedim> FVB;
+    Assert (shape_function < fe_values.fe->dofs_per_cell,
+            ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell));
+    Assert (fe_values.update_flags & update_values,
+            typename FVB::ExcAccessToUninitializedField());
+
+    // similar to the vector case where we
+    // have more then one index and we need
+    // to convert between unrolled and
+    // component indexing for tensors
+    const int snc
+      = shape_function_data[shape_function].single_nonzero_component;
+
+    if (snc == -2)
+      {
+        // shape function is zero for the
+        // selected components
+        return value_type();
+
+      }
+    else if (snc != -1)
+      {
+        value_type return_value;
+        const unsigned int comp =
+          shape_function_data[shape_function].single_nonzero_component_index;
+        const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(comp);
+        return_value[indices] = fe_values.shape_values(snc,q_point);//last index first [jj][ii]
+        return return_value;
+      }
+    else
+      {
+        value_type return_value;
+        for (unsigned int d = 0; d < dim*dim; ++d)
+          if (shape_function_data[shape_function].is_nonzero_shape_function_component[d]) {
+            const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(d);
+            return_value[indices]
+              = fe_values.shape_values(shape_function_data[shape_function].row_index[d],q_point);//last index first [jj][ii]
+          }
+        return return_value;
+      }
+  }
+
+
+  template <int dim, int spacedim>
+  inline
+  typename Tensor<2, dim, spacedim>::divergence_type
+  Tensor<2, dim, spacedim>::divergence(const unsigned int shape_function,
+                                       const unsigned int q_point) const
+  {
+    typedef FEValuesBase<dim,spacedim> FVB;
+    Assert (shape_function < fe_values.fe->dofs_per_cell,
+            ExcIndexRange (shape_function, 0, fe_values.fe->dofs_per_cell));
+    Assert (fe_values.update_flags & update_gradients,
+            typename FVB::ExcAccessToUninitializedField());
+
+    const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+    if (snc == -2)
+      {
+        // shape function is zero for the
+        // selected components
+        return divergence_type();
+      }
+    else if (snc != -1)
+      {
+        // we have a single non-zero component
+        // when the tensor is
+        // represented in unrolled form.
+        //
+        // the divergence of a second-order tensor
+        // is a first order tensor.
+        //
+        // assume the second-order tensor is
+        // A with components A_{ij}.  
+        // divergence as:
+        // b_j := \dfrac{\partial phi_{ij}}{\partial x_i}.
+        //
+        // Now, we know the nonzero component
+        // in unrolled form: it is indicated
+        // by 'snc'. we can figure out which
+        // tensor components belong to this:
+        const unsigned int comp =
+          shape_function_data[shape_function].single_nonzero_component_index;
+        const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(comp);
+        const unsigned int ii = indices[0];
+        const unsigned int jj = indices[1];
+
+        const dealii::Tensor<1, spacedim> phi_grad = fe_values.shape_gradients[snc][q_point];
+
+        divergence_type return_value;
+        return_value[jj] = phi_grad[ii];
+
+        return return_value;
+
+      }
+    else
+      {
+        Assert (false, ExcNotImplemented());
+        divergence_type return_value;
+        return return_value;
+      }
+  }
 }
 
 
diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc
index 7957a09afc..23b1b8c226 100644
--- a/deal.II/source/fe/fe_values.cc
+++ b/deal.II/source/fe/fe_values.cc
@@ -1170,12 +1170,10 @@ namespace FEValuesViews
                     shape_function_data[shape_function].single_nonzero_component_index;
 
                   const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(comp);
-                  const unsigned int i = indices[0],
-                                     j = indices[1];
 
                   const double *shape_value_ptr = &shape_values(snc,0);
                   for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
-                    values[q_point][i][j] += value * *shape_value_ptr++;
+                    values[q_point][indices] += value * *shape_value_ptr++;//last index first [j][i]
                 }
               else
                 for (unsigned int d=0;
@@ -1183,13 +1181,11 @@ namespace FEValuesViews
                   if (shape_function_data[shape_function].is_nonzero_shape_function_component[d])
                     {
                       const TableIndices<2> indices = dealii::Tensor<2,spacedim>::unrolled_to_component_indices(d);
-                      const unsigned int i = indices[0],
-                                         j = indices[1];
 
                       const double *shape_value_ptr =
                         &shape_values(shape_function_data[shape_function].row_index[d],0);
                       for (unsigned int q_point=0; q_point<n_quadrature_points; ++q_point)
-                        values[q_point][i][j] += value * *shape_value_ptr++;
+                        values[q_point][indices] += value * *shape_value_ptr++;//last index first [j][i]
                     }
             }
         }
@@ -1239,10 +1235,7 @@ namespace FEValuesViews
                   for (unsigned int q_point = 0; q_point < n_quadrature_points;
                        ++q_point, ++shape_gradient_ptr)
                     {
-                      divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj];
-
-                      if (ii != jj)
-                        divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii];
+                      divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii];
                     }
                 }
               else