From: Reza Rastak <rastak@stanford.edu>
Date: Wed, 8 May 2019 02:05:33 +0000 (-0700)
Subject: Tensor<rank,dim,Number> is now constexpr
X-Git-Tag: v9.2.0-rc1~1402^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F8064%2Fhead;p=dealii.git

Tensor<rank,dim,Number> is now constexpr

CONSEXPR macro added to Tensor<rank,dim,Number> functions. New constructors added to VectorizedArray<Number>

Tensor free functions all made constexpr

constexpr added to tensor_accessors.h

default initialization for VectorizedArray<Number>

constexpr feature added to Tensor<rank,dim,Number> class

adding DEAL_II_ALWAYS_INLINE back to functions in tensor_accessors.h

Add assignement operator ans default constructor back

Apply suggestions from code review

Typos fixed

Co-Authored-By: Daniel Arndt <arndtd@ornl.gov>

Indentation

Improve check_01_cxx_features.cmake

Fix changelog entry

static constexpr converted to back static const for VectorizedArray

one more constexpr tensor test

test failure fixed

added DEAL_II_CONSTEXPR to three more functions

removed member initialization in VectorizedArray

Avoid warning

constructor added to VectorizedArray<Number>

reverted the changes in vectorization.h

asserting VectorizedArray is POD
---

diff --git a/cmake/checks/check_01_cxx_features.cmake b/cmake/checks/check_01_cxx_features.cmake
index edd3f7996c..2285e206fa 100644
--- a/cmake/checks/check_01_cxx_features.cmake
+++ b/cmake/checks/check_01_cxx_features.cmake
@@ -28,6 +28,7 @@
 #   DEAL_II_HAVE_FP_EXCEPTIONS
 #   DEAL_II_HAVE_COMPLEX_OPERATOR_OVERLOADS
 #
+#   DEAL_II_CONSTEXPR
 #   DEAL_II_FALLTHROUGH
 #
 
@@ -753,6 +754,18 @@ CHECK_CXX_SOURCE_COMPILES(
   "
   DEAL_II_HAVE_CXX14_CONSTEXPR_CAN_CALL_NONCONSTEXPR)
 
+#
+# The macro DEAL_II_CONSTEXPR allows using c++ constexpr features in a portable way.
+# Here we enable it only when a constexpr function can call simple non-constexpr
+# functions. This requirement is probabely very conservative in most cases, but
+# it will prevent breaking builds with certain compilers.
+#
+IF (DEAL_II_HAVE_CXX14_CONSTEXPR_CAN_CALL_NONCONSTEXPR)
+  SET(DEAL_II_CONSTEXPR "constexpr")
+ELSE()
+  SET(DEAL_II_CONSTEXPR " ")
+ENDIF()
+
 #
 # Not all compilers with C++17 support include the new special math
 # functions. Check this separately so that we can use C++17 compilers that don't
diff --git a/doc/news/changes/minor/20190526RezaRastak b/doc/news/changes/minor/20190526RezaRastak
new file mode 100644
index 0000000000..475e038763
--- /dev/null
+++ b/doc/news/changes/minor/20190526RezaRastak
@@ -0,0 +1,4 @@
+New: Tensor<rank, dim, Number> can be constructed and manipulated in constexpr
+setting (if supported by the compiler).
+<br>
+(Reza Rastak, 2019/05/26)
diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in
index d4e79ef91c..a86966b272 100644
--- a/include/deal.II/base/config.h.in
+++ b/include/deal.II/base/config.h.in
@@ -99,6 +99,7 @@
 #cmakedefine DEAL_II_ALWAYS_INLINE @DEAL_II_ALWAYS_INLINE@
 #cmakedefine DEAL_II_RESTRICT @DEAL_II_RESTRICT@
 #cmakedefine DEAL_II_COMPILER_HAS_DIAGNOSTIC_PRAGMA
+#cmakedefine DEAL_II_CONSTEXPR @DEAL_II_CONSTEXPR@
 
 /*
  * A variable to tell if the compiler used in the current compilation process
diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index d0a350fc8d..d9f8f4c676 100644
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -20,6 +20,7 @@
 
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/numbers.h>
+#include <deal.II/base/std_cxx14/utility.h>
 #include <deal.II/base/table_indices.h>
 #include <deal.II/base/template_constraints.h>
 #include <deal.II/base/tensor_accessors.h>
@@ -142,7 +143,7 @@ public:
    *
    * @note This function can also be used in CUDA device code.
    */
-  DEAL_II_CUDA_HOST_DEV
+  constexpr DEAL_II_CUDA_HOST_DEV
   Tensor();
 
   /**
@@ -151,19 +152,19 @@ public:
    * Number.
    */
   template <typename OtherNumber>
-  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
+  constexpr Tensor(const Tensor<0, dim, OtherNumber> &initializer);
 
   /**
    * Constructor, where the data is copied from a C-style array.
    */
   template <typename OtherNumber>
-  Tensor(const OtherNumber &initializer);
+  constexpr Tensor(const OtherNumber &initializer);
 
   /**
    * Return a pointer to the first element of the underlying storage.
    */
-  Number *
-  begin_raw();
+  DEAL_II_CONSTEXPR Number *
+                    begin_raw();
 
   /**
    * Return a const pointer to the first element of the underlying storage.
@@ -174,8 +175,8 @@ public:
   /**
    * Return a pointer to the element past the end of the underlying storage.
    */
-  Number *
-  end_raw();
+  DEAL_II_CONSTEXPR Number *
+                    end_raw();
 
   /**
    * Return a const pointer to the element past the end of the underlying
@@ -193,7 +194,7 @@ public:
    *
    * @note This function can also be used in CUDA device code.
    */
-  DEAL_II_CUDA_HOST_DEV operator Number &();
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator Number &();
 
   /**
    * Return a reference to the encapsulated Number object. Since rank-0
@@ -203,7 +204,7 @@ public:
    *
    * @note This function can also be used in CUDA device code.
    */
-  DEAL_II_CUDA_HOST_DEV operator const Number &() const;
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV operator const Number &() const;
 
   /**
    * Assignment from tensors with different underlying scalar type. This
@@ -211,8 +212,8 @@ public:
    * Number.
    */
   template <typename OtherNumber>
-  Tensor &
-  operator=(const Tensor<0, dim, OtherNumber> &rhs);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator=(const Tensor<0, dim, OtherNumber> &rhs);
 
 #ifdef __INTEL_COMPILER
   /**
@@ -221,8 +222,8 @@ public:
    * copy constructor for Sacado::Rad::ADvar types automatically.
    * See https://github.com/dealii/dealii/pull/5865.
    */
-  Tensor &
-  operator=(const Tensor<0, dim, Number> &rhs);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator=(const Tensor<0, dim, Number> &rhs);
 #endif
 
   /**
@@ -230,14 +231,14 @@ public:
    * that the @p OtherNumber type is convertible to @p Number.
    */
   template <typename OtherNumber>
-  Tensor &
-  operator=(const OtherNumber &d);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator=(const OtherNumber &d);
 
   /**
    * Test for equality of two tensors.
    */
   template <typename OtherNumber>
-  bool
+  DEAL_II_CONSTEXPR bool
   operator==(const Tensor<0, dim, OtherNumber> &rhs) const;
 
   /**
@@ -251,15 +252,15 @@ public:
    * Add another scalar
    */
   template <typename OtherNumber>
-  Tensor &
-  operator+=(const Tensor<0, dim, OtherNumber> &rhs);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator+=(const Tensor<0, dim, OtherNumber> &rhs);
 
   /**
    * Subtract another scalar.
    */
   template <typename OtherNumber>
-  Tensor &
-  operator-=(const Tensor<0, dim, OtherNumber> &rhs);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator-=(const Tensor<0, dim, OtherNumber> &rhs);
 
   /**
    * Multiply the scalar with a <tt>factor</tt>.
@@ -267,15 +268,15 @@ public:
    * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CUDA_HOST_DEV Tensor &
-                        operator*=(const OtherNumber &factor);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator*=(const OtherNumber &factor);
 
   /**
    * Divide the scalar by <tt>factor</tt>.
    */
   template <typename OtherNumber>
-  Tensor &
-  operator/=(const OtherNumber &factor);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator/=(const OtherNumber &factor);
 
   /**
    * Tensor with inverted entries.
@@ -295,7 +296,7 @@ public:
    * and indeed the state where all elements have a zero value is the state
    * right after construction of such an object.
    */
-  void
+  DEAL_II_CONSTEXPR void
   clear();
 
   /**
@@ -303,8 +304,8 @@ public:
    * the absolute squares of all entries. For the present case of rank-1
    * tensors, this equals the usual <tt>l<sub>2</sub></tt> norm of the vector.
    */
-  real_type
-  norm() const;
+  DEAL_II_CONSTEXPR real_type
+                    norm() const;
 
   /**
    * Return the square of the Frobenius-norm of a tensor, i.e. the sum of the
@@ -339,7 +340,7 @@ private:
    * Internal helper function for unroll.
    */
   template <typename OtherNumber>
-  void
+  DEAL_II_CONSTEXPR void
   unroll_recursion(Vector<OtherNumber> &result,
                    unsigned int &       start_index) const;
 
@@ -444,7 +445,7 @@ public:
   /**
    * Constructor, where the data is copied from a C-style array.
    */
-  explicit Tensor(const array_type &initializer);
+  constexpr explicit Tensor(const array_type &initializer);
 
   /**
    * Constructor from tensors with different underlying scalar type. This
@@ -452,13 +453,13 @@ public:
    * Number.
    */
   template <typename OtherNumber>
-  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
+  constexpr Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
 
   /**
    * Constructor that converts from a "tensor of tensors".
    */
   template <typename OtherNumber>
-  Tensor(
+  constexpr Tensor(
     const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer);
 
   /**
@@ -473,7 +474,8 @@ public:
    *
    * @note This function can also be used in CUDA device code.
    */
-  DEAL_II_CUDA_HOST_DEV value_type &operator[](const unsigned int i);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV value_type &
+                                          operator[](const unsigned int i);
 
   /**
    * Read-only access operator.
@@ -486,18 +488,19 @@ public:
   /**
    * Read access using TableIndices <tt>indices</tt>
    */
-  const Number &operator[](const TableIndices<rank_> &indices) const;
+  DEAL_II_CONSTEXPR const Number &
+                          operator[](const TableIndices<rank_> &indices) const;
 
   /**
    * Read and write access using TableIndices <tt>indices</tt>
    */
-  Number &operator[](const TableIndices<rank_> &indices);
+  DEAL_II_CONSTEXPR Number &operator[](const TableIndices<rank_> &indices);
 
   /**
    * Return a pointer to the first element of the underlying storage.
    */
-  Number *
-  begin_raw();
+  DEAL_II_CONSTEXPR Number *
+                    begin_raw();
 
   /**
    * Return a const pointer to the first element of the underlying storage.
@@ -508,8 +511,8 @@ public:
   /**
    * Return a pointer to the element past the end of the underlying storage.
    */
-  Number *
-  end_raw();
+  DEAL_II_CONSTEXPR Number *
+                    end_raw();
 
   /**
    * Return a pointer to the element past the end of the underlying storage.
@@ -523,8 +526,8 @@ public:
    * Number.
    */
   template <typename OtherNumber>
-  Tensor &
-  operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
 
   /**
    * This operator assigns a scalar to a tensor. To avoid confusion with what
@@ -532,14 +535,14 @@ public:
    * value allowed for <tt>d</tt>, allowing the intuitive notation
    * <tt>t=0</tt> to reset all elements of the tensor to zero.
    */
-  Tensor &
-  operator=(const Number &d);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator=(const Number &d);
 
   /**
    * Test for equality of two tensors.
    */
   template <typename OtherNumber>
-  bool
+  DEAL_II_CONSTEXPR bool
   operator==(const Tensor<rank_, dim, OtherNumber> &) const;
 
   /**
@@ -553,15 +556,15 @@ public:
    * Add another tensor.
    */
   template <typename OtherNumber>
-  Tensor &
-  operator+=(const Tensor<rank_, dim, OtherNumber> &);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator+=(const Tensor<rank_, dim, OtherNumber> &);
 
   /**
    * Subtract another tensor.
    */
   template <typename OtherNumber>
-  Tensor &
-  operator-=(const Tensor<rank_, dim, OtherNumber> &);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator-=(const Tensor<rank_, dim, OtherNumber> &);
 
   /**
    * Scale the tensor by <tt>factor</tt>, i.e. multiply all components by
@@ -570,21 +573,21 @@ public:
    * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CUDA_HOST_DEV Tensor &
-                        operator*=(const OtherNumber &factor);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator*=(const OtherNumber &factor);
 
   /**
    * Scale the vector by <tt>1/factor</tt>.
    */
   template <typename OtherNumber>
-  Tensor &
-  operator/=(const OtherNumber &factor);
+  DEAL_II_CONSTEXPR Tensor &
+                    operator/=(const OtherNumber &factor);
 
   /**
    * Unary minus operator. Negate all entries of a tensor.
    */
-  Tensor
-  operator-() const;
+  DEAL_II_CONSTEXPR Tensor
+                    operator-() const;
 
   /**
    * Reset all values to zero.
@@ -598,7 +601,7 @@ public:
    * and indeed the state where all elements have a zero value is the state
    * right after construction of such an object.
    */
-  void
+  DEAL_II_CONSTEXPR void
   clear();
 
   /**
@@ -607,7 +610,7 @@ public:
    * tensors, this equals the usual <tt>l<sub>2</sub></tt> norm of the vector.
    */
 
-  typename numbers::NumberTraits<Number>::real_type
+  DEAL_II_CONSTEXPR typename numbers::NumberTraits<Number>::real_type
   norm() const;
 
   /**
@@ -616,8 +619,9 @@ public:
    *
    * @note This function can also be used in CUDA device code.
    */
-  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
-  norm_square() const;
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV
+    typename numbers::NumberTraits<Number>::real_type
+    norm_square() const;
 
   /**
    * Fill a vector with all tensor elements.
@@ -627,22 +631,22 @@ public:
    * fastest.
    */
   template <typename OtherNumber>
-  void
+  DEAL_II_CONSTEXPR void
   unroll(Vector<OtherNumber> &result) const;
 
   /**
    * Return an unrolled index in the range [0,dim^rank-1] for the element of
    * the tensor indexed by the argument to the function.
    */
-  static unsigned int
+  static DEAL_II_CONSTEXPR unsigned int
   component_to_unrolled_index(const TableIndices<rank_> &indices);
 
   /**
    * Opposite of  component_to_unrolled_index: For an index in the range
    * [0,dim^rank-1], return which set of indices it would correspond to.
    */
-  static TableIndices<rank_>
-  unrolled_to_component_indices(const unsigned int i);
+  static DEAL_II_CONSTEXPR TableIndices<rank_>
+                           unrolled_to_component_indices(const unsigned int i);
 
   /**
    * Determine an estimate for the memory consumption (in bytes) of this
@@ -677,10 +681,18 @@ private:
    * Internal helper function for unroll.
    */
   template <typename OtherNumber>
-  void
+  DEAL_II_CONSTEXPR void
   unroll_recursion(Vector<OtherNumber> &result,
                    unsigned int &       start_index) const;
 
+  /**
+   * This constructor is for internal use. It provides a way
+   *  to create constexpr constructors for Tensor<rank, dim, Number>
+   */
+  template <typename ArrayLike, std::size_t... Indices>
+  constexpr Tensor(const ArrayLike &initializer,
+                   std_cxx14::index_sequence<Indices...>);
+
   /**
    * Allow an arbitrary Tensor to access the underlying values.
    */
@@ -719,8 +731,8 @@ namespace internal
       return t;
     }
 
-    static Tensor<rank, dim, T>
-    value(const T &t)
+    static DEAL_II_CONSTEXPR Tensor<rank, dim, T>
+                             value(const T &t)
     {
       Tensor<rank, dim, T> tmp;
       tmp = t;
@@ -737,16 +749,16 @@ namespace internal
       return t;
     }
 
-    static Tensor<rank, dim, VectorizedArray<T>>
-    value(const T &t)
+    static DEAL_II_CONSTEXPR Tensor<rank, dim, VectorizedArray<T>>
+                             value(const T &t)
     {
       Tensor<rank, dim, VectorizedArray<T>> tmp;
       tmp = internal::NumberType<VectorizedArray<T>>::value(t);
       return tmp;
     }
 
-    static Tensor<rank, dim, VectorizedArray<T>>
-    value(const VectorizedArray<T> &t)
+    static DEAL_II_CONSTEXPR Tensor<rank, dim, VectorizedArray<T>>
+                             value(const VectorizedArray<T> &t)
     {
       Tensor<rank, dim, VectorizedArray<T>> tmp;
       tmp = t;
@@ -760,34 +772,33 @@ namespace internal
 
 
 template <int dim, typename Number>
-DEAL_II_CUDA_HOST_DEV inline Tensor<0, dim, Number>::Tensor()
+constexpr DEAL_II_CUDA_HOST_DEV
+Tensor<0, dim, Number>::Tensor()
   // Some auto-differentiable numbers need explicit
-  // zero initialization.
-  : value(internal::NumberType<Number>::value(0.0))
+  // zero initialization such as adtl::adouble.
+  : Tensor{0.0}
 {}
 
 
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
-{
-  value = internal::NumberType<Number>::value(initializer);
-}
+constexpr Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
+  : value(internal::NumberType<Number>::value(initializer))
+{}
 
 
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<0, dim, Number>::Tensor(const Tensor<0, dim, OtherNumber> &p)
-{
-  value = p.value;
-}
+constexpr Tensor<0, dim, Number>::Tensor(const Tensor<0, dim, OtherNumber> &p)
+  : Tensor{p.value}
+{}
 
 
 
 template <int dim, typename Number>
-inline Number *
+DEAL_II_CONSTEXPR inline Number *
 Tensor<0, dim, Number>::begin_raw()
 {
   return std::addressof(value);
@@ -805,7 +816,7 @@ Tensor<0, dim, Number>::begin_raw() const
 
 
 template <int dim, typename Number>
-inline Number *
+DEAL_II_CONSTEXPR inline Number *
 Tensor<0, dim, Number>::end_raw()
 {
   return begin_raw() + n_independent_components;
@@ -823,7 +834,8 @@ Tensor<0, dim, Number>::end_raw() const
 
 
 template <int dim, typename Number>
-inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &()
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::
+                                               operator Number &()
 {
   // We cannot use Assert inside a CUDA kernel
 #ifndef __CUDA_ARCH__
@@ -835,8 +847,8 @@ inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::operator Number &()
 
 
 template <int dim, typename Number>
-inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::
-                             operator const Number &() const
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::
+                                               operator const Number &() const
 {
   // We cannot use Assert inside a CUDA kernel
 #ifndef __CUDA_ARCH__
@@ -849,7 +861,7 @@ inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p)
 {
   value = internal::NumberType<Number>::value(p);
@@ -859,7 +871,7 @@ Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p)
 
 #ifdef __INTEL_COMPILER
 template <int dim, typename Number>
-inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator=(const Tensor<0, dim, Number> &p)
 {
   value = p.value;
@@ -870,7 +882,7 @@ Tensor<0, dim, Number>::operator=(const Tensor<0, dim, Number> &p)
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator=(const OtherNumber &d)
 {
   value = internal::NumberType<Number>::value(d);
@@ -880,7 +892,7 @@ Tensor<0, dim, Number>::operator=(const OtherNumber &d)
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline bool
+DEAL_II_CONSTEXPR inline bool
 Tensor<0, dim, Number>::operator==(const Tensor<0, dim, OtherNumber> &p) const
 {
 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
@@ -906,7 +918,7 @@ Tensor<0, dim, Number>::operator!=(const Tensor<0, dim, OtherNumber> &p) const
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator+=(const Tensor<0, dim, OtherNumber> &p)
 {
   value += p.value;
@@ -916,7 +928,7 @@ Tensor<0, dim, Number>::operator+=(const Tensor<0, dim, OtherNumber> &p)
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator-=(const Tensor<0, dim, OtherNumber> &p)
 {
   value -= p.value;
@@ -930,7 +942,7 @@ namespace internal
   namespace ComplexWorkaround
   {
     template <typename Number, typename OtherNumber>
-    inline DEAL_II_CUDA_HOST_DEV void
+    DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV void
     multiply_assign_scalar(Number &val, const OtherNumber &s)
     {
       val *= s;
@@ -938,7 +950,7 @@ namespace internal
 
 #ifdef __CUDA_ARCH__
     template <typename Number, typename OtherNumber>
-    inline DEAL_II_CUDA_HOST_DEV void
+    DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV void
     multiply_assign_scalar(std::complex<Number> &, const OtherNumber &)
     {
       printf("This function is not implemented for std::complex<Number>!\n");
@@ -951,7 +963,7 @@ namespace internal
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
 {
   internal::ComplexWorkaround::multiply_assign_scalar(value, s);
@@ -962,7 +974,7 @@ Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
 {
   value /= s;
@@ -979,7 +991,7 @@ Tensor<0, dim, Number>::operator-() const
 
 
 template <int dim, typename Number>
-inline typename Tensor<0, dim, Number>::real_type
+DEAL_II_CONSTEXPR inline typename Tensor<0, dim, Number>::real_type
 Tensor<0, dim, Number>::norm() const
 {
   Assert(dim != 0,
@@ -1003,7 +1015,7 @@ inline typename Tensor<0, dim, Number>::real_type DEAL_II_CUDA_HOST_DEV
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline void
+DEAL_II_CONSTEXPR inline void
 Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
                                          unsigned int &       index) const
 {
@@ -1015,7 +1027,7 @@ Tensor<0, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
 
 
 template <int dim, typename Number>
-inline void
+DEAL_II_CONSTEXPR inline void
 Tensor<0, dim, Number>::clear()
 {
   // Some auto-differentiable numbers need explicit
@@ -1035,42 +1047,45 @@ Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
 
 /*-------------------- Inline functions: Tensor<rank,dim> --------------------*/
 
-
 template <int rank_, int dim, typename Number>
-inline DEAL_II_ALWAYS_INLINE
-Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
+template <typename ArrayLike, std::size_t... indices>
+DEAL_II_ALWAYS_INLINE constexpr Tensor<rank_, dim, Number>::Tensor(
+  const ArrayLike &initializer,
+  std_cxx14::index_sequence<indices...>)
+  : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
 {
-  for (unsigned int i = 0; i < dim; ++i)
-    values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
+  static_assert(sizeof...(indices) == dim,
+                "dim should match the number of indices");
 }
 
 
+template <int rank_, int dim, typename Number>
+DEAL_II_ALWAYS_INLINE constexpr Tensor<rank_, dim, Number>::Tensor(
+  const array_type &initializer)
+  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
+{}
+
+
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE
-Tensor<rank_, dim, Number>::Tensor(
+DEAL_II_ALWAYS_INLINE constexpr Tensor<rank_, dim, Number>::Tensor(
   const Tensor<rank_, dim, OtherNumber> &initializer)
-{
-  for (unsigned int i = 0; i != dim; ++i)
-    values[i] = Tensor<rank_ - 1, dim, Number>(initializer[i]);
-}
+  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
+{}
 
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE
-Tensor<rank_, dim, Number>::Tensor(
+DEAL_II_ALWAYS_INLINE constexpr Tensor<rank_, dim, Number>::Tensor(
   const Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>> &initializer)
-{
-  for (unsigned int i = 0; i < dim; ++i)
-    values[i] = initializer[i];
-}
+  : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
+{}
 
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-constexpr DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number>::
-                                operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
+DEAL_II_ALWAYS_INLINE constexpr Tensor<rank_, dim, Number>::
+operator Tensor<1, dim, Tensor<rank_ - 1, dim, OtherNumber>>() const
 {
   return Tensor<1, dim, Tensor<rank_ - 1, dim, Number>>(values);
 }
@@ -1082,10 +1097,11 @@ namespace internal
   namespace TensorSubscriptor
   {
     template <typename ArrayElementType, int dim>
-    inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV ArrayElementType &
-                                                       subscript(ArrayElementType * values,
-                                                                 const unsigned int i,
-                                                                 std::integral_constant<int, dim>)
+    DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
+      DEAL_II_CUDA_HOST_DEV ArrayElementType &
+                            subscript(ArrayElementType * values,
+                                      const unsigned int i,
+                                      std::integral_constant<int, dim>)
     {
       // We cannot use Assert in a CUDA kernel
 #ifndef __CUDA_ARCH__
@@ -1094,9 +1110,19 @@ namespace internal
       return values[i];
     }
 
+    // The variables within this struct will be referenced in the next function.
+    // It is a workaround that allows returning a reference to a static variable
+    // while allowing constexpr evaluation of the function.
+    // It has to be defined outside the function because constexpr functions
+    // cannot define static variables
+    template <typename ArrayElementType>
+    struct Uninitialized
+    {
+      static ArrayElementType value;
+    };
 
     template <typename ArrayElementType>
-    ArrayElementType &
+    DEAL_II_CONSTEXPR inline ArrayElementType &
     subscript(ArrayElementType *,
               const unsigned int,
               std::integral_constant<int, 0>)
@@ -1105,15 +1131,14 @@ namespace internal
         false,
         ExcMessage(
           "Cannot access elements of an object of type Tensor<rank,0,Number>."));
-      static ArrayElementType t;
-      return t;
+      return Uninitialized<ArrayElementType>::value;
     }
   } // namespace TensorSubscriptor
 } // namespace internal
 
 
 template <int rank_, int dim, typename Number>
-inline DEAL_II_ALWAYS_INLINE                       DEAL_II_CUDA_HOST_DEV
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE     DEAL_II_CUDA_HOST_DEV
   typename Tensor<rank_, dim, Number>::value_type &Tensor<rank_, dim, Number>::
                                                    operator[](const unsigned int i)
 {
@@ -1123,9 +1148,11 @@ inline DEAL_II_ALWAYS_INLINE                       DEAL_II_CUDA_HOST_DEV
 
 
 template <int rank_, int dim, typename Number>
-constexpr DEAL_II_ALWAYS_INLINE
-    DEAL_II_CUDA_HOST_DEV const typename Tensor<rank_, dim, Number>::value_type &
-    Tensor<rank_, dim, Number>::operator[](const unsigned int i) const
+DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV const typename Tensor<
+  rank_,
+  dim,
+  Number>::value_type &Tensor<rank_, dim, Number>::
+                       operator[](const unsigned int i) const
 {
   return dealii::internal::TensorSubscriptor::subscript(
     values, i, std::integral_constant<int, dim>());
@@ -1133,8 +1160,8 @@ constexpr DEAL_II_ALWAYS_INLINE
 
 
 template <int rank_, int dim, typename Number>
-inline const Number &Tensor<rank_, dim, Number>::
-                     operator[](const TableIndices<rank_> &indices) const
+DEAL_II_CONSTEXPR inline const Number &Tensor<rank_, dim, Number>::
+                                       operator[](const TableIndices<rank_> &indices) const
 {
   Assert(dim != 0,
          ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
@@ -1145,8 +1172,8 @@ inline const Number &Tensor<rank_, dim, Number>::
 
 
 template <int rank_, int dim, typename Number>
-inline Number &Tensor<rank_, dim, Number>::
-               operator[](const TableIndices<rank_> &indices)
+DEAL_II_CONSTEXPR inline Number &Tensor<rank_, dim, Number>::
+                                 operator[](const TableIndices<rank_> &indices)
 {
   Assert(dim != 0,
          ExcMessage("Cannot access an object of type Tensor<rank_,0,Number>"));
@@ -1157,7 +1184,7 @@ inline Number &Tensor<rank_, dim, Number>::
 
 
 template <int rank_, int dim, typename Number>
-inline Number *
+DEAL_II_CONSTEXPR inline Number *
 Tensor<rank_, dim, Number>::begin_raw()
 {
   return std::addressof(
@@ -1177,7 +1204,7 @@ Tensor<rank_, dim, Number>::begin_raw() const
 
 
 template <int rank_, int dim, typename Number>
-inline Number *
+DEAL_II_CONSTEXPR inline Number *
 Tensor<rank_, dim, Number>::end_raw()
 {
   return begin_raw() + n_independent_components;
@@ -1196,7 +1223,7 @@ Tensor<rank_, dim, Number>::end_raw() const
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
 Tensor<rank_, dim, Number>::operator=(const Tensor<rank_, dim, OtherNumber> &t)
 {
   if (dim > 0)
@@ -1206,7 +1233,7 @@ Tensor<rank_, dim, Number>::operator=(const Tensor<rank_, dim, OtherNumber> &t)
 
 
 template <int rank_, int dim, typename Number>
-inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<rank_, dim, Number> &
 Tensor<rank_, dim, Number>::operator=(const Number &d)
 {
   Assert(numbers::value_is_zero(d),
@@ -1221,7 +1248,7 @@ Tensor<rank_, dim, Number>::operator=(const Number &d)
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline bool
+DEAL_II_CONSTEXPR inline bool
 Tensor<rank_, dim, Number>::
 operator==(const Tensor<rank_, dim, OtherNumber> &p) const
 {
@@ -1239,7 +1266,7 @@ operator==(const Tensor<rank_, dim, OtherNumber> &p) const
 // implement this function here
 template <>
 template <>
-inline bool
+DEAL_II_CONSTEXPR inline bool
 Tensor<1, 0, double>::operator==(const Tensor<1, 0, double> &) const
 {
   return true;
@@ -1258,7 +1285,7 @@ operator!=(const Tensor<rank_, dim, OtherNumber> &p) const
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<rank_, dim, Number> &
+DEAL_II_CONSTEXPR inline Tensor<rank_, dim, Number> &
 Tensor<rank_, dim, Number>::operator+=(const Tensor<rank_, dim, OtherNumber> &p)
 {
   for (unsigned int i = 0; i < dim; ++i)
@@ -1269,7 +1296,7 @@ Tensor<rank_, dim, Number>::operator+=(const Tensor<rank_, dim, OtherNumber> &p)
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<rank_, dim, Number> &
+DEAL_II_CONSTEXPR inline Tensor<rank_, dim, Number> &
 Tensor<rank_, dim, Number>::operator-=(const Tensor<rank_, dim, OtherNumber> &p)
 {
   for (unsigned int i = 0; i < dim; ++i)
@@ -1280,7 +1307,7 @@ Tensor<rank_, dim, Number>::operator-=(const Tensor<rank_, dim, OtherNumber> &p)
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
 Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
 {
   for (unsigned int i = 0; i < dim; ++i)
@@ -1291,7 +1318,7 @@ Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline Tensor<rank_, dim, Number> &
+DEAL_II_CONSTEXPR inline Tensor<rank_, dim, Number> &
 Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
 {
   for (unsigned int i = 0; i < dim; ++i)
@@ -1301,7 +1328,7 @@ Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
 
 
 template <int rank_, int dim, typename Number>
-inline Tensor<rank_, dim, Number>
+DEAL_II_CONSTEXPR inline Tensor<rank_, dim, Number>
 Tensor<rank_, dim, Number>::operator-() const
 {
   Tensor<rank_, dim, Number> tmp;
@@ -1314,7 +1341,7 @@ Tensor<rank_, dim, Number>::operator-() const
 
 
 template <int rank_, int dim, typename Number>
-inline typename numbers::NumberTraits<Number>::real_type
+DEAL_II_CONSTEXPR inline typename numbers::NumberTraits<Number>::real_type
 Tensor<rank_, dim, Number>::norm() const
 {
   return std::sqrt(norm_square());
@@ -1322,8 +1349,9 @@ Tensor<rank_, dim, Number>::norm() const
 
 
 template <int rank_, int dim, typename Number>
-inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
-Tensor<rank_, dim, Number>::norm_square() const
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV
+  typename numbers::NumberTraits<Number>::real_type
+  Tensor<rank_, dim, Number>::norm_square() const
 {
   typename numbers::NumberTraits<Number>::real_type s = internal::NumberType<
     typename numbers::NumberTraits<Number>::real_type>::value(0.0);
@@ -1336,7 +1364,7 @@ Tensor<rank_, dim, Number>::norm_square() const
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline void
+DEAL_II_CONSTEXPR inline void
 Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
 {
   AssertDimension(result.size(),
@@ -1349,7 +1377,7 @@ Tensor<rank_, dim, Number>::unroll(Vector<OtherNumber> &result) const
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-inline void
+DEAL_II_CONSTEXPR inline void
 Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
                                              unsigned int &       index) const
 {
@@ -1359,7 +1387,7 @@ Tensor<rank_, dim, Number>::unroll_recursion(Vector<OtherNumber> &result,
 
 
 template <int rank_, int dim, typename Number>
-inline unsigned int
+DEAL_II_CONSTEXPR inline unsigned int
 Tensor<rank_, dim, Number>::component_to_unrolled_index(
   const TableIndices<rank_> &indices)
 {
@@ -1412,7 +1440,7 @@ namespace internal
 
 
 template <int rank_, int dim, typename Number>
-inline TableIndices<rank_>
+DEAL_II_CONSTEXPR inline TableIndices<rank_>
 Tensor<rank_, dim, Number>::unrolled_to_component_indices(const unsigned int i)
 {
   Assert(i < n_independent_components,
@@ -1433,7 +1461,7 @@ Tensor<rank_, dim, Number>::unrolled_to_component_indices(const unsigned int i)
 
 
 template <int rank_, int dim, typename Number>
-inline void
+DEAL_II_CONSTEXPR inline void
 Tensor<rank_, dim, Number>::clear()
 {
   for (unsigned int i = 0; i < dim; ++i)
@@ -1619,7 +1647,7 @@ constexpr DEAL_II_ALWAYS_INLINE
  * @relatesalso Tensor
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
   Tensor<rank,
          dim,
          typename ProductType<Number,
@@ -1665,12 +1693,12 @@ constexpr DEAL_II_ALWAYS_INLINE
  * @relatesalso Tensor
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-inline Tensor<
-  rank,
-  dim,
-  typename ProductType<Number,
-                       typename EnableIfScalar<OtherNumber>::type>::type>
-operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
+  Tensor<rank,
+         dim,
+         typename ProductType<Number,
+                              typename EnableIfScalar<OtherNumber>::type>::type>
+  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
 {
   // recurse over the base objects
   Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tt;
@@ -1688,7 +1716,7 @@ operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
  * @relatesalso Tensor
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
   Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
   operator+(const Tensor<rank, dim, Number> &     p,
             const Tensor<rank, dim, OtherNumber> &q)
@@ -1710,7 +1738,7 @@ inline DEAL_II_ALWAYS_INLINE
  * @relatesalso Tensor
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
   Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
   operator-(const Tensor<rank, dim, Number> &     p,
             const Tensor<rank, dim, OtherNumber> &q)
@@ -1759,7 +1787,7 @@ template <int rank_1,
           int dim,
           typename Number,
           typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
   typename Tensor<rank_1 + rank_2 - 2,
                   dim,
                   typename ProductType<Number, OtherNumber>::type>::tensor_type
@@ -1816,7 +1844,7 @@ template <int index_1,
           int dim,
           typename Number,
           typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
   typename Tensor<rank_1 + rank_2 - 2,
                   dim,
                   typename ProductType<Number, OtherNumber>::type>::tensor_type
@@ -1890,7 +1918,7 @@ template <int index_1,
           int dim,
           typename Number,
           typename OtherNumber>
-inline
+DEAL_II_CONSTEXPR inline
   typename Tensor<rank_1 + rank_2 - 4,
                   dim,
                   typename ProductType<Number, OtherNumber>::type>::tensor_type
@@ -1972,9 +2000,10 @@ inline
  * @author Matthias Maier, 2015
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
-scalar_product(const Tensor<rank, dim, Number> &     left,
-               const Tensor<rank, dim, OtherNumber> &right)
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
+  typename ProductType<Number, OtherNumber>::type
+  scalar_product(const Tensor<rank, dim, Number> &     left,
+                 const Tensor<rank, dim, OtherNumber> &right)
 {
   typename ProductType<Number, OtherNumber>::type result;
   TensorAccessors::contract<rank, rank, rank, dim>(result, left, right);
@@ -2009,10 +2038,11 @@ template <template <int, int, typename> class TensorT1,
           typename T1,
           typename T2,
           typename T3>
-typename ProductType<T1, typename ProductType<T2, T3>::type>::type
-contract3(const TensorT1<rank_1, dim, T1> &         left,
-          const TensorT2<rank_1 + rank_2, dim, T2> &middle,
-          const TensorT3<rank_2, dim, T3> &         right)
+DEAL_II_CONSTEXPR
+  typename ProductType<T1, typename ProductType<T2, T3>::type>::type
+  contract3(const TensorT1<rank_1, dim, T1> &         left,
+            const TensorT2<rank_1 + rank_2, dim, T2> &middle,
+            const TensorT3<rank_2, dim, T3> &         right)
 {
   using return_type =
     typename ProductType<T1, typename ProductType<T2, T3>::type>::type;
@@ -2038,7 +2068,7 @@ template <int rank_1,
           int dim,
           typename Number,
           typename OtherNumber>
-inline DEAL_II_ALWAYS_INLINE
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
   Tensor<rank_1 + rank_2, dim, typename ProductType<Number, OtherNumber>::type>
   outer_product(const Tensor<rank_1, dim, Number> &     src1,
                 const Tensor<rank_2, dim, OtherNumber> &src2)
@@ -2071,8 +2101,8 @@ inline DEAL_II_ALWAYS_INLINE
  * @author Guido Kanschat, 2001
  */
 template <int dim, typename Number>
-inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
-                             cross_product_2d(const Tensor<1, dim, Number> &src)
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
+                                               cross_product_2d(const Tensor<1, dim, Number> &src)
 {
   Assert(dim == 2, ExcInternalError());
 
@@ -2096,7 +2126,7 @@ inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
  * @author Guido Kanschat, 2001
  */
 template <int dim, typename Number1, typename Number2>
-inline DEAL_II_ALWAYS_INLINE
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
   Tensor<1, dim, typename ProductType<Number1, Number2>::type>
   cross_product_3d(const Tensor<1, dim, Number1> &src1,
                    const Tensor<1, dim, Number2> &src2)
@@ -2127,7 +2157,7 @@ inline DEAL_II_ALWAYS_INLINE
  * @author Wolfgang Bangerth, 2009
  */
 template <int dim, typename Number>
-inline Number
+DEAL_II_CONSTEXPR inline Number
 determinant(const Tensor<2, dim, Number> &t)
 {
   // Compute the determinant using the Laplace expansion of the
@@ -2170,8 +2200,8 @@ determinant(const Tensor<2, 1, Number> &t)
  * @author Wolfgang Bangerth, 2001
  */
 template <int dim, typename Number>
-inline DEAL_II_ALWAYS_INLINE Number
-                             trace(const Tensor<2, dim, Number> &d)
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Number
+                                               trace(const Tensor<2, dim, Number> &d)
 {
   Number t = d[0][0];
   for (unsigned int i = 1; i < dim; ++i)
@@ -2190,7 +2220,7 @@ inline DEAL_II_ALWAYS_INLINE Number
  * @author Wolfgang Bangerth, 2000
  */
 template <int dim, typename Number>
-inline Tensor<2, dim, Number>
+DEAL_II_CONSTEXPR inline Tensor<2, dim, Number>
 invert(const Tensor<2, dim, Number> &)
 {
   Number return_tensor[dim][dim];
@@ -2207,19 +2237,19 @@ invert(const Tensor<2, dim, Number> &)
 #ifndef DOXYGEN
 
 template <typename Number>
-inline Tensor<2, 1, Number>
+DEAL_II_CONSTEXPR inline Tensor<2, 1, Number>
 invert(const Tensor<2, 1, Number> &t)
 {
-  Number return_tensor[1][1];
+  Tensor<2, 1, Number> return_tensor;
 
   return_tensor[0][0] = internal::NumberType<Number>::value(1.0 / t[0][0]);
 
-  return Tensor<2, 1, Number>(return_tensor);
+  return return_tensor;
 }
 
 
 template <typename Number>
-inline Tensor<2, 2, Number>
+DEAL_II_CONSTEXPR inline Tensor<2, 2, Number>
 invert(const Tensor<2, 2, Number> &t)
 {
   Tensor<2, 2, Number> return_tensor;
@@ -2239,7 +2269,7 @@ invert(const Tensor<2, 2, Number> &t)
 
 
 template <typename Number>
-inline Tensor<2, 3, Number>
+DEAL_II_CONSTEXPR inline Tensor<2, 3, Number>
 invert(const Tensor<2, 3, Number> &t)
 {
   Tensor<2, 3, Number> return_tensor;
@@ -2284,8 +2314,8 @@ invert(const Tensor<2, 3, Number> &t)
  * @author Wolfgang Bangerth, 2002
  */
 template <int dim, typename Number>
-inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
-                             transpose(const Tensor<2, dim, Number> &t)
+DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
+                                               transpose(const Tensor<2, dim, Number> &t)
 {
   Tensor<2, dim, Number> tt;
   for (unsigned int i = 0; i < dim; ++i)
diff --git a/include/deal.II/base/tensor_accessors.h b/include/deal.II/base/tensor_accessors.h
index df64e1ac21..b465113e1c 100644
--- a/include/deal.II/base/tensor_accessors.h
+++ b/include/deal.II/base/tensor_accessors.h
@@ -188,8 +188,8 @@ namespace TensorAccessors
    * @author Matthias Maier, 2015
    */
   template <int index, int rank, typename T>
-  inline DEAL_II_ALWAYS_INLINE internal::ReorderedIndexView<index, rank, T>
-                               reordered_index_view(T &t)
+  constexpr DEAL_II_ALWAYS_INLINE internal::ReorderedIndexView<index, rank, T>
+                                  reordered_index_view(T &t)
   {
     static_assert(0 <= index && index < rank,
                   "The specified index must lie within the range [0,rank)");
@@ -222,7 +222,7 @@ namespace TensorAccessors
    * @author Matthias Maier, 2015
    */
   template <int rank, typename T, typename ArrayType>
-  typename ReturnType<rank, T>::value_type &
+  constexpr typename ReturnType<rank, T>::value_type &
   extract(T &t, const ArrayType &indices)
   {
     return internal::ExtractHelper<0, rank>::template extract<T, ArrayType>(
@@ -277,7 +277,7 @@ namespace TensorAccessors
             typename T1,
             typename T2,
             typename T3>
-  inline DEAL_II_ALWAYS_INLINE void
+  DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE void
   contract(T1 &result, const T2 &left, const T3 &right)
   {
     static_assert(rank_1 >= no_contr,
@@ -331,7 +331,7 @@ namespace TensorAccessors
             typename T2,
             typename T3,
             typename T4>
-  T1
+  constexpr T1
   contract3(const T2 &left, const T3 &middle, const T4 &right)
   {
     return internal::Contract3<rank_1, rank_2, dim>::
@@ -412,7 +412,7 @@ namespace TensorAccessors
     class ReorderedIndexView
     {
     public:
-      ReorderedIndexView(typename ReferenceType<T>::type t)
+      constexpr ReorderedIndexView(typename ReferenceType<T>::type t)
         : t_(t)
       {}
 
@@ -421,7 +421,8 @@ namespace TensorAccessors
                                             typename ValueType<T>::value_type>;
 
       // Recurse by applying index j directly:
-      inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
+      constexpr DEAL_II_ALWAYS_INLINE value_type
+                                      operator[](unsigned int j) const
       {
         return value_type(t_[j]);
       }
@@ -445,13 +446,14 @@ namespace TensorAccessors
     class ReorderedIndexView<0, rank, T>
     {
     public:
-      ReorderedIndexView(typename ReferenceType<T>::type t)
+      constexpr ReorderedIndexView(typename ReferenceType<T>::type t)
         : t_(t)
       {}
 
       using value_type = StoreIndex<rank - 1, internal::Identity<T>>;
 
-      inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
+      constexpr DEAL_II_ALWAYS_INLINE value_type
+                                      operator[](unsigned int j) const
       {
         return value_type(Identity<T>(t_), j);
       }
@@ -467,14 +469,15 @@ namespace TensorAccessors
     class ReorderedIndexView<0, 1, T>
     {
     public:
-      ReorderedIndexView(typename ReferenceType<T>::type t)
+      constexpr ReorderedIndexView(typename ReferenceType<T>::type t)
         : t_(t)
       {}
 
       using value_type =
         typename ReferenceType<typename ValueType<T>::value_type>::type;
 
-      inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
+      constexpr DEAL_II_ALWAYS_INLINE value_type
+                                      operator[](unsigned int j) const
       {
         return t_[j];
       }
@@ -492,13 +495,13 @@ namespace TensorAccessors
     class Identity
     {
     public:
-      Identity(typename ReferenceType<T>::type t)
+      constexpr Identity(typename ReferenceType<T>::type t)
         : t_(t)
       {}
 
       using return_type = typename ValueType<T>::value_type;
 
-      inline DEAL_II_ALWAYS_INLINE typename ReferenceType<return_type>::type
+      constexpr DEAL_II_ALWAYS_INLINE typename ReferenceType<return_type>::type
       apply(unsigned int j) const
       {
         return t_[j];
@@ -520,14 +523,15 @@ namespace TensorAccessors
     class StoreIndex
     {
     public:
-      StoreIndex(S s, int i)
+      constexpr StoreIndex(S s, int i)
         : s_(s)
         , i_(i)
       {}
 
       using value_type = StoreIndex<rank - 1, StoreIndex<rank, S>>;
 
-      inline DEAL_II_ALWAYS_INLINE value_type operator[](unsigned int j) const
+      constexpr DEAL_II_ALWAYS_INLINE value_type
+                                      operator[](unsigned int j) const
       {
         return value_type(*this, j);
       }
@@ -535,7 +539,7 @@ namespace TensorAccessors
       using return_type =
         typename ValueType<typename S::return_type>::value_type;
 
-      inline typename ReferenceType<return_type>::type
+      constexpr typename ReferenceType<return_type>::type
       apply(unsigned int j) const
       {
         return s_.apply(j)[i_];
@@ -554,7 +558,7 @@ namespace TensorAccessors
     class StoreIndex<1, S>
     {
     public:
-      StoreIndex(S s, int i)
+      constexpr StoreIndex(S s, int i)
         : s_(s)
         , i_(i)
       {}
@@ -563,7 +567,8 @@ namespace TensorAccessors
         typename ValueType<typename S::return_type>::value_type;
       using value_type = return_type;
 
-      inline DEAL_II_ALWAYS_INLINE return_type &operator[](unsigned int j) const
+      constexpr DEAL_II_ALWAYS_INLINE return_type &
+                                      operator[](unsigned int j) const
       {
         return s_.apply(j)[i_];
       }
@@ -585,7 +590,7 @@ namespace TensorAccessors
     struct ExtractHelper
     {
       template <typename T, typename ArrayType>
-      inline static typename ReturnType<rank - position, T>::value_type &
+      constexpr static typename ReturnType<rank - position, T>::value_type &
       extract(T &t, const ArrayType &indices)
       {
         return ExtractHelper<position + 1, rank>::template extract<
@@ -600,7 +605,7 @@ namespace TensorAccessors
     struct ExtractHelper<rank, rank>
     {
       template <typename T, typename ArrayType>
-      inline static T &
+      constexpr static T &
       extract(T &t, const ArrayType &)
       {
         return t;
@@ -628,7 +633,7 @@ namespace TensorAccessors
     {
     public:
       template <typename T1, typename T2, typename T3>
-      inline DEAL_II_ALWAYS_INLINE static void
+      DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE static void
       contract(T1 &result, const T2 &left, const T3 &right)
       {
         for (unsigned int i = 0; i < dim; ++i)
@@ -658,7 +663,7 @@ namespace TensorAccessors
     {
     public:
       template <typename T1, typename T2, typename T3>
-      inline DEAL_II_ALWAYS_INLINE static void
+      DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE static void
       contract(T1 &result, const T2 &left, const T3 &right)
       {
         for (unsigned int i = 0; i < dim; ++i)
@@ -688,7 +693,7 @@ namespace TensorAccessors
     {
     public:
       template <typename T1, typename T2, typename T3>
-      inline DEAL_II_ALWAYS_INLINE static void
+      DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE static void
       contract(T1 &result, const T2 &left, const T3 &right)
       {
         result = Contract2<no_contr, dim>::template contract2<T1>(left, right);
@@ -704,7 +709,7 @@ namespace TensorAccessors
     {
     public:
       template <typename T1, typename T2, typename T3>
-      inline DEAL_II_ALWAYS_INLINE static T1
+      DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE static T1
       contract2(const T2 &left, const T3 &right)
       {
         // Some auto-differentiable numbers need explicit
@@ -726,7 +731,7 @@ namespace TensorAccessors
     {
     public:
       template <typename T1, typename T2, typename T3>
-      inline DEAL_II_ALWAYS_INLINE static T1
+      constexpr DEAL_II_ALWAYS_INLINE static T1
       contract2(const T2 &left, const T3 &right)
       {
         return left * right;
@@ -753,7 +758,7 @@ namespace TensorAccessors
     {
     public:
       template <typename T1, typename T2, typename T3, typename T4>
-      static inline T1
+      DEAL_II_CONSTEXPR static inline T1
       contract3(const T2 &left, const T3 &middle, const T4 &right)
       {
         // Some auto-differentiable numbers need explicit
@@ -783,7 +788,7 @@ namespace TensorAccessors
     {
     public:
       template <typename T1, typename T2, typename T3, typename T4>
-      static inline T1
+      DEAL_II_CONSTEXPR static inline T1
       contract3(const T2 &left, const T3 &middle, const T4 &right)
       {
         // Some auto-differentiable numbers need explicit
@@ -806,7 +811,7 @@ namespace TensorAccessors
     {
     public:
       template <typename T1, typename T2, typename T3, typename T4>
-      static inline T1
+      constexpr static T1
       contract3(const T2 &left, const T3 &middle, const T4 &right)
       {
         return left * middle * right;
diff --git a/source/base/vectorization.cc b/source/base/vectorization.cc
index 9c0774e5b7..ebd3c66961 100644
--- a/source/base/vectorization.cc
+++ b/source/base/vectorization.cc
@@ -20,6 +20,19 @@ DEAL_II_NAMESPACE_OPEN
 #if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && !defined(DEAL_II_MSVC)
 const unsigned int VectorizedArray<double>::n_array_elements;
 const unsigned int VectorizedArray<float>::n_array_elements;
+
+// VectorizedArray must be a POD (plain old data) type to make sure it
+// can use maximum level of compiler optimization.
+// A type is POD if it has standard layout (similar to a C struct)
+// and it is trivial (can be statically default initialized)
+// Here, the trait std::is_pod cannot be used because it is deprecated
+// in C++20.
+static_assert(std::is_standard_layout<VectorizedArray<double>>::value &&
+                std::is_trivial<VectorizedArray<double>>::value,
+              "VectorizedArray<double> must be a POD type");
+static_assert(std::is_standard_layout<VectorizedArray<float>>::value &&
+                std::is_trivial<VectorizedArray<float>>::value,
+              "VectorizedArray<float> must be a POD type");
 #endif
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/base/vectorization_07.cc b/tests/base/vectorization_07.cc
deleted file mode 100644
index 17ae0138a3..0000000000
--- a/tests/base/vectorization_07.cc
+++ /dev/null
@@ -1,36 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// Copyright (C) 2012 - 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.
-//
-// ---------------------------------------------------------------------
-
-#include <deal.II/base/vectorization.h>
-
-#include <type_traits>
-
-#include "../tests.h"
-
-int
-main()
-{
-  initlog();
-
-  if (!std::is_pod<VectorizedArray<double>>::value)
-    deallog.get_file_stream()
-      << "Not OK because VectorizedArray<double> is not POD!" << std::endl;
-
-  if (!std::is_trivial<VectorizedArray<double>>::value)
-    deallog.get_file_stream()
-      << "Not OK because VectorizedArray<double> is not trivial!" << std::endl;
-
-  deallog.get_file_stream() << "OK" << std::endl;
-}
diff --git a/tests/base/vectorization_07.output b/tests/base/vectorization_07.output
deleted file mode 100644
index 6acbd104b1..0000000000
--- a/tests/base/vectorization_07.output
+++ /dev/null
@@ -1,2 +0,0 @@
-
-OK
diff --git a/tests/tensors/constexpr_tensor.cc b/tests/tensors/constexpr_tensor.cc
new file mode 100644
index 0000000000..1afa257a92
--- /dev/null
+++ b/tests/tensors/constexpr_tensor.cc
@@ -0,0 +1,95 @@
+// ---------------------------------------------------------------------
+//
+// 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Create and copy tensors in constexpr setting
+
+#include <deal.II/base/tensor.h>
+
+#include "../tests.h"
+
+template <int rank, int dim, typename Number>
+void
+test_consexpr_tensor_constructors()
+{
+  constexpr dealii::Tensor<rank, dim, Number> A;
+  constexpr dealii::Tensor<rank, dim, Number> B(A);
+  constexpr dealii::Tensor<rank, dim, Number> C = A;
+  deallog << " Tensor<" << rank << "," << dim << ">" << std::endl;
+  deallog << A << std::endl;
+  deallog << B << std::endl;
+  deallog << C << std::endl;
+}
+
+DEAL_II_CONSTEXPR double
+tensor_rank0_constexpr(const double value)
+{
+  const Tensor<0, 2> a{value};
+  return 2. * a;
+}
+
+DEAL_II_CONSTEXPR double
+tensor_rank2_constexpr()
+{
+  constexpr double       a[3][3] = {{4., 2., 0.}, {2., 3., 0.}, {0., 0., 5.}};
+  constexpr Tensor<2, 3> A(a);
+  DEAL_II_CONSTEXPR const auto det = determinant(A);
+  DEAL_II_CONSTEXPR const auto tr  = trace(transpose(A) * A);
+  return det + tr;
+}
+
+int
+main()
+{
+  initlog();
+
+  deallog << "Cheching constexpr default constructor of Tensor<rank,dim,Number>"
+          << std::endl;
+  {
+    dealii::LogStream::Prefix p("float");
+    test_consexpr_tensor_constructors<0, 1, float>();
+    test_consexpr_tensor_constructors<0, 2, float>();
+    test_consexpr_tensor_constructors<0, 3, float>();
+    test_consexpr_tensor_constructors<1, 1, float>();
+    test_consexpr_tensor_constructors<1, 2, float>();
+    test_consexpr_tensor_constructors<1, 3, float>();
+    test_consexpr_tensor_constructors<2, 1, float>();
+    test_consexpr_tensor_constructors<2, 2, float>();
+    test_consexpr_tensor_constructors<2, 3, float>();
+  }
+
+  {
+    dealii::LogStream::Prefix p("double");
+    test_consexpr_tensor_constructors<0, 1, double>();
+    test_consexpr_tensor_constructors<0, 2, double>();
+    test_consexpr_tensor_constructors<0, 3, double>();
+    test_consexpr_tensor_constructors<1, 1, double>();
+    test_consexpr_tensor_constructors<1, 2, double>();
+    test_consexpr_tensor_constructors<1, 3, double>();
+    test_consexpr_tensor_constructors<2, 1, double>();
+    test_consexpr_tensor_constructors<2, 2, double>();
+    test_consexpr_tensor_constructors<2, 3, double>();
+  }
+
+  deallog << "Using Tensor within constexpr functions" << std::endl;
+  constexpr double             number = 7.6;
+  DEAL_II_CONSTEXPR const auto tensor_rank0_result =
+    tensor_rank0_constexpr(number);
+  deallog << tensor_rank0_result << std::endl;
+  DEAL_II_CONSTEXPR const auto tensor_rank2_result = tensor_rank2_constexpr();
+  deallog << tensor_rank2_result << std::endl;
+
+  deallog << "OK" << std::endl;
+  return 0;
+}
diff --git a/tests/tensors/constexpr_tensor.output b/tests/tensors/constexpr_tensor.output
new file mode 100644
index 0000000000..7c10bec9a8
--- /dev/null
+++ b/tests/tensors/constexpr_tensor.output
@@ -0,0 +1,78 @@
+
+DEAL::Cheching constexpr default constructor of Tensor<rank,dim,Number>
+DEAL:float:: Tensor<0,1>
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float:: Tensor<0,2>
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float:: Tensor<0,3>
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float:: Tensor<1,1>
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float:: Tensor<1,2>
+DEAL:float::0.00000 0.00000
+DEAL:float::0.00000 0.00000
+DEAL:float::0.00000 0.00000
+DEAL:float:: Tensor<1,3>
+DEAL:float::0.00000 0.00000 0.00000
+DEAL:float::0.00000 0.00000 0.00000
+DEAL:float::0.00000 0.00000 0.00000
+DEAL:float:: Tensor<2,1>
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float::0.00000
+DEAL:float:: Tensor<2,2>
+DEAL:float::0.00000 0.00000 0.00000 0.00000
+DEAL:float::0.00000 0.00000 0.00000 0.00000
+DEAL:float::0.00000 0.00000 0.00000 0.00000
+DEAL:float:: Tensor<2,3>
+DEAL:float::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:float::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:float::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:double:: Tensor<0,1>
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double:: Tensor<0,2>
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double:: Tensor<0,3>
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double:: Tensor<1,1>
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double:: Tensor<1,2>
+DEAL:double::0.00000 0.00000
+DEAL:double::0.00000 0.00000
+DEAL:double::0.00000 0.00000
+DEAL:double:: Tensor<1,3>
+DEAL:double::0.00000 0.00000 0.00000
+DEAL:double::0.00000 0.00000 0.00000
+DEAL:double::0.00000 0.00000 0.00000
+DEAL:double:: Tensor<2,1>
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double::0.00000
+DEAL:double:: Tensor<2,2>
+DEAL:double::0.00000 0.00000 0.00000 0.00000
+DEAL:double::0.00000 0.00000 0.00000 0.00000
+DEAL:double::0.00000 0.00000 0.00000 0.00000
+DEAL:double:: Tensor<2,3>
+DEAL:double::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:double::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:double::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL::Using Tensor within constexpr functions
+DEAL::15.2000
+DEAL::98.0000
+DEAL::OK