From: Jean-Paul Pelteret Date: Fri, 27 Apr 2018 00:47:02 +0000 (+0200) Subject: Adol-C advanced branching: Fix SymmetricTensor X-Git-Tag: v9.0.0-rc1~35^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=47a70e19d08c683fd0f58a8a77406047b8e32d50;p=dealii.git Adol-C advanced branching: Fix SymmetricTensor --- diff --git a/include/deal.II/base/symmetric_tensor.h b/include/deal.II/base/symmetric_tensor.h index d9611b895b..877bcc7cd5 100644 --- a/include/deal.II/base/symmetric_tensor.h +++ b/include/deal.II/base/symmetric_tensor.h @@ -1073,7 +1073,7 @@ inline SymmetricTensor & SymmetricTensor::operator = (const Number &d) { - Assert (d==internal::NumberType::value(0.0), ExcMessage ("Only assignment with zero is allowed")); + Assert (numbers::value_is_zero(d), ExcMessage ("Only assignment with zero is allowed")); (void) d; data = internal::NumberType::value(0.0); diff --git a/include/deal.II/base/symmetric_tensor.templates.h b/include/deal.II/base/symmetric_tensor.templates.h index 5d38eadbb2..80d5c7760d 100644 --- a/include/deal.II/base/symmetric_tensor.templates.h +++ b/include/deal.II/base/symmetric_tensor.templates.h @@ -42,7 +42,7 @@ std::array eigenvalues (const SymmetricTensor<2,2,Number> &T) { const Number upp_tri_sq = T[0][1]*T[0][1]; - if (upp_tri_sq == Number(0.0)) + if (upp_tri_sq == internal::NumberType::value(0.0)) { // The tensor is diagonal std::array eig_vals = @@ -59,7 +59,7 @@ eigenvalues (const SymmetricTensor<2,2,Number> &T) const Number tr_T = trace(T); const Number det_T = determinant(T); const Number descrim = tr_T*tr_T - 4.0*det_T; - Assert(descrim > Number(0.0), ExcMessage("The roots of the characteristic polynomial are complex valued.")); + Assert(descrim > internal::NumberType::value(0.0), ExcMessage("The roots of the characteristic polynomial are complex valued.")); const Number sqrt_desc = std::sqrt(descrim); const std::array eig_vals = @@ -81,7 +81,7 @@ std::array eigenvalues (const SymmetricTensor<2,3,Number> &T) { const Number upp_tri_sq = T[0][1]*T[0][1] + T[0][2]*T[0][2] + T[1][2]*T[1][2]; - if (upp_tri_sq == Number(0.0)) + if (upp_tri_sq == internal::NumberType::value(0.0)) { // The tensor is diagonal std::array eig_vals @@ -711,7 +711,6 @@ namespace internal } - template Tensor<2,1,Number> dediagonalize_tensor (const dealii::SymmetricTensor<2,1,Number> &T, @@ -959,6 +958,78 @@ eigenvectors (const SymmetricTensor<2,dim,Number> &T, +#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING + +// Specializations of the above functions for taped Adol-C numbers +// when the advanced branching feature is activated. +// We could copy-paste all of these functions and add the appropriate +// conditional assignments (see the Adol-C manual, section 1.8). +// However, some of the conditions are quite complicated (with possibly +// no 1-1 correspondence for the operations along each branch) so this +// needs some careful attention. For the sake of simplicity and until +// this can be rigourously checked, at this point in time we simply +// disable these functions. + +namespace internal +{ + namespace SymmetricTensorImplementation + { + template<> + struct Inverse<4,3,adouble> + { + static dealii::SymmetricTensor<4,3,adouble> + value (const dealii::SymmetricTensor<4,3,adouble> &/*t*/) + { + AssertThrow(false, ExcADOLCAdvancedBranching()); + return dealii::SymmetricTensor<4,3,adouble>(); + } + }; + } +} + +template<> +std::array +eigenvalues (const SymmetricTensor<2,1,adouble> &/*T*/) +{ + AssertThrow(false, ExcADOLCAdvancedBranching()); + return std::array (); +} + + + +template<> +std::array +eigenvalues (const SymmetricTensor<2,2,adouble> &/*T*/) +{ + AssertThrow(false, ExcADOLCAdvancedBranching()); + return std::array (); +} + + + +template<> +std::array +eigenvalues (const SymmetricTensor<2,3,adouble> &/*T*/) +{ + AssertThrow(false, ExcADOLCAdvancedBranching()); + return std::array (); +} + + + +template +std::array >,std::integral_constant::value> +eigenvectors (const SymmetricTensor<2,dim,adouble> &/*T*/, + const SymmetricTensorEigenvectorMethod /*method*/) +{ + AssertThrow(false, ExcADOLCAdvancedBranching()); + return std::array >,std::integral_constant::value> (); +} + +#endif + + + DEAL_II_NAMESPACE_CLOSE #endif