Add support for Nth order tensor computations

include/numerics/type_n_tensor.h

@@ -37,9 +37,9 @@ namespace libMesh
  * This class will eventually define a rank-N tensor in \p LIBMESH_DIM
  * dimensional space of type T.
  *
- * Right now it defines a shim to allow for rank-independent code to
- * compile (but not give correct results) in the case of vector-valued
- * elements and second derivatives.
+ * For routines that are too hard to implement with a general N-dimensional,
+ * it defines a shim to allow for rank-independent code to
+ * compile (but error on unimplemented).
  *
  * \author Roy Stogner
  * \date 2012
@@ -55,18 +55,44 @@ class TypeNTensor
 
   TypeNTensor () : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N))) {}
 
-  TypeNTensor (const T &) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N))) {}
+  TypeNTensor(const T &) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N)))
+  {
+    libmesh_not_implemented();
+  }
 
-  TypeNTensor (const TypeVector<T> &) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N))) {}
+  TypeNTensor(const TypeVector<T> &) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N)))
+  {
+    libmesh_not_implemented();
+  }
 
-  TypeNTensor (const TypeTensor<T> &) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N))) {}
+  TypeNTensor(const TypeTensor<T> &) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N)))
+  {
+    libmesh_not_implemented();
+  }
 
-  TypeNTensor (const TypeNTensor<N,T> &) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N))) {}
+  TypeNTensor(const TypeNTensor<N, T> & TN) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N)))
+  {
+    _coords = TN._coords;
+  }
 
-  TypeNTensor & operator = (const TypeNTensor<N,T> &) { libmesh_not_implemented(); return *this; }
+  TypeNTensor(const std::vector<T> & vec) : _coords(std::vector<T>(int_pow(LIBMESH_DIM, N)))
+  {
+    libmesh_assert(int_pow(LIBMESH_DIM, N) == vec.size());
+    _coords = vec;
+  }
 
-  operator TypeVector<T> () const { libmesh_not_implemented(); return 0; }
-  operator VectorValue<T> () const { libmesh_not_implemented(); return 0; }
+  TypeNTensor & operator=(const TypeNTensor<N, T> &) { return *this; }
+
+  operator TypeVector<T>() const
+  {
+    libmesh_assert(N == 1);
+    return TypeVector<T>(_coords[0], _coords[1], _coords[2]);
+  }
+  operator VectorValue<T>() const
+  {
+    libmesh_assert(N == 1);
+    return VectorValue<T>(_coords[0], _coords[1], _coords[2]);
+  }
 
   operator TypeTensor<T> () const { libmesh_not_implemented(); return 0; }
   operator TensorValue<T> () const { libmesh_not_implemented(); return 0; }
@@ -79,21 +105,28 @@ class TypeNTensor
   /**
    * \returns A proxy for the \f$ i^{th} \f$ slice of the tensor.
    */
-  const TypeNTensor<N-1,T> slice (const unsigned int /*i*/) const
+  const TypeNTensor<N - 1, T> slice(const unsigned int i) const
   {
-    libmesh_not_implemented();
-    return TypeNTensor<N-1,T>();
+    return const_cast<TypeNTensor<N, T> *>(this)->slice(i);
   }
 
   /**
    * \returns A writable proxy for the \f$ i^{th} \f$ slice of the tensor.
    */
-  TypeNTensor<N-1,T> slice (const unsigned int /*i*/)
+  TypeNTensor<N - 1, T> slice(const unsigned int i)
   {
-    libmesh_not_implemented();
-    return TypeNTensor<N-1,T>();
+    libmesh_assert(i <= N);
+    TypeNTensor<N - 1, T> slice;
+    unsigned int slice_size = int_pow(LIBMESH_DIM, N - 1);
+    slice._coords.resize(slice_size);
+    for (unsigned int j = 0; j < slice_size; j++)
+      slice._coords[j] = _coords[i * slice_size + j];
+    return slice;
   }
 
+  /**
+   * Zero out a tensor by setting it with = 0
+   */
   template <typename Scalar>
   typename boostcopy::enable_if_c<
     ScalarTraits<Scalar>::value,
@@ -104,42 +137,54 @@ class TypeNTensor
   /**
    * Add two tensors.
    */
-  template<typename T2>
-  TypeNTensor<N,typename CompareTypes<T, T2>::supertype>
-  operator + (const TypeNTensor<N,T2> &) const
+  template <typename T2>
+  TypeNTensor<N, typename CompareTypes<T, T2>::supertype>
+  operator+(const TypeNTensor<N, T2> & other) const
   {
-    libmesh_not_implemented();
-    return TypeNTensor<N,typename CompareTypes<T,T2>::supertype>();
+    TypeNTensor<N, typename CompareTypes<T, T2>::supertype> sum;
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    sum._coords.resize(size);
+    for (unsigned int i = 0; i < size; i++)
+      sum._coords[i] = _coords[i] + other._coords[i];
+    return sum;
   }
 
   /**
    * Add to this tensor.
    */
-  template<typename T2>
-  const TypeNTensor<N,T> & operator += (const TypeNTensor<N,T2> &/*rhs*/)
+  template <typename T2>
+  const TypeNTensor<N, T> & operator+=(const TypeNTensor<N, T2> & other)
   {
-    libmesh_not_implemented();
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    for (unsigned int i = 0; i < size; i++)
+      _coords[i] += other._coords[i];
     return *this;
   }
 
   /**
    * Subtract two tensors.
    */
-  template<typename T2>
-  TypeNTensor<N,typename CompareTypes<T, T2>::supertype>
-  operator - (const TypeNTensor<N,T2> &) const
+  template <typename T2>
+  TypeNTensor<N, typename CompareTypes<T, T2>::supertype>
+  operator-(const TypeNTensor<N, T2> & other) const
   {
-    libmesh_not_implemented();
-    return TypeNTensor<N,typename CompareTypes<T,T2>::supertype>();
+    TypeNTensor<N, typename CompareTypes<T, T2>::supertype> subtract;
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    subtract._coords.resize(size);
+    for (unsigned int i = 0; i < size; i++)
+      subtract._coords[i] = _coords[i] - other._coords[i];
+    return subtract;
   }
 
   /**
    * Subtract from this tensor.
    */
-  template<typename T2>
-  const TypeNTensor<N,T> & operator -= (const TypeNTensor<N,T2> &)
+  template <typename T2>
+  const TypeNTensor<N, T> & operator-=(const TypeNTensor<N, T2> & other)
   {
-    libmesh_not_implemented();
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    for (unsigned int i = 0; i < size; i++)
+      _coords[i] -= other._coords[i];
     return *this;
   }
 
@@ -148,52 +193,64 @@ class TypeNTensor
    */
   TypeNTensor<N,T> operator - () const
   {
-    libmesh_not_implemented();
-    return *this;
+    TypeNTensor<N, T> minus;
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    minus._coords.resize(size);
+    for (unsigned int i = 0; i < size; i++)
+      minus._coords[i] = -_coords[i];
+    return minus;
   }
 
   /**
    * Multiply every entry of a tensor by a number.
    */
   template <typename Scalar>
-  typename boostcopy::enable_if_c<
-    ScalarTraits<Scalar>::value,
-    TypeNTensor<N,typename CompareTypes<T, Scalar>::supertype>>::type
-  operator * (const Scalar) const
+  typename boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
+                                  TypeNTensor<N, typename CompareTypes<T, Scalar>::supertype>>::type
+  operator*(const Scalar factor) const
   {
-    libmesh_not_implemented();
-    return TypeNTensor<N,typename CompareTypes<T, Scalar>::supertype>();
+    TypeNTensor<N, T> multiplied;
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    multiplied._coords.resize(size);
+    for (unsigned int i = 0; i < size; i++)
+      multiplied._coords[i] = factor * _coords[i];
+    return multiplied;
   }
 
   /**
    * Multiply every entry of this tensor by a number.
    */
   template <typename Scalar>
-  const TypeNTensor<N,T> & operator *= (const Scalar)
+  const TypeNTensor<N, T> & operator*=(const Scalar factor)
   {
-    libmesh_not_implemented();
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    for (unsigned int i = 0; i < size; i++)
+      _coords[i] *= factor;
     return *this;
   }
 
   /**
    * Divide every entry of a tensor by a number.
    */
   template <typename Scalar>
-  typename boostcopy::enable_if_c<
-    ScalarTraits<Scalar>::value,
-    TypeNTensor<N,typename CompareTypes<T, Scalar>::supertype>>::type
-  operator / (const Scalar) const
+  typename boostcopy::enable_if_c<ScalarTraits<Scalar>::value,
+                                  TypeNTensor<N, typename CompareTypes<T, Scalar>::supertype>>::type
+  operator/(const Scalar factor) const
   {
-    libmesh_not_implemented();
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    for (unsigned int i = 0; i < size; i++)
+      _coords[i] /= factor;
     return *this;
   }
 
   /**
    * Divide every entry of this tensor by a number.
    */
-  const TypeNTensor<N,T> & operator /= (const T)
+  const TypeNTensor<N, T> & operator/=(const T factor)
   {
-    libmesh_not_implemented();
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    for (unsigned int i = 0; i < size; i++)
+      _coords[i] /= factor;
     return *this;
   }
 
@@ -218,26 +275,30 @@ class TypeNTensor
    * \returns The Frobenius norm of the tensor, i.e. the square-root of
    * the sum of the elements squared.
    */
-  auto norm() const
-  {
-    libmesh_not_implemented();
-    return 0.;
-  }
+  auto norm() const -> decltype(std::norm(T())) { return std::sqrt(norm_sq()); }
 
   /**
    * \returns The Frobenius norm of the tensor squared, i.e. the sum of the
    * entry magnitudes squared.
    */
   auto norm_sq() const
   {
-    libmesh_not_implemented();
-    return 0.;
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    auto norm = 0.;
+    for (unsigned int i = 0; i < size; i++)
+      norm += _coords[i] * _coords[i];
+    return norm;
   }
 
   /**
    * Set all entries of the tensor to 0.
    */
-  void zero() { libmesh_not_implemented(); }
+  void zero()
+  {
+    unsigned int size = int_pow(LIBMESH_DIM, N);
+    for (unsigned int i = 0; i < size; i++)
+      _coords[i] = T(0);
+  }
 
   /**
    * \returns \p true if two tensors are equal, \p false otherwise.
@@ -294,7 +355,10 @@ class TypeNTensor
   void add_scaled (const TypeNTensor<N, T2> &, const T &);
 
   /**
-   * The coordinates of the \p TypeNTensor
+   * The coordinates of the \p TypeNTensor.
+   * The vector is indexed such that slice i is contiguous.
+   * In dimension 2, this would mean that slice i is the ith row,
+   * and _coords has row-major ordering and column minor-ordering
    */
   std::vector<T> _coords;
 

tests/Makefile.am

@@ -108,6 +108,7 @@ unit_tests_sources = \
   numerics/type_vector_test.h \
   numerics/vector_value_test.C \
   numerics/type_tensor_test.C \
+  numerics/type_N_tensor_test.C \
   numerics/sparse_matrix_test.h \
   numerics/dense_matrix_test.C \
   numerics/petsc_matrix_test.C \