From 9bac38c6ecafe3da017b9b49b013bd603866dc94 Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Mon, 11 Feb 2019 18:46:21 +0100 Subject: [PATCH] AD Helpers: Add some internal functions to assist when using Extractors --- .../deal.II/differentiation/ad/ad_helpers.h | 970 ++++++++++++++++++ 1 file changed, 970 insertions(+) diff --git a/include/deal.II/differentiation/ad/ad_helpers.h b/include/deal.II/differentiation/ad/ad_helpers.h index bf3ea90c62..5460362a72 100644 --- a/include/deal.II/differentiation/ad/ad_helpers.h +++ b/include/deal.II/differentiation/ad/ad_helpers.h @@ -21,6 +21,8 @@ #if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO) # include +# include +# include # include # include @@ -29,6 +31,8 @@ # include # include +# include + # include # include @@ -1647,6 +1651,972 @@ namespace Differentiation }; // class ADHelperResidualLinearization + + namespace internal + { + /** + * A helper struct that assists with the extraction of data associated + * with fields that are defined by FEExtractors. + */ + template + struct Extractor; + + + /** + * A helper struct that assists with the extraction of data associated + * with fields that are defined by FEExtractors. + * This particular specialization is for scalar fields. + */ + template + struct Extractor + { + /** + * The number of components of the field. + */ + static const unsigned int n_components = 1; + + /** + * The tensor rank of the field. + */ + static const unsigned int rank = 0; + + /** + * The tensor type associated with this field. + */ + template + using tensor_type = Tensor; + + static_assert( + n_components == tensor_type::n_independent_components, + "The number of components doesn't match that of the corresponding tensor type."); + static_assert( + rank == tensor_type::rank, + "The rank doesn't match that of the corresponding tensor type."); + + /** + * The value type associated with this field. + */ + // Note: FEValuesViews::Scalar::tensor_type is double, so we can't + // use it (FEValuesViews) in this context. + // In fact, sadly, all FEValuesViews objects expect doubles as value + // types. + template + using value_type = NumberType; + + /** + * The gradient type associated with this field. + */ + template + using gradient_type = Tensor; // NumberType; + + /** + * Return the first global component of this field. + */ + static inline unsigned int + first_component(const FEValuesExtractors::Scalar &extractor) + { + return extractor.component; + } + + /** + * Return a flag that indicates if the input @p unrolled_index + * corresponds to a symmetric component of the field. + * + * For a scalar field, the single component is defined + * to not have a symmetric counterpart. + */ + static bool + symmetric_component(const unsigned int unrolled_index) + { + (void)unrolled_index; + return false; + } + + /** + * Return the local unrolled component corresponding to + * @p column_offset entry of the @p table_indices. + * + * For a scalar field, the local component is always + * equal to zero. + */ + template + static IndexType + local_component(const TableIndices &table_indices, + const unsigned int column_offset) + { + Assert(column_offset <= rank_in, ExcInternalError()); + (void)table_indices; + (void)column_offset; + return 0; + } + }; + + + /** + * A helper struct that assists with the extraction of data associated + * with fields that are defined by FEExtractors. + * This particular specialization is for vector fields. + */ + template + struct Extractor + { + /** + * The number of components of the field. + */ + static const unsigned int n_components = dim; + + /** + * The tensor rank of the field. + */ + static const unsigned int rank = 1; + + /** + * The tensor type associated with this field. + */ + template + using tensor_type = Tensor; + + static_assert( + n_components == tensor_type::n_independent_components, + "The number of components doesn't match that of the corresponding tensor type."); + static_assert( + rank == tensor_type::rank, + "The rank doesn't match that of the corresponding tensor type."); + + /** + * The value type associated with this field. + */ + template + using value_type = tensor_type; + + /** + * The gradient type associated with this field. + */ + template + using gradient_type = Tensor; + + /** + * Return the first global component of this field. + */ + static inline unsigned int + first_component(const FEValuesExtractors::Vector &extractor) + { + return extractor.first_vector_component; + } + + /** + * + * Return a flag that indicates if the input @p unrolled_index + * corresponds to a symmetric component of the field. + * + * For a vector field, the none of the vector components + * have a symmetric counterpart. + */ + static bool + symmetric_component(const unsigned int unrolled_index) + { + (void)unrolled_index; + return false; + } + + /** + * Return the table index corresponding to + * @p column_offset entry of the input @p table_indices. + */ + template + static TableIndices + table_index_view(const TableIndices &table_indices, + const unsigned int column_offset) + { + Assert(0 + column_offset < rank_in, ExcInternalError()); + return TableIndices(table_indices[column_offset]); + } + + /** + * Return the local unrolled component corresponding to + * @p column_offset entry of the @p table_indices. + * + * This function computes and returns a local component + * associated with the extractor's @p tensor_type from a + * set of @p table_indices that are generally associated + * with a tensor of equal or greater order. In effect, it + * creates a view of a selected number of indices of the + * input table, and interprets that subtable's indices as + * the local index to be returned. Since the @p table_indices + * may be of size greater than the extractor's @p rank, + * the @p column_offset specifies the first index of the + * input table to create the view from. + */ + template + static IndexType + local_component(const TableIndices &table_indices, + const unsigned int column_offset) + { + static_assert( + rank_in >= rank, + "Cannot extract more table indices than the input table has!"); + using TensorType = tensor_type; + return TensorType::component_to_unrolled_index( + table_index_view(table_indices, column_offset)); + } + }; + + + /** + * A helper struct that assists with the extraction of data associated + * with fields that are defined by FEExtractors. + * This particular specialization is for rank-1 tensor fields. + */ + template + struct Extractor> + { + /** + * The number of components of the field. + */ + static const unsigned int n_components = + Tensor<1, dim>::n_independent_components; + + /** + * The tensor rank of the field. + */ + static const unsigned int rank = 1; + + /** + * The tensor type associated with this field. + */ + template + using tensor_type = Tensor; + + /** + * The value type associated with this field. + */ + template + using value_type = tensor_type; + + /** + * The gradient type associated with this field. + */ + template + using gradient_type = Tensor; + + /** + * Return the first global component of this field. + */ + static inline unsigned int + first_component(const FEValuesExtractors::Tensor<1> &extractor) + { + return extractor.first_tensor_component; + } + + /** + * Return a flag that indicates if the input @p unrolled_index + * corresponds to a symmetric component of the field. + * + * For a vector field, the none of the vector components + * have a symmetric counterpart. + */ + static bool + symmetric_component(const unsigned int unrolled_index) + { + (void)unrolled_index; + return false; + } + + /** + * Return the table index corresponding to + * @p column_offset entry of the input @p table_indices. + */ + template + static TableIndices + table_index_view(const TableIndices &table_indices, + const unsigned int column_offset) + { + Assert(column_offset < rank_in, ExcInternalError()); + return TableIndices(table_indices[column_offset]); + } + + /** + * Return the local unrolled component corresponding to + * a subset of table indices from the input @p table_indices. + * + * This function computes and returns a local component + * associated with the extractor's @p tensor_type from a + * set of @p table_indices that are generally associated + * with a tensor of equal or greater order. In effect, it + * creates a view of a selected number of indices of the + * input table, and interprets that subtable's indices as + * the local index to be returned. Since the @p table_indices + * may be of size greater than the extractor's @p rank, + * the @p column_offset specifies the first index of the + * input table to create the view from. + */ + template + static IndexType + local_component(const TableIndices &table_indices, + const unsigned int column_offset) + { + static_assert( + rank_in >= rank, + "Cannot extract more table indices than the input table has!"); + using TensorType = tensor_type; + return TensorType::component_to_unrolled_index( + table_index_view(table_indices, column_offset)); + } + }; + + + /** + * A helper struct that assists with the extraction of data associated + * with fields that are defined by FEExtractors. + * This particular specialization is for rank-2 tensor fields. + */ + template + struct Extractor> + { + /** + * The number of components of the field. + */ + static const unsigned int n_components = + Tensor<2, dim>::n_independent_components; + + /** + * The tensor rank of the field. + */ + static const unsigned int rank = Tensor<2, dim>::rank; + + /** + * The tensor type associated with this field. + */ + template + using tensor_type = Tensor; + + /** + * The value type associated with this field. + */ + template + using value_type = tensor_type; + + /** + * The gradient type associated with this field. + */ + template + using gradient_type = Tensor; + + /** + * Return the first global component of this field. + */ + static inline unsigned int + first_component(const FEValuesExtractors::Tensor<2> &extractor) + { + return extractor.first_tensor_component; + } + + /** + * Return a flag that indicates if the input @p unrolled_index + * corresponds to a symmetric component of the field. + * + * For a rank-2 tensor field, the none of the tensor + * components have a symmetric counterpart. + */ + static bool + symmetric_component(const unsigned int unrolled_index) + { + (void)unrolled_index; + return false; + } + + /** + * Return the table indices corresponding to + * @p column_offset entry of the input @p table_indices. + */ + template + static TableIndices + table_index_view(const TableIndices &table_indices, + const unsigned int column_offset) + { + Assert(column_offset < rank_in, ExcInternalError()); + Assert(column_offset + 1 < rank_in, ExcInternalError()); + return TableIndices(table_indices[column_offset], + table_indices[column_offset + 1]); + } + + /** + * Return the local unrolled component corresponding to + * @p column_offset entry of the @p table_indices. + * + * This function computes and returns a local component + * associated with the extractor's @p tensor_type from a + * set of @p table_indices that are generally associated + * with a tensor of equal or greater order. In effect, it + * creates a view of a selected number of indices of the + * input table, and interprets that subtable's indices as + * the local index to be returned. Since the @p table_indices + * may be of size greater than the extractor's @p rank, + * the @p column_offset specifies the first index of the + * input table to create the view from. + */ + template + static IndexType + local_component(const TableIndices &table_indices, + const unsigned int column_offset) + { + static_assert( + rank_in >= rank, + "Cannot extract more table indices than the input table has!"); + using TensorType = tensor_type; + return TensorType::component_to_unrolled_index( + table_index_view(table_indices, column_offset)); + } + }; + + + /** + * A helper struct that assists with the extraction of data associated + * with fields that are defined by FEExtractors. + * This particular specialization is for rank-2 symmetric tensor fields. + */ + template + struct Extractor> + { + /** + * The number of components of the field. + */ + static const unsigned int n_components = + SymmetricTensor<2, dim>::n_independent_components; + + /** + * The tensor rank of the field. + */ + static const unsigned int rank = SymmetricTensor<2, dim>::rank; + + /** + * The tensor type associated with this field. + */ + template + using tensor_type = SymmetricTensor; + + /** + * The value type associated with this field. + */ + template + using value_type = tensor_type; + + /** + * The gradient type associated with this field. + */ + template + using gradient_type = Tensor; + + /** + * Return the first global component of this field. + */ + static inline unsigned int + first_component(const FEValuesExtractors::SymmetricTensor<2> &extractor) + { + return extractor.first_tensor_component; + } + + /** + * Return a flag that indicates if the input @p unrolled_index + * corresponds to a symmetric component of the field. + * + * For a rank-2 symmetric tensor field, each of the + * off-diagonal components have a symmetric counterpart, + * while the diagonal components do not. + */ + static bool + symmetric_component(const unsigned int unrolled_index) + { + const TableIndices<2> table_indices = + tensor_type::unrolled_to_component_indices(unrolled_index); + return table_indices[0] != table_indices[1]; + } + + /** + * Return the table indices corresponding to + * @p column_offset entry of the input @p table_indices. + */ + template + static TableIndices + table_index_view(const TableIndices &table_indices, + const unsigned int column_offset) + { + Assert(column_offset < rank_in, ExcInternalError()); + Assert(column_offset + 1 < rank_in, ExcInternalError()); + return TableIndices(table_indices[column_offset], + table_indices[column_offset + 1]); + } + + /** + * Return the local unrolled component corresponding to + * @p column_offset entry of the @p table_indices. + * + * This function computes and returns a local component + * associated with the extractor's @p tensor_type from a + * set of @p table_indices that are generally associated + * with a tensor of equal or greater order. In effect, it + * creates a view of a selected number of indices of the + * input table, and interprets that subtable's indices as + * the local index to be returned. Since the @p table_indices + * may be of size greater than the extractor's @p rank, + * the @p column_offset specifies the first index of the + * input table to create the view from. + */ + template + static IndexType + local_component(const TableIndices &table_indices, + const unsigned int column_offset) + { + static_assert( + rank_in >= rank, + "Cannot extract more table indices than the input table has!"); + using TensorType = tensor_type; + return TensorType::component_to_unrolled_index( + table_index_view(table_indices, column_offset)); + } + }; + + + /** + * A helper struct that defines the return type of gradient (first + * derivative) calculations of scalar fields with respect to a field + * defined by the @p ExtractorType template parameter. + */ + template + struct ScalarFieldGradient + { + /** + * The type associated with computing the gradient of a scalar + * field with respect to the given @p ExtractorType. + */ + using type = + typename Extractor::template tensor_type; + }; + + + /** + * An intermediate helper struct that defines the return type of Hessian + * (second derivative) calculations of scalar fields with respect to + * fields defined by the two extractor-type template parameters. + * The first, @p ExtractorType_Row, defines the field that the first + * derivatives are taken with respect to while the second, + * @p ExtractorType_Col, defines the field that the second derivatives + * are taken with respect to. + */ + template + struct HessianType + { + /** + * The type associated with computing the gradient of a scalar + * field with respect to the given @p ExtractorType_Row + * followed by the @p ExtractorType_Col. + * + * @note We set the return type for + * HessianType + * as a normal Tensor. This is because if one has two vector components, + * the coupling tensor (i.e. Hessian component) is in + * general not symmetric. + */ + template + using type = Tensor; + }; + + + /** + * An intermediate helper struct that defines the return type of Hessian + * (second derivative) calculations of scalar fields with respect to + * fields defined by the two extractor-type template parameters. This + * particular specialization is for taking the first derivative with + * respect to a symmetric tensor field, and the second with respect to a + * scalar field. + */ + template <> + struct HessianType, + FEValuesExtractors::Scalar> + { + /** + * The type associated with computing the gradient of a scalar + * field with respect to the given + * ExtractorType_Row = + * FEValuesExtractors::SymmetricTensor<2> followed by the + * ExtractorType_Col = FEValuesExtractors::Scalar. + */ + template + using type = SymmetricTensor<2 /*rank*/, dim, NumberType>; + }; + + + /** + * An intermediate helper struct that defines the return type of Hessian + * (second derivative) calculations of scalar fields with respect to + * fields defined by the two extractor-type template parameters. This + * particular specialization is for taking the first derivative with + * respect to a scalar field, and the second with respect to a symmetric + * tensor field. + */ + template <> + struct HessianType> + { + /** + * The type associated with computing the gradient of a scalar + * field with respect to the given + * ExtractorType_Row = FEValuesExtractors::Scalar + * followed by the + * ExtractorType_Col = + * FEValuesExtractors::SymmetricTensor<2>. + */ + template + using type = SymmetricTensor<2 /*rank*/, dim, NumberType>; + }; + + + /** + * An intermediate helper struct that defines the return type of Hessian + * (second derivative) calculations of scalar fields with respect to + * fields defined by the two extractor-type template parameters. This + * particular specialization is for taking both the first and second + * derivatives with respect to symmetric tensor fields. + */ + template <> + struct HessianType, + FEValuesExtractors::SymmetricTensor<2>> + { + /** + * The type associated with computing the gradient of a scalar + * field with respect to the given + * ExtractorType_Row = + * FEValuesExtractors::SymmetricTensor<2> followed by the + * ExtractorType_Col = + * FEValuesExtractors::SymmetricTensor<2>. + */ + template + using type = SymmetricTensor<4 /*rank*/, dim, NumberType>; + }; + + + /** + * A helper struct that defines the final return type of Hessian + * (second derivative) calculations of scalar fields with respect to + * fields defined by the two extractor-type + * template parameters. The first, @p ExtractorType_Row, defines the field + * that the first derivatives are taken with respect to while the second, + * @p ExtractorType_Col, defines the field that the second derivatives + * are taken with respect to. + */ + template + struct ScalarFieldHessian + { + /** + * The tensor rank of the resulting derivative computation. + */ + static const int rank = Extractor::rank + + Extractor::rank; + + /** + * The type associated with computing the Hessian of a scalar + * field with first respect to the field defined by the + * @p ExtractorType_Row and then with respect to the field defined by + * the @p ExtractorType_Col. + */ + using type = + typename HessianType:: + template type; + }; + + + /** + * A helper struct that defines the return type of value computations + * of vector fields the @p ExtractorType_Field template parameter. + */ + template + using VectorFieldValue = + ScalarFieldGradient; + + + /** + * A helper struct that defines the final return type of Jacobian + * (first derivative) calculations of vector fields with respect to + * fields as defined by the two extractor-type template parameters. + * The first, @p ExtractorType_Field, defines the field from which + * the initial field values are computed while the second, + * @p ExtractorType_Derivative, defines the field that the derivatives + * are taken with respect to. + */ + template + using VectorFieldJacobian = ScalarFieldHessian; + + + /** + * Return a global view of the field component indices that correspond to + * the input @p extractor. For this general function the + * @p ignore_symmetries flag has no effect. + */ + template + std::vector + extract_field_component_indices(const ExtractorType &extractor, + const bool ignore_symmetries = true) + { + (void)ignore_symmetries; + const IndexType n_components = + internal::Extractor::n_components; + const IndexType comp_first = + internal::Extractor::first_component(extractor); + std::vector indices(n_components); + std::iota(indices.begin(), indices.end(), comp_first); + return indices; + } + + + /** + * Return a global view of the field component indices that correspond to + * the input FEValuesExtractors::SymmetricTensor @p extractor_symm_tensor. + * If the @p ignore_symmetries is set true, then all + * component of the tensor are considered to be independent. If set to + * code>false, then the set of returned indices will contain + * duplicate entries for components that are symmetric. + */ + template + std::vector + extract_field_component_indices( + const FEValuesExtractors::SymmetricTensor<2> &extractor_symm_tensor, + const bool ignore_symmetries = true) + { + using ExtractorType = FEValuesExtractors::SymmetricTensor<2>; + const IndexType n_components = + internal::Extractor::n_components; + if (ignore_symmetries == true) + { + const IndexType comp_first = + internal::Extractor::first_component( + extractor_symm_tensor); + std::vector indices(n_components); + std::iota(indices.begin(), indices.end(), comp_first); + return indices; + } + else + { + // First get all of the indices of the non-symmetric tensor + const FEValuesExtractors::Tensor<2> extractor_tensor( + extractor_symm_tensor.first_tensor_component); + std::vector indices = + extract_field_component_indices(extractor_tensor, true); + + // Then we overwrite any illegal entries with the equivalent indices + // from the symmetric tensor + for (unsigned int i = 0; i < indices.size(); ++i) + { + if (indices[i] >= n_components) + { + const TableIndices<2> ti_tensor = + Tensor<2, dim>::unrolled_to_component_indices(indices[i]); + const IndexType sti_new_index = + SymmetricTensor<2, dim>::component_to_unrolled_index( + ti_tensor); + indices[i] = sti_new_index; + } + } + + return indices; + } + } + + + /** + * Set the unrolled component given by @p index in the generic tensor + * @p t to the given @p value. + */ + template + inline void + set_tensor_entry(TensorType & t, + const unsigned int &unrolled_index, + const NumberType & value) + { + // Where possible, set values using TableIndices + Assert(unrolled_index < t.n_independent_components, + ExcIndexRange(unrolled_index, 0, t.n_independent_components)); + t[TensorType::unrolled_to_component_indices(unrolled_index)] = value; + } + + + /** + * Set the unrolled component given by @p index in the rank-0 tensor + * @p t to the given @p value. + */ + template + inline void set_tensor_entry(Tensor<0, dim, NumberType> &t, + const unsigned int & unrolled_index, + const NumberType & value) + { + Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1)); + (void)unrolled_index; + t = value; + } + + + /** + * Set the value of @p t to the given @p value. + * This function exists to provide compatibility with similar functions + * that exist for use with the tensor classes. + */ + template + inline void + set_tensor_entry(NumberType & t, + const unsigned int &unrolled_index, + const NumberType & value) + { + Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1)); + (void)unrolled_index; + t = value; + } + + + /** + * Set the unrolled component given by the @p index_row and + * the @p index_col in the fourth-order symmetric tensor + * @p t to the given @p value. + */ + template + inline void set_tensor_entry(SymmetricTensor<4, dim, NumberType> &t, + const unsigned int &unrolled_index_row, + const unsigned int &unrolled_index_col, + const NumberType & value) + { + // Fourth order symmetric tensors require a specialized interface + // to extract values. + using SubTensorType = SymmetricTensor<2, dim, NumberType>; + Assert(unrolled_index_row < SubTensorType::n_independent_components, + ExcIndexRange(unrolled_index_row, + 0, + SubTensorType::n_independent_components)); + Assert(unrolled_index_col < SubTensorType::n_independent_components, + ExcIndexRange(unrolled_index_col, + 0, + SubTensorType::n_independent_components)); + const TableIndices<2> indices_row = + SubTensorType::unrolled_to_component_indices(unrolled_index_row); + const TableIndices<2> indices_col = + SubTensorType::unrolled_to_component_indices(unrolled_index_col); + t[indices_row[0]][indices_row[1]][indices_col[0]][indices_col[1]] = + value; + } + + + /** + * Return the value of the @p index'th unrolled component of the + * generic tensor @p t. + */ + template class TensorType> + inline NumberType + get_tensor_entry(const TensorType &t, + const unsigned int & unrolled_index) + { + // Where possible, get values using TableIndices + Assert(unrolled_index < t.n_independent_components, + ExcIndexRange(unrolled_index, 0, t.n_independent_components)); + return t[TensorType:: + unrolled_to_component_indices(unrolled_index)]; + } + + + /** + * Return the value of the @p index'th unrolled component of the + * rank-0 tensor @p t. + */ + template class TensorType> + inline NumberType + get_tensor_entry(const TensorType<0, dim, NumberType> &t, + const unsigned int & unrolled_index) + { + Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1)); + (void)unrolled_index; + return t; + } + + + /** + * Return the value of @p t. + * This function exists to provide compatibility with similar functions + * that exist for use with the tensor classes. + */ + template + inline const NumberType & + get_tensor_entry(const NumberType &t, const unsigned int &unrolled_index) + { + Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1)); + (void)unrolled_index; + return t; + } + + + /** + * Return a reference to the entry stored in the @p index'th unrolled + * component of the generic tensor @p t. + */ + template class TensorType> + inline NumberType & + get_tensor_entry(TensorType &t, + const unsigned int & unrolled_index) + { + // Where possible, get values using TableIndices + Assert(unrolled_index < t.n_independent_components, + ExcIndexRange(unrolled_index, 0, t.n_independent_components)); + return t[TensorType:: + unrolled_to_component_indices(unrolled_index)]; + } + + + /** + * Return a reference to the entry stored in the @p index'th unrolled + * component of the rank-0 tensor @p t. + */ + template class TensorType> + NumberType &get_tensor_entry(TensorType<0, dim, NumberType> &t, + const unsigned int & index) + { + Assert(index == 0, ExcIndexRange(index, 0, 1)); + (void)index; + return t; + } + + + /** + * Return a reference to @p t. + * This function exists to provide compatibility with similar functions + * that exist for use with the tensor classes. + */ + template + inline NumberType & + get_tensor_entry(NumberType &t, const unsigned int &index) + { + Assert(index == 0, ExcIndexRange(index, 0, 1)); + (void)index; + return t; + } + + } // namespace internal } // namespace AD } // namespace Differentiation -- 2.39.5