Reference documentation for deal.II version GIT d77e5ebb0a 2023-01-27 22:35:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
ad_helpers.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_differentiation_ad_ad_helpers_h
17 #define dealii_differentiation_ad_ad_helpers_h
18 
19 #include <deal.II/base/config.h>
20 
21 #if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
22 
23 # include <deal.II/base/numbers.h>
25 # include <deal.II/base/tensor.h>
26 
33 
35 
36 # include <deal.II/lac/full_matrix.h>
37 # include <deal.II/lac/vector.h>
38 
39 # include <algorithm>
40 # include <iostream>
41 # include <iterator>
42 # include <numeric>
43 
45 
46 namespace Differentiation
47 {
48  namespace AD
49  {
167  template <enum AD::NumberTypes ADNumberTypeCode,
168  typename ScalarType = double>
170  {
171  public:
176  using scalar_type =
178 
183  using ad_type =
185 
205  HelperBase(const unsigned int n_independent_variables,
206  const unsigned int n_dependent_variables);
207 
211  virtual ~HelperBase() = default;
212 
224  std::size_t
225  n_independent_variables() const;
226 
231  std::size_t
232  n_dependent_variables() const;
233 
242  void
243  print(std::ostream &stream) const;
244 
251  void
252  print_values(std::ostream &stream) const;
253 
264  void
265  print_tape_stats(const typename Types<ad_type>::tape_index tape_index,
266  std::ostream & stream) const;
267 
293  static void
295  const bool ensure_persistent_setting = true);
296 
339  virtual void
340  reset(const unsigned int n_independent_variables =
342  const unsigned int n_dependent_variables =
344  const bool clear_registered_tapes = true);
345 
350  bool
351  is_recording() const;
352 
358  active_tape_index() const;
359 
364  bool
366  const typename Types<ad_type>::tape_index tape_index) const;
367 
393  void
395  const typename Types<ad_type>::tape_buffer_sizes obufsize = 64 * 1024 *
396  1024,
397  const typename Types<ad_type>::tape_buffer_sizes lbufsize = 64 * 1024 *
398  1024,
399  const typename Types<ad_type>::tape_buffer_sizes vbufsize = 64 * 1024 *
400  1024,
401  const typename Types<ad_type>::tape_buffer_sizes tbufsize = 64 * 1024 *
402  1024);
403 
449  bool
451  const typename Types<ad_type>::tape_index tape_index,
452  const bool overwrite_tape = false,
453  const bool keep_independent_values = true);
454 
466  void
467  stop_recording_operations(const bool write_tapes_to_file = false);
468 
478  void
480  const typename Types<ad_type>::tape_index tape_index);
481 
521  bool
523  const typename Types<ad_type>::tape_index tape_index) const;
524 
564  bool
566 
574  void
576 
579  protected:
592 
600 
617  void
618  activate_tape(const typename Types<ad_type>::tape_index tape_index,
619  const bool read_mode);
620 
635  mutable std::vector<scalar_type> independent_variable_values;
636 
645  mutable std::vector<ad_type> independent_variables;
646 
652 
657  mutable std::vector<bool> registered_marked_independent_variables;
658 
664  void
666 
673  void
674  set_sensitivity_value(const unsigned int index, const scalar_type &value);
675 
694  void
695  mark_independent_variable(const unsigned int index, ad_type &out) const;
696 
707  void
709 
719  void
720  initialize_non_sensitive_independent_variable(const unsigned int index,
721  ad_type &out) const;
722 
727  unsigned int
729 
749  std::vector<ad_type> dependent_variables;
750 
755 
762  void
763  reset_registered_dependent_variables(const bool flag = false);
764 
768  unsigned int
770 
782  void
783  register_dependent_variable(const unsigned int index,
784  const ad_type & func);
785 
788  }; // class HelperBase
789 
790 
791 
833  template <enum AD::NumberTypes ADNumberTypeCode,
834  typename ScalarType = double>
835  class CellLevelBase : public HelperBase<ADNumberTypeCode, ScalarType>
836  {
837  public:
842  using scalar_type =
844 
849  using ad_type =
851 
871  CellLevelBase(const unsigned int n_independent_variables,
872  const unsigned int n_dependent_variables);
873 
877  virtual ~CellLevelBase() = default;
878 
903  void
904  register_dof_values(const std::vector<scalar_type> &dof_values);
905 
924  template <typename VectorType>
925  void
927  const VectorType & values,
928  const std::vector<::types::global_dof_index> &local_dof_indices);
929 
950  const std::vector<ad_type> &
951  get_sensitive_dof_values() const;
952 
979  void
980  set_dof_values(const std::vector<scalar_type> &dof_values);
981 
999  template <typename VectorType>
1000  void
1002  const VectorType & values,
1003  const std::vector<::types::global_dof_index> &local_dof_indices);
1004 
1026  virtual void
1028 
1046  virtual void
1048 
1051  }; // class CellLevelBase
1052 
1053 
1054 
1215  template <enum AD::NumberTypes ADNumberTypeCode,
1216  typename ScalarType = double>
1217  class EnergyFunctional : public CellLevelBase<ADNumberTypeCode, ScalarType>
1218  {
1219  public:
1224  using scalar_type =
1226 
1231  using ad_type =
1233 
1254  EnergyFunctional(const unsigned int n_independent_variables);
1255 
1259  virtual ~EnergyFunctional() = default;
1260 
1281  void
1282  register_energy_functional(const ad_type &energy);
1283 
1297  scalar_type
1298  compute_energy() const;
1299 
1319  void
1320  compute_residual(Vector<scalar_type> &residual) const override;
1321 
1342  virtual void
1344  FullMatrix<scalar_type> &linearization) const override;
1345 
1348  }; // class EnergyFunctional
1349 
1350 
1351 
1526  template <enum AD::NumberTypes ADNumberTypeCode,
1527  typename ScalarType = double>
1529  : public CellLevelBase<ADNumberTypeCode, ScalarType>
1530  {
1531  public:
1536  using scalar_type =
1538 
1543  using ad_type =
1545 
1566  const unsigned int n_dependent_variables);
1567 
1571  virtual ~ResidualLinearization() = default;
1572 
1593  void
1594  register_residual_vector(const std::vector<ad_type> &residual);
1595 
1612  virtual void
1613  compute_residual(Vector<scalar_type> &residual) const override;
1614 
1632  virtual void
1634  FullMatrix<scalar_type> &linearization) const override;
1635 
1638  }; // class ResidualLinearization
1639 
1640 
1641 
1642  namespace internal
1643  {
1648  template <int dim, typename ExtractorType>
1649  struct Extractor;
1650 
1651 
1657  template <int dim>
1658  struct Extractor<dim, FEValuesExtractors::Scalar>
1659  {
1663  static const unsigned int n_components = 1;
1664 
1668  static const unsigned int rank = 0;
1669 
1673  template <typename NumberType>
1675 
1676  static_assert(
1678  "The number of components doesn't match that of the corresponding tensor type.");
1679  static_assert(
1680  rank == tensor_type<double>::rank,
1681  "The rank doesn't match that of the corresponding tensor type.");
1682 
1686  // Note: FEValuesViews::Scalar::tensor_type is double, so we can't
1687  // use it (FEValuesViews) in this context.
1688  // In fact, sadly, all FEValuesViews objects expect doubles as value
1689  // types.
1690  template <typename NumberType>
1692 
1696  template <typename NumberType>
1698 
1702  static inline unsigned int
1704  {
1705  return extractor.component;
1706  }
1707 
1715  static bool
1716  symmetric_component(const unsigned int unrolled_index)
1717  {
1718  (void)unrolled_index;
1719  return false;
1720  }
1721 
1729  template <typename IndexType = unsigned int, int rank_in>
1730  static IndexType
1732  const unsigned int column_offset)
1733  {
1734  Assert(column_offset <= rank_in, ExcInternalError());
1735  (void)table_indices;
1736  (void)column_offset;
1737  return 0;
1738  }
1739  };
1740 
1741 
1747  template <int dim>
1749  {
1753  static const unsigned int n_components = dim;
1754 
1758  static const unsigned int rank = 1;
1759 
1763  template <typename NumberType>
1765 
1766  static_assert(
1768  "The number of components doesn't match that of the corresponding tensor type.");
1769  static_assert(
1770  rank == tensor_type<double>::rank,
1771  "The rank doesn't match that of the corresponding tensor type.");
1772 
1776  template <typename NumberType>
1778 
1782  template <typename NumberType>
1784 
1788  static inline unsigned int
1790  {
1791  return extractor.first_vector_component;
1792  }
1793 
1802  static bool
1803  symmetric_component(const unsigned int unrolled_index)
1804  {
1805  (void)unrolled_index;
1806  return false;
1807  }
1808 
1813  template <int rank_in>
1814  static TableIndices<rank>
1816  const unsigned int column_offset)
1817  {
1818  Assert(0 + column_offset < rank_in, ExcInternalError());
1819  return TableIndices<rank>(table_indices[column_offset]);
1820  }
1821 
1837  template <typename IndexType = unsigned int, int rank_in>
1838  static IndexType
1840  const unsigned int column_offset)
1841  {
1842  static_assert(
1843  rank_in >= rank,
1844  "Cannot extract more table indices than the input table has!");
1845  using TensorType = tensor_type<double>;
1846  return TensorType::component_to_unrolled_index(
1847  table_index_view(table_indices, column_offset));
1848  }
1849  };
1850 
1851 
1857  template <int dim>
1859  {
1863  static const unsigned int n_components =
1865 
1869  static const unsigned int rank = 1;
1870 
1874  template <typename NumberType>
1876 
1880  template <typename NumberType>
1882 
1886  template <typename NumberType>
1888 
1892  static inline unsigned int
1894  {
1895  return extractor.first_tensor_component;
1896  }
1897 
1905  static bool
1906  symmetric_component(const unsigned int unrolled_index)
1907  {
1908  (void)unrolled_index;
1909  return false;
1910  }
1911 
1916  template <int rank_in>
1917  static TableIndices<rank>
1919  const unsigned int column_offset)
1920  {
1921  Assert(column_offset < rank_in, ExcInternalError());
1922  return TableIndices<rank>(table_indices[column_offset]);
1923  }
1924 
1940  template <typename IndexType = unsigned int, int rank_in>
1941  static IndexType
1943  const unsigned int column_offset)
1944  {
1945  static_assert(
1946  rank_in >= rank,
1947  "Cannot extract more table indices than the input table has!");
1948  using TensorType = tensor_type<double>;
1949  return TensorType::component_to_unrolled_index(
1950  table_index_view(table_indices, column_offset));
1951  }
1952  };
1953 
1954 
1960  template <int dim>
1962  {
1966  static const unsigned int n_components =
1968 
1972  static const unsigned int rank = Tensor<2, dim>::rank;
1973 
1977  template <typename NumberType>
1979 
1983  template <typename NumberType>
1985 
1989  template <typename NumberType>
1991 
1995  static inline unsigned int
1997  {
1998  return extractor.first_tensor_component;
1999  }
2000 
2008  static bool
2009  symmetric_component(const unsigned int unrolled_index)
2010  {
2011  (void)unrolled_index;
2012  return false;
2013  }
2014 
2019  template <int rank_in>
2020  static TableIndices<rank>
2022  const unsigned int column_offset)
2023  {
2024  Assert(column_offset < rank_in, ExcInternalError());
2025  Assert(column_offset + 1 < rank_in, ExcInternalError());
2026  return TableIndices<rank>(table_indices[column_offset],
2027  table_indices[column_offset + 1]);
2028  }
2029 
2045  template <typename IndexType = unsigned int, int rank_in>
2046  static IndexType
2048  const unsigned int column_offset)
2049  {
2050  static_assert(
2051  rank_in >= rank,
2052  "Cannot extract more table indices than the input table has!");
2053  using TensorType = tensor_type<double>;
2054  return TensorType::component_to_unrolled_index(
2055  table_index_view(table_indices, column_offset));
2056  }
2057  };
2058 
2059 
2065  template <int dim>
2067  {
2071  static const unsigned int n_components =
2073 
2077  static const unsigned int rank = SymmetricTensor<2, dim>::rank;
2078 
2082  template <typename NumberType>
2084 
2088  template <typename NumberType>
2090 
2094  template <typename NumberType>
2096 
2100  static inline unsigned int
2102  {
2103  return extractor.first_tensor_component;
2104  }
2105 
2114  static bool
2115  symmetric_component(const unsigned int unrolled_index)
2116  {
2117  const TableIndices<2> table_indices =
2119  return table_indices[0] != table_indices[1];
2120  }
2121 
2126  template <int rank_in>
2127  static TableIndices<rank>
2129  const unsigned int column_offset)
2130  {
2131  Assert(column_offset < rank_in, ExcInternalError());
2132  Assert(column_offset + 1 < rank_in, ExcInternalError());
2133  return TableIndices<rank>(table_indices[column_offset],
2134  table_indices[column_offset + 1]);
2135  }
2136 
2152  template <typename IndexType = unsigned int, int rank_in>
2153  static IndexType
2155  const unsigned int column_offset)
2156  {
2157  static_assert(
2158  rank_in >= rank,
2159  "Cannot extract more table indices than the input table has!");
2160  using TensorType = tensor_type<double>;
2161  return TensorType::component_to_unrolled_index(
2162  table_index_view(table_indices, column_offset));
2163  }
2164  };
2165 
2166 
2172  template <int dim, typename NumberType, typename ExtractorType>
2174  {
2179  using type =
2180  typename Extractor<dim,
2181  ExtractorType>::template tensor_type<NumberType>;
2182  };
2183 
2184 
2194  template <typename ExtractorType_Row, typename ExtractorType_Col>
2196  {
2208  template <int rank, int dim, typename NumberType>
2210  };
2211 
2212 
2221  template <>
2224  {
2232  template <int /*rank*/, int dim, typename NumberType>
2233  using type = SymmetricTensor<2 /*rank*/, dim, NumberType>;
2234  };
2235 
2236 
2245  template <>
2248  {
2257  template <int /*rank*/, int dim, typename NumberType>
2258  using type = SymmetricTensor<2 /*rank*/, dim, NumberType>;
2259  };
2260 
2261 
2269  template <>
2272  {
2281  template <int /*rank*/, int dim, typename NumberType>
2282  using type = SymmetricTensor<4 /*rank*/, dim, NumberType>;
2283  };
2284 
2285 
2295  template <int dim,
2296  typename NumberType,
2297  typename ExtractorType_Row,
2298  typename ExtractorType_Col>
2300  {
2306 
2313  using type =
2315  template type<rank, dim, NumberType>;
2316  };
2317 
2318 
2323  template <int dim, typename NumberType, typename ExtractorType_Field>
2326 
2327 
2337  template <int dim,
2338  typename NumberType,
2339  typename ExtractorType_Field,
2340  typename ExtractorType_Derivative>
2342  NumberType,
2343  ExtractorType_Field,
2344  ExtractorType_Derivative>;
2345 
2346 
2352  template <int dim,
2353  typename IndexType = unsigned int,
2354  typename ExtractorType>
2355  std::vector<IndexType>
2356  extract_field_component_indices(const ExtractorType &extractor,
2357  const bool ignore_symmetries = true)
2358  {
2359  (void)ignore_symmetries;
2360  const IndexType n_components =
2362  const IndexType comp_first =
2364  std::vector<IndexType> indices(n_components);
2365  std::iota(indices.begin(), indices.end(), comp_first);
2366  return indices;
2367  }
2368 
2369 
2378  template <int dim, typename IndexType = unsigned int>
2379  std::vector<IndexType>
2381  const FEValuesExtractors::SymmetricTensor<2> &extractor_symm_tensor,
2382  const bool ignore_symmetries = true)
2383  {
2384  using ExtractorType = FEValuesExtractors::SymmetricTensor<2>;
2385  const IndexType n_components =
2387  const IndexType comp_first =
2389  extractor_symm_tensor);
2390 
2391  if (ignore_symmetries == true)
2392  {
2393  std::vector<IndexType> indices(n_components);
2394  std::iota(indices.begin(), indices.end(), comp_first);
2395  return indices;
2396  }
2397  else
2398  {
2399  // First get all of the indices of the non-symmetric tensor
2400  const FEValuesExtractors::Tensor<2> extractor_tensor(
2401  extractor_symm_tensor.first_tensor_component);
2402  std::vector<IndexType> indices =
2403  extract_field_component_indices<dim>(extractor_tensor, true);
2404 
2405  // Then we overwrite any illegal entries with the equivalent indices
2406  // from the symmetric tensor
2407  for (unsigned int i = 0; i < indices.size(); ++i)
2408  {
2409  // The indices stored in the vector start with the extractor's
2410  // first_component_index. We need to account for this when
2411  // retrieving the tensor (local) index.
2412  const IndexType local_index_i = indices[i] - comp_first;
2413  if (local_index_i >= n_components)
2414  {
2415  const TableIndices<2> ti_tensor =
2417  local_index_i);
2418  const IndexType sti_new_index =
2420  ti_tensor);
2421  indices[i] = comp_first + sti_new_index;
2422  }
2423  }
2424 
2425  return indices;
2426  }
2427  }
2428 
2429 
2434  template <typename TensorType, typename NumberType>
2435  inline void
2436  set_tensor_entry(TensorType & t,
2437  const unsigned int unrolled_index,
2438  const NumberType & value)
2439  {
2440  // Where possible, set values using TableIndices
2441  AssertIndexRange(unrolled_index, t.n_independent_components);
2442  t[TensorType::unrolled_to_component_indices(unrolled_index)] = value;
2443  }
2444 
2445 
2450  template <int dim, typename NumberType>
2451  inline void
2453  const unsigned int unrolled_index,
2454  const NumberType & value)
2455  {
2456  AssertIndexRange(unrolled_index, 1);
2457  (void)unrolled_index;
2458  t = value;
2459  }
2460 
2461 
2467  template <typename NumberType>
2468  inline void
2470  const unsigned int unrolled_index,
2471  const NumberType & value)
2472  {
2473  AssertIndexRange(unrolled_index, 1);
2474  (void)unrolled_index;
2475  t = value;
2476  }
2477 
2478 
2484  template <int dim, typename NumberType>
2485  inline void
2487  const unsigned int unrolled_index_row,
2488  const unsigned int unrolled_index_col,
2489  const NumberType & value)
2490  {
2491  // Fourth order symmetric tensors require a specialized interface
2492  // to extract values.
2493  using SubTensorType = SymmetricTensor<2, dim, NumberType>;
2494  AssertIndexRange(unrolled_index_row,
2495  SubTensorType::n_independent_components);
2496  AssertIndexRange(unrolled_index_col,
2497  SubTensorType::n_independent_components);
2498  const TableIndices<2> indices_row =
2499  SubTensorType::unrolled_to_component_indices(unrolled_index_row);
2500  const TableIndices<2> indices_col =
2501  SubTensorType::unrolled_to_component_indices(unrolled_index_col);
2502  t[indices_row[0]][indices_row[1]][indices_col[0]][indices_col[1]] =
2503  value;
2504  }
2505 
2506 
2511  template <int rank,
2512  int dim,
2513  typename NumberType,
2514  template <int, int, typename>
2515  class TensorType>
2516  inline NumberType
2517  get_tensor_entry(const TensorType<rank, dim, NumberType> &t,
2518  const unsigned int unrolled_index)
2519  {
2520  // Where possible, get values using TableIndices
2521  AssertIndexRange(unrolled_index, t.n_independent_components);
2522  return t[TensorType<rank, dim, NumberType>::
2523  unrolled_to_component_indices(unrolled_index)];
2524  }
2525 
2526 
2531  template <int dim,
2532  typename NumberType,
2533  template <int, int, typename>
2534  class TensorType>
2535  inline NumberType
2536  get_tensor_entry(const TensorType<0, dim, NumberType> &t,
2537  const unsigned int unrolled_index)
2538  {
2539  AssertIndexRange(unrolled_index, 1);
2540  (void)unrolled_index;
2541  return t;
2542  }
2543 
2544 
2550  template <typename NumberType>
2551  inline const NumberType &
2552  get_tensor_entry(const NumberType &t, const unsigned int unrolled_index)
2553  {
2554  AssertIndexRange(unrolled_index, 1);
2555  (void)unrolled_index;
2556  return t;
2557  }
2558 
2559 
2564  template <int rank,
2565  int dim,
2566  typename NumberType,
2567  template <int, int, typename>
2568  class TensorType>
2569  inline NumberType &
2570  get_tensor_entry(TensorType<rank, dim, NumberType> &t,
2571  const unsigned int unrolled_index)
2572  {
2573  // Where possible, get values using TableIndices
2574  AssertIndexRange(unrolled_index, t.n_independent_components);
2575  return t[TensorType<rank, dim, NumberType>::
2576  unrolled_to_component_indices(unrolled_index)];
2577  }
2578 
2579 
2584  template <int dim,
2585  typename NumberType,
2586  template <int, int, typename>
2587  class TensorType>
2588  NumberType &
2589  get_tensor_entry(TensorType<0, dim, NumberType> &t,
2590  const unsigned int index)
2591  {
2592  AssertIndexRange(index, 1);
2593  (void)index;
2594  return t;
2595  }
2596 
2597 
2603  template <typename NumberType>
2604  inline NumberType &
2605  get_tensor_entry(NumberType &t, const unsigned int index)
2606  {
2607  AssertIndexRange(index, 1);
2608  (void)index;
2609  return t;
2610  }
2611 
2612  } // namespace internal
2613 
2614 
2615 
2637  template <int dim,
2638  enum AD::NumberTypes ADNumberTypeCode,
2639  typename ScalarType = double>
2641  : public HelperBase<ADNumberTypeCode, ScalarType>
2642  {
2643  public:
2648  static constexpr unsigned int dimension = dim;
2649 
2654  using scalar_type =
2656 
2661  using ad_type =
2663 
2684  const unsigned int n_dependent_variables);
2685 
2689  virtual ~PointLevelFunctionsBase() = default;
2690 
2733  virtual void
2734  reset(const unsigned int n_independent_variables =
2736  const unsigned int n_dependent_variables =
2738  const bool clear_registered_tapes = true) override;
2739 
2755  void
2756  register_independent_variables(const std::vector<scalar_type> &values);
2757 
2786  template <typename ValueType, typename ExtractorType>
2787  void
2788  register_independent_variable(const ValueType & value,
2789  const ExtractorType &extractor);
2790 
2811  const std::vector<ad_type> &
2812  get_sensitive_variables() const;
2813 
2814  /*
2815  * Extract a subset of the independent variables as represented by
2816  * auto-differentiable numbers. These are the independent
2817  * variables @f$\mathbf{A} \subset \mathbf{X}@f$ at which the dependent values
2818  * are evaluated and differentiated.
2819  *
2820  * This function indicates to the AD library that implements the
2821  * auto-differentiable number type that operations performed on these
2822  * numbers are to be tracked so they are considered "sensitive"
2823  * variables. This is, therefore, the set of variables with which one
2824  * would then perform computations, and based on which one can then
2825  * extract both the value of the function and its derivatives with the
2826  * member functions below. The values of the components of the returned
2827  * object are initialized to the values set with
2828  * register_independent_variable().
2829  *
2830  * @param[in] extractor An extractor associated with the input field
2831  * variables. This effectively defines which components of the global set
2832  * of independent variables this field is associated with.
2833  * @return An object of auto-differentiable type numbers. The return type is
2834  * based on the input extractor, and will be either a scalar,
2835  * Tensor<1,dim>, Tensor<2,dim>, or SymmetricTensor<2,dim>.
2836  *
2837  * @note For taped AD numbers, this operation is only valid in recording mode.
2838  */
2839  template <typename ExtractorType>
2840  typename internal::Extractor<dim,
2841  ExtractorType>::template tensor_type<ad_type>
2842  get_sensitive_variables(const ExtractorType &extractor) const;
2843 
2868  void
2869  set_independent_variables(const std::vector<scalar_type> &values);
2870 
2897  template <typename ValueType, typename ExtractorType>
2898  void
2899  set_independent_variable(const ValueType & value,
2900  const ExtractorType &extractor);
2901 
2904  protected:
2920  void
2921  set_sensitivity_value(const unsigned int index,
2922  const bool symmetric_component,
2923  const scalar_type &value);
2924 
2930  bool
2931  is_symmetric_independent_variable(const unsigned int index) const;
2932 
2937  unsigned int
2939 
2942  private:
2953 
2956  }; // class PointLevelFunctionsBase
2957 
2958 
2959 
3114  template <int dim,
3115  enum AD::NumberTypes ADNumberTypeCode,
3116  typename ScalarType = double>
3118  : public PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>
3119  {
3120  public:
3125  using scalar_type =
3127 
3132  using ad_type =
3134 
3149  ScalarFunction(const unsigned int n_independent_variables);
3150 
3154  virtual ~ScalarFunction() = default;
3155 
3174  void
3175  register_dependent_variable(const ad_type &func);
3176 
3184  scalar_type
3185  compute_value() const;
3186 
3199  void
3200  compute_gradient(Vector<scalar_type> &gradient) const;
3201 
3216  void
3217  compute_hessian(FullMatrix<scalar_type> &hessian) const;
3218 
3253  template <typename ExtractorType_Row>
3254  static typename internal::
3257  const ExtractorType_Row & extractor_row);
3258 
3297  template <typename ExtractorType_Row, typename ExtractorType_Col>
3298  static typename internal::ScalarFieldHessian<dim,
3299  scalar_type,
3300  ExtractorType_Row,
3301  ExtractorType_Col>::type
3303  const ExtractorType_Row & extractor_row,
3304  const ExtractorType_Col & extractor_col);
3305 
3321  const FullMatrix<scalar_type> & hessian,
3322  const FEValuesExtractors::Scalar &extractor_row,
3323  const FEValuesExtractors::Scalar &extractor_col);
3324 
3338  const FullMatrix<scalar_type> & hessian,
3339  const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
3340  const FEValuesExtractors::SymmetricTensor<2> &extractor_col);
3341 
3344  }; // class ScalarFunction
3345 
3346 
3347 
3505  template <int dim,
3506  enum AD::NumberTypes ADNumberTypeCode,
3507  typename ScalarType = double>
3509  : public PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>
3510  {
3511  public:
3516  using scalar_type =
3518 
3523  using ad_type =
3525 
3545  VectorFunction(const unsigned int n_independent_variables,
3546  const unsigned int n_dependent_variables);
3547 
3551  virtual ~VectorFunction() = default;
3552 
3572  void
3573  register_dependent_variables(const std::vector<ad_type> &funcs);
3574 
3593  template <typename ValueType, typename ExtractorType>
3594  void
3595  register_dependent_variable(const ValueType & funcs,
3596  const ExtractorType &extractor);
3597 
3606  void
3608 
3623  void
3624  compute_jacobian(FullMatrix<scalar_type> &jacobian) const;
3625 
3626 
3639  template <typename ExtractorType_Row>
3640  static typename internal::
3643  const ExtractorType_Row & extractor_row);
3644 
3692  template <typename ExtractorType_Row, typename ExtractorType_Col>
3693  static typename internal::VectorFieldJacobian<dim,
3694  scalar_type,
3695  ExtractorType_Row,
3696  ExtractorType_Col>::type
3698  const ExtractorType_Row & extractor_row,
3699  const ExtractorType_Col & extractor_col);
3700 
3718  const FullMatrix<scalar_type> & jacobian,
3719  const FEValuesExtractors::Scalar &extractor_row,
3720  const FEValuesExtractors::Scalar &extractor_col);
3721 
3737  const FullMatrix<scalar_type> & jacobian,
3738  const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
3739  const FEValuesExtractors::SymmetricTensor<2> &extractor_col);
3740 
3743  }; // class VectorFunction
3744 
3745 
3746  } // namespace AD
3747 } // namespace Differentiation
3748 
3749 
3750 /* ----------------- inline and template functions ----------------- */
3751 
3752 
3753 # ifndef DOXYGEN
3754 
3755 namespace Differentiation
3756 {
3757  namespace AD
3758  {
3759  /* ----------------- CellLevelBase ----------------- */
3760 
3761 
3762 
3763  template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
3764  template <typename VectorType>
3765  void
3767  const VectorType & values,
3768  const std::vector<::types::global_dof_index> &local_dof_indices)
3769  {
3770  // This is actually the same thing the set_dof_values() function,
3771  // in the sense that we simply populate our array of independent values
3772  // with a meaningful number. However, in this case we need to double check
3773  // that we're not registering these variables twice
3774  Assert(
3775  local_dof_indices.size() == this->n_independent_variables(),
3776  ExcMessage(
3777  "Degree of freedom index vector size does not match number of independent variables"));
3778  for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
3779  {
3780  Assert(this->registered_independent_variable_values[i] == false,
3781  ExcMessage("Independent variables already registered."));
3782  }
3783  set_dof_values(values, local_dof_indices);
3784  }
3785 
3786 
3787 
3788  template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
3789  template <typename VectorType>
3790  void
3792  const VectorType & values,
3793  const std::vector<::types::global_dof_index> &local_dof_indices)
3794  {
3795  Assert(local_dof_indices.size() == this->n_independent_variables(),
3796  ExcMessage(
3797  "Vector size does not match number of independent variables"));
3798  for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
3800  i, values[local_dof_indices[i]]);
3801  }
3802 
3803 
3804 
3805  /* ----------------- PointLevelFunctionsBase ----------------- */
3806 
3807 
3808 
3809  template <int dim,
3810  enum AD::NumberTypes ADNumberTypeCode,
3811  typename ScalarType>
3812  template <typename ValueType, typename ExtractorType>
3813  void
3815  register_independent_variable(const ValueType & value,
3816  const ExtractorType &extractor)
3817  {
3818  // This is actually the same thing as the set_independent_variable
3819  // function, in the sense that we simply populate our array of independent
3820  // values with a meaningful number. However, in this case we need to
3821  // double check that we're not registering these variables twice
3822 # ifdef DEBUG
3823  const std::vector<unsigned int> index_set(
3824  internal::extract_field_component_indices<dim>(extractor));
3825  for (const unsigned int index : index_set)
3826  {
3827  Assert(
3828  this->registered_independent_variable_values[index] == false,
3829  ExcMessage(
3830  "Overlapping indices for independent variables. "
3831  "One or more indices associated with the field that "
3832  "is being registered as an independent variable have "
3833  "already been associated with another field. This suggests "
3834  "that the component offsets used to construct their counterpart "
3835  "extractors are incompatible with one another. Make sure that "
3836  "the first component for each extractor properly takes into "
3837  "account the dimensionality of the preceding fields."));
3838  }
3839 # endif
3840  set_independent_variable(value, extractor);
3841  }
3842 
3843 
3844 
3845  template <int dim,
3846  enum AD::NumberTypes ADNumberTypeCode,
3847  typename ScalarType>
3848  template <typename ValueType, typename ExtractorType>
3849  void
3851  set_independent_variable(const ValueType & value,
3852  const ExtractorType &extractor)
3853  {
3854  const std::vector<unsigned int> index_set(
3855  internal::extract_field_component_indices<dim>(extractor));
3856  for (unsigned int i = 0; i < index_set.size(); ++i)
3857  {
3858  set_sensitivity_value(
3859  index_set[i],
3860  internal::Extractor<dim, ExtractorType>::symmetric_component(i),
3861  internal::get_tensor_entry(value, i));
3862  }
3863  }
3864 
3865 
3866 
3867  template <int dim,
3868  enum AD::NumberTypes ADNumberTypeCode,
3869  typename ScalarType>
3870  template <typename ExtractorType>
3871  typename internal::Extractor<dim, ExtractorType>::template tensor_type<
3872  typename PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>::
3873  ad_type>
3875  get_sensitive_variables(const ExtractorType &extractor) const
3876  {
3877  if (ADNumberTraits<ad_type>::is_taped == true)
3878  {
3879  Assert(this->active_tape_index() !=
3881  ExcMessage("Invalid tape index"));
3882  }
3883 
3884  // If necessary, finalize the internally stored vector of
3885  // AD numbers that represents the independent variables
3886  this->finalize_sensitive_independent_variables();
3887  Assert(this->independent_variables.size() ==
3888  this->n_independent_variables(),
3889  ExcDimensionMismatch(this->independent_variables.size(),
3890  this->n_independent_variables()));
3891 
3892  const std::vector<unsigned int> index_set(
3893  internal::extract_field_component_indices<dim>(extractor));
3894  typename internal::Extractor<dim,
3895  ExtractorType>::template tensor_type<ad_type>
3896  out;
3897 
3898  for (unsigned int i = 0; i < index_set.size(); ++i)
3899  {
3900  const unsigned int index = index_set[i];
3901  Assert(index < this->n_independent_variables(), ExcInternalError());
3902  Assert(this->registered_independent_variable_values[index] == true,
3903  ExcInternalError());
3904  internal::get_tensor_entry(out, i) =
3905  this->independent_variables[index];
3906  }
3907 
3908  return out;
3909  }
3910 
3911 
3912 
3913  /* ----------------- ScalarFunction ----------------- */
3914 
3915 
3916 
3917  template <int dim,
3918  enum AD::NumberTypes ADNumberTypeCode,
3919  typename ScalarType>
3920  template <typename ExtractorType_Row>
3921  typename internal::ScalarFieldGradient<
3922  dim,
3924  ExtractorType_Row>::type
3927  const ExtractorType_Row & extractor_row)
3928  {
3929  // NOTE: The order of components must be consistently defined throughout
3930  // this class.
3931  typename internal::
3933 
3934  // Get indexsets for the subblock from which we wish to extract the
3935  // gradient values
3936  const std::vector<unsigned int> row_index_set(
3937  internal::extract_field_component_indices<dim>(extractor_row));
3938  Assert(out.n_independent_components == row_index_set.size(),
3939  ExcMessage("Not all tensor components have been extracted!"));
3940  for (unsigned int r = 0; r < row_index_set.size(); ++r)
3941  internal::set_tensor_entry(out, r, gradient[row_index_set[r]]);
3942 
3943  return out;
3944  }
3945 
3946 
3947 
3948  template <int dim,
3949  enum AD::NumberTypes ADNumberTypeCode,
3950  typename ScalarType>
3951  template <typename ExtractorType_Row, typename ExtractorType_Col>
3952  typename internal::ScalarFieldHessian<
3953  dim,
3955  ExtractorType_Row,
3956  ExtractorType_Col>::type
3959  const ExtractorType_Row & extractor_row,
3960  const ExtractorType_Col & extractor_col)
3961  {
3962  using InternalHessian = internal::ScalarFieldHessian<dim,
3963  scalar_type,
3964  ExtractorType_Row,
3965  ExtractorType_Col>;
3966  using InternalExtractorRow = internal::Extractor<dim, ExtractorType_Row>;
3967  using InternalExtractorCol = internal::Extractor<dim, ExtractorType_Col>;
3968  using HessianType = typename InternalHessian::type;
3969 
3970  // NOTE: The order of components must be consistently defined throughout
3971  // this class.
3972  HessianType out;
3973 
3974  // Get indexsets for the subblocks from which we wish to extract the
3975  // Hessian values
3976  // NOTE: Here we have to do some clever accounting when the
3977  // one extractor is a symmetric Tensor and the other is not, e.g.
3978  // <SymmTensor,Vector>. In this scenario the return type is a
3979  // non-symmetric Tensor<3,dim> but we have to fetch information from a
3980  // SymmTensor row/column that has too few entries to fill the output
3981  // tensor. So we must duplicate the relevant entries in the row/column
3982  // indexset to fetch off-diagonal components that are otherwise
3983  // non-existent in a SymmTensor.
3984  const std::vector<unsigned int> row_index_set(
3985  internal::extract_field_component_indices<dim>(
3986  extractor_row, false /*ignore_symmetries*/));
3987  const std::vector<unsigned int> col_index_set(
3988  internal::extract_field_component_indices<dim>(
3989  extractor_col, false /*ignore_symmetries*/));
3990 
3991  for (unsigned int index = 0;
3992  index < HessianType::n_independent_components;
3993  ++index)
3994  {
3995  const TableIndices<HessianType::rank> ti_out =
3996  HessianType::unrolled_to_component_indices(index);
3997  const unsigned int r =
3998  InternalExtractorRow::local_component(ti_out, 0);
3999  const unsigned int c =
4000  InternalExtractorCol::local_component(ti_out,
4001  InternalExtractorRow::rank);
4002 
4004  out, index, hessian[row_index_set[r]][col_index_set[c]]);
4005  }
4006 
4007  return out;
4008  }
4009 
4010 
4011 
4012  /* ----------------- VectorFunction ----------------- */
4013 
4014 
4015 
4016  template <int dim,
4017  enum AD::NumberTypes ADNumberTypeCode,
4018  typename ScalarType>
4019  template <typename ValueType, typename ExtractorType>
4020  void
4022  register_dependent_variable(const ValueType & funcs,
4023  const ExtractorType &extractor)
4024  {
4025  const std::vector<unsigned int> index_set(
4026  internal::extract_field_component_indices<dim>(extractor));
4027  for (unsigned int i = 0; i < index_set.size(); ++i)
4028  {
4029  Assert(this->registered_marked_dependent_variables[index_set[i]] ==
4030  false,
4031  ExcMessage("Overlapping indices for dependent variables."));
4033  index_set[i], internal::get_tensor_entry(funcs, i));
4034  }
4035  }
4036 
4037 
4038 
4039  template <int dim,
4040  enum AD::NumberTypes ADNumberTypeCode,
4041  typename ScalarType>
4042  template <typename ExtractorType_Row>
4043  typename internal::VectorFieldValue<
4044  dim,
4046  ExtractorType_Row>::type
4048  const Vector<scalar_type> &values,
4049  const ExtractorType_Row & extractor_row)
4050  {
4051  // NOTE: The order of components must be consistently defined throughout
4052  // this class.
4053  typename internal::VectorFieldValue<dim, scalar_type, ExtractorType_Row>::
4054  type out;
4055 
4056  // Get indexsets for the subblock from which we wish to extract the
4057  // gradient values
4058  const std::vector<unsigned int> row_index_set(
4059  internal::extract_field_component_indices<dim>(extractor_row));
4060  Assert(out.n_independent_components == row_index_set.size(),
4061  ExcMessage("Not all tensor components have been extracted!"));
4062  for (unsigned int r = 0; r < row_index_set.size(); ++r)
4063  internal::set_tensor_entry(out, r, values[row_index_set[r]]);
4064 
4065  return out;
4066  }
4067 
4068 
4069 
4070  template <int dim,
4071  enum AD::NumberTypes ADNumberTypeCode,
4072  typename ScalarType>
4073  template <typename ExtractorType_Row, typename ExtractorType_Col>
4075  dim,
4077  ExtractorType_Row,
4078  ExtractorType_Col>::type
4081  const ExtractorType_Row & extractor_row,
4082  const ExtractorType_Col & extractor_col)
4083  {
4084  using InternalJacobian = internal::VectorFieldJacobian<dim,
4085  scalar_type,
4086  ExtractorType_Row,
4087  ExtractorType_Col>;
4088  using InternalExtractorRow = internal::Extractor<dim, ExtractorType_Row>;
4089  using InternalExtractorCol = internal::Extractor<dim, ExtractorType_Col>;
4090  using JacobianType = typename InternalJacobian::type;
4091 
4092  // NOTE: The order of components must be consistently defined throughout
4093  // this class.
4094  JacobianType out;
4095 
4096  // Get indexsets for the subblocks from which we wish to extract the
4097  // Hessian values.
4098  // NOTE: Here we have to do some clever accounting when the
4099  // one extractor is a symmetric Tensor and the other is not, e.g.
4100  // <SymmTensor,Vector>. In this scenario the return type is a
4101  // non-symmetric Tensor<3,dim> but we have to fetch information from a
4102  // SymmTensor row/column that has too few entries to fill the output
4103  // tensor. So we must duplicate the relevant entries in the row/column
4104  // indexset to fetch off-diagonal components that are otherwise
4105  // non-existent in a SymmTensor.
4106  const std::vector<unsigned int> row_index_set(
4107  internal::extract_field_component_indices<dim>(
4108  extractor_row, false /*ignore_symmetries*/));
4109  const std::vector<unsigned int> col_index_set(
4110  internal::extract_field_component_indices<dim>(
4111  extractor_col, false /*ignore_symmetries*/));
4112 
4113  for (unsigned int index = 0;
4114  index < JacobianType::n_independent_components;
4115  ++index)
4116  {
4117  const TableIndices<JacobianType::rank> ti_out =
4118  JacobianType::unrolled_to_component_indices(index);
4119  const unsigned int r =
4120  InternalExtractorRow::local_component(ti_out, 0);
4121  const unsigned int c =
4122  InternalExtractorCol::local_component(ti_out,
4123  InternalExtractorRow::rank);
4124 
4126  out, index, jacobian[row_index_set[r]][col_index_set[c]]);
4127  }
4128 
4129  return out;
4130  }
4131 
4132 
4133  } // namespace AD
4134 } // namespace Differentiation
4135 
4136 
4137 # endif // DOXYGEN
4138 
4139 
4141 
4142 #endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
4143 
4144 #endif // dealii_differentiation_ad_ad_helpers_h
void register_dof_values(const VectorType &values, const std::vector<::types::global_dof_index > &local_dof_indices)
const std::vector< ad_type > & get_sensitive_dof_values() const
Definition: ad_helpers.cc:744
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:850
void set_dof_values(const std::vector< scalar_type > &dof_values)
Definition: ad_helpers.cc:769
virtual void compute_residual(Vector< scalar_type > &residual) const =0
void register_dof_values(const std::vector< scalar_type > &dof_values)
Definition: ad_helpers.cc:721
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const =0
CellLevelBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:710
void set_dof_values(const VectorType &values, const std::vector<::types::global_dof_index > &local_dof_indices)
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:843
void register_energy_functional(const ad_type &energy)
Definition: ad_helpers.cc:802
EnergyFunctional(const unsigned int n_independent_variables)
Definition: ad_helpers.cc:793
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:1232
void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:868
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:933
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:1225
unsigned int n_registered_dependent_variables() const
Definition: ad_helpers.cc:250
TapedDrivers< ad_type, scalar_type > taped_driver
Definition: ad_helpers.h:591
void stop_recording_operations(const bool write_tapes_to_file=false)
Definition: ad_helpers.cc:643
std::vector< bool > registered_marked_independent_variables
Definition: ad_helpers.h:657
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::scalar_type scalar_type
Definition: ad_helpers.h:177
std::vector< bool > registered_independent_variable_values
Definition: ad_helpers.h:651
void print_tape_stats(const typename Types< ad_type >::tape_index tape_index, std::ostream &stream) const
Definition: ad_helpers.cc:373
std::size_t n_independent_variables() const
Definition: ad_helpers.cc:241
bool start_recording_operations(const typename Types< ad_type >::tape_index tape_index, const bool overwrite_tape=false, const bool keep_independent_values=true)
Definition: ad_helpers.cc:578
std::vector< scalar_type > independent_variable_values
Definition: ad_helpers.h:635
void mark_independent_variable(const unsigned int index, ad_type &out) const
Definition: ad_helpers.cc:142
bool is_registered_tape(const typename Types< ad_type >::tape_index tape_index) const
Definition: ad_helpers.cc:296
std::vector< bool > registered_marked_dependent_variables
Definition: ad_helpers.h:754
void reset_registered_dependent_variables(const bool flag=false)
Definition: ad_helpers.cc:94
void activate_tape(const typename Types< ad_type >::tape_index tape_index, const bool read_mode)
Definition: ad_helpers.cc:522
void set_tape_buffer_sizes(const typename Types< ad_type >::tape_buffer_sizes obufsize=64 *1024 *1024, const typename Types< ad_type >::tape_buffer_sizes lbufsize=64 *1024 *1024, const typename Types< ad_type >::tape_buffer_sizes vbufsize=64 *1024 *1024, const typename Types< ad_type >::tape_buffer_sizes tbufsize=64 *1024 *1024)
Definition: ad_helpers.cc:559
unsigned int n_registered_independent_variables() const
Definition: ad_helpers.cc:230
void print(std::ostream &stream) const
Definition: ad_helpers.cc:309
void set_sensitivity_value(const unsigned int index, const scalar_type &value)
Definition: ad_helpers.cc:105
bool active_tape_requires_retaping() const
Definition: ad_helpers.cc:497
std::vector< ad_type > independent_variables
Definition: ad_helpers.h:645
Types< ad_type >::tape_index active_tape_index() const
Definition: ad_helpers.cc:284
HelperBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:38
virtual void reset(const unsigned int n_independent_variables=::numbers::invalid_unsigned_int, const unsigned int n_dependent_variables=::numbers::invalid_unsigned_int, const bool clear_registered_tapes=true)
Definition: ad_helpers.cc:390
TapelessDrivers< ad_type, scalar_type > tapeless_driver
Definition: ad_helpers.h:599
bool recorded_tape_requires_retaping(const typename Types< ad_type >::tape_index tape_index) const
Definition: ad_helpers.cc:484
std::vector< ad_type > dependent_variables
Definition: ad_helpers.h:749
std::size_t n_dependent_variables() const
Definition: ad_helpers.cc:262
void print_values(std::ostream &stream) const
Definition: ad_helpers.cc:359
void register_dependent_variable(const unsigned int index, const ad_type &func)
Definition: ad_helpers.cc:679
void finalize_sensitive_independent_variables() const
Definition: ad_helpers.cc:181
static void configure_tapeless_mode(const unsigned int n_independent_variables, const bool ensure_persistent_setting=true)
Definition: ad_helpers.cc:448
void initialize_non_sensitive_independent_variable(const unsigned int index, ad_type &out) const
Definition: ad_helpers.cc:205
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::ad_type ad_type
Definition: ad_helpers.h:184
void activate_recorded_tape(const typename Types< ad_type >::tape_index tape_index)
Definition: ad_helpers.cc:474
void set_independent_variable(const ValueType &value, const ExtractorType &extractor)
bool is_symmetric_independent_variable(const unsigned int index) const
Definition: ad_helpers.cc:1192
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:2662
PointLevelFunctionsBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1155
static constexpr unsigned int dimension
Definition: ad_helpers.h:2648
void set_independent_variables(const std::vector< scalar_type > &values)
Definition: ad_helpers.cc:1296
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:2655
std::vector< bool > symmetric_independent_variables
Definition: ad_helpers.h:2952
unsigned int n_symmetric_independent_variables() const
Definition: ad_helpers.cc:1206
void register_independent_variables(const std::vector< scalar_type > &values)
Definition: ad_helpers.cc:1220
void set_sensitivity_value(const unsigned int index, const bool symmetric_component, const scalar_type &value)
Definition: ad_helpers.cc:1274
void register_independent_variable(const ValueType &value, const ExtractorType &extractor)
const std::vector< ad_type > & get_sensitive_variables() const
Definition: ad_helpers.cc:1246
virtual void reset(const unsigned int n_independent_variables=::numbers::invalid_unsigned_int, const unsigned int n_dependent_variables=::numbers::invalid_unsigned_int, const bool clear_registered_tapes=true) override
Definition: ad_helpers.cc:1168
internal::Extractor< dim, ExtractorType >::template tensor_type< ad_type > get_sensitive_variables(const ExtractorType &extractor) const
ResidualLinearization(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1006
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:1537
virtual void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:1032
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:1086
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:1544
void register_residual_vector(const std::vector< ad_type > &residual)
Definition: ad_helpers.cc:1018
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:3133
void register_dependent_variable(const ad_type &func)
Definition: ad_helpers.cc:1335
static internal::ScalarFieldGradient< dim, scalar_type, ExtractorType_Row >::type extract_gradient_component(const Vector< scalar_type > &gradient, const ExtractorType_Row &extractor_row)
ScalarFunction(const unsigned int n_independent_variables)
Definition: ad_helpers.cc:1321
static internal::ScalarFieldHessian< dim, scalar_type, ExtractorType_Row, ExtractorType_Col >::type extract_hessian_component(const FullMatrix< scalar_type > &hessian, const ExtractorType_Row &extractor_row, const ExtractorType_Col &extractor_col)
void compute_gradient(Vector< scalar_type > &gradient) const
Definition: ad_helpers.cc:1398
void compute_hessian(FullMatrix< scalar_type > &hessian) const
Definition: ad_helpers.cc:1472
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3126
typename HelperBase< ADNumberTypeCode, ScalarType >::ad_type ad_type
Definition: ad_helpers.h:3524
void compute_values(Vector< scalar_type > &values) const
Definition: ad_helpers.cc:1677
void register_dependent_variables(const std::vector< ad_type > &funcs)
Definition: ad_helpers.cc:1661
static internal::VectorFieldValue< dim, scalar_type, ExtractorType_Row >::type extract_value_component(const Vector< scalar_type > &values, const ExtractorType_Row &extractor_row)
static 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)
void compute_jacobian(FullMatrix< scalar_type > &jacobian) const
Definition: ad_helpers.cc:1733
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3517
void register_dependent_variable(const ValueType &funcs, const ExtractorType &extractor)
VectorFunction(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1646
static constexpr unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
Definition: tensor.h:503
static constexpr TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:461
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:462
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1509
#define AssertIndexRange(index, range)
Definition: exceptions.h:1768
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< IndexType > extract_field_component_indices(const ExtractorType &extractor, const bool ignore_symmetries=true)
Definition: ad_helpers.h:2356
ScalarFieldHessian< dim, NumberType, ExtractorType_Field, ExtractorType_Derivative > VectorFieldJacobian
Definition: ad_helpers.h:2344
NumberType get_tensor_entry(const TensorType< rank, dim, NumberType > &t, const unsigned int unrolled_index)
Definition: ad_helpers.h:2517
void set_tensor_entry(TensorType &t, const unsigned int unrolled_index, const NumberType &value)
Definition: ad_helpers.h:2436
ScalarFieldGradient< dim, NumberType, ExtractorType_Field > VectorFieldValue
Definition: ad_helpers.h:2325
void set_dof_values(const DoFCellAccessor< dim, spacedim, lda > &cell, const Vector< number > &local_values, OutputVector &values, const bool perform_check)
static const unsigned int invalid_unsigned_int
Definition: types.h:206
static const Types< ADNumberType >::tape_index invalid_tape_index
Definition: ad_drivers.h:122
unsigned int tape_buffer_sizes
Definition: ad_drivers.h:105
static unsigned int first_component(const FEValuesExtractors::Scalar &extractor)
Definition: ad_helpers.h:1703
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1731
static bool symmetric_component(const unsigned int unrolled_index)
Definition: ad_helpers.h:1716
static TableIndices< rank > table_index_view(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:2128
static unsigned int first_component(const FEValuesExtractors::SymmetricTensor< 2 > &extractor)
Definition: ad_helpers.h:2101
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:2154
static unsigned int first_component(const FEValuesExtractors::Tensor< 1 > &extractor)
Definition: ad_helpers.h:1893
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1942
static TableIndices< rank > table_index_view(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1918
static unsigned int first_component(const FEValuesExtractors::Tensor< 2 > &extractor)
Definition: ad_helpers.h:1996
static TableIndices< rank > table_index_view(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:2021
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:2047
static bool symmetric_component(const unsigned int unrolled_index)
Definition: ad_helpers.h:1803
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1839
static unsigned int first_component(const FEValuesExtractors::Vector &extractor)
Definition: ad_helpers.h:1789
static TableIndices< rank > table_index_view(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1815
typename Extractor< dim, ExtractorType >::template tensor_type< NumberType > type
Definition: ad_helpers.h:2181
typename HessianType< ExtractorType_Row, ExtractorType_Col >::template type< rank, dim, NumberType > type
Definition: ad_helpers.h:2315