From: Jean-Paul Pelteret Date: Mon, 22 Apr 2019 07:59:18 +0000 (+0200) Subject: Add functions to perform SD tensor differentation X-Git-Tag: v9.1.0-rc1~183^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=25d10474931dadd4809720aa49a4aae7554c0246;p=dealii.git Add functions to perform SD tensor differentation --- diff --git a/include/deal.II/differentiation/sd/symengine_tensor_operations.h b/include/deal.II/differentiation/sd/symengine_tensor_operations.h index d8ba720180..b14f120eaf 100644 --- a/include/deal.II/differentiation/sd/symengine_tensor_operations.h +++ b/include/deal.II/differentiation/sd/symengine_tensor_operations.h @@ -185,9 +185,574 @@ namespace Differentiation //@} + /** + * @name Symbolic differentiation + */ + //@{ + + /** + * Return the symbolic result of computing the partial derivative of the + * scalar @p f with respect to the tensor @p T. + * + * @param[in] f A scalar symbolic function or (dependent) expression. + * @param[in] T A tensor of symbolic (independent) variables. + * @return The tensor of symbolic functions or expressions representing + * the result $\frac{\partial f}{\partial \mathbf{T}}$. + */ + template + Tensor + differentiate(const Expression &f, const Tensor &T); + + /** + * Return the symbolic result of computing the partial derivative of the + * scalar @p f with respect to the symmetric tensor @p S. + * + * @param[in] f A scalar symbolic function or (dependent) expression. + * @param[in] S A symmetric tensor of symbolic (independent) variables. + * @return The symmetric tensor of symbolic functions or expressions representing + * the result $\frac{\partial f}{\partial \mathbf{S}}$. + */ + template + SymmetricTensor + differentiate(const Expression & f, + const SymmetricTensor &S); + + /** + * Return the symbolic result of computing the partial derivative of the + * rank-0 tensor (or scalar) @p f with respect to the tensor @p T. + * + * @param[in] f A rank-0 tensor symbolic function or (dependent) expression. + * @param[in] T A tensor of symbolic (independent) variables. + * @return The tensor of symbolic functions or expressions representing + * the result $\frac{\partial f}{\partial \mathbf{T}}$. + */ + template + Tensor + differentiate(const Tensor<0, dim, Expression> & f, + const Tensor &T); + + /** + * Return the symbolic result of computing the partial derivative of the + * rank-0 tensor (or scalar) @p f with respect to the symmetric tensor @p S. + * + * @param[in] f A rank-0 tensor symbolic function or (dependent) expression. + * @param[in] S A symmetric tensor of symbolic (independent) variables. + * @return The symmetric tensor of symbolic functions or expressions representing + * the result $\frac{\partial f}{\partial \mathbf{S}}$. + */ + template + SymmetricTensor + differentiate(const Tensor<0, dim, Expression> & f, + const SymmetricTensor &S); + + /** + * Return the symbolic result of computing the partial derivative of the + * tensor @p T with respect to the scalar @p x. + * + * @param[in] T A tensor of symbolic functions or (dependent) expressions. + * @param[in] x A scalar symbolic (independent) variable. + * @return The tensor of symbolic functions or expressions representing + * the result $\frac{\partial \mathbf{T}}{\partial x}$. + */ + template + Tensor + differentiate(const Tensor &T, const Expression &x); + + /** + * Return the symbolic result of computing the partial derivative of the + * symmetric tensor @p S with respect to the scalar @p x. + * + * @param[in] S A symmetric tensor of symbolic functions or (dependent) + * expressions. + * @param[in] x A scalar symbolic (independent) variable. + * @return The symmetric tensor of symbolic functions or expressions representing + * the result $\frac{\partial \mathbf{S}}{\partial x}$. + */ + template + SymmetricTensor + differentiate(const SymmetricTensor &S, + const Expression & x); + + /** + * Return the symbolic result of computing the partial derivative of the + * tensor @p T with respect to the rank-0 tensor @p x. + * + * @param[in] T A tensor of symbolic functions or (dependent) expressions. + * @param[in] x A rank-0 tensor symbolic symbolic (independent) variable. + * @return The tensor of symbolic functions or expressions representing + * the result $\frac{\partial \mathbf{T}}{\partial x}$. + */ + template + Tensor + differentiate(const Tensor &T, + const Tensor<0, dim, Expression> & x); + + /** + * Return the symbolic result of computing the partial derivative of the + * symmetric tensor @p S with respect to the rank-0 tensor @p x. + * + * @param[in] S A symmetric tensor of symbolic functions or (dependent) + * expressions. + * @param[in] x A rank-0 tensor symbolic symbolic (independent) variable. + * @return The symmetric tensor of symbolic functions or expressions representing + * the result $\frac{\partial \mathbf{S}}{\partial x}$. + */ + template + SymmetricTensor + differentiate(const SymmetricTensor &S, + const Tensor<0, dim, Expression> & x); + + /** + * Return the symbolic result of computing the partial derivative of the + * tensor @p T1 with respect to the tensor @p T2. + * + * @param[in] T1 A tensor of symbolic functions or (dependent) expressions. + * @param[in] T2 A tensor symbolic symbolic (independent) variables. + * @return The tensor of symbolic functions or variables representing + * the result $\frac{\partial \mathbf{T}_{1}}{\partial + * \mathbf{T}_{2}}$. + */ + template + Tensor + differentiate(const Tensor &T1, + const Tensor &T2); + + /** + * Return the symbolic result of computing the partial derivative of the + * symmetric tensor @p S1 with respect to the symmetric tensor @p S2. + * + * @param[in] S1 A symmetric tensor of symbolic functions or (dependent) + * expressions. + * @param[in] S2 A symmetric tensor symbolic symbolic (independent) + * variables. + * @return The symmetric tensor of symbolic functions or variables representing + * the result $\frac{\partial \mathbf{S}_{1}}{\partial + * \mathbf{S}_{2}}$. + */ + template + SymmetricTensor + differentiate(const SymmetricTensor &S1, + const SymmetricTensor &S2); + + /** + * Return the symbolic result of computing the partial derivative of the + * tensor @p T with respect to the symmetric tensor @p S. + * + * @param[in] T A tensor of symbolic functions or (dependent) expressions. + * @param[in] S A symmetric tensor symbolic symbolic (independent) + * variables. + * @return The tensor of symbolic functions or variables representing + * the result $\frac{\partial \mathbf{T}}{\partial \mathbf{S}}$. + */ + template + Tensor + differentiate(const Tensor & T, + const SymmetricTensor &S); + + /** + * Return the symbolic result of computing the partial derivative of the + * symmetric tensor @p S with respect to the tensor @p T. + * + * @param[in] S A symmetric tensor of symbolic functions or (dependent) + * expressions. + * @param[in] T A tensor symbolic symbolic (independent) variables. + * @return The tensor of symbolic functions or variables representing + * the result $\frac{\partial \mathbf{S}}{\partial \mathbf{T}}$. + */ + template + Tensor + differentiate(const SymmetricTensor &S, + const Tensor & T); + + //@} + + } // namespace SD +} // namespace Differentiation + + +/* -------------------- inline and template functions ------------------ */ + + +# ifndef DOXYGEN + +namespace Differentiation +{ + namespace SD + { + // --- Symbolic differentiation --- + + + namespace internal + { + template + TableIndices + concatenate_indices(const TableIndices &indices_1, + const TableIndices &indices_2) + { + TableIndices indices_out; + for (unsigned int i = 0; i < rank_1; ++i) + indices_out[i] = indices_1[i]; + for (unsigned int j = 0; j < rank_2; ++j) + indices_out[rank_1 + j] = indices_2[j]; + return indices_out; + } + + + template + TableIndices + transpose_indices(const TableIndices &indices) + { + return indices; + } + + + template <> + inline TableIndices<2> + transpose_indices(const TableIndices<2> &indices) + { + return TableIndices<2>(indices[1], indices[0]); + } + + + template + bool + is_symmetric_component(const TableIndices &, + const Tensor &) + { + return false; + } + + + template + bool + is_symmetric_component(const TableIndices &, + const SymmetricTensor &) + { + static_assert( + rank == 0 || rank == 2, + "Querying symmetric component for non rank-2 symmetric tensor index is not allowed."); + return false; + } + + + template + bool + is_symmetric_component(const TableIndices<2> &table_indices, + const SymmetricTensor<2, dim, ValueType> &) + { + return table_indices[0] != table_indices[1]; + } + + + template class TensorType> + TensorType<0, dim, ValueType> + scalar_diff_tensor(const ValueType & func, + const TensorType<0, dim, ValueType> &op) + { + return differentiate(func, op); + } + + + template class TensorType> + TensorType + scalar_diff_tensor(const ValueType & func, + const TensorType &op) + { + TensorType out; + for (unsigned int i = 0; i < out.n_independent_components; ++i) + { + const TableIndices indices( + out.unrolled_to_component_indices(i)); + out[indices] = differentiate(func, op[indices]); + + if (is_symmetric_component(indices, op)) + out[indices] *= 0.5; + } + return out; + } + + + // Specialization for rank-0 tensor + template class TensorType> + TensorType + tensor_diff_tensor(const TensorType<0, dim, ValueType> & func, + const TensorType &op) + { + TensorType out; + for (unsigned int i = 0; i < out.n_independent_components; ++i) + { + const TableIndices indices( + out.unrolled_to_component_indices(i)); + out[indices] = differentiate(func, op[indices]); + + if (is_symmetric_component(indices, op)) + out[indices] *= 0.5; + } + return out; + } + + + template class TensorType> + TensorType + tensor_diff_scalar(const TensorType &funcs, + const ValueType & op) + { + TensorType out; + for (unsigned int i = 0; i < out.n_independent_components; ++i) + { + const TableIndices indices( + out.unrolled_to_component_indices(i)); + out[indices] = differentiate(funcs[indices], op); + } + return out; + } + + + // Specialization for rank-0 tensor + template class TensorType> + TensorType + tensor_diff_tensor(const TensorType &funcs, + const TensorType<0, dim, ValueType> & op) + { + TensorType out; + for (unsigned int i = 0; i < out.n_independent_components; ++i) + { + const TableIndices indices( + out.unrolled_to_component_indices(i)); + out[indices] = differentiate(funcs[indices], op); + } + return out; + } + + + // For either symmetric or normal tensors + template class TensorType> + TensorType + tensor_diff_tensor(const TensorType &funcs, + const TensorType &op) + { + TensorType out; + for (unsigned int i = 0; i < funcs.n_independent_components; ++i) + { + const TableIndices indices_i( + funcs.unrolled_to_component_indices(i)); + for (unsigned int j = 0; j < op.n_independent_components; ++j) + { + const TableIndices indices_j( + op.unrolled_to_component_indices(j)); + const TableIndices indices_out = + concatenate_indices(indices_i, indices_j); + + out[indices_out] = + differentiate(funcs[indices_i], op[indices_j]); + + if (is_symmetric_component(indices_j, op)) + out[indices_out] *= 0.5; + } + } + return out; + } + + + // For mixed symmetric/standard tensors + // The return type is always a standard tensor, since we cannot be sure + // that any symmetries exist in either the function tensor or the + // differential operator. + template class TensorType_1, + template class TensorType_2> + Tensor + tensor_diff_tensor(const TensorType_1 &funcs, + const TensorType_2 &op) + { + Tensor out; + for (unsigned int i = 0; i < funcs.n_independent_components; ++i) + { + const TableIndices indices_i( + funcs.unrolled_to_component_indices(i)); + for (unsigned int j = 0; j < op.n_independent_components; ++j) + { + const TableIndices indices_j( + op.unrolled_to_component_indices(j)); + const TableIndices indices_out = + concatenate_indices(indices_i, indices_j); + + out[indices_out] = + differentiate(funcs[indices_i], op[indices_j]); + + if (is_symmetric_component(indices_j, op)) + out[indices_out] *= 0.5; + + // TODO: Implement for SymmetricTensor<4,dim,...> + if (std::is_same, + SymmetricTensor<2, dim, ValueType>>:: + value) // Symmetric function + { + const TableIndices indices_out_t = + concatenate_indices(transpose_indices(indices_i), + indices_j); + out[indices_out_t] = out[indices_out]; + } + else if (std::is_same, + SymmetricTensor<2, dim, ValueType>>:: + value) // Symmetric operator + { + const TableIndices indices_out_t = + concatenate_indices(indices_i, + transpose_indices(indices_j)); + out[indices_out_t] = out[indices_out]; + } + else + { + Assert( + false, + ExcMessage( + "Expect mixed tensor differentiation to have at least " + "one component stemming from a symmetric tensor.")); + } + } + } + return out; + } + + } // namespace internal + + + template + Tensor + differentiate(const Expression & func, + const Tensor &op) + { + return internal::scalar_diff_tensor(func, op); + } + + + template + SymmetricTensor + differentiate(const Expression & func, + const SymmetricTensor &op) + { + return internal::scalar_diff_tensor(func, op); + } + + + template + Tensor + differentiate(const Tensor<0, dim, Expression> & func, + const Tensor &op) + { + return internal::scalar_diff_tensor(func, op); + } + + + template + SymmetricTensor + differentiate(const Tensor<0, dim, Expression> & func, + const SymmetricTensor &op) + { + // Ensure that this returns a symmetric tensor by + // invoking the scalar value function + const Expression tmp = func; + return internal::scalar_diff_tensor(tmp, op); + } + + + template + Tensor + differentiate(const Tensor &symbol_tensor, + const Expression & op) + { + return internal::tensor_diff_scalar(symbol_tensor, op); + } + + + template + Tensor + differentiate(const Tensor &symbol_tensor, + const Tensor<0, dim, Expression> & op) + { + return internal::tensor_diff_scalar(symbol_tensor, op); + } + + + template + SymmetricTensor + differentiate(const SymmetricTensor &symbol_tensor, + const Expression & op) + { + return internal::tensor_diff_scalar(symbol_tensor, op); + } + + + template + SymmetricTensor + differentiate(const SymmetricTensor &symbol_tensor, + const Tensor<0, dim, Expression> & op) + { + return internal::tensor_diff_scalar(symbol_tensor, op); + } + + + template + Tensor + differentiate(const Tensor &symbol_tensor, + const Tensor &op) + { + return internal::tensor_diff_tensor(symbol_tensor, op); + } + + + template + SymmetricTensor + differentiate(const SymmetricTensor &symbol_tensor, + const SymmetricTensor &op) + { + return internal::tensor_diff_tensor(symbol_tensor, op); + } + + + template + Tensor + differentiate(const Tensor & symbol_tensor, + const SymmetricTensor &op) + { + return internal::tensor_diff_tensor(symbol_tensor, op); + } + + + template + Tensor + differentiate(const SymmetricTensor &symbol_tensor, + const Tensor & op) + { + return internal::tensor_diff_tensor(symbol_tensor, op); + } + } // namespace SD } // namespace Differentiation +# endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE