From: Denis Davydov Date: Thu, 12 Apr 2018 18:36:29 +0000 (+0200) Subject: implement gradient of FEValuesViews::Tensor X-Git-Tag: v9.0.0-rc1~185^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=770c819616db1a0063dc70f53aba3d47f9647d5e;p=dealii.git implement gradient of FEValuesViews::Tensor Also change meaning of divergence so that Grad(T) : I = Div(T). --- diff --git a/include/deal.II/fe/fe_values.h b/include/deal.II/fe/fe_values.h index 37878b2a7d..406404588f 100644 --- a/include/deal.II/fe/fe_values.h +++ b/include/deal.II/fe/fe_values.h @@ -1376,17 +1376,18 @@ namespace FEValuesViews * @ref vector_valued * module. * - * This class allows to query the value and divergence of (components of) + * This class allows to query the value, gradient and divergence of (components of) * shape functions and solutions representing tensors. The divergence of a - * tensor $T_{ij}, 0\le i,j<\text{dim}$ is defined as $d_i = \sum_j - * \frac{\partial T_{ji}}{\partial x_j}, 0\le i<\text{dim}$. + * tensor $T_{ij},\, 0\le i,j<\text{dim}$ is defined as $d_i = \sum_j + * \frac{\partial T_{ij}}{\partial x_j}, \, 0\le i<\text{dim}$, + * whereas its gradient is $G_{ijk} = \frac{\partial T_{ij}}{\partial x_k}$. * * You get an object of this type if you apply a FEValuesExtractors::Tensor * to an FEValues, FEFaceValues or FESubfaceValues object. * * @ingroup feaccess vector_valued * - * @author Denis Davydov, 2013 + * @author Denis Davydov, 2013, 2018 */ template class Tensor<2,dim,spacedim> @@ -1404,6 +1405,11 @@ namespace FEValuesViews */ typedef dealii::Tensor<1, spacedim> divergence_type; + /** + * Data type for taking the gradient of a second order tensor: a third order tensor. + */ + typedef dealii::Tensor<3, spacedim> gradient_type; + /** * A struct that provides the output type for the product of the value * and derivatives of basis functions of the Tensor view and any @p Number type. @@ -1422,6 +1428,12 @@ namespace FEValuesViews * divergences of the view the Tensor class. */ typedef typename ProductType::divergence_type>::type divergence_type; + + /** + * A typedef for the data type of the product of a @p Number and the + * gradient of the view the Tensor class. + */ + typedef typename ProductType::gradient_type>::type gradient_type; }; /** @@ -1528,6 +1540,23 @@ namespace FEValuesViews divergence (const unsigned int shape_function, const unsigned int q_point) const; + /** + * Return the gradient (3-rd order tensor) of the vector components selected by this + * view, for the shape function and quadrature point selected by the + * arguments. + * + * See the general discussion of this class for a definition of the + * gradient. + * + * @note The meaning of the arguments is as documented for the value() + * function. + * + * @dealiiRequiresUpdateFlags{update_gradients} + */ + gradient_type + gradient (const unsigned int shape_function, + const unsigned int q_point) const; + /** * Return the values of the selected vector components of the finite * element function characterized by fe_function at the @@ -1603,8 +1632,34 @@ namespace FEValuesViews */ template void get_function_divergences_from_local_dof_values (const InputVector &dof_values, - std::vector::divergence_type> &values) const; + std::vector::divergence_type> &divergences) const; + /** + * Return the gradient of the selected vector components of the finite + * element function characterized by fe_function at the + * quadrature points of the cell, face or subface selected the last time + * the reinit function of the FEValues object was called. + * + * See the general discussion of this class for a definition of the + * gradient. + * + * The data type stored by the output vector must be what you get when you + * multiply the gradients of shape functions (i.e., @p gradient_type) + * times the type used to store the values of the unknowns $U_j$ of your + * finite element vector $U$ (represented by the @p fe_function argument). + * + * @dealiiRequiresUpdateFlags{update_gradients} + */ + template + void get_function_gradients (const InputVector &fe_function, + std::vector::type> &gradients) const; + + /** + * @copydoc FEValuesViews::Tensor<2,dim,spacedim>::get_function_values_from_local_dof_values() + */ + template + void get_function_gradients_from_local_dof_values (const InputVector &dof_values, + std::vector::gradient_type> &gradients) const; private: /** @@ -4227,7 +4282,8 @@ namespace FEValuesViews const dealii::Tensor<1, spacedim> &phi_grad = fe_values->finite_element_output.shape_gradients[snc][q_point]; divergence_type return_value; - return_value[jj] = phi_grad[ii]; + // note that we contract \nabla from the right + return_value[ii] = phi_grad[jj]; return return_value; } @@ -4238,6 +4294,58 @@ namespace FEValuesViews return return_value; } } + + template + inline + typename Tensor<2, dim, spacedim>::gradient_type + Tensor<2, dim, spacedim>::gradient(const unsigned int shape_function, + const unsigned int q_point) const + { + 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 FEValuesBase::ExcAccessToUninitializedField("update_gradients"))); + + const int snc = shape_function_data[shape_function].single_nonzero_component; + + if (snc == -2) + { + // shape function is zero for the selected components + return gradient_type(); + } + else if (snc != -1) + { + // we have a single non-zero component when the tensor is + // represented in unrolled form. + // + // the gradient of a second-order tensor is a third order tensor. + // + // assume the second-order tensor is A with components A_{ij}, + // then gradient is B_{ijk} := \frac{\partial A_{ij}}{\partial x_k} + // + // 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->finite_element_output.shape_gradients[snc][q_point]; + + gradient_type return_value; + return_value[ii][jj] = phi_grad; + + return return_value; + } + else + { + Assert (false, ExcNotImplemented()); + gradient_type return_value; + return return_value; + } + } + } diff --git a/include/deal.II/fe/fe_values_extractors.h b/include/deal.II/fe/fe_values_extractors.h index c086a94d1d..0e10b6569d 100644 --- a/include/deal.II/fe/fe_values_extractors.h +++ b/include/deal.II/fe/fe_values_extractors.h @@ -168,7 +168,7 @@ namespace FEValuesExtractors /** - * Extractor for a (possible non-)symmetric tensor of a rank specified by + * Extractor for a general tensor of a given rank specified by * the template argument. For a second order tensor, this represents a * collection of (dim*dim) components of a vector-valued * element. The value of dim is defined by the FEValues object diff --git a/source/fe/fe_values.cc b/source/fe/fe_values.cc index b4518e7031..ddbfb4888c 100644 --- a/source/fe/fe_values.cc +++ b/source/fe/fe_values.cc @@ -1308,7 +1308,69 @@ namespace FEValuesViews for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point, ++shape_gradient_ptr) { - divergences[q_point][jj] += value * (*shape_gradient_ptr)[ii]; + divergences[q_point][ii] += value * (*shape_gradient_ptr)[jj]; + } + } + else + { + for (unsigned int d = 0; + d < dim*dim; ++d) + if (shape_function_data[shape_function].is_nonzero_shape_function_component[d]) + { + Assert (false, ExcNotImplemented()); + } + } + } + } + + template + void + do_function_gradients (const ArrayView &dof_values, + const Table<2,dealii::Tensor<1,spacedim> > &shape_gradients, + const std::vector::ShapeFunctionData> &shape_function_data, + std::vector::template OutputType::gradient_type> &gradients) + { + const unsigned int dofs_per_cell = dof_values.size(); + const unsigned int n_quadrature_points = dofs_per_cell > 0 ? + shape_gradients[0].size() : gradients.size(); + AssertDimension (gradients.size(), n_quadrature_points); + + std::fill (gradients.begin(), gradients.end(), + typename Tensor<2,dim,spacedim>::template OutputType::gradient_type()); + + for (unsigned int shape_function=0; + shape_function::value) + if (value == dealii::internal::NumberType::value(0.0)) + continue; + + if (snc != -1) + { + const unsigned int comp = + shape_function_data[shape_function].single_nonzero_component_index; + + const dealii::Tensor < 1, spacedim> *shape_gradient_ptr = + &shape_gradients[snc][0]; + + 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]; + + for (unsigned int q_point = 0; q_point < n_quadrature_points; + ++q_point, ++shape_gradient_ptr) + { + gradients[q_point][ii][jj] += value * (*shape_gradient_ptr); } } else @@ -2080,6 +2142,52 @@ namespace FEValuesViews (make_array_view(dof_values.begin(), dof_values.end()), fe_values->finite_element_output.shape_gradients, shape_function_data, divergences); } + + + + template + template + void + Tensor<2, dim, spacedim>:: + get_function_gradients(const InputVector &fe_function, + std::vector::type> &gradients) const + { + Assert(fe_values->update_flags & update_gradients, + (typename FEValuesBase::ExcAccessToUninitializedField("update_gradients"))); + Assert(fe_values->present_cell.get() != nullptr, + ExcMessage("FEValues object is not reinit'ed to any cell")); + AssertDimension(fe_function.size(), + fe_values->present_cell->n_dofs_for_dof_handler()); + + // get function values of dofs + // on this cell + dealii::Vector dof_values(fe_values->dofs_per_cell); + fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values); + internal::do_function_gradients + (make_array_view(dof_values.begin(), dof_values.end()), + fe_values->finite_element_output.shape_gradients, shape_function_data, gradients); + } + + + + template + template + void + Tensor<2, dim, spacedim>:: + get_function_gradients_from_local_dof_values (const InputVector &dof_values, + std::vector::gradient_type> &gradients) const + { + Assert(fe_values->update_flags & update_gradients, + (typename FEValuesBase::ExcAccessToUninitializedField("update_gradients"))); + Assert(fe_values->present_cell.get() != nullptr, + ExcMessage("FEValues object is not reinit'ed to any cell")); + AssertDimension (dof_values.size(), fe_values->dofs_per_cell); + + internal::do_function_gradients + (make_array_view(dof_values.begin(), dof_values.end()), + fe_values->finite_element_output.shape_gradients, shape_function_data, gradients); + } + } diff --git a/source/fe/fe_values.inst.in b/source/fe/fe_values.inst.in index 7a78ca1539..de759298b6 100644 --- a/source/fe/fe_values.inst.in +++ b/source/fe/fe_values.inst.in @@ -135,6 +135,10 @@ for (VEC : VECTOR_TYPES; deal_II_dimension : DIMENSIONS; deal_II_space_dimension void FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension> ::get_function_divergences (const dealii::VEC&, std::vector >::type>&) const; + template + void FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension> + ::get_function_gradients + (const dealii::VEC&, std::vector >::type>&) const; #endif } @@ -235,6 +239,11 @@ for (VEC : GENERAL_CONTAINER_TYPES; void FEValuesViews::Tensor<2,deal_II_dimension, deal_II_space_dimension> ::get_function_divergences_from_local_dof_values> (const VEC&, std::vector::divergence_type>&) const; + + template + void FEValuesViews::Tensor<2,deal_II_dimension, deal_II_space_dimension> + ::get_function_gradients_from_local_dof_values> + (const VEC&, std::vector::gradient_type>&) const; #endif } @@ -396,6 +405,9 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : DIMENSIONS) template void FEValuesViews::Tensor<2,deal_II_dimension,deal_II_space_dimension>::get_function_divergences (const dealii::IndexSet&, std::vector >::type>&) const; + template + void FEValuesViews::Tensor<2,deal_II_dimension,deal_II_space_dimension>::get_function_gradients + (const dealii::IndexSet&, std::vector >::type>&) const; #endif }