Fix Tensor assignment of nested Tensors and provide a function overload of apply_tran...

diff --git a/include/deal.II/base/derivative_form.h b/include/deal.II/base/derivative_form.h
index bb4fa64e3c0ea2eb177ff2b7a03d5797add866d8..99ffd053a548572e56ed57288833b3a7a87305f8 100644 (file)
--- a/include/deal.II/base/derivative_form.h
+++ b/include/deal.II/base/derivative_form.h
@@ -108,14 +108,12 @@ public:
    * Number>. In particular, if order == 1 and the derivative is the Jacobian of
    * $\mathbf F(\mathbf x)$, then Tensor[i] = $\nabla F_i(\mathbf x)$.
    */
-  template <typename Number2>
-  operator Tensor<order + 1, dim, Number2>() const;
+  operator Tensor<order + 1, dim, Number>() const;
 
   /**
    * Converts a DerivativeForm<1, dim, 1, Number> to Tensor<1, dim, Number>.
    */
-  template <typename Number2>
-  operator Tensor<1, dim, Number2>() const;
+  operator Tensor<1, dim, Number>() const;
 
   /**
    * Return the transpose of a rectangular DerivativeForm,
@@ -275,9 +273,8 @@ DerivativeForm<order, dim, spacedim, Number>::operator[](
 
 
 template <int order, int dim, int spacedim, typename Number>
-template <typename Number2>
 inline DerivativeForm<order, dim, spacedim, Number>::
-operator Tensor<1, dim, Number2>() const
+operator Tensor<1, dim, Number>() const
 {
   Assert((1 == spacedim) && (order == 1),
          ExcMessage("Only allowed for spacedim==1."));
@@ -288,13 +285,12 @@ operator Tensor<1, dim, Number2>() const
 
 
 template <int order, int dim, int spacedim, typename Number>
-template <typename Number2>
 inline DerivativeForm<order, dim, spacedim, Number>::
-operator Tensor<order + 1, dim, Number2>() const
+operator Tensor<order + 1, dim, Number>() const
 {
   Assert((dim == spacedim), ExcMessage("Only allowed when dim==spacedim."));
 
-  Tensor<order + 1, dim, Number2> t;
+  Tensor<order + 1, dim, Number> t;
 
   if (dim == spacedim)
     for (unsigned int j = 0; j < dim; ++j)
@@ -470,6 +466,31 @@ apply_transformation(const DerivativeForm<1, dim, spacedim, Number1> &grad_F,
 
 
 
+/**
+ * Similar to the previous apply_transformation(), specialized for the case `dim
+ * == spacedim` where we can return a rank-2 tensor instead of the more general
+ * `DerivativeForm`.
+ * Each row of the result corresponds to one of the rows of @p D_X transformed
+ * by @p grad_F, equivalent to $\mathrm{D\_X} \, \mathrm{grad\_F}^T$ in matrix
+ * notation.
+ *
+ * @relatesalso DerivativeForm
+ */
+// rank=2
+template <int dim, typename Number1, typename Number2>
+inline Tensor<2, dim, typename ProductType<Number1, Number2>::type>
+apply_transformation(const DerivativeForm<1, dim, dim, Number1> &grad_F,
+                     const Tensor<2, dim, Number2> &             D_X)
+{
+  Tensor<2, dim, typename ProductType<Number1, Number2>::type> dest;
+  for (unsigned int i = 0; i < dim; ++i)
+    dest[i] = apply_transformation(grad_F, D_X[i]);
+
+  return dest;
+}
+
+
+
 /**
  * Similar to the previous apply_transformation().
  * Each row of the result corresponds to one of the rows of @p D_X transformed

diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index c9798ba3834d3f5c48acbe6aef6ae98ba94261f6..d14ab9bf5d27b8fa0be9958573bb0fa5dc66bc48 100644 (file)
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -1069,7 +1069,7 @@ constexpr inline DEAL_II_ALWAYS_INLINE
   DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
   Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p)
 {
-  value = internal::NumberType<Number>::value(p);
+  value = internal::NumberType<Number>::value(p.value);
   return *this;
 }
 
@@ -1361,7 +1361,7 @@ template <typename OtherNumber>
 constexpr DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number>::
 operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
 {
-  return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
+  return Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>(values);
 }
 
 

diff --git a/include/deal.II/matrix_free/fe_point_evaluation.h b/include/deal.II/matrix_free/fe_point_evaluation.h
index 4e64cfc2df80ec449ffec76b5f43e08dabea68fb..21969e79c118518b9b8dadf434c20b2015a20e80 100644 (file)
--- a/include/deal.II/matrix_free/fe_point_evaluation.h
+++ b/include/deal.II/matrix_free/fe_point_evaluation.h
@@ -947,11 +947,11 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate(
                     Number>::set_gradient(val_and_grad.second,
                                           j,
                                           unit_gradients[i + j]);
-                  gradients[i + j] = static_cast<gradient_type>(
+                  gradients[i + j] =
                     apply_transformation(mapping_info->get_mapping_data()
                                            .inverse_jacobians[i + j]
                                            .transpose(),
-                                         unit_gradients[i + j]));
+                                         unit_gradients[i + j]);
                 }
             }
         }
@@ -1087,10 +1087,9 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate(
           if (integration_flags & EvaluationFlags::gradients)
             for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
               {
-                gradients[i + j] =
-                  static_cast<gradient_type>(apply_transformation(
-                    mapping_info->get_mapping_data().inverse_jacobians[i + j],
-                    gradients[i + j]));
+                gradients[i + j] = apply_transformation(
+                  mapping_info->get_mapping_data().inverse_jacobians[i + j],
+                  gradients[i + j]);
                 internal::FEPointEvaluation::
                   EvaluatorTypeTraits<dim, n_components, Number>::get_gradient(
                     gradient, j, gradients[i + j]);

diff --git a/tests/tensors/tensor_mixing_types_03.cc b/tests/tensors/tensor_mixing_types_03.cc
new file mode 100644 (file)
index 0000000..69ce653
--- /dev/null
+++ b/tests/tensors/tensor_mixing_types_03.cc
@@ -0,0 +1,40 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Test for mixed Number type assignments and casts of nested Tensors.
+
+#include <deal.II/base/tensor.h>
+
+#include <complex>
+
+#include "../tests.h"
+
+
+int
+main()
+{
+  initlog();
+
+  Tensor<1, 1, Tensor<1, 2, float>>  nested_f;
+  Tensor<1, 1, Tensor<1, 2, double>> nested_d;
+
+  nested_d = nested_f;
+  nested_f = nested_d;
+
+  nested_d = static_cast<Tensor<1, 1, Tensor<1, 2, double>>>(nested_f);
+  nested_f = static_cast<Tensor<1, 1, Tensor<1, 2, float>>>(nested_d);
+
+  deallog << "OK!" << std::endl;
+}

diff --git a/tests/tensors/tensor_mixing_types_03.output b/tests/tensors/tensor_mixing_types_03.output
new file mode 100644 (file)
index 0000000..5cfb783
--- /dev/null
+++ b/tests/tensors/tensor_mixing_types_03.output
@@ -0,0 +1,2 @@
+
+DEAL::OK!