* @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 <int dim, int spacedim>
class Tensor<2,dim,spacedim>
*/
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.
* divergences of the view the Tensor class.
*/
typedef typename ProductType<Number, typename Tensor<2,dim,spacedim>::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<Number, typename Tensor<2,dim,spacedim>::gradient_type>::type gradient_type;
};
/**
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 <tt>fe_function</tt> at the
*/
template <class InputVector>
void get_function_divergences_from_local_dof_values (const InputVector &dof_values,
- std::vector<typename OutputType<typename InputVector::value_type>::divergence_type> &values) const;
+ std::vector<typename OutputType<typename InputVector::value_type>::divergence_type> &divergences) const;
+ /**
+ * Return the gradient of the selected vector components of the finite
+ * element function characterized by <tt>fe_function</tt> at the
+ * quadrature points of the cell, face or subface selected the last time
+ * the <tt>reinit</tt> 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 <class InputVector>
+ void get_function_gradients (const InputVector &fe_function,
+ std::vector<typename ProductType<gradient_type,typename InputVector::value_type>::type> &gradients) const;
+
+ /**
+ * @copydoc FEValuesViews::Tensor<2,dim,spacedim>::get_function_values_from_local_dof_values()
+ */
+ template <class InputVector>
+ void get_function_gradients_from_local_dof_values (const InputVector &dof_values,
+ std::vector<typename OutputType<typename InputVector::value_type>::gradient_type> &gradients) const;
private:
/**
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;
}
return return_value;
}
}
+
+ template <int dim, int spacedim>
+ 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<dim,spacedim>::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;
+ }
+ }
+
}
/**
- * 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 <code>(dim*dim)</code> components of a vector-valued
* element. The value of <code>dim</code> is defined by the FEValues object
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 <int dim, int spacedim, typename Number>
+ void
+ do_function_gradients (const ArrayView<Number> &dof_values,
+ const Table<2,dealii::Tensor<1,spacedim> > &shape_gradients,
+ const std::vector<typename Tensor<2,dim,spacedim>::ShapeFunctionData> &shape_function_data,
+ std::vector<typename Tensor<2,dim,spacedim>::template OutputType<Number>::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<Number>::gradient_type());
+
+ for (unsigned int shape_function=0;
+ shape_function<dofs_per_cell; ++shape_function)
+ {
+ const int snc = shape_function_data[shape_function].single_nonzero_component;
+
+ if (snc == -2)
+ // shape function is zero for the selected components
+ continue;
+
+ const Number &value = dof_values[shape_function];
+ // For auto-differentiable numbers, the fact that a DoF value is zero
+ // does not imply that its derivatives are zero as well. So we
+ // can't filter by value for these number types.
+ if (!Differentiation::AD::is_ad_number<Number>::value)
+ if (value == dealii::internal::NumberType<Number>::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
(make_array_view(dof_values.begin(), dof_values.end()),
fe_values->finite_element_output.shape_gradients, shape_function_data, divergences);
}
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Tensor<2, dim, spacedim>::
+ get_function_gradients(const InputVector &fe_function,
+ std::vector<typename ProductType<gradient_type,typename InputVector::value_type>::type> &gradients) const
+ {
+ Assert(fe_values->update_flags & update_gradients,
+ (typename FEValuesBase<dim,spacedim>::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<typename InputVector::value_type> dof_values(fe_values->dofs_per_cell);
+ fe_values->present_cell->get_interpolated_dof_values(fe_function, dof_values);
+ internal::do_function_gradients<dim,spacedim>
+ (make_array_view(dof_values.begin(), dof_values.end()),
+ fe_values->finite_element_output.shape_gradients, shape_function_data, gradients);
+ }
+
+
+
+ template <int dim, int spacedim>
+ template <class InputVector>
+ void
+ Tensor<2, dim, spacedim>::
+ get_function_gradients_from_local_dof_values (const InputVector &dof_values,
+ std::vector<typename OutputType<typename InputVector::value_type>::gradient_type> &gradients) const
+ {
+ Assert(fe_values->update_flags & update_gradients,
+ (typename FEValuesBase<dim,spacedim>::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<dim,spacedim>
+ (make_array_view(dof_values.begin(), dof_values.end()),
+ fe_values->finite_element_output.shape_gradients, shape_function_data, gradients);
+ }
+
}
void FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension>
::get_function_divergences<dealii::VEC>
(const dealii::VEC&, std::vector<ProductType<dealii::VEC::value_type,dealii::Tensor<1,deal_II_space_dimension> >::type>&) const;
+ template
+ void FEValuesViews::Tensor<2, deal_II_dimension, deal_II_space_dimension>
+ ::get_function_gradients<dealii::VEC>
+ (const dealii::VEC&, std::vector<ProductType<dealii::VEC::value_type,dealii::Tensor<3,deal_II_space_dimension> >::type>&) const;
#endif
}
void FEValuesViews::Tensor<2,deal_II_dimension, deal_II_space_dimension>
::get_function_divergences_from_local_dof_values<VEC<Number>>
(const VEC<Number>&, std::vector<typename OutputType<Number>::divergence_type>&) const;
+
+ template
+ void FEValuesViews::Tensor<2,deal_II_dimension, deal_II_space_dimension>
+ ::get_function_gradients_from_local_dof_values<VEC<Number>>
+ (const VEC<Number>&, std::vector<typename OutputType<Number>::gradient_type>&) const;
#endif
}
template
void FEValuesViews::Tensor<2,deal_II_dimension,deal_II_space_dimension>::get_function_divergences<dealii::IndexSet>
(const dealii::IndexSet&, std::vector<ProductType<IndexSet::value_type,dealii::Tensor<1,deal_II_space_dimension> >::type>&) const;
+ template
+ void FEValuesViews::Tensor<2,deal_II_dimension,deal_II_space_dimension>::get_function_gradients<dealii::IndexSet>
+ (const dealii::IndexSet&, std::vector<ProductType<IndexSet::value_type,dealii::Tensor<3,deal_II_space_dimension> >::type>&) const;
#endif
}
projects / dealii.git / commitdiff