Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
ad_helpers.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2017 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 
27 # include <deal.II/differentiation/ad/ad_drivers.h>
28 # include <deal.II/differentiation/ad/ad_number_traits.h>
29 # include <deal.II/differentiation/ad/adolc_number_types.h>
30 # include <deal.II/differentiation/ad/adolc_product_types.h>
31 # include <deal.II/differentiation/ad/sacado_number_types.h>
32 # include <deal.II/differentiation/ad/sacado_product_types.h>
33 
34 # include <deal.II/fe/fe_values_extractors.h>
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 # include <set>
44 
45 DEAL_II_NAMESPACE_OPEN
46 
47 namespace Differentiation
48 {
49  namespace AD
50  {
170  template <enum AD::NumberTypes ADNumberTypeCode,
171  typename ScalarType = double>
173  {
174  public:
179  using scalar_type =
181 
186  using ad_type =
188 
193 
208  HelperBase(const unsigned int n_independent_variables,
209  const unsigned int n_dependent_variables);
210 
214  virtual ~HelperBase() = default;
215 
217 
222 
227  std::size_t
228  n_independent_variables() const;
229 
234  std::size_t
235  n_dependent_variables() const;
236 
245  void
246  print(std::ostream &stream) const;
247 
254  void
255  print_values(std::ostream &stream) const;
256 
267  void
268  print_tape_stats(const typename Types<ad_type>::tape_index tape_index,
269  std::ostream & stream) const;
270 
272 
277 
296  static void
298  const bool ensure_persistent_setting = true);
299 
301 
306 
342  virtual void
343  reset(const unsigned int n_independent_variables =
345  const unsigned int n_dependent_variables =
347  const bool clear_registered_tapes = true);
348 
353  bool
354  is_recording() const;
355 
361  active_tape_index() const;
362 
367  bool
369  const typename Types<ad_type>::tape_index tape_index) const;
370 
396  void
398  const typename Types<ad_type>::tape_buffer_sizes obufsize = 64 * 1024 *
399  1024,
400  const typename Types<ad_type>::tape_buffer_sizes lbufsize = 64 * 1024 *
401  1024,
402  const typename Types<ad_type>::tape_buffer_sizes vbufsize = 64 * 1024 *
403  1024,
404  const typename Types<ad_type>::tape_buffer_sizes tbufsize = 64 * 1024 *
405  1024);
406 
452  bool
454  const typename Types<ad_type>::tape_index tape_index,
455  const bool overwrite_tape = false,
456  const bool keep_independent_values = true);
457 
469  void
470  stop_recording_operations(const bool write_tapes_to_file = false);
471 
481  void
483  const typename Types<ad_type>::tape_index tape_index);
484 
524  bool
526  const typename Types<ad_type>::tape_index tape_index) const;
527 
567  bool
569 
577  void
579 
581 
582  protected:
587 
595 
603 
620  void
621  activate_tape(const typename Types<ad_type>::tape_index tape_index,
622  const bool read_mode);
623 
625 
630 
638  mutable std::vector<scalar_type> independent_variable_values;
639 
648  mutable std::vector<ad_type> independent_variables;
649 
655 
660  mutable std::vector<bool> registered_marked_independent_variables;
661 
667  void
669 
676  void
677  set_sensitivity_value(const unsigned int index, const scalar_type &value);
678 
697  void
698  mark_independent_variable(const unsigned int index, ad_type &out) const;
699 
710  void
712 
722  void
723  initialize_non_sensitive_independent_variable(const unsigned int index,
724  ad_type &out) const;
725 
730  unsigned int
732 
734 
739 
752  std::vector<ad_type> dependent_variables;
753 
758 
765  void
766  reset_registered_dependent_variables(const bool flag = false);
767 
771  unsigned int
773 
785  void
786  register_dependent_variable(const unsigned int index,
787  const ad_type & func);
788 
790 
791  }; // class HelperBase
792 
793 
794 
838  template <enum AD::NumberTypes ADNumberTypeCode,
839  typename ScalarType = double>
840  class CellLevelBase : public HelperBase<ADNumberTypeCode, ScalarType>
841  {
842  public:
847  using scalar_type =
849 
854  using ad_type =
856 
861 
876  CellLevelBase(const unsigned int n_independent_variables,
877  const unsigned int n_dependent_variables);
878 
882  virtual ~CellLevelBase() = default;
883 
885 
890 
908  void
909  register_dof_values(const std::vector<scalar_type> &dof_values);
910 
929  template <typename VectorType>
930  void
932  const VectorType & values,
933  const std::vector<::types::global_dof_index> &local_dof_indices);
934 
955  const std::vector<ad_type> &
956  get_sensitive_dof_values() const;
957 
959 
964 
984  void
985  set_dof_values(const std::vector<scalar_type> &dof_values);
986 
1004  template <typename VectorType>
1005  void
1007  const VectorType & values,
1008  const std::vector<::types::global_dof_index> &local_dof_indices);
1009 
1011 
1016 
1031  virtual void
1032  compute_residual(Vector<scalar_type> &residual) const = 0;
1033 
1051  virtual void
1052  compute_linearization(FullMatrix<scalar_type> &linearization) const = 0;
1053 
1055 
1056  }; // class CellLevelBase
1057 
1058 
1059 
1222  template <enum AD::NumberTypes ADNumberTypeCode,
1223  typename ScalarType = double>
1224  class EnergyFunctional : public CellLevelBase<ADNumberTypeCode, ScalarType>
1225  {
1226  public:
1231  using scalar_type =
1233 
1238  using ad_type =
1240 
1245 
1261  EnergyFunctional(const unsigned int n_independent_variables);
1262 
1266  virtual ~EnergyFunctional() = default;
1267 
1269 
1274 
1288  void
1289  register_energy_functional(const ad_type &energy);
1290 
1304  scalar_type
1305  compute_energy() const;
1306 
1326  void
1327  compute_residual(Vector<scalar_type> &residual) const override;
1328 
1349  virtual void
1351  FullMatrix<scalar_type> &linearization) const override;
1352 
1354 
1355  }; // class EnergyFunctional
1356 
1357 
1358 
1537  template <enum AD::NumberTypes ADNumberTypeCode,
1538  typename ScalarType = double>
1540  : public CellLevelBase<ADNumberTypeCode, ScalarType>
1541  {
1542  public:
1547  using scalar_type =
1549 
1554  using ad_type =
1556 
1561 
1577  const unsigned int n_dependent_variables);
1578 
1582  virtual ~ResidualLinearization() = default;
1583 
1585 
1590 
1604  void
1605  register_residual_vector(const std::vector<ad_type> &residual);
1606 
1623  virtual void
1624  compute_residual(Vector<scalar_type> &residual) const override;
1625 
1643  virtual void
1645  FullMatrix<scalar_type> &linearization) const override;
1646 
1648 
1649  }; // class ResidualLinearization
1650 
1651 
1652 
1653  namespace internal
1654  {
1659  template <int dim, typename ExtractorType>
1660  struct Extractor;
1661 
1662 
1668  template <int dim>
1669  struct Extractor<dim, FEValuesExtractors::Scalar>
1670  {
1674  static const unsigned int n_components = 1;
1675 
1679  static const unsigned int rank = 0;
1680 
1684  template <typename NumberType>
1686 
1687  static_assert(
1689  "The number of components doesn't match that of the corresponding tensor type.");
1690  static_assert(
1691  rank == tensor_type<double>::rank,
1692  "The rank doesn't match that of the corresponding tensor type.");
1693 
1697  // Note: FEValuesViews::Scalar::tensor_type is double, so we can't
1698  // use it (FEValuesViews) in this context.
1699  // In fact, sadly, all FEValuesViews objects expect doubles as value
1700  // types.
1701  template <typename NumberType>
1703 
1707  template <typename NumberType>
1709 
1713  static inline unsigned int
1715  {
1716  return extractor.component;
1717  }
1718 
1726  static bool
1727  symmetric_component(const unsigned int unrolled_index)
1728  {
1729  (void)unrolled_index;
1730  return false;
1731  }
1732 
1740  template <typename IndexType = unsigned int, int rank_in>
1741  static IndexType
1743  const unsigned int column_offset)
1744  {
1745  Assert(column_offset <= rank_in, ExcInternalError());
1746  (void)table_indices;
1747  (void)column_offset;
1748  return 0;
1749  }
1750  };
1751 
1752 
1758  template <int dim>
1760  {
1764  static const unsigned int n_components = dim;
1765 
1769  static const unsigned int rank = 1;
1770 
1774  template <typename NumberType>
1776 
1777  static_assert(
1779  "The number of components doesn't match that of the corresponding tensor type.");
1780  static_assert(
1781  rank == tensor_type<double>::rank,
1782  "The rank doesn't match that of the corresponding tensor type.");
1783 
1787  template <typename NumberType>
1789 
1793  template <typename NumberType>
1795 
1799  static inline unsigned int
1801  {
1802  return extractor.first_vector_component;
1803  }
1804 
1813  static bool
1814  symmetric_component(const unsigned int unrolled_index)
1815  {
1816  (void)unrolled_index;
1817  return false;
1818  }
1819 
1824  template <int rank_in>
1825  static TableIndices<rank>
1827  const unsigned int column_offset)
1828  {
1829  Assert(0 + column_offset < rank_in, ExcInternalError());
1830  return TableIndices<rank>(table_indices[column_offset]);
1831  }
1832 
1848  template <typename IndexType = unsigned int, int rank_in>
1849  static IndexType
1851  const unsigned int column_offset)
1852  {
1853  static_assert(
1854  rank_in >= rank,
1855  "Cannot extract more table indices than the input table has!");
1856  using TensorType = tensor_type<double>;
1857  return TensorType::component_to_unrolled_index(
1858  table_index_view(table_indices, column_offset));
1859  }
1860  };
1861 
1862 
1868  template <int dim>
1870  {
1874  static const unsigned int n_components =
1876 
1880  static const unsigned int rank = 1;
1881 
1885  template <typename NumberType>
1887 
1891  template <typename NumberType>
1893 
1897  template <typename NumberType>
1899 
1903  static inline unsigned int
1905  {
1906  return extractor.first_tensor_component;
1907  }
1908 
1916  static bool
1917  symmetric_component(const unsigned int unrolled_index)
1918  {
1919  (void)unrolled_index;
1920  return false;
1921  }
1922 
1927  template <int rank_in>
1928  static TableIndices<rank>
1930  const unsigned int column_offset)
1931  {
1932  Assert(column_offset < rank_in, ExcInternalError());
1933  return TableIndices<rank>(table_indices[column_offset]);
1934  }
1935 
1951  template <typename IndexType = unsigned int, int rank_in>
1952  static IndexType
1954  const unsigned int column_offset)
1955  {
1956  static_assert(
1957  rank_in >= rank,
1958  "Cannot extract more table indices than the input table has!");
1959  using TensorType = tensor_type<double>;
1960  return TensorType::component_to_unrolled_index(
1961  table_index_view(table_indices, column_offset));
1962  }
1963  };
1964 
1965 
1971  template <int dim>
1973  {
1977  static const unsigned int n_components =
1979 
1983  static const unsigned int rank = Tensor<2, dim>::rank;
1984 
1988  template <typename NumberType>
1990 
1994  template <typename NumberType>
1996 
2000  template <typename NumberType>
2002 
2006  static inline unsigned int
2008  {
2009  return extractor.first_tensor_component;
2010  }
2011 
2019  static bool
2020  symmetric_component(const unsigned int unrolled_index)
2021  {
2022  (void)unrolled_index;
2023  return false;
2024  }
2025 
2030  template <int rank_in>
2031  static TableIndices<rank>
2033  const unsigned int column_offset)
2034  {
2035  Assert(column_offset < rank_in, ExcInternalError());
2036  Assert(column_offset + 1 < rank_in, ExcInternalError());
2037  return TableIndices<rank>(table_indices[column_offset],
2038  table_indices[column_offset + 1]);
2039  }
2040 
2056  template <typename IndexType = unsigned int, int rank_in>
2057  static IndexType
2059  const unsigned int column_offset)
2060  {
2061  static_assert(
2062  rank_in >= rank,
2063  "Cannot extract more table indices than the input table has!");
2064  using TensorType = tensor_type<double>;
2065  return TensorType::component_to_unrolled_index(
2066  table_index_view(table_indices, column_offset));
2067  }
2068  };
2069 
2070 
2076  template <int dim>
2078  {
2082  static const unsigned int n_components =
2084 
2088  static const unsigned int rank = SymmetricTensor<2, dim>::rank;
2089 
2093  template <typename NumberType>
2095 
2099  template <typename NumberType>
2101 
2105  template <typename NumberType>
2107 
2111  static inline unsigned int
2113  {
2114  return extractor.first_tensor_component;
2115  }
2116 
2125  static bool
2126  symmetric_component(const unsigned int unrolled_index)
2127  {
2128  const TableIndices<2> table_indices =
2130  return table_indices[0] != table_indices[1];
2131  }
2132 
2137  template <int rank_in>
2138  static TableIndices<rank>
2140  const unsigned int column_offset)
2141  {
2142  Assert(column_offset < rank_in, ExcInternalError());
2143  Assert(column_offset + 1 < rank_in, ExcInternalError());
2144  return TableIndices<rank>(table_indices[column_offset],
2145  table_indices[column_offset + 1]);
2146  }
2147 
2163  template <typename IndexType = unsigned int, int rank_in>
2164  static IndexType
2166  const unsigned int column_offset)
2167  {
2168  static_assert(
2169  rank_in >= rank,
2170  "Cannot extract more table indices than the input table has!");
2171  using TensorType = tensor_type<double>;
2172  return TensorType::component_to_unrolled_index(
2173  table_index_view(table_indices, column_offset));
2174  }
2175  };
2176 
2177 
2183  template <int dim, typename NumberType, typename ExtractorType>
2185  {
2190  using type =
2191  typename Extractor<dim,
2192  ExtractorType>::template tensor_type<NumberType>;
2193  };
2194 
2195 
2205  template <typename ExtractorType_Row, typename ExtractorType_Col>
2207  {
2219  template <int rank, int dim, typename NumberType>
2221  };
2222 
2223 
2232  template <>
2235  {
2243  template <int /*rank*/, int dim, typename NumberType>
2244  using type = SymmetricTensor<2 /*rank*/, dim, NumberType>;
2245  };
2246 
2247 
2256  template <>
2259  {
2268  template <int /*rank*/, int dim, typename NumberType>
2269  using type = SymmetricTensor<2 /*rank*/, dim, NumberType>;
2270  };
2271 
2272 
2280  template <>
2283  {
2292  template <int /*rank*/, int dim, typename NumberType>
2293  using type = SymmetricTensor<4 /*rank*/, dim, NumberType>;
2294  };
2295 
2296 
2306  template <int dim,
2307  typename NumberType,
2308  typename ExtractorType_Row,
2309  typename ExtractorType_Col>
2311  {
2317 
2324  using type =
2327  };
2328 
2329 
2334  template <int dim, typename NumberType, typename ExtractorType_Field>
2335  using VectorFieldValue =
2337 
2338 
2348  template <int dim,
2349  typename NumberType,
2350  typename ExtractorType_Field,
2351  typename ExtractorType_Derivative>
2353  NumberType,
2354  ExtractorType_Field,
2355  ExtractorType_Derivative>;
2356 
2357 
2363  template <int dim,
2364  typename IndexType = unsigned int,
2365  typename ExtractorType>
2366  std::vector<IndexType>
2367  extract_field_component_indices(const ExtractorType &extractor,
2368  const bool ignore_symmetries = true)
2369  {
2370  (void)ignore_symmetries;
2371  const IndexType n_components =
2373  const IndexType comp_first =
2375  std::vector<IndexType> indices(n_components);
2376  std::iota(indices.begin(), indices.end(), comp_first);
2377  return indices;
2378  }
2379 
2380 
2389  template <int dim, typename IndexType = unsigned int>
2390  std::vector<IndexType>
2392  const FEValuesExtractors::SymmetricTensor<2> &extractor_symm_tensor,
2393  const bool ignore_symmetries = true)
2394  {
2395  using ExtractorType = FEValuesExtractors::SymmetricTensor<2>;
2396  const IndexType n_components =
2398  if (ignore_symmetries == true)
2399  {
2400  const IndexType comp_first =
2402  extractor_symm_tensor);
2403  std::vector<IndexType> indices(n_components);
2404  std::iota(indices.begin(), indices.end(), comp_first);
2405  return indices;
2406  }
2407  else
2408  {
2409  // First get all of the indices of the non-symmetric tensor
2410  const FEValuesExtractors::Tensor<2> extractor_tensor(
2411  extractor_symm_tensor.first_tensor_component);
2412  std::vector<IndexType> indices =
2413  extract_field_component_indices<dim>(extractor_tensor, true);
2414 
2415  // Then we overwrite any illegal entries with the equivalent indices
2416  // from the symmetric tensor
2417  for (unsigned int i = 0; i < indices.size(); ++i)
2418  {
2419  if (indices[i] >= n_components)
2420  {
2421  const TableIndices<2> ti_tensor =
2423  const IndexType sti_new_index =
2425  ti_tensor);
2426  indices[i] = sti_new_index;
2427  }
2428  }
2429 
2430  return indices;
2431  }
2432  }
2433 
2434 
2439  template <typename TensorType, typename NumberType>
2440  inline void
2441  set_tensor_entry(TensorType & t,
2442  const unsigned int &unrolled_index,
2443  const NumberType & value)
2444  {
2445  // Where possible, set values using TableIndices
2446  Assert(unrolled_index < t.n_independent_components,
2447  ExcIndexRange(unrolled_index, 0, t.n_independent_components));
2448  t[TensorType::unrolled_to_component_indices(unrolled_index)] = value;
2449  }
2450 
2451 
2456  template <int dim, typename NumberType>
2458  const unsigned int & unrolled_index,
2459  const NumberType & value)
2460  {
2461  Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1));
2462  (void)unrolled_index;
2463  t = value;
2464  }
2465 
2466 
2472  template <typename NumberType>
2473  inline void
2475  const unsigned int &unrolled_index,
2476  const NumberType & value)
2477  {
2478  Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1));
2479  (void)unrolled_index;
2480  t = value;
2481  }
2482 
2483 
2489  template <int dim, typename NumberType>
2491  const unsigned int &unrolled_index_row,
2492  const unsigned int &unrolled_index_col,
2493  const NumberType & value)
2494  {
2495  // Fourth order symmetric tensors require a specialized interface
2496  // to extract values.
2497  using SubTensorType = SymmetricTensor<2, dim, NumberType>;
2498  Assert(unrolled_index_row < SubTensorType::n_independent_components,
2499  ExcIndexRange(unrolled_index_row,
2500  0,
2501  SubTensorType::n_independent_components));
2502  Assert(unrolled_index_col < SubTensorType::n_independent_components,
2503  ExcIndexRange(unrolled_index_col,
2504  0,
2505  SubTensorType::n_independent_components));
2506  const TableIndices<2> indices_row =
2507  SubTensorType::unrolled_to_component_indices(unrolled_index_row);
2508  const TableIndices<2> indices_col =
2509  SubTensorType::unrolled_to_component_indices(unrolled_index_col);
2510  t[indices_row[0]][indices_row[1]][indices_col[0]][indices_col[1]] =
2511  value;
2512  }
2513 
2514 
2519  template <int rank,
2520  int dim,
2521  typename NumberType,
2522  template <int, int, typename> class TensorType>
2523  inline NumberType
2524  get_tensor_entry(const TensorType<rank, dim, NumberType> &t,
2525  const unsigned int & unrolled_index)
2526  {
2527  // Where possible, get values using TableIndices
2528  Assert(unrolled_index < t.n_independent_components,
2529  ExcIndexRange(unrolled_index, 0, t.n_independent_components));
2530  return t[TensorType<rank, dim, NumberType>::
2531  unrolled_to_component_indices(unrolled_index)];
2532  }
2533 
2534 
2539  template <int dim,
2540  typename NumberType,
2541  template <int, int, typename> class TensorType>
2542  inline NumberType
2543  get_tensor_entry(const TensorType<0, dim, NumberType> &t,
2544  const unsigned int & unrolled_index)
2545  {
2546  Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1));
2547  (void)unrolled_index;
2548  return t;
2549  }
2550 
2551 
2557  template <typename NumberType>
2558  inline const NumberType &
2559  get_tensor_entry(const NumberType &t, const unsigned int &unrolled_index)
2560  {
2561  Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1));
2562  (void)unrolled_index;
2563  return t;
2564  }
2565 
2566 
2571  template <int rank,
2572  int dim,
2573  typename NumberType,
2574  template <int, int, typename> class TensorType>
2575  inline NumberType &
2576  get_tensor_entry(TensorType<rank, dim, NumberType> &t,
2577  const unsigned int & unrolled_index)
2578  {
2579  // Where possible, get values using TableIndices
2580  Assert(unrolled_index < t.n_independent_components,
2581  ExcIndexRange(unrolled_index, 0, t.n_independent_components));
2582  return t[TensorType<rank, dim, NumberType>::
2583  unrolled_to_component_indices(unrolled_index)];
2584  }
2585 
2586 
2591  template <int dim,
2592  typename NumberType,
2593  template <int, int, typename> class TensorType>
2594  NumberType &get_tensor_entry(TensorType<0, dim, NumberType> &t,
2595  const unsigned int & index)
2596  {
2597  Assert(index == 0, ExcIndexRange(index, 0, 1));
2598  (void)index;
2599  return t;
2600  }
2601 
2602 
2608  template <typename NumberType>
2609  inline NumberType &
2610  get_tensor_entry(NumberType &t, const unsigned int &index)
2611  {
2612  Assert(index == 0, ExcIndexRange(index, 0, 1));
2613  (void)index;
2614  return t;
2615  }
2616 
2617  } // namespace internal
2618 
2619 
2620 
2644  template <int dim,
2645  enum AD::NumberTypes ADNumberTypeCode,
2646  typename ScalarType = double>
2648  : public HelperBase<ADNumberTypeCode, ScalarType>
2649  {
2650  public:
2655  static const unsigned int dimension = dim;
2656 
2661  using scalar_type =
2663 
2668  using ad_type =
2670 
2675 
2691  const unsigned int n_dependent_variables);
2692 
2696  virtual ~PointLevelFunctionsBase() = default;
2697 
2699 
2704 
2708  virtual void
2709  reset(const unsigned int n_independent_variables =
2711  const unsigned int n_dependent_variables =
2713  const bool clear_registered_tapes = true) override;
2714 
2730  void
2731  register_independent_variables(const std::vector<scalar_type> &values);
2732 
2761  template <typename ValueType, typename ExtractorType>
2762  void
2763  register_independent_variable(const ValueType & value,
2764  const ExtractorType &extractor);
2765 
2786  const std::vector<ad_type> &
2787  get_sensitive_variables() const;
2788 
2789  /*
2790  * Extract a subset of the independent variables as represented by
2791  * auto-differentiable numbers. These are the independent
2792  * variables @f$\mathbf{A} \subset \mathbf{X}@f$ at which the dependent values
2793  * are evaluated and differentiated.
2794  *
2795  * This function indicates to the AD library that implements the
2796  * auto-differentiable number type that operations performed on these
2797  * numbers are to be tracked so they are considered "sensitive"
2798  * variables. This is, therefore, the set of variables with which one
2799  * would then perform computations, and based on which one can then
2800  * extract both the value of the function and its derivatives with the
2801  * member functions below. The values of the components of the returned
2802  * object are initialized to the values set with
2803  * register_independent_variable().
2804  *
2805  * @param[in] extractor An extractor associated with the input field
2806  * variables. This effectively defines which components of the global set
2807  * of independent variables this field is associated with.
2808  * @return An object of auto-differentiable type numbers. The return type is
2809  * based on the input extractor, and will be either a scalar,
2810  * Tensor<1,dim>, Tensor<2,dim>, or SymmetricTensor<2,dim>.
2811  *
2812  * @note For taped AD numbers, this operation is only valid in recording mode.
2813  */
2814  template <typename ExtractorType>
2815  typename internal::Extractor<dim,
2816  ExtractorType>::template tensor_type<ad_type>
2817  get_sensitive_variables(const ExtractorType &extractor) const;
2818 
2820 
2825 
2843  void
2844  set_independent_variables(const std::vector<scalar_type> &values);
2845 
2872  template <typename ValueType, typename ExtractorType>
2873  void
2874  set_independent_variable(const ValueType & value,
2875  const ExtractorType &extractor);
2876 
2878 
2879  protected:
2884 
2895  void
2896  set_sensitivity_value(const unsigned int index,
2897  const bool symmetric_component,
2898  const scalar_type &value);
2899 
2905  bool
2906  is_symmetric_independent_variable(const unsigned int index) const;
2907 
2912  unsigned int
2914 
2916 
2917  private:
2922 
2928 
2930 
2931  }; // class PointLevelFunctionsBase
2932 
2933 
2934 
3091  template <int dim,
3092  enum AD::NumberTypes ADNumberTypeCode,
3093  typename ScalarType = double>
3095  : public PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>
3096  {
3097  public:
3102  using scalar_type =
3104 
3109  using ad_type =
3111 
3116 
3126  ScalarFunction(const unsigned int n_independent_variables);
3127 
3131  virtual ~ScalarFunction() = default;
3132 
3134 
3139 
3151  void
3152  register_dependent_variable(const ad_type &func);
3153 
3161  scalar_type
3162  compute_value() const;
3163 
3176  void
3177  compute_gradient(Vector<scalar_type> &gradient) const;
3178 
3193  void
3194  compute_hessian(FullMatrix<scalar_type> &hessian) const;
3195 
3230  template <typename ExtractorType_Row>
3231  static typename internal::
3233  extract_gradient_component(const Vector<scalar_type> &gradient,
3234  const ExtractorType_Row & extractor_row);
3235 
3274  template <typename ExtractorType_Row, typename ExtractorType_Col>
3275  static typename internal::ScalarFieldHessian<dim,
3276  scalar_type,
3277  ExtractorType_Row,
3278  ExtractorType_Col>::type
3280  const ExtractorType_Row & extractor_row,
3281  const ExtractorType_Col & extractor_col);
3282 
3298  const FullMatrix<scalar_type> & hessian,
3299  const FEValuesExtractors::Scalar &extractor_row,
3300  const FEValuesExtractors::Scalar &extractor_col);
3301 
3315  const FullMatrix<scalar_type> & hessian,
3316  const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
3317  const FEValuesExtractors::SymmetricTensor<2> &extractor_col);
3318 
3320 
3321  }; // class ScalarFunction
3322 
3323 
3324 
3484  template <int dim,
3485  enum AD::NumberTypes ADNumberTypeCode,
3486  typename ScalarType = double>
3488  : public PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>
3489  {
3490  public:
3495  using scalar_type =
3497 
3502  using ad_type =
3504 
3509 
3524  VectorFunction(const unsigned int n_independent_variables,
3525  const unsigned int n_dependent_variables);
3526 
3530  virtual ~VectorFunction() = default;
3531 
3533 
3538 
3551  void
3552  register_dependent_variables(const std::vector<ad_type> &funcs);
3553 
3572  template <typename ValueType, typename ExtractorType>
3573  void
3574  register_dependent_variable(const ValueType & funcs,
3575  const ExtractorType &extractor);
3576 
3585  void
3586  compute_values(Vector<scalar_type> &values) const;
3587 
3602  void
3603  compute_jacobian(FullMatrix<scalar_type> &jacobian) const;
3604 
3605 
3618  template <typename ExtractorType_Row>
3619  static typename internal::
3621  extract_value_component(const Vector<scalar_type> &values,
3622  const ExtractorType_Row & extractor_row);
3623 
3671  template <typename ExtractorType_Row, typename ExtractorType_Col>
3672  static typename internal::VectorFieldJacobian<dim,
3673  scalar_type,
3674  ExtractorType_Row,
3675  ExtractorType_Col>::type
3677  const ExtractorType_Row & extractor_row,
3678  const ExtractorType_Col & extractor_col);
3679 
3697  const FullMatrix<scalar_type> & jacobian,
3698  const FEValuesExtractors::Scalar &extractor_row,
3699  const FEValuesExtractors::Scalar &extractor_col);
3700 
3716  const FullMatrix<scalar_type> & jacobian,
3717  const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
3718  const FEValuesExtractors::SymmetricTensor<2> &extractor_col);
3719 
3721 
3722  }; // class VectorFunction
3723 
3724 
3725  } // namespace AD
3726 } // namespace Differentiation
3727 
3728 
3729 /* ----------------- inline and template functions ----------------- */
3730 
3731 
3732 # ifndef DOXYGEN
3733 
3734 namespace Differentiation
3735 {
3736  namespace AD
3737  {
3738  /* ----------------- CellLevelBase ----------------- */
3739 
3740 
3741 
3742  template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
3743  template <typename VectorType>
3744  void
3746  const VectorType & values,
3747  const std::vector<::types::global_dof_index> &local_dof_indices)
3748  {
3749  // This is actually the same thing the set_dof_values() function,
3750  // in the sense that we simply populate our array of independent values
3751  // with a meaningful number. However, in this case we need to double check
3752  // that we're not registering these variables twice
3753  Assert(
3754  local_dof_indices.size() == this->n_independent_variables(),
3755  ExcMessage(
3756  "Degree of freedom index vector size does not match number of independent variables"));
3757  for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
3758  {
3759  Assert(this->registered_independent_variable_values[i] == false,
3760  ExcMessage("Independent variables already registered."));
3761  }
3762  set_dof_values(values, local_dof_indices);
3763  }
3764 
3765 
3766 
3767  template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
3768  template <typename VectorType>
3769  void
3771  const VectorType & values,
3772  const std::vector<::types::global_dof_index> &local_dof_indices)
3773  {
3774  Assert(local_dof_indices.size() == this->n_independent_variables(),
3775  ExcMessage(
3776  "Vector size does not match number of independent variables"));
3777  for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
3779  i, values[local_dof_indices[i]]);
3780  }
3781 
3782 
3783 
3784  /* ----------------- PointLevelFunctionsBase ----------------- */
3785 
3786 
3787 
3788  template <int dim,
3789  enum AD::NumberTypes ADNumberTypeCode,
3790  typename ScalarType>
3791  template <typename ValueType, typename ExtractorType>
3792  void
3794  register_independent_variable(const ValueType & value,
3795  const ExtractorType &extractor)
3796  {
3797  // This is actually the same thing as the set_independent_variable
3798  // function, in the sense that we simply populate our array of independent
3799  // values with a meaningful number. However, in this case we need to
3800  // double check that we're not registering these variables twice
3801 # ifdef DEBUG
3802  const std::vector<unsigned int> index_set(
3803  internal::extract_field_component_indices<dim>(extractor));
3804  for (unsigned int i = 0; i < index_set.size(); ++i)
3805  {
3806  Assert(
3807  this->registered_independent_variable_values[index_set[i]] == false,
3808  ExcMessage(
3809  "Overlapping indices for independent variables. "
3810  "One or more indices associated with the field that "
3811  "is being registered as an independent variable have "
3812  "already been associated with another field. This suggests "
3813  "that the component offsets used to construct their counterpart "
3814  "extractors are incompatible with one another. Make sure that "
3815  "the first component for each extractor properly takes into "
3816  "account the dimensionality of the preceding fields."));
3817  }
3818 # endif
3819  set_independent_variable(value, extractor);
3820  }
3821 
3822 
3823 
3824  template <int dim,
3825  enum AD::NumberTypes ADNumberTypeCode,
3826  typename ScalarType>
3827  template <typename ValueType, typename ExtractorType>
3828  void
3830  set_independent_variable(const ValueType & value,
3831  const ExtractorType &extractor)
3832  {
3833  const std::vector<unsigned int> index_set(
3834  internal::extract_field_component_indices<dim>(extractor));
3835  for (unsigned int i = 0; i < index_set.size(); ++i)
3836  {
3837  set_sensitivity_value(
3838  index_set[i],
3839  internal::Extractor<dim, ExtractorType>::symmetric_component(i),
3840  internal::get_tensor_entry(value, i));
3841  }
3842  }
3843 
3844 
3845 
3846  template <int dim,
3847  enum AD::NumberTypes ADNumberTypeCode,
3848  typename ScalarType>
3849  template <typename ExtractorType>
3850  typename internal::Extractor<dim, ExtractorType>::template tensor_type<
3853  get_sensitive_variables(const ExtractorType &extractor) const
3854  {
3855  if (ADNumberTraits<ad_type>::is_taped == true)
3856  {
3857  Assert(this->active_tape_index() !=
3859  ExcMessage("Invalid tape index"));
3860  }
3861 
3862  // If necessary, finalize the internally stored vector of
3863  // AD numbers that represents the independent variables
3864  this->finalize_sensitive_independent_variables();
3865  Assert(this->independent_variables.size() ==
3866  this->n_independent_variables(),
3867  ExcDimensionMismatch(this->independent_variables.size(),
3868  this->n_independent_variables()));
3869 
3870  const std::vector<unsigned int> index_set(
3871  internal::extract_field_component_indices<dim>(extractor));
3872  typename internal::Extractor<dim,
3873  ExtractorType>::template tensor_type<ad_type>
3874  out;
3875 
3876  for (unsigned int i = 0; i < index_set.size(); ++i)
3877  {
3878  const unsigned int index = index_set[i];
3879  Assert(index < this->n_independent_variables(), ExcInternalError());
3880  Assert(this->registered_independent_variable_values[index] == true,
3881  ExcInternalError());
3882  internal::get_tensor_entry(out, i) =
3883  this->independent_variables[index];
3884  }
3885 
3886  return out;
3887  }
3888 
3889 
3890 
3891  /* ----------------- ScalarFunction ----------------- */
3892 
3893 
3894 
3895  template <int dim,
3896  enum AD::NumberTypes ADNumberTypeCode,
3897  typename ScalarType>
3898  template <typename ExtractorType_Row>
3899  typename internal::ScalarFieldGradient<
3900  dim,
3902  ExtractorType_Row>::type
3904  extract_gradient_component(const Vector<scalar_type> &gradient,
3905  const ExtractorType_Row & extractor_row)
3906  {
3907  // NOTE: The order of components must be consistently defined throughout
3908  // this class.
3909  typename internal::
3911 
3912  // Get indexsets for the subblock from which we wish to extract the
3913  // gradient values
3914  const std::vector<unsigned int> row_index_set(
3915  internal::extract_field_component_indices<dim>(extractor_row));
3916  Assert(out.n_independent_components == row_index_set.size(),
3917  ExcMessage("Not all tensor components have been extracted!"));
3918  for (unsigned int r = 0; r < row_index_set.size(); ++r)
3919  internal::set_tensor_entry(out, r, gradient[row_index_set[r]]);
3920 
3921  return out;
3922  }
3923 
3924 
3925 
3926  template <int dim,
3927  enum AD::NumberTypes ADNumberTypeCode,
3928  typename ScalarType>
3929  template <typename ExtractorType_Row, typename ExtractorType_Col>
3930  typename internal::ScalarFieldHessian<
3931  dim,
3933  ExtractorType_Row,
3934  ExtractorType_Col>::type
3937  const ExtractorType_Row & extractor_row,
3938  const ExtractorType_Col & extractor_col)
3939  {
3940  using InternalHessian = internal::ScalarFieldHessian<dim,
3941  scalar_type,
3942  ExtractorType_Row,
3943  ExtractorType_Col>;
3944  using InternalExtractorRow = internal::Extractor<dim, ExtractorType_Row>;
3945  using InternalExtractorCol = internal::Extractor<dim, ExtractorType_Col>;
3946  using HessianType = typename InternalHessian::type;
3947 
3948  // NOTE: The order of components must be consistently defined throughout
3949  // this class.
3950  HessianType out;
3951 
3952  // Get indexsets for the subblocks from which we wish to extract the
3953  // Hessian values
3954  // NOTE: Here we have to do some clever accounting when the
3955  // one extractor is a symmetric Tensor and the other is not, e.g.
3956  // <SymmTensor,Vector>. In this scenario the return type is a
3957  // non-symmetric Tensor<3,dim> but we have to fetch information from a
3958  // SymmTensor row/column that has too few entries to fill the output
3959  // tensor. So we must duplicate the relevant entries in the row/column
3960  // indexset to fetch off-diagonal components that are Otherwise
3961  // non-existent in a SymmTensor.
3962  const std::vector<unsigned int> row_index_set(
3963  internal::extract_field_component_indices<dim>(
3964  extractor_row, false /*ignore_symmetries*/));
3965  const std::vector<unsigned int> col_index_set(
3966  internal::extract_field_component_indices<dim>(
3967  extractor_col, false /*ignore_symmetries*/));
3968 
3969  for (unsigned int index = 0;
3970  index < HessianType::n_independent_components;
3971  ++index)
3972  {
3973  const TableIndices<HessianType::rank> ti_out =
3974  HessianType::unrolled_to_component_indices(index);
3975  const unsigned int r =
3976  InternalExtractorRow::local_component(ti_out, 0);
3977  const unsigned int c =
3978  InternalExtractorCol::local_component(ti_out,
3979  InternalExtractorRow::rank);
3980 
3982  out, index, hessian[row_index_set[r]][col_index_set[c]]);
3983  }
3984 
3985  return out;
3986  }
3987 
3988 
3989 
3990  /* ----------------- VectorFunction ----------------- */
3991 
3992 
3993 
3994  template <int dim,
3995  enum AD::NumberTypes ADNumberTypeCode,
3996  typename ScalarType>
3997  template <typename ValueType, typename ExtractorType>
3998  void
4000  register_dependent_variable(const ValueType & funcs,
4001  const ExtractorType &extractor)
4002  {
4003  const std::vector<unsigned int> index_set(
4004  internal::extract_field_component_indices<dim>(extractor));
4005  for (unsigned int i = 0; i < index_set.size(); ++i)
4006  {
4007  Assert(this->registered_marked_dependent_variables[index_set[i]] ==
4008  false,
4009  ExcMessage("Overlapping indices for dependent variables."));
4011  index_set[i], internal::get_tensor_entry(funcs, i));
4012  }
4013  }
4014 
4015 
4016 
4017  template <int dim,
4018  enum AD::NumberTypes ADNumberTypeCode,
4019  typename ScalarType>
4020  template <typename ExtractorType_Row>
4021  typename internal::VectorFieldValue<
4022  dim,
4024  ExtractorType_Row>::type
4026  const Vector<scalar_type> &values,
4027  const ExtractorType_Row & extractor_row)
4028  {
4029  // NOTE: The order of components must be consistently defined throughout
4030  // this class.
4031  typename internal::VectorFieldValue<dim, scalar_type, ExtractorType_Row>::
4032  type out;
4033 
4034  // Get indexsets for the subblock from which we wish to extract the
4035  // gradient values
4036  const std::vector<unsigned int> row_index_set(
4037  internal::extract_field_component_indices<dim>(extractor_row));
4038  Assert(out.n_independent_components == row_index_set.size(),
4039  ExcMessage("Not all tensor components have been extracted!"));
4040  for (unsigned int r = 0; r < row_index_set.size(); ++r)
4041  internal::set_tensor_entry(out, r, values[row_index_set[r]]);
4042 
4043  return out;
4044  }
4045 
4046 
4047 
4048  template <int dim,
4049  enum AD::NumberTypes ADNumberTypeCode,
4050  typename ScalarType>
4051  template <typename ExtractorType_Row, typename ExtractorType_Col>
4053  dim,
4055  ExtractorType_Row,
4056  ExtractorType_Col>::type
4059  const ExtractorType_Row & extractor_row,
4060  const ExtractorType_Col & extractor_col)
4061  {
4062  using InternalJacobian = internal::VectorFieldJacobian<dim,
4063  scalar_type,
4064  ExtractorType_Row,
4065  ExtractorType_Col>;
4066  using InternalExtractorRow = internal::Extractor<dim, ExtractorType_Row>;
4067  using InternalExtractorCol = internal::Extractor<dim, ExtractorType_Col>;
4068  using JacobianType = typename InternalJacobian::type;
4069 
4070  // NOTE: The order of components must be consistently defined throughout
4071  // this class.
4072  JacobianType out;
4073 
4074  // Get indexsets for the subblocks from which we wish to extract the
4075  // Hessian values.
4076  // NOTE: Here we have to do some clever accounting when the
4077  // one extractor is a symmetric Tensor and the other is not, e.g.
4078  // <SymmTensor,Vector>. In this scenario the return type is a
4079  // non-symmetric Tensor<3,dim> but we have to fetch information from a
4080  // SymmTensor row/column that has too few entries to fill the output
4081  // tensor. So we must duplicate the relevant entries in the row/column
4082  // indexset to fetch off-diagonal components that are Otherwise
4083  // non-existent in a SymmTensor.
4084  const std::vector<unsigned int> row_index_set(
4085  internal::extract_field_component_indices<dim>(
4086  extractor_row, false /*ignore_symmetries*/));
4087  const std::vector<unsigned int> col_index_set(
4088  internal::extract_field_component_indices<dim>(
4089  extractor_col, false /*ignore_symmetries*/));
4090 
4091  for (unsigned int index = 0;
4092  index < JacobianType::n_independent_components;
4093  ++index)
4094  {
4095  const TableIndices<JacobianType::rank> ti_out =
4096  JacobianType::unrolled_to_component_indices(index);
4097  const unsigned int r =
4098  InternalExtractorRow::local_component(ti_out, 0);
4099  const unsigned int c =
4100  InternalExtractorCol::local_component(ti_out,
4101  InternalExtractorRow::rank);
4102 
4104  out, index, jacobian[row_index_set[r]][col_index_set[c]]);
4105  }
4106 
4107  return out;
4108  }
4109 
4110 
4111  } // namespace AD
4112 } // namespace Differentiation
4113 
4114 
4115 # endif // DOXYGEN
4116 
4117 
4118 DEAL_II_NAMESPACE_CLOSE
4119 
4120 #endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)
4121 
4122 #endif // dealii_differentiation_ad_ad_helpers_h
TapedDrivers< ad_type, scalar_type > taped_driver
Definition: ad_helpers.h:594
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:399
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:2165
static const unsigned int invalid_unsigned_int
Definition: types.h:173
HelperBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:38
static unsigned int first_component(const FEValuesExtractors::SymmetricTensor< 2 > &extractor)
Definition: ad_helpers.h:2112
PointLevelFunctionsBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1164
static unsigned int component_to_unrolled_index(const TableIndices< rank_ > &indices)
void compute_jacobian(FullMatrix< scalar_type > &jacobian) const
Definition: ad_helpers.cc:1742
VectorFunction(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1655
static unsigned int first_component(const FEValuesExtractors::Tensor< 1 > &extractor)
Definition: ad_helpers.h:1904
void set_sensitivity_value(const unsigned int index, const scalar_type &value)
Definition: ad_helpers.cc:109
std::vector< IndexType > extract_field_component_indices(const ExtractorType &extractor, const bool ignore_symmetries=true)
Definition: ad_helpers.h:2367
static TableIndices< rank > table_index_view(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:2139
void mark_independent_variable(const unsigned int index, ad_type &out) const
Definition: ad_helpers.cc:146
void set_tensor_entry(TensorType &t, const unsigned int &unrolled_index, const NumberType &value)
Definition: ad_helpers.h:2441
void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:877
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1850
CellLevelBase(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:719
void register_independent_variables(const std::vector< scalar_type > &values)
Definition: ad_helpers.cc:1229
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3103
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)
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)
static TableIndices< rank > table_index_view(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:2032
virtual void compute_residual(Vector< scalar_type > &residual) const override
Definition: ad_helpers.cc:1041
typename Extractor< dim, ExtractorType >::template tensor_type< NumberType > type
Definition: ad_helpers.h:2192
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:568
void register_independent_variable(const ValueType &value, const ExtractorType &extractor)
const std::vector< ad_type > & get_sensitive_dof_values() const
Definition: ad_helpers.cc:753
static unsigned int first_component(const FEValuesExtractors::Tensor< 2 > &extractor)
Definition: ad_helpers.h:2007
unsigned int tape_buffer_sizes
Definition: ad_drivers.h:104
void register_dependent_variable(const ad_type &func)
Definition: ad_helpers.cc:1344
void stop_recording_operations(const bool write_tapes_to_file=false)
Definition: ad_helpers.cc:652
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const =0
static TableIndices< rank > table_index_view(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1929
void compute_gradient(Vector< scalar_type > &gradient) const
Definition: ad_helpers.cc:1407
static void configure_tapeless_mode(const unsigned int n_independent_variables, const bool ensure_persistent_setting=true)
Definition: ad_helpers.cc:457
std::vector< scalar_type > independent_variable_values
Definition: ad_helpers.h:638
Types< ad_type >::tape_index active_tape_index() const
Definition: ad_helpers.cc:288
LinearAlgebra::distributed::Vector< Number > Vector
std::vector< bool > symmetric_independent_variables
Definition: ad_helpers.h:2927
void set_sensitivity_value(const unsigned int index, const bool symmetric_component, const scalar_type &value)
Definition: ad_helpers.cc:1283
std::vector< bool > registered_marked_independent_variables
Definition: ad_helpers.h:660
ResidualLinearization(const unsigned int n_independent_variables, const unsigned int n_dependent_variables)
Definition: ad_helpers.cc:1015
static internal::ScalarFieldGradient< dim, scalar_type, ExtractorType_Row >::type extract_gradient_component(const Vector< scalar_type > &gradient, const ExtractorType_Row &extractor_row)
static const Types< ADNumberType >::tape_index invalid_tape_index
Definition: ad_drivers.h:121
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:587
void compute_values(Vector< scalar_type > &values) const
Definition: ad_helpers.cc:1686
unsigned int n_symmetric_independent_variables() const
Definition: ad_helpers.cc:1215
void reset_registered_dependent_variables(const bool flag=false)
Definition: ad_helpers.cc:96
void register_dependent_variable(const unsigned int index, const ad_type &func)
Definition: ad_helpers.cc:688
void print_values(std::ostream &stream) const
Definition: ad_helpers.cc:368
static ::ExceptionBase & ExcMessage(std::string arg1)
typename HessianType< ExtractorType_Row, ExtractorType_Col >::template type< rank, dim, NumberType > type
Definition: ad_helpers.h:2326
std::vector< ad_type > dependent_variables
Definition: ad_helpers.h:752
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::ad_type ad_type
Definition: ad_helpers.h:187
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:1095
static unsigned int first_component(const FEValuesExtractors::Vector &extractor)
Definition: ad_helpers.h:1800
void compute_hessian(FullMatrix< scalar_type > &hessian) const
Definition: ad_helpers.cc:1481
std::size_t n_independent_variables() const
Definition: ad_helpers.cc:245
typename AD::NumberTraits< ScalarType, ADNumberTypeCode >::scalar_type scalar_type
Definition: ad_helpers.h:180
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< bool > registered_marked_dependent_variables
Definition: ad_helpers.h:757
ScalarFunction(const unsigned int n_independent_variables)
Definition: ad_helpers.cc:1330
void register_dependent_variable(const ValueType &funcs, const ExtractorType &extractor)
unsigned int n_registered_dependent_variables() const
Definition: ad_helpers.cc:254
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:1177
void register_energy_functional(const ad_type &energy)
Definition: ad_helpers.cc:811
static TableIndices< rank > table_index_view(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1826
static bool symmetric_component(const unsigned int unrolled_index)
Definition: ad_helpers.h:1727
ScalarFieldHessian< dim, NumberType, ExtractorType_Field, ExtractorType_Derivative > VectorFieldJacobian
Definition: ad_helpers.h:2355
void print(std::ostream &stream) const
Definition: ad_helpers.cc:313
static TableIndices< rank_ > unrolled_to_component_indices(const unsigned int i)
Definition: tensor.h:1416
void finalize_sensitive_independent_variables() const
Definition: ad_helpers.cc:185
ScalarFieldGradient< dim, NumberType, ExtractorType_Field > VectorFieldValue
Definition: ad_helpers.h:2336
virtual void compute_linearization(FullMatrix< scalar_type > &linearization) const override
Definition: ad_helpers.cc:942
static unsigned int first_component(const FEValuesExtractors::Scalar &extractor)
Definition: ad_helpers.h:1714
bool is_symmetric_independent_variable(const unsigned int index) const
Definition: ad_helpers.cc:1201
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1742
void register_dof_values(const std::vector< scalar_type > &dof_values)
Definition: ad_helpers.cc:730
NumberType get_tensor_entry(const TensorType< rank, dim, NumberType > &t, const unsigned int &unrolled_index)
Definition: ad_helpers.h:2524
bool active_tape_requires_retaping() const
Definition: ad_helpers.cc:506
void register_dependent_variables(const std::vector< ad_type > &funcs)
Definition: ad_helpers.cc:1670
Definition: mpi.h:90
std::vector< bool > registered_independent_variable_values
Definition: ad_helpers.h:654
static internal::VectorFieldValue< dim, scalar_type, ExtractorType_Row >::type extract_value_component(const Vector< scalar_type > &values, const ExtractorType_Row &extractor_row)
static bool symmetric_component(const unsigned int unrolled_index)
Definition: ad_helpers.h:1814
virtual void compute_residual(Vector< scalar_type > &residual) const =0
void initialize_non_sensitive_independent_variable(const unsigned int index, ad_type &out) const
Definition: ad_helpers.cc:209
bool recorded_tape_requires_retaping(const typename Types< ad_type >::tape_index tape_index) const
Definition: ad_helpers.cc:493
void set_independent_variable(const ValueType &value, const ExtractorType &extractor)
void register_residual_vector(const std::vector< ad_type > &residual)
Definition: ad_helpers.cc:1027
void print_tape_stats(const typename Types< ad_type >::tape_index tape_index, std::ostream &stream) const
Definition: ad_helpers.cc:382
TapelessDrivers< ad_type, scalar_type > tapeless_driver
Definition: ad_helpers.h:602
bool is_registered_tape(const typename Types< ad_type >::tape_index tape_index) const
Definition: ad_helpers.cc:300
std::vector< ad_type > independent_variables
Definition: ad_helpers.h:648
std::size_t n_dependent_variables() const
Definition: ad_helpers.cc:266
void activate_recorded_tape(const typename Types< ad_type >::tape_index tape_index)
Definition: ad_helpers.cc:483
typename HelperBase< ADNumberTypeCode, ScalarType >::scalar_type scalar_type
Definition: ad_helpers.h:3496
void set_independent_variables(const std::vector< scalar_type > &values)
Definition: ad_helpers.cc:1305
EnergyFunctional(const unsigned int n_independent_variables)
Definition: ad_helpers.cc:802
const std::vector< ad_type > & get_sensitive_variables() const
Definition: ad_helpers.cc:1255
void set_dof_values(const std::vector< scalar_type > &dof_values)
Definition: ad_helpers.cc:778
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:1953
unsigned int n_registered_independent_variables() const
Definition: ad_helpers.cc:234
static IndexType local_component(const TableIndices< rank_in > &table_indices, const unsigned int column_offset)
Definition: ad_helpers.h:2058
static ::ExceptionBase & ExcInternalError()
void activate_tape(const typename Types< ad_type >::tape_index tape_index, const bool read_mode)
Definition: ad_helpers.cc:531