From 8b3eeb3caa6fc270a41c143b89d1be29bcf4a27b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 8 Jul 2013 19:43:27 +0000 Subject: [PATCH] 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 --- deal.II/include/deal.II/fe/fe_values.h | 108 +++++++++++++++++++++++++ deal.II/source/fe/fe_values.cc | 13 +-- 2 files changed, 111 insertions(+), 10 deletions(-) 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 + 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 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 + 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 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 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