]> https://gitweb.dealii.org/ - dealii.git/commitdiff
AD Helpers: Add helper for vector functions (QP-level)
authorJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 21 Mar 2019 14:26:06 +0000 (15:26 +0100)
committerJean-Paul Pelteret <jppelteret@gmail.com>
Thu, 21 Mar 2019 14:26:50 +0000 (15:26 +0100)
include/deal.II/differentiation/ad/ad_helpers.h
source/differentiation/ad/ad_helpers.cc
source/differentiation/ad/ad_helpers.inst1.in
source/differentiation/ad/ad_helpers.inst2.in

index 828ad50c932695d99b9e5bb27a5340e13fcbd10a..7619901a82bc00961d3d89892b6c93abcbf6d7ef 100644 (file)
@@ -3326,6 +3326,409 @@ namespace Differentiation
     }; // class ADHelperScalarFunction
 
 
+
+    /**
+     * A helper class that facilitates the evaluation of a vector of functions,
+     * typically one that represents a collection of coupled, multi-dimensional
+     * fields. This class would typically be used to compute the linearization
+     * of a set of kinetic field variables defined at the quadrature point
+     * level.
+     *
+     * An example of its usage in the case of linearizing the kinetic variables
+     * derived from 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;
+     *
+     *   // Declare how many dependent variables we expect to compute.
+     *   // In this case, we will be computing a stress field (a symmetric
+     *   // rank-2 tensor) and the magnetic induction (a vector field).
+     *   // At the same time we define some additional extractors associated
+     *   // with these kinetic fields. In general, these need not be of the same
+     *   // layout as the independent variables.
+     *   const FEValuesExtractors::SymmetricTensor<2> S_dofs (0);
+     *   const FEValuesExtractors::Vector             B_dofs
+     *     (dealii::SymmetricTensor<2,dim>::n_independent_components);
+     *   const unsigned int n_dependent_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.
+     *   ADHelperVectorFunction<dim,double> ad_helper (n_independent_variables,
+     *                                                 n_dependent_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<ADNumberType>::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 Tensor<1,dim,ADNumberType>          H_AD =
+     *       ad_helper.get_sensitive_variables(H_dofs);
+     *
+     *     // Here we define the stress and magnetic induction in terms
+     *     // of the independent values C_AD and H_AD.
+     *     const SymmetricTensor<2, dim, ad_type> S_AD = ...;
+     *     const Tensor<1, dim, ad_type>          B_AD = ...;
+     *
+     *     // Register the definition of the kinetic fields. The second
+     *     // argument to the function provides a non-overlapping ordering
+     *     // of the
+     *     ad_helper.register_dependent_variable(S_AD, S_dofs);
+     *     ad_helper.register_dependent_variable(B_AD, B_dofs);
+     *
+     *     // 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<double> values (ad_helper.n_dependent_variables());
+     *   FullMatrix<double> jacobian (ad_helper.n_dependent_variables(),
+     *                                ad_helper.n_independent_variables());
+     *   ad_helper.compute_values(values);
+     *   ad_helper.compute_jacobian(jacobian);
+     *
+     *   // Extract the desired components of the value vector and its Jacobian
+     *   // matrix. In this example, we use them to compute the Piola-Kirchhoff
+     *   // stress tensor S and its associated tangent, defined
+     *   // as HH = 2*dS/dC...
+     *   const SymmetricTensor<2,dim> S =
+     *     ad_helper.extract_value_component(values,S_dofs);
+     *   const SymmetricTensor<4,dim> HH =
+     *     2.0*ad_helper.extract_jacobian_component(jacobian,S_dofs,C_dofs);
+     *
+     *   // ... the magnetic induction B and its associated tangent defined
+     *   // as BB = dB/dH...
+     *   const Tensor<1,dim> B =
+     *     ad_helper.extract_value_component(values,H_dofs);
+     *   const SymmetricTensor<2,dim> BB =
+     *     symmetrize(ad_helper.extract_jacobian_component(jacobian,B_dofs,H_dofs));
+     *
+     *   // ... and finally the magnetoelastic coupling tangent, defined
+     *   // as PP = -dS/dH. Here the order of the extractor arguments is
+     *   // especially important, as it dictates the field that is being
+     *   // differentiated, and which the directional derivatives are being
+     *   // computed.
+     *   const Tensor<3,dim,double> PP =
+     *     -ad_helper.extract_jacobian_component(jacobian,S_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 <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType = double>
+    class ADHelperVectorFunction
+      : public ADHelperPointLevelFunctionsBase<dim,
+                                               ADNumberTypeCode,
+                                               ScalarType>
+    {
+    public:
+      /**
+       * Type definition for the floating point number type that is used in,
+       * and results from, all computations.
+       */
+      using scalar_type =
+        typename ADHelperBase<ADNumberTypeCode, ScalarType>::scalar_type;
+
+      /**
+       * Type definition for the auto-differentiation number type that is used
+       * in all computations.
+       */
+      using ad_type =
+        typename ADHelperBase<ADNumberTypeCode, ScalarType>::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.
+       * @param[in] n_dependent_variables The number of scalar functions to be
+       * defined that will have a sensitivity to the given independent
+       * variables. In the computation of $\mathbf{f}(\mathbf{X})$, this will
+       * be the number of outputs $\mathbf{f}$, i.e., the dimension of the
+       * image space.
+       */
+      ADHelperVectorFunction(const unsigned int n_independent_variables,
+                             const unsigned int n_dependent_variables);
+
+      /**
+       * Destructor.
+       */
+      virtual ~ADHelperVectorFunction() = default;
+
+      //@}
+
+      /**
+       * @name Dependent variables
+       */
+      //@{
+
+      /**
+       * Register the definition of the vector field
+       * $\boldsymbol{\Psi}(\mathbf{X})$.
+       *
+       * @param[in] funcs A vector of recorded functions that defines the
+       * dependent variables.
+       *
+       * @note For this class that expects only vector field of dependent
+       * variables, 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_variables(const std::vector<ad_type> &funcs);
+
+      /**
+       * Register the definition of the vector field
+       * $\hat{\mathbf{g}}(\mathbf{X}) \subset \boldsymbol{\Psi}(\mathbf{X})$
+       * that may represent a subset of the dependent variables.
+       *
+       * @param[in] funcs The recorded functions that define a set of dependent
+       * variables.
+       * @param[in] extractor An extractor associated with the input field
+       * variables. This effectively defines which components of the global set
+       * of dependent variables this field is associated with.
+       *
+       * @note The input extractor must correspond to the input @p ValueType.
+       * So, for example, if a value is a rank-1 tensor
+       * (i.e. of type Tensor<1,dim,scalar_type>), then the extractor must
+       * be an FEValuesExtractors::Vector or FEValuesExtractors::Tensor<1>.
+       *
+       * @note For taped AD numbers, this operation is only valid in recording mode.
+       */
+      template <typename ValueType, typename ExtractorType>
+      void
+      register_dependent_variable(const ValueType &    funcs,
+                                  const ExtractorType &extractor);
+
+      /**
+       * Compute the value of the vector field $\boldsymbol{\Psi}(\mathbf{X})$.
+       *
+       * @param[out] values A Vector object with the value for each component
+       * of the vector field evaluated at the point defined by the independent
+       * variable values. The output @p values vector has a length
+       * corresponding to @p n_dependent_variables.
+       */
+      void
+      compute_values(Vector<scalar_type> &values) const;
+
+      /**
+       * Compute the Jacobian (first derivative) of the vector field with
+       * respect to all independent variables, i.e.
+       * @f[
+       *   \mathbf{J}(\boldsymbol{\Psi})
+       *      = \frac{\partial\boldsymbol{\Psi}(\mathbf{X})}{\partial\mathbf{X}}
+       * @f]
+       *
+       * @param[out] jacobian A FullMatrix with the gradient of each component
+       * of the vector field evaluated at the point defined by the independent
+       * variable values. The output @p jacobian matrix has
+       * dimensions corresponding to
+       * <code>n_dependent_variables</code>$\times$<code>n_independent_variables</code>.
+       */
+      void
+      compute_jacobian(FullMatrix<scalar_type> &jacobian) const;
+
+
+      /**
+       * Extract the set of functions' values for a subset of dependent
+       * variables
+       * $\mathbf{g} \subset \boldsymbol{\Psi}(\mathbf{X})$.
+       *
+       * @param[in] values A Vector object with the value for each component of
+       * the vector field evaluated at the point defined by the independent
+       * variable values.
+       * @param[in] extractor_row An extractor associated with the input field
+       * variables. This effectively defines which components of the global set
+       * of dependent variables this field is associated with.
+       */
+      template <typename ExtractorType_Row>
+      static typename internal::
+        VectorFieldValue<dim, scalar_type, ExtractorType_Row>::type
+        extract_value_component(const Vector<scalar_type> &values,
+                                const ExtractorType_Row &  extractor_row);
+
+      /**
+       * Extract the Jacobian of the subset of dependent functions
+       * $\mathbf{g} \subset \boldsymbol{\Psi}(\mathbf{X})$
+       * for a subset of independent variables
+       * $\mathbf{A} \subset \mathbf{X}$, i.e.
+       * @f[
+       *   \mathbf{J}(\mathbf{g})
+       *      = \frac{\partial\mathbf{g}(\mathbf{X})}{\partial\mathbf{A}}
+       * @f]
+       * The first index of the Jacobian matrix $\mathbf{J}(\mathbf{g})$
+       * relates to the dependent variables, while the second index relates
+       * to the independent variables.
+       *
+       * @param[in] jacobian The Jacobian of the vector function with respect to
+       * all independent variables, i.e., that returned by compute_jacobian().
+       * @param[in] extractor_row An extractor associated with the input field
+       * variables for which the first index of the Jacobian is extracted.
+       * This effectively defines the correspondence between components of the
+       * global set of dependent variables and the field (representing a
+       * subset of dependent functions) associated with the extractor.
+       * @param[in] extractor_col An extractor associated with the input field
+       * variables for which the second index of the Jacobian is extracted.
+       * This effectively defines the correspondence between components of the
+       * global set of independent variables and the field (representing a
+       * subset of independent variables) associated with the extractor.
+       *
+       * @return A Tensor or SymmetricTensor with its rank and symmetries
+       * determined by the @p extractor_row and @p extractor_col .
+       * This corresponds to subsetting a whole set of rows and columns of the
+       * Jacobian matrix, scaling those entries to take account of tensor
+       * symmetries, and then reshaping the (sub-)matrix so obtained into a
+       * tensor, the final result.
+       * For example, if
+       * @p extractor_row is a FEValuesExtractors::Vector and
+       * @p extractor_col is a FEValuesExtractors::Tensor,
+       * then the returned object is a Tensor of rank 3, with its first
+       * index associated with the field corresponding to the row extractor and
+       * the second and third indices associated with the field corresponding to
+       * the column extractor.
+       * Similarly, if
+       * @p extractor_row is a FEValuesExtractors::SymmetricTensor and
+       * @p extractor_col is a FEValuesExtractors::SymmetricTensor,
+       * then the returned object is a SymmetricTensor of rank 4, with its first
+       * two indices associated with the field corresponding to the row
+       * extractor and the last two indices associated with the field
+       * corresponding to the column extractor.
+       */
+      template <typename ExtractorType_Row, typename ExtractorType_Col>
+      static typename internal::VectorFieldJacobian<dim,
+                                                    scalar_type,
+                                                    ExtractorType_Row,
+                                                    ExtractorType_Col>::type
+      extract_jacobian_component(const FullMatrix<scalar_type> &jacobian,
+                                 const ExtractorType_Row &      extractor_row,
+                                 const ExtractorType_Col &      extractor_col);
+
+      /**
+       * Extract the Jacobian of the subset of dependent functions
+       * $\mathbf{g} \subset \boldsymbol{\Psi}(\mathbf{X})$
+       * for a subset of independent variables
+       * $\mathbf{A} \subset \mathbf{X}$, i.e.
+       * @f[
+       *   \mathbf{J}(\mathbf{g})
+       *      = \frac{\partial\mathbf{g}(\mathbf{X})}{\partial\mathbf{A}}
+       * @f]
+       *
+       * This function is a specialization of the above for rank-0 tensors
+       * (scalars). This corresponds to extracting a single entry of the
+       * Jacobian matrix because both extractors imply selection of just a
+       * single row or column of the matrix.
+       */
+      static Tensor<0, dim, scalar_type>
+      extract_jacobian_component(
+        const FullMatrix<scalar_type> &   jacobian,
+        const FEValuesExtractors::Scalar &extractor_row,
+        const FEValuesExtractors::Scalar &extractor_col);
+
+      /**
+       * Extract the Jacobian of the subset of dependent functions
+       * $\mathbf{g} \subset \boldsymbol{\Psi}(\mathbf{X})$
+       * for a subset of independent variables
+       * $\mathbf{A} \subset \mathbf{X}$, i.e.
+       * @f[
+       *   \mathbf{J}(\mathbf{g})
+       *      = \frac{\partial\mathbf{g}(\mathbf{X})}{\partial\mathbf{A}}
+       * @f]
+       *
+       * This function is a specialization of the above for rank-4 symmetric
+       * tensors.
+       */
+      static SymmetricTensor<4, dim, scalar_type>
+      extract_jacobian_component(
+        const FullMatrix<scalar_type> &               jacobian,
+        const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
+        const FEValuesExtractors::SymmetricTensor<2> &extractor_col);
+
+      //@}
+
+    }; // class ADHelperVectorFunction
+
+
   } // namespace AD
 } // namespace Differentiation
 
@@ -3592,6 +3995,131 @@ namespace Differentiation
     }
 
 
+
+    /* ----------------- ADHelperVectorFunction ----------------- */
+
+
+
+    template <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType>
+    template <typename ValueType, typename ExtractorType>
+    void
+    ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+      register_dependent_variable(const ValueType &    funcs,
+                                  const ExtractorType &extractor)
+    {
+      const std::vector<unsigned int> index_set(
+        internal::extract_field_component_indices<dim>(extractor));
+      for (unsigned int i = 0; i < index_set.size(); ++i)
+        {
+          Assert(this->registered_marked_dependent_variables[index_set[i]] ==
+                   false,
+                 ExcMessage("Overlapping indices for dependent variables."));
+          ADHelperBase<ADNumberTypeCode, ScalarType>::
+            register_dependent_variable(index_set[i],
+                                        internal::get_tensor_entry(funcs, i));
+        }
+    }
+
+
+
+    template <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType>
+    template <typename ExtractorType_Row>
+    typename internal::VectorFieldValue<
+      dim,
+      typename ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+        scalar_type,
+      ExtractorType_Row>::type
+    ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+      extract_value_component(const Vector<scalar_type> &values,
+                              const ExtractorType_Row &  extractor_row)
+    {
+      // NOTE: The order of components must be consistently defined throughout
+      // this class.
+      typename internal::VectorFieldValue<dim, scalar_type, ExtractorType_Row>::
+        type out;
+
+      // Get indexsets for the subblock from which we wish to extract the
+      // gradient values
+      const std::vector<unsigned int> row_index_set(
+        internal::extract_field_component_indices<dim>(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, values[row_index_set[r]]);
+
+      return out;
+    }
+
+
+
+    template <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType>
+    template <typename ExtractorType_Row, typename ExtractorType_Col>
+    typename internal::VectorFieldJacobian<
+      dim,
+      typename ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+        scalar_type,
+      ExtractorType_Row,
+      ExtractorType_Col>::type
+    ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+      extract_jacobian_component(const FullMatrix<scalar_type> &jacobian,
+                                 const ExtractorType_Row &      extractor_row,
+                                 const ExtractorType_Col &      extractor_col)
+    {
+      using InternalJacobian     = internal::VectorFieldJacobian<dim,
+                                                             scalar_type,
+                                                             ExtractorType_Row,
+                                                             ExtractorType_Col>;
+      using InternalExtractorRow = internal::Extractor<dim, ExtractorType_Row>;
+      using InternalExtractorCol = internal::Extractor<dim, ExtractorType_Col>;
+      using JacobianType         = typename InternalJacobian::type;
+
+      // NOTE: The order of components must be consistently defined throughout
+      // this class.
+      JacobianType 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.
+      // <SymmTensor,Vector>. 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<unsigned int> row_index_set(
+        internal::extract_field_component_indices<dim>(
+          extractor_row, false /*ignore_symmetries*/));
+      const std::vector<unsigned int> col_index_set(
+        internal::extract_field_component_indices<dim>(
+          extractor_col, false /*ignore_symmetries*/));
+
+      for (unsigned int index = 0;
+           index < JacobianType::n_independent_components;
+           ++index)
+        {
+          const TableIndices<JacobianType::rank> ti_out =
+            JacobianType::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, jacobian[row_index_set[r]][col_index_set[c]]);
+        }
+
+      return out;
+    }
+
+
   } // namespace AD
 } // namespace Differentiation
 
index 1722791301b5304e1fbefbd01b7d5dd5554c65e2..4d67b34ad7f905fa17a36fe8941b41f16dac8017 100644 (file)
@@ -1640,6 +1640,249 @@ namespace Differentiation
     }
 
 
+
+    /* -------------------- ADHelperVectorFunction -------------------- */
+
+
+
+    template <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType>
+    ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+      ADHelperVectorFunction(const unsigned int n_independent_variables,
+                             const unsigned int n_dependent_variables)
+      : ADHelperPointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>(
+          n_independent_variables,
+          n_dependent_variables)
+    {}
+
+
+
+    template <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType>
+    void
+    ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+      register_dependent_variables(const std::vector<ad_type> &funcs)
+    {
+      Assert(funcs.size() == this->n_dependent_variables(),
+             ExcMessage(
+               "Vector size does not match number of dependent variables"));
+      for (unsigned int i = 0; i < this->n_dependent_variables(); ++i)
+        ADHelperBase<ADNumberTypeCode, ScalarType>::register_dependent_variable(
+          i, funcs[i]);
+    }
+
+
+
+    template <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType>
+    void
+    ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::compute_values(
+      Vector<scalar_type> &values) const
+    {
+      if ((ADNumberTraits<ad_type>::is_taped == true &&
+           this->taped_driver.keep_independent_values() == false) ||
+          ADNumberTraits<ad_type>::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."));
+
+      // We can neglect correctly initializing the entries as
+      // we'll be overwriting them immediately in the succeeding call to
+      // Drivers::values().
+      if (values.size() != this->n_dependent_variables())
+        values.reinit(this->n_dependent_variables(),
+                      true /*omit_zeroing_entries*/);
+
+      if (ADNumberTraits<ad_type>::is_taped == true)
+        {
+          Assert(this->active_tape_index() !=
+                   Numbers<ad_type>::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()));
+
+          this->taped_driver.values(this->active_tape_index(),
+                                    this->n_dependent_variables(),
+                                    this->independent_variable_values,
+                                    values);
+        }
+      else
+        {
+          Assert(ADNumberTraits<ad_type>::is_tapeless == true,
+                 ExcInternalError());
+          this->tapeless_driver.values(this->dependent_variables, values);
+        }
+    }
+
+
+
+    template <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType>
+    void
+    ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::compute_jacobian(
+      FullMatrix<scalar_type> &jacobian) const
+    {
+      if ((ADNumberTraits<ad_type>::is_taped == true &&
+           this->taped_driver.keep_independent_values() == false) ||
+          ADNumberTraits<ad_type>::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."));
+
+      // We can neglect correctly initializing the entries as
+      // we'll be overwriting them immediately in the succeeding call to
+      // Drivers::jacobian().
+      if (jacobian.m() != this->n_dependent_variables() ||
+          jacobian.n() != this->n_independent_variables())
+        jacobian.reinit({this->n_dependent_variables(),
+                         this->n_independent_variables()},
+                        true /*omit_default_initialization*/);
+
+      if (ADNumberTraits<ad_type>::is_taped == true)
+        {
+          Assert(this->active_tape_index() !=
+                   Numbers<ad_type>::invalid_tape_index,
+                 ExcMessage("Invalid tape index"));
+          Assert(this->is_recording() == false,
+                 ExcMessage(
+                   "Cannot compute Jacobian 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.jacobian(this->active_tape_index(),
+                                      this->n_dependent_variables(),
+                                      this->independent_variable_values,
+                                      jacobian);
+        }
+      else
+        {
+          Assert(ADNumberTraits<ad_type>::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.jacobian(this->independent_variables,
+                                         this->dependent_variables,
+                                         jacobian);
+        }
+
+      for (unsigned int j = 0; j < this->n_independent_variables(); j++)
+        {
+          // Because we perform just a single differentiation
+          // operation with respect to the "column" variables,
+          // we only need to consider them for symmetry conditions.
+          if (this->is_symmetric_independent_variable(j) == true)
+            for (unsigned int i = 0; i < this->n_dependent_variables(); i++)
+              jacobian[i][j] *= 0.5;
+        }
+    }
+
+
+
+    template <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType>
+    Tensor<0,
+           dim,
+           typename ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+             scalar_type>
+    ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+      extract_jacobian_component(
+        const FullMatrix<scalar_type> &   jacobian,
+        const FEValuesExtractors::Scalar &extractor_row,
+        const FEValuesExtractors::Scalar &extractor_col)
+    {
+      // 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<unsigned int> row_index_set(
+        internal::extract_field_component_indices<dim>(extractor_row));
+      const std::vector<unsigned int> col_index_set(
+        internal::extract_field_component_indices<dim>(extractor_col));
+      Assert(row_index_set.size() == 1, ExcInternalError());
+      Assert(col_index_set.size() == 1, ExcInternalError());
+
+      internal::set_tensor_entry(out,
+                                 0,
+                                 jacobian[row_index_set[0]][col_index_set[0]]);
+
+      return out;
+    }
+
+
+
+    template <int                  dim,
+              enum AD::NumberTypes ADNumberTypeCode,
+              typename ScalarType>
+    SymmetricTensor<4,
+                    dim,
+                    typename ADHelperVectorFunction<dim,
+                                                    ADNumberTypeCode,
+                                                    ScalarType>::scalar_type>
+    ADHelperVectorFunction<dim, ADNumberTypeCode, ScalarType>::
+      extract_jacobian_component(
+        const FullMatrix<scalar_type> &               jacobian,
+        const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
+        const FEValuesExtractors::SymmetricTensor<2> &extractor_col)
+    {
+      // 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<unsigned int> row_index_set(
+        internal::extract_field_component_indices<dim>(extractor_row));
+      const std::vector<unsigned int> col_index_set(
+        internal::extract_field_component_indices<dim>(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, jacobian[row_index_set[r]][col_index_set[c]]);
+          }
+
+      return out;
+    }
+
+
   } // namespace AD
 } // namespace Differentiation
 
index a5c9b87da3bfe3c4104c25dec92d30d8f104168f..1c1ff5a534b8f6d3de756a582ae88f2e541de4ef 100644 (file)
@@ -125,6 +125,14 @@ for (deal_II_dimension : DIMENSIONS ; number : REAL_SCALARS)
     template
     class ADHelperScalarFunction<deal_II_dimension,NumberTypes::adolc_tapeless,number>;
 
+    // -------------------------- ADHelperVectorFunction ----------------------
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::adolc_taped,number>;
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::adolc_tapeless,number>;
+
     \}
     \}
 }
@@ -153,6 +161,14 @@ for (deal_II_dimension : DIMENSIONS)
     template
     class ADHelperScalarFunction<deal_II_dimension,NumberTypes::adolc_tapeless,typename NumberTraits<double,NumberTypes::adolc_tapeless>::ad_type>;
 
+    // -------------------------- ADHelperVectorFunction ----------------------
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::adolc_taped,typename NumberTraits<double,NumberTypes::adolc_taped>::ad_type>;
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::adolc_tapeless,typename NumberTraits<double,NumberTypes::adolc_tapeless>::ad_type>;
+
     \}
     \}
 }
index e886a8575ad0ab436a8ac31df608218d1af1ef2a..800b697f99aea1e145a3a565881afefa5731d912 100644 (file)
@@ -186,6 +186,20 @@ for (deal_II_dimension : DIMENSIONS ; number : REAL_SCALARS)
     template
     class ADHelperScalarFunction<deal_II_dimension,NumberTypes::sacado_rad_dfad,number>;
 
+    // -------------------------- ADHelperVectorFunction ----------------------
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::sacado_dfad_dfad,number>;
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::sacado_dfad,number>;
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::sacado_rad,number>;
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::sacado_rad_dfad,number>;
+
     \}
     \}
 }
@@ -225,6 +239,20 @@ for (deal_II_dimension : DIMENSIONS)
 
     template
     class ADHelperScalarFunction<deal_II_dimension,NumberTypes::sacado_rad_dfad,typename NumberTraits<double,NumberTypes::sacado_rad_dfad>::ad_type>;
+
+    // -------------------------- ADHelperVectorFunction ----------------------
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::sacado_dfad_dfad,typename NumberTraits<double,NumberTypes::sacado_dfad_dfad>::ad_type>;
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::sacado_dfad,typename NumberTraits<double,NumberTypes::sacado_dfad>::ad_type>;
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::sacado_rad,typename NumberTraits<double,NumberTypes::sacado_rad>::ad_type>;
+
+    template
+    class ADHelperVectorFunction<deal_II_dimension,NumberTypes::sacado_rad_dfad,typename NumberTraits<double,NumberTypes::sacado_rad_dfad>::ad_type>;
     
     \}
     \}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.