Add functions to add tensor types to SD subsitution maps

diff --git a/include/deal.II/differentiation/sd/symengine_tensor_operations.h b/include/deal.II/differentiation/sd/symengine_tensor_operations.h
index 22ab4c7dde04637aca208f4eb955b516b67d1a11..c325bf9cc89c412c8b3a0bbe709e939b872fa597 100644 (file)
--- a/include/deal.II/differentiation/sd/symengine_tensor_operations.h
+++ b/include/deal.II/differentiation/sd/symengine_tensor_operations.h
@@ -366,6 +366,61 @@ namespace Differentiation
 
     //@}
 
+    /**
+     * @name Symbol substitution map enlargement
+     */
+    //@{
+
+    /**
+     * A convenience function for adding an entry to the @p substitution_map.
+     * The new entries will have the keys given in the @p symbol_tensor with
+     * their paired values extracted from the corresponding elements of the
+     * @p value_tensor.
+     *
+     * For more context which this function is used, see the other
+     * \ref add_to_substitution_map(types::substitution_map &,const
+     * Expression &, const Expression &) function.
+     *
+     * @tparam ignore_invalid_symbols See the other
+     * \ref add_to_substitution_map(types::substitution_map &,const
+     * Expression &,const Expression &) function for a detailed
+     * discussion on the role of this template argument.
+     */
+    template <bool ignore_invalid_symbols = false,
+              int  rank,
+              int  dim,
+              typename ValueType>
+    void
+    add_to_substitution_map(types::substitution_map &substitution_map,
+                            const Tensor<rank, dim, Expression> &symbol_tensor,
+                            const Tensor<rank, dim, ValueType> & value_tensor);
+
+    /**
+     * A convenience function for adding an entry to the @p substitution_map.
+     * The new entries will have the keys given in the @p symbol_tensor with
+     * their paired values extracted from the corresponding elements of the @p value_tensor.
+     *
+     * For more context which this function is used, see the other
+     * \ref add_to_substitution_map(types::substitution_map &,const
+     * Expression &, const Expression &) function.
+     *
+     * @tparam ignore_invalid_symbols See the other
+     * \ref add_to_substitution_map(types::substitution_map &,const
+     * Expression &,const Expression &) function for a detailed
+     * discussion on the role of this template argument.
+     */
+    template <bool ignore_invalid_symbols = false,
+              int  rank,
+              int  dim,
+              typename ValueType>
+    void
+    add_to_substitution_map(
+      types::substitution_map &                     substitution_map,
+      const SymmetricTensor<rank, dim, Expression> &symbol_tensor,
+      const SymmetricTensor<rank, dim, ValueType> & value_tensor);
+
+    //@}
+
   } // namespace SD
 } // namespace Differentiation
 
@@ -384,6 +439,12 @@ namespace Differentiation
 
     namespace internal
     {
+      template <int dim>
+      TableIndices<4>
+      make_rank_4_tensor_indices(const unsigned int &idx_i,
+                                 const unsigned int &idx_j);
+
+
       template <int rank_1, int rank_2>
       TableIndices<rank_1 + rank_2>
       concatenate_indices(const TableIndices<rank_1> &indices_1,
@@ -749,6 +810,88 @@ namespace Differentiation
       return internal::tensor_diff_tensor(symbol_tensor, op);
     }
 
+
+    /* ---------------- Symbolic substitution map enlargement --------------*/
+
+
+    namespace internal
+    {
+      template <int rank,
+                int dim,
+                typename ValueType,
+                template <int, int, typename> class TensorType>
+      std::vector<std::pair<Expression, ValueType>>
+      tensor_substitution_map(
+        const TensorType<rank, dim, Expression> &symbol_tensor,
+        const TensorType<rank, dim, ValueType> & value_tensor)
+      {
+        std::vector<std::pair<Expression, ValueType>> symbol_values;
+        for (unsigned int i = 0; i < symbol_tensor.n_independent_components;
+             ++i)
+          {
+            const TableIndices<rank> indices(
+              symbol_tensor.unrolled_to_component_indices(i));
+            symbol_values.push_back(
+              std::make_pair(symbol_tensor[indices], value_tensor[indices]));
+          }
+        return symbol_values;
+      }
+
+
+      template <int dim, typename ValueType>
+      std::vector<std::pair<Expression, ValueType>>
+      tensor_substitution_map(
+        const SymmetricTensor<4, dim, Expression> &symbol_tensor,
+        const SymmetricTensor<4, dim, ValueType> & value_tensor)
+      {
+        std::vector<std::pair<Expression, ValueType>> symbol_values;
+        for (unsigned int i = 0;
+             i < SymmetricTensor<2, dim>::n_independent_components;
+             ++i)
+          for (unsigned int j = 0;
+               j < SymmetricTensor<2, dim>::n_independent_components;
+               ++j)
+            {
+              const TableIndices<4> indices =
+                make_rank_4_tensor_indices<dim>(i, j);
+              symbol_values.push_back(
+                std::make_pair(symbol_tensor[indices], value_tensor[indices]));
+            }
+        return symbol_values;
+      }
+    } // namespace internal
+
+
+    template <bool ignore_invalid_symbols,
+              int  rank,
+              int  dim,
+              typename ValueType>
+    void
+    add_to_substitution_map(types::substitution_map &substitution_map,
+                            const Tensor<rank, dim, Expression> &symbol_tensor,
+                            const Tensor<rank, dim, ValueType> & value_tensor)
+    {
+      add_to_substitution_map<ignore_invalid_symbols>(
+        substitution_map,
+        internal::tensor_substitution_map(symbol_tensor, value_tensor));
+    }
+
+
+    template <bool ignore_invalid_symbols,
+              int  rank,
+              int  dim,
+              typename ValueType>
+    void
+    add_to_substitution_map(
+      types::substitution_map &                     substitution_map,
+      const SymmetricTensor<rank, dim, Expression> &symbol_tensor,
+      const SymmetricTensor<rank, dim, ValueType> & value_tensor)
+    {
+      add_to_substitution_map<ignore_invalid_symbols>(
+        substitution_map,
+        internal::tensor_substitution_map(symbol_tensor, value_tensor));
+    }
+
   } // namespace SD
 } // namespace Differentiation
 

diff --git a/source/differentiation/sd/symengine_tensor_operations.cc b/source/differentiation/sd/symengine_tensor_operations.cc
index bd2e743623f92947614f641f6b236e11f63359de..e07a3b49113dde3b225d9cfcce2715758b0ada09 100644 (file)
--- a/source/differentiation/sd/symengine_tensor_operations.cc
+++ b/source/differentiation/sd/symengine_tensor_operations.cc
@@ -38,24 +38,23 @@ namespace Differentiation
 
     namespace internal
     {
-      namespace
+      template <int dim>
+      TableIndices<4>
+      make_rank_4_tensor_indices(const unsigned int &idx_i,
+                                 const unsigned int &idx_j)
       {
-        template <int dim>
-        TableIndices<4>
-        make_rank_4_tensor_indices(const unsigned int &idx_i,
-                                   const unsigned int &idx_j)
-        {
-          const TableIndices<2> indices_i(
-            SymmetricTensor<2, dim>::unrolled_to_component_indices(idx_i));
-          const TableIndices<2> indices_j(
-            SymmetricTensor<2, dim>::unrolled_to_component_indices(idx_j));
-          return TableIndices<4>(indices_i[0],
-                                 indices_i[1],
-                                 indices_j[0],
-                                 indices_j[1]);
-        }
-
+        const TableIndices<2> indices_i(
+          SymmetricTensor<2, dim>::unrolled_to_component_indices(idx_i));
+        const TableIndices<2> indices_j(
+          SymmetricTensor<2, dim>::unrolled_to_component_indices(idx_j));
+        return TableIndices<4>(indices_i[0],
+                               indices_i[1],
+                               indices_j[0],
+                               indices_j[1]);
+      }
 
+      namespace
+      {
         template <int rank>
         std::string
         make_index_string(const TableIndices<rank> &indices)

diff --git a/tests/symengine/substitution_maps_tensor_01.cc b/tests/symengine/substitution_maps_tensor_01.cc
new file mode 100644 (file)
index 0000000..c043a6b
--- /dev/null
+++ b/tests/symengine/substitution_maps_tensor_01.cc
@@ -0,0 +1,62 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Check that functions to add tensor variables to substitution maps work
+// correctly.
+
+#include <deal.II/differentiation/sd.h>
+
+#include <string>
+
+#include "../tests.h"
+
+using namespace dealii;
+namespace SD = Differentiation::SD;
+namespace SE = ::SymEngine;
+
+int
+main()
+{
+  initlog();
+  const unsigned int dim = 2;
+
+  Tensor<2, dim>          t;
+  SymmetricTensor<2, dim> st;
+  for (unsigned int i = 0; i < t.n_independent_components; ++i)
+    t[t.unrolled_to_component_indices(i)] = i;
+  for (unsigned int i = 0; i < st.n_independent_components; ++i)
+    st[st.unrolled_to_component_indices(i)] = i;
+
+
+  SD::types::substitution_map substitution_map;
+  SD::add_to_substitution_map(substitution_map,
+                              SD::Expression("x1"),
+                              SD::Expression(1));
+  SD::add_to_substitution_map(substitution_map,
+                              SD::make_tensor_of_symbols<2, dim>("t1"),
+                              t);
+  SD::add_to_substitution_map(
+    substitution_map, SD::make_symmetric_tensor_of_symbols<2, dim>("st1"), st);
+
+  SD::add_to_substitution_map(
+    substitution_map,
+    std::make_pair(SD::make_tensor_of_symbols<2, dim>("t2"), t),
+    std::make_pair(SD::make_symmetric_tensor_of_symbols<2, dim>("st2"), st));
+
+  SD::Utilities::print_substitution_map(deallog, substitution_map);
+
+  deallog << "OK" << std::endl;
+}

diff --git a/tests/symengine/substitution_maps_tensor_01.output b/tests/symengine/substitution_maps_tensor_01.output
new file mode 100644 (file)
index 0000000..c562f76
--- /dev/null
+++ b/tests/symengine/substitution_maps_tensor_01.output
@@ -0,0 +1,17 @@
+
+DEAL::x1 = 1
+t2_01 = 1.0
+t2_00 = 0.0
+t2_10 = 2.0
+t2_11 = 3.0
+t1_00 = 0.0
+t1_01 = 1.0
+t1_11 = 3.0
+t1_10 = 2.0
+st2_11 = 1.0
+st2_00 = 0.0
+st2_01 = 2.0
+st1_11 = 1.0
+st1_00 = 0.0
+st1_01 = 2.0
+OK