From: Jean-Paul Pelteret Date: Fri, 15 Feb 2019 10:52:49 +0000 (+0100) Subject: AD Helpers: Add helper for scalar functions (QP-level) X-Git-Tag: v9.1.0-rc1~327^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cb435e24a1f9c72adc9759363d4e845ec0b23970;p=dealii.git AD Helpers: Add helper for scalar functions (QP-level) --- diff --git a/doc/news/changes/major/20190220Jean-PaulPelteret b/doc/news/changes/major/20190220Jean-PaulPelteret new file mode 100644 index 0000000000..c9df4e7cb7 --- /dev/null +++ b/doc/news/changes/major/20190220Jean-PaulPelteret @@ -0,0 +1,9 @@ +New: A new class Differentiation::AD::ADHelperScalarFunction has been +added to help implement point-wise scalar functions using automatic +differentiation. In particular, this class is designed to compute the +gradient and Hessian of a scalar function that is parameterized in +terms of scalar, vector, and tensor arguments. One example of its use +would be to compute the derivatives of a multi-field constitutive law +that is expressed in terms of an energy function. +
+(Jean-Paul Pelteret, 2019/02/20) diff --git a/include/deal.II/differentiation/ad/ad_helpers.h b/include/deal.II/differentiation/ad/ad_helpers.h index 4f36c1dbfe..f884ee4694 100644 --- a/include/deal.II/differentiation/ad/ad_helpers.h +++ b/include/deal.II/differentiation/ad/ad_helpers.h @@ -2934,6 +2934,352 @@ namespace Differentiation }; // class ADHelperPointLevelFunctionsBase + + /** + * A helper class that facilitates the evaluation of a scalar function, + * its first derivatives (gradient), and its second derivatives (Hessian). + * This class would typically be used to compute the first and second + * derivatives of a stored energy function defined at a quadrature + * point. It can also be used to compute derivatives of any other scalar + * field so long as all its dependencies on the independent variables are + * explicit (that is to say, no independent variables may have some implicit + * dependence on one another). + * + * An example of its usage in the case of a multi-field constitutive law + * might be as follows: + * @code + * // Define some extractors that will help us set independent variables + * // and later get the computed values related to the dependent + * // variables. Each of these extractors is related to the gradient of a + * // component of the solution field (in this case, displacement and + * // magnetic scalar potential). Here "C" is the right Cauchy-Green + * // tensor and "H" is the magnetic field. + * const FEValuesExtractors::SymmetricTensor<2> C_dofs (0); + * const FEValuesExtractors::Vector H_dofs + * (dealii::SymmetricTensor<2,dim>::n_independent_components); + * const unsigned int n_independent_variables = + * SymmetricTensor<2,dim>::n_independent_components + + * Tensor<1,dim>::n_independent_components; + * + * // Define the helper that we will use in the AD computations for our + * // scalar energy function. Note that we expect it to return values of + * // type double. + * ADHelperScalarFunction ad_helper (n_independent_variables); + * using ADNumberType = typename ADHelper::ad_type; + * + * // Compute the fields that provide the independent values. + * // When the tape is being replayed, these should be set to something + * // meaningful. + * const Tensor<1,dim> H = ...; + * const SymmetricTensor<2,dim> C = ...; + * + * // If using a taped AD number, then at this point we would initiate + * // taping of the expression for the material stored energy function + * // for this particular set of material parameters: + * + * // Select a tape number to record to + * const typename Types::tape_index tape_index = ...; + * + * // Indicate that we are about to start tracing the operations for + * // function evaluation on the tape. If this tape has already been + * // used (i.e., the operations are already recorded) then we + * // (optionally) load the tape and reuse this data. + * const bool is_recording + * = ad_helper.start_recording_operations(tape_index); + * + * // The steps that follow in the recording phase are required for + * // tapeless methods as well. + * if (is_recording == true) + * { + * // This is the "recording" phase of the operations. + * + * // First, we set the values for all fields. + * // These could happily be set to anything, unless the function will + * // be evaluated along a branch not otherwise traversed during later + * // use. For this reason, in this example instead of using some dummy + * // values, we'll actually map out the function at the same point + * // around which we'll later linearize it. + * ad_helper.register_independent_variable(H, H_dofs); + * ad_helper.register_independent_variable(C, C_dofs); + * + * // NOTE: We have to extract the sensitivities in the order we wish to + * // introduce them. So this means we have to do it by logical order + * // of the extractors that we've created. + * const SymmetricTensor<2,dim,ADNumberType> C_AD = + * ad_helper.get_sensitive_variables(C_dofs); const + * const Tensor<1,dim,ADNumberType> H_AD = + * ad_helper.get_sensitive_variables(H_dofs); + * + * // Here we define the material stored energy function. + * // This example is sufficiently complex to warrant the use of AD to, + * // at the very least, verify an unassisted implementation. + * const double mu_e = 10; // Shear modulus + * const double lambda_e = 15; // Lame parameter + * const double mu_0 = 4*M_PI*1e-7; // Magnetic permeability constant + * const double mu_r = 5; // Relative magnetic permeability + * + * const ADNumberType J = std::sqrt(determinant(C_AD)); + * const SymmetricTensor<2,dim,ADNumberType> C_inv_AD = invert(C_AD); + * const ADNumberType psi = + * 0.5*mu_e*(1.0+std::tanh((H_AD*H_AD)/100.0))* + * (trace(C_AD) - dim - 2*std::log(J)) + + * lambda_e*std::log(J)*std::log(J) - + * 0.5*mu_0*mu_r*J*H_AD*C_inv_AD*H_AD; + * + * // Register the definition of the total stored energy + * ad_helper.register_dependent_variable(psi_CH); + * + * // Indicate that we have completed tracing the operations onto + * // the tape. + * ad_helper.stop_recording_operations(false); // write_tapes_to_file + * } + * else + * { + * // This is the "tape reuse" phase of the operations. + * // Here we will leverage the already traced operations that reside + * // on a tape, and simply re-evaluate the tape at a different point + * // to get the function values and their derivatives. + * + * // Load the existing tape to be reused + * ad_helper.activate_recorded_tape(tape_index); + * + * // Set the new values of the independent variables where the + * // recorded dependent functions are to be evaluated (and + * // differentiated around). + * ad_helper.set_independent_variable(C, C_dofs); + * ad_helper.set_independent_variable(H, H_dofs); + * } + * + * // Play the tape and store the output function value, its gradient and + * // linearization. These are expensive to compute, so we'll do this once + * // and extract the desired values from these intermediate outputs. + * Vector Dpsi (ad_helper.n_dependent_variables()); + * FullMatrix D2psi (ad_helper.n_dependent_variables(), + * ad_helper.n_dependent_variables()); + * const double psi = ad_helper.compute_value(); + * ad_helper.compute_gradient(Dpsi); + * ad_helper.compute_hessian(D2psi); + * + * // Extract the desired components of the gradient vector and Hessian + * // matrix. In this example, we use them to compute the Piola-Kirchhoff + * // stress tensor and its associated tangent, defined by thermodynamic + * // arguments as S = 2*dpsi/dC and HH = 2*dS/dC... + * const SymmetricTensor<2,dim> S = + * 2.0*ad_helper.extract_gradient_component(Dpsi,C_dofs); + * const SymmetricTensor<4,dim> HH = + * 4.0*ad_helper.extract_hessian_component(D2psi,C_dofs,C_dofs); + * + * // ... the magnetic induction and its associated tangent defined + * // as B = -dpsi/dH and BB = dB/dH... + * const Tensor<1,dim> B = + * -ad_helper.extract_gradient_component(Dpsi,H_dofs); + * const SymmetricTensor<2,dim> BB = + * -symmetrize(ad_helper.extract_hessian_component(D2psi,H_dofs,H_dofs)); + * + * // ... and finally the magnetoelastic coupling tangent, defined + * // as PP = -dS/dH = -d/dH(2*dpsi/dC). Here the order of the extractor + * // arguments is especially important, as it dictates the order in which + * // the directional derivatives are taken. + * const Tensor<3,dim,double> PP = + * -2.0*ad_helper.extract_hessian_component(D2psi,C_dofs,H_dofs) + * @endcode + * + * @warning ADOL-C does not support the standard threading models used by + * deal.II, so this class should @b not be embedded within a multithreaded + * function when using ADOL-C number types. It is, however, suitable for use + * in both serial and MPI routines. + * + * @author Jean-Paul Pelteret, 2016, 2017, 2018 + */ + template + class ADHelperScalarFunction + : public ADHelperPointLevelFunctionsBase + { + public: + /** + * Type definition for the floating point number type that is used in, + * and results from, all computations. + */ + using scalar_type = + typename ADHelperBase::scalar_type; + + /** + * Type definition for the auto-differentiation number type that is used + * in all computations. + */ + using ad_type = + typename ADHelperBase::ad_type; + + /** + * @name Constructor / destructor + */ + //@{ + + /** + * The constructor for the class. + * + * @param[in] n_independent_variables The number of independent variables + * that will be used in the definition of the functions that it is + * desired to compute the sensitivities of. In the computation of + * $\mathbf{f}(\mathbf{X})$, this will be the number of inputs + * $\mathbf{X}$, i.e., the dimension of the domain space. + */ + ADHelperScalarFunction(const unsigned int n_independent_variables); + + /** + * Destructor. + */ + virtual ~ADHelperScalarFunction() = default; + + //@} + + /** + * @name Dependent variables + */ + //@{ + + /** + * Register the definition of the scalar field $\Psi(\mathbf{X})$. + * + * @param[in] func The recorded function that defines a dependent + * variable. + * + * @note For this class that expects only one dependent variable, this + * function must only be called once per tape. + * + * @note For taped AD numbers, this operation is only valid in recording mode. + */ + void + register_dependent_variable(const ad_type &func); + + /** + * Compute the value of the scalar field $\Psi(\mathbf{X})$ using the + * tape as opposed to executing the source code. + * + * @return A scalar object with the value for the scalar field evaluated + * at the point defined by the independent variable values. + */ + scalar_type + compute_value() const; + + /** + * Compute the gradient (first derivative) of the scalar field with + * respect to all independent variables, i.e. + * @f[ + * \frac{\partial\Psi(\mathbf{X})}{\partial\mathbf{X}} + * @f] + * + * @param[out] gradient A Vector with the values for the scalar field + * gradient (first derivatives) evaluated at the point defined by the + * independent variable values. + */ + void + compute_gradient(Vector &gradient) const; + + /** + * Compute the Hessian (second derivative) of the scalar field with + * respect to all independent variables, i.e. + * @f[ + * \frac{\partial^{2}\Psi(\mathbf{X})}{\partial\mathbf{X} \otimes + * \partial\mathbf{X}} + * @f] + * + * @param[out] hessian A FullMatrix with the values for the scalar field + * Hessian (second derivatives) evaluated at the point defined by the + * independent variable values. + */ + void + compute_hessian(FullMatrix &hessian) const; + + /** + * Extract the function gradient for a subset of independent variables + * $\mathbf{A} \subset \mathbf{X}$, i.e. + * @f[ + * \frac{\partial\Psi(\mathbf{X})}{\partial\mathbf{A}} + * @f] + * + * @param[in] gradient The gradient of the scalar function with respect to + * all independent variables, i.e. that returned by compute_gradient(). + * @param[in] extractor_row An extractor associated with the input field + * variables. This effectively defines which components of the global set + * of independent variables this field is associated with. + */ + template + typename internal::ScalarFieldGradient::type + extract_gradient_component(const Vector &gradient, + const ExtractorType_Row &extractor_row) const; + + /** + * Extract the function Hessian for a subset of independent variables + * $\mathbf{A},\mathbf{B} \subset \mathbf{X}$, i.e. + * @f[ + * \frac{}{\partial\mathbf{B}} \left[ + * \frac{\partial\Psi(\mathbf{X})}{\partial\mathbf{A}} \right] = + * \frac{\partial^{2}\Psi(\mathbf{X})}{\partial\mathbf{B} \otimes + * \partial\mathbf{A}} + * @f] + * + * @param[in] hessian The Hessian of the scalar function with respect to + * all independent variables, i.e. that returned by compute_hessian(). + * @param[in] extractor_row An extractor associated with the input field + * variables for which the first index of the Hessian is extracted. + * @param[in] extractor_col An extractor associated with the input field + * variables for which the second index of the Hessian is extracted. + */ + template + typename internal::ScalarFieldHessian::type + extract_hessian_component(const FullMatrix &hessian, + const ExtractorType_Row & extractor_row, + const ExtractorType_Col &extractor_col) const; + + /** + * Extract the function Hessian for a subset of independent variables + * $\mathbf{A},\mathbf{B} \subset \mathbf{X}$, i.e. + * @f[ + * \frac{}{\partial\mathbf{B}} \left[ + * \frac{\partial\Psi(\mathbf{X})}{\partial\mathbf{A}} \right] + * @f] + * + * This function is a specialization of the above for rank-0 tensors + * (scalars) + */ + Tensor<0, dim, scalar_type> + extract_hessian_component( + const FullMatrix & hessian, + const FEValuesExtractors::Scalar &extractor_row, + const FEValuesExtractors::Scalar &extractor_col) const; + + /** + * Extract the function Hessian for a subset of independent variables + * $\mathbf{A},\mathbf{B} \subset \mathbf{X}$, i.e. + * @f[ + * \frac{}{\partial\mathbf{B}} \left[ + * \frac{\partial\Psi(\mathbf{X})}{\partial\mathbf{A}} \right] + * @f] + * + * This function is a specialization of the above for rank-4 symmetric + * tensors + */ + SymmetricTensor<4, dim, scalar_type> + extract_hessian_component( + const FullMatrix & hessian, + const FEValuesExtractors::SymmetricTensor<2> &extractor_row, + const FEValuesExtractors::SymmetricTensor<2> &extractor_col) const; + + //@} + + }; // class ADHelperScalarFunction + + } // namespace AD } // namespace Differentiation @@ -3099,6 +3445,107 @@ namespace Differentiation } + + /* ----------------- ADHelperScalarFunction ----------------- */ + + + + template + template + typename internal::ScalarFieldGradient< + dim, + typename ADHelperScalarFunction:: + scalar_type, + ExtractorType_Row>::type + ADHelperScalarFunction:: + extract_gradient_component(const Vector &gradient, + const ExtractorType_Row & extractor_row) const + { + // NOTE: The order of components must be consistently defined throughout + // this class. + typename internal:: + ScalarFieldGradient::type out; + + // Get indexsets for the subblock from which we wish to extract the + // gradient values + const std::vector row_index_set( + internal::extract_field_component_indices(extractor_row)); + Assert(out.n_independent_components == row_index_set.size(), + ExcMessage("Not all tensor components have been extracted!")); + for (unsigned int r = 0; r < row_index_set.size(); ++r) + internal::set_tensor_entry(out, r, gradient[row_index_set[r]]); + + return out; + } + + + + template + template + typename internal::ScalarFieldHessian< + dim, + typename ADHelperScalarFunction:: + scalar_type, + ExtractorType_Row, + ExtractorType_Col>::type + ADHelperScalarFunction:: + extract_hessian_component(const FullMatrix &hessian, + const ExtractorType_Row & extractor_row, + const ExtractorType_Col &extractor_col) const + { + using InternalHessian = internal::ScalarFieldHessian; + using InternalExtractorRow = internal::Extractor; + using InternalExtractorCol = internal::Extractor; + using HessianType = typename InternalHessian::type; + + // NOTE: The order of components must be consistently defined throughout + // this class. + HessianType out; + + // Get indexsets for the subblocks from which we wish to extract the + // Hessian values + // NOTE: Here we have to do some clever accounting when the + // one extractor is a symmetric Tensor and the other is not, e.g. + // . In this scenario the return type is a + // non-symmetric Tensor<3,dim> but we have to fetch information from a + // SymmTensor row/column that has too few entries to fill the output + // tensor. So we must duplicate the relevant entries in the row/column + // indexset to fetch off-diagonal components that are Otherwise + // non-existent in a SymmTensor. + const std::vector row_index_set( + internal::extract_field_component_indices( + extractor_row, false /*ignore_symmetries*/)); + const std::vector col_index_set( + internal::extract_field_component_indices( + extractor_col, false /*ignore_symmetries*/)); + + for (unsigned int index = 0; + index < HessianType::n_independent_components; + ++index) + { + const TableIndices ti_out = + HessianType::unrolled_to_component_indices(index); + const unsigned int r = + InternalExtractorRow::local_component(ti_out, 0); + const unsigned int c = + InternalExtractorCol::local_component(ti_out, + InternalExtractorRow::rank); + + internal::set_tensor_entry( + out, index, hessian[row_index_set[r]][col_index_set[c]]); + } + + return out; + } + + } // namespace AD } // namespace Differentiation diff --git a/source/differentiation/ad/ad_helpers.cc b/source/differentiation/ad/ad_helpers.cc index 7789a47543..a1e0809b5b 100644 --- a/source/differentiation/ad/ad_helpers.cc +++ b/source/differentiation/ad/ad_helpers.cc @@ -1312,6 +1312,335 @@ namespace Differentiation } + + /* -------------------- ADHelperScalarFunction -------------------- */ + + + + template + ADHelperScalarFunction:: + ADHelperScalarFunction(const unsigned int n_independent_variables) + : ADHelperPointLevelFunctionsBase( + n_independent_variables, + 1) + {} + + + + template + void + ADHelperScalarFunction:: + register_dependent_variable(const ad_type &func) + { + Assert(this->n_dependent_variables() == 1, ExcInternalError()); + ADHelperBase::register_dependent_variable( + 0, func); + } + + + + template + typename ADHelperScalarFunction:: + scalar_type + ADHelperScalarFunction::compute_value() + const + { + if ((ADNumberTraits::is_taped == true && + this->taped_driver.keep_independent_values() == false) || + ADNumberTraits::is_tapeless == true) + { + Assert( + this->n_registered_independent_variables() == + this->n_independent_variables(), + ExcMessage( + "Not all values of sensitivities have been registered or subsequently set!")); + } + Assert(this->n_registered_dependent_variables() == + this->n_dependent_variables(), + ExcMessage("Not all dependent variables have been registered.")); + + Assert( + this->n_dependent_variables() == 1, + ExcMessage( + "The ADHelperScalarFunction class expects there to be only one dependent variable.")); + + if (ADNumberTraits::is_taped == true) + { + Assert(this->active_tape_index() != + Numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(this->is_recording() == false, + ExcMessage( + "Cannot compute values while tape is being recorded.")); + Assert(this->independent_variable_values.size() == + this->n_independent_variables(), + ExcDimensionMismatch(this->independent_variable_values.size(), + this->n_independent_variables())); + + return this->taped_driver.value(this->active_tape_index(), + this->independent_variable_values); + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + return this->tapeless_driver.value(this->dependent_variables); + } + } + + + template + void + ADHelperScalarFunction::compute_gradient( + Vector &gradient) const + { + if ((ADNumberTraits::is_taped == true && + this->taped_driver.keep_independent_values() == false) || + ADNumberTraits::is_tapeless == true) + { + Assert( + this->n_registered_independent_variables() == + this->n_independent_variables(), + ExcMessage( + "Not all values of sensitivities have been registered or subsequently set!")); + } + Assert(this->n_registered_dependent_variables() == + this->n_dependent_variables(), + ExcMessage("Not all dependent variables have been registered.")); + + Assert( + this->n_dependent_variables() == 1, + ExcMessage( + "The ADHelperScalarFunction class expects there to be only one dependent variable.")); + + // We can neglect correctly initializing the entries as + // we'll be overwriting them immediately in the succeeding call to + // Drivers::gradient(). + if (gradient.size() != this->n_independent_variables()) + gradient.reinit(this->n_independent_variables(), + true /*omit_zeroing_entries*/); + + if (ADNumberTraits::is_taped == true) + { + Assert(this->active_tape_index() != + Numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(this->is_recording() == false, + ExcMessage( + "Cannot compute gradient while tape is being recorded.")); + Assert(this->independent_variable_values.size() == + this->n_independent_variables(), + ExcDimensionMismatch(this->independent_variable_values.size(), + this->n_independent_variables())); + + this->taped_driver.gradient(this->active_tape_index(), + this->independent_variable_values, + gradient); + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + Assert(this->independent_variables.size() == + this->n_independent_variables(), + ExcDimensionMismatch(this->independent_variables.size(), + this->n_independent_variables())); + + this->tapeless_driver.gradient(this->independent_variables, + this->dependent_variables, + gradient); + } + + // Account for symmetries of tensor components + for (unsigned int i = 0; i < this->n_independent_variables(); i++) + { + if (this->is_symmetric_independent_variable(i) == true) + gradient[i] *= 0.5; + } + } + + + + template + void + ADHelperScalarFunction::compute_hessian( + FullMatrix &hessian) const + { + Assert(AD::ADNumberTraits::n_supported_derivative_levels >= 2, + ExcMessage( + "Cannot computed function Hessian: AD number type does" + "not support the calculation of second order derivatives.")); + + if ((ADNumberTraits::is_taped == true && + this->taped_driver.keep_independent_values() == false)) + { + Assert( + this->n_registered_independent_variables() == + this->n_independent_variables(), + ExcMessage( + "Not all values of sensitivities have been registered or subsequently set!")); + } + Assert(this->n_registered_dependent_variables() == + this->n_dependent_variables(), + ExcMessage("Not all dependent variables have been registered.")); + + Assert( + this->n_dependent_variables() == 1, + ExcMessage( + "The ADHelperScalarFunction class expects there to be only one dependent variable.")); + + // We can neglect correctly initializing the entries as + // we'll be overwriting them immediately in the succeeding call to + // Drivers::hessian(). + if (hessian.m() != this->n_independent_variables() || + hessian.n() != this->n_independent_variables()) + hessian.reinit({this->n_independent_variables(), + this->n_independent_variables()}, + true /*omit_default_initialization*/); + + if (ADNumberTraits::is_taped == true) + { + Assert(this->active_tape_index() != + Numbers::invalid_tape_index, + ExcMessage("Invalid tape index")); + Assert(this->is_recording() == false, + ExcMessage( + "Cannot compute Hessian while tape is being recorded.")); + Assert(this->independent_variable_values.size() == + this->n_independent_variables(), + ExcDimensionMismatch(this->independent_variable_values.size(), + this->n_independent_variables())); + + this->taped_driver.hessian(this->active_tape_index(), + this->independent_variable_values, + hessian); + } + else + { + Assert(ADNumberTraits::is_tapeless == true, + ExcInternalError()); + Assert(this->independent_variables.size() == + this->n_independent_variables(), + ExcDimensionMismatch(this->independent_variables.size(), + this->n_independent_variables())); + + this->tapeless_driver.hessian(this->independent_variables, + this->dependent_variables, + hessian); + } + + // Account for symmetries of tensor components + for (unsigned int i = 0; i < this->n_independent_variables(); i++) + for (unsigned int j = 0; j < i + 1; j++) + { + if (this->is_symmetric_independent_variable(i) == true && + this->is_symmetric_independent_variable(j) == true) + { + hessian[i][j] *= 0.25; + if (i != j) + hessian[j][i] *= 0.25; + } + else if ((this->is_symmetric_independent_variable(i) == true && + this->is_symmetric_independent_variable(j) == false) || + (this->is_symmetric_independent_variable(j) == true && + this->is_symmetric_independent_variable(i) == false)) + { + hessian[i][j] *= 0.5; + if (i != j) + hessian[j][i] *= 0.5; + } + } + } + + + + template + Tensor<0, + dim, + typename ADHelperScalarFunction:: + scalar_type> + ADHelperScalarFunction:: + extract_hessian_component( + const FullMatrix & hessian, + const FEValuesExtractors::Scalar &extractor_row, + const FEValuesExtractors::Scalar &extractor_col) const + { + // NOTE: It is necessary to make special provision for the case when the + // HessianType is scalar. Unfortunately Tensor<0,dim> does not provide + // the function unrolled_to_component_indices! + // NOTE: The order of components must be consistently defined throughout + // this class. + Tensor<0, dim, scalar_type> out; + + // Get indexsets for the subblocks from which we wish to extract the + // matrix values + const std::vector row_index_set( + internal::extract_field_component_indices(extractor_row)); + const std::vector col_index_set( + internal::extract_field_component_indices(extractor_col)); + Assert(row_index_set.size() == 1, ExcInternalError()); + Assert(col_index_set.size() == 1, ExcInternalError()); + + internal::set_tensor_entry(out, + 0, + hessian[row_index_set[0]][col_index_set[0]]); + + return out; + } + + + + template + SymmetricTensor<4, + dim, + typename ADHelperScalarFunction::scalar_type> + ADHelperScalarFunction:: + extract_hessian_component( + const FullMatrix & hessian, + const FEValuesExtractors::SymmetricTensor<2> &extractor_row, + const FEValuesExtractors::SymmetricTensor<2> &extractor_col) const + { + // NOTE: The order of components must be consistently defined throughout + // this class. NOTE: We require a specialisation for rank-4 symmetric + // tensors because they + // do not define their rank, and setting data using TableIndices is + // somewhat specialised as well. + SymmetricTensor<4, dim, scalar_type> out; + + // Get indexsets for the subblocks from which we wish to extract the + // matrix values + const std::vector row_index_set( + internal::extract_field_component_indices(extractor_row)); + const std::vector col_index_set( + internal::extract_field_component_indices(extractor_col)); + + for (unsigned int r = 0; r < row_index_set.size(); ++r) + for (unsigned int c = 0; c < col_index_set.size(); ++c) + { + internal::set_tensor_entry( + out, r, c, hessian[row_index_set[r]][col_index_set[c]]); + } + + return out; + } + + } // namespace AD } // namespace Differentiation diff --git a/source/differentiation/ad/ad_helpers.inst1.in b/source/differentiation/ad/ad_helpers.inst1.in index ef44c93984..a5c9b87da3 100644 --- a/source/differentiation/ad/ad_helpers.inst1.in +++ b/source/differentiation/ad/ad_helpers.inst1.in @@ -117,6 +117,14 @@ for (deal_II_dimension : DIMENSIONS ; number : REAL_SCALARS) template class ADHelperPointLevelFunctionsBase; + // -------------------------- ADHelperScalarFunction ---------------------- + + template + class ADHelperScalarFunction; + + template + class ADHelperScalarFunction; + \} \} } @@ -137,6 +145,14 @@ for (deal_II_dimension : DIMENSIONS) template class ADHelperPointLevelFunctionsBase::ad_type>; + // -------------------------- ADHelperScalarFunction ---------------------- + + template + class ADHelperScalarFunction::ad_type>; + + template + class ADHelperScalarFunction::ad_type>; + \} \} } diff --git a/source/differentiation/ad/ad_helpers.inst2.in b/source/differentiation/ad/ad_helpers.inst2.in index 5fa4b2511e..e886a8575a 100644 --- a/source/differentiation/ad/ad_helpers.inst2.in +++ b/source/differentiation/ad/ad_helpers.inst2.in @@ -172,6 +172,20 @@ for (deal_II_dimension : DIMENSIONS ; number : REAL_SCALARS) template class ADHelperPointLevelFunctionsBase; + // -------------------------- ADHelperScalarFunction ---------------------- + + template + class ADHelperScalarFunction; + + template + class ADHelperScalarFunction; + + template + class ADHelperScalarFunction; + + template + class ADHelperScalarFunction; + \} \} } @@ -198,6 +212,20 @@ for (deal_II_dimension : DIMENSIONS) template class ADHelperPointLevelFunctionsBase::ad_type>; + // -------------------------- ADHelperScalarFunction ---------------------- + + template + class ADHelperScalarFunction::ad_type>; + + template + class ADHelperScalarFunction::ad_type>; + + template + class ADHelperScalarFunction::ad_type>; + + template + class ADHelperScalarFunction::ad_type>; + \} \} }