/usr/src/RPM/BUILD/trilinos-11.12.1/packages/intrepid/src/Shared/MiniTensor/Intrepid_MiniTensor_Tensor.i.h

// @HEADER
// ************************************************************************
//
//                           Intrepid Package
//                 Copyright (2007) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions: Alejandro Mota (amota@sandia.gov)
//
// ************************************************************************
// @HEADER

#if !defined(Intrepid_MiniTensor_Tensor_i_h)
#define Intrepid_MiniTensor_Tensor_i_h

namespace Intrepid {

//
// Constructor that initializes to NaNs
//
template<typename T, Index N>
inline
Tensor<T, N>::Tensor() :
TensorBase<T, Store>::TensorBase()
{
  return;
}

template<typename T, Index N>
inline
Tensor<T, N>::Tensor(Index const dimension) :
TensorBase<T, Store>::TensorBase(dimension, ORDER)
{
  return;
}

template<typename T, Index N>
inline
Tensor<T, N>::Tensor(ComponentValue const value) :
TensorBase<T, Store>::TensorBase(N, ORDER, value)
{
  return;
}

template<typename T, Index N>
inline
Tensor<T, N>::Tensor(Index const dimension, ComponentValue const value) :
TensorBase<T, Store>::TensorBase(dimension, ORDER, value)
{
  return;
}

//
//  Create tensor from array
//
template<typename T, Index N>
inline
Tensor<T, N>::Tensor(T const * data_ptr) :
TensorBase<T, Store>::TensorBase(N, ORDER, data_ptr)
{
  return;
}

template<typename T, Index N>
inline
Tensor<T, N>::Tensor(Index const dimension, T const * data_ptr) :
TensorBase<T, Store>::TensorBase(dimension, ORDER, data_ptr)
{
  return;
}

//
// Copy constructor
//
template<typename T, Index N>
inline
Tensor<T, N>::Tensor(Tensor<T, N> const & A) :
TensorBase<T, Store>::TensorBase(A)
{
  return;
}

//
// Create tensor specifying components
// \param  s00, s01, ... components in the R^2 canonical basis
//
template<typename T, Index N>
inline
Tensor<T, N>::Tensor(
    T const & s00, T const & s01,
    T const & s10, T const & s11)
{
  Tensor<T, N> &
  self = (*this);

  self.set_dimension(2);

  self[0] = s00;
  self[1] = s01;

  self[2] = s10;
  self[3] = s11;

  return;
}

//
// Create tensor specifying components
// \param  s00, s01, ... components in the R^3 canonical basis
//
template<typename T, Index N>
inline
Tensor<T, N>::Tensor(
    T const & s00, T const & s01, T const & s02,
    T const & s10, T const & s11, T const & s12,
    T const & s20, T const & s21, T const & s22)
{
  Tensor<T, N> &
  self = (*this);

  self.set_dimension(3);

  self[0] = s00;
  self[1] = s01;
  self[2] = s02;

  self[3] = s10;
  self[4] = s11;
  self[5] = s12;

  self[6] = s20;
  self[7] = s21;
  self[8] = s22;

  return;
}

//
//  Create tensor from array with component order
//
template<typename T, Index N>
inline
Tensor<T, N>::Tensor(T const * data_ptr, ComponentOrder const component_order)
{
  assert(data_ptr != NULL);

  fill(data_ptr, component_order);

  return;
}

template<typename T, Index N>
inline
Tensor<T, N>::Tensor(
    Index const dimension,
    T const * data_ptr,
    ComponentOrder const component_order)
{
  assert(data_ptr != NULL);

  Tensor<T, N> &
  self = (*this);

  self.set_dimension(dimension);

  fill(data_ptr, component_order);

  return;
}

//
// 2nd-order tensor from 4th-order tensor
//
template<typename T, Index N>
inline
Tensor<T, N>::Tensor(Tensor4<T, dimension_sqrt<N>::value> const & A)
{
  Index const
  dimension_4th = A.get_dimension();

  Index const
  dimension_2nd = dimension_4th * dimension_4th;

  Tensor<T, N> &
  self = (*this);

  self.set_dimension(dimension_2nd);

  Index const
  number_components = dimension_2nd * dimension_2nd;

  for (Index i = 0; i < number_components; ++i) {
    self[i] = A[i];
  }

  return;
}

//
// Simple destructor
//
template<typename T, Index N>
inline
Tensor<T, N>::~Tensor()
{
  return;
}

//
// Get dimension
//
template<typename T, Index N>
inline
Index
Tensor<T, N>::get_dimension() const
{
  return IS_DYNAMIC == true ? TensorBase<T, Store>::get_dimension() : N;
}

//
// Set dimension
//
template<typename T, Index N>
inline
void
Tensor<T, N>::set_dimension(Index const dimension)
{
  if (IS_DYNAMIC == true) {
    TensorBase<T, Store>::set_dimension(dimension, ORDER);
  }
  else {
    assert(dimension == N);
  }

  return;
}

//
// Indexing for constant tensor
//
template<typename T, Index N>
inline
T const &
Tensor<T, N>::operator()(Index const i, Index const j) const
{
  Tensor<T, N> const &
  self = (*this);

  Index const
  dimension = self.get_dimension();

  return self[i * dimension + j];
}

//
//Tensor indexing
//
template<typename T, Index N>
inline
T &
Tensor<T, N>::operator()(Index const i, Index const j)
{
  Tensor<T, N> &
  self = (*this);

  Index const
  dimension = self.get_dimension();

  return self[i * dimension + j];
}

//
// Fill components with value specification
//
template<typename T, Index N>
inline
void
Tensor<T, N>::fill(ComponentValue const value)
{
  TensorBase<T, Store>::fill(value);
  return;
}

//
// Fill components with value as parameter
//
template<typename T, Index N>
inline
void
Tensor<T, N>::fill(T const & s)
{
  TensorBase<T, Store>::fill(s);
  return;
}

//
// Fill components from array defined by pointer.
//
template<typename T, Index N>
inline
void
Tensor<T, N>::fill(T const * data_ptr)
{
  TensorBase<T, Store>::fill(data_ptr);
  return;
}

//
// Fill components from array defined by pointer.
//
template<typename T, Index N>
inline
void
Tensor<T, N>::fill(T const * data_ptr, ComponentOrder const component_order)
{
  assert(data_ptr != NULL);

  Tensor<T, N> &
  self = (*this);

  Index const
  dimension = self.get_dimension();

  switch (dimension) {

    default:
      TensorBase<T, Store>::fill(data_ptr);
      break;

    case 3:

      switch (component_order) {

        case CANONICAL:
          TensorBase<T, Store>::fill(data_ptr);
          break;

        case SIERRA_FULL:

          //  0  1  2  3  4  5  6  7  8
          // XX YY ZZ XY YZ ZX YX ZY XZ
          //  0  4  8  1  5  6  3  7  2
          self[0] = data_ptr[0];
          self[4] = data_ptr[1];
          self[8] = data_ptr[2];

          self[1] = data_ptr[3];
          self[5] = data_ptr[4];
          self[6] = data_ptr[5];

          self[3] = data_ptr[6];
          self[7] = data_ptr[7];
          self[2] = data_ptr[8];
          break;

        case SIERRA_SYMMETRIC:
          self[0] = data_ptr[0];
          self[4] = data_ptr[1];
          self[8] = data_ptr[2];

          self[1] = data_ptr[3];
          self[5] = data_ptr[4];
          self[6] = data_ptr[5];

          self[3] = data_ptr[3];
          self[7] = data_ptr[4];
          self[2] = data_ptr[5];
          break;

        default:
          std::cerr << "ERROR: " << __PRETTY_FUNCTION__;
          std::cerr << std::endl;
          std::cerr << "Unknown component order.";
          std::cerr << std::endl;
          exit(1);
          break;

      }

      break;
  }

  return;
}

//
// Tensor addition
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
operator+(Tensor<S, N> const & A, Tensor<T, N> const & B)
{
  Tensor<typename Promote<S, T>::type, N>
  C(A.get_dimension());

  add(A, B, C);

  return C;
}

//
// Tensor subtraction
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
operator-(Tensor<S, N> const & A, Tensor<T, N> const & B)
{
  Tensor<typename Promote<S, T>::type, N>
  C(A.get_dimension());

  subtract(A, B, C);

  return C;
}

//
// Tensor minus
//
template<typename T, Index N>
inline
Tensor<T, N>
operator-(Tensor<T, N> const & A)
{
  Tensor<T, N>
  B(A.get_dimension());

  minus(A, B);

  return B;
}

//
// Tensor equality
//
template<typename T, Index N>
inline
bool
operator==(Tensor<T, N> const & A, Tensor<T, N> const & B)
{
  return equal(A, B);
}

//
// Tensor inequality
//
template<typename T, Index N>
inline
bool
operator!=(Tensor<T, N> const & A, Tensor<T, N> const & B)
{
  return not_equal(A, B);
}

//
// Scalar tensor product
//
template<typename S, typename T, Index N>
inline
typename lazy_disable_if< order_1234<S>, apply_tensor< Promote<S,T>, N> >::type
operator*(S const & s, Tensor<T, N> const & A)
{
  Tensor<typename Promote<S, T>::type, N>
  B(A.get_dimension());

  scale(A, s, B);

  return B;
}

//
// Tensor scalar product
//
template<typename S, typename T, Index N>
inline
typename lazy_disable_if< order_1234<S>, apply_tensor< Promote<S,T>, N> >::type
operator*(Tensor<T, N> const & A, S const & s)
{
  Tensor<typename Promote<S, T>::type, N>
  B(A.get_dimension());

  scale(A, s, B);

  return B;
}

//
// Tensor scalar division
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
operator/(Tensor<T, N> const & A, S const & s)
{
  Tensor<typename Promote<S, T>::type, N>
  B(A.get_dimension());

  divide(A, s, B);

  return B;
}

//
// Scalar tensor division
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
operator/(S const & s, Tensor<T, N> const & A)
{
  Tensor<typename Promote<S, T>::type, N>
  B(A.get_dimension());

  split(A, s, B);

  return B;
}

//
// Tensor vector product v = A u
//
template<typename S, typename T, Index N>
inline
Vector<typename Promote<S, T>::type, N>
operator*(Tensor<T, N> const & A, Vector<S, N> const & u)
{
  return dot(A, u);
}

//
// Vector tensor product v = u A
//
template<typename S, typename T, Index N>
inline
Vector<typename Promote<S, T>::type, N>
operator*(Vector<S, N> const & u, Tensor<T, N> const & A)
{
  return dot(u, A);
}

//
// Tensor dot product C = A B
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
operator*(Tensor<S, N> const & A, Tensor<T, N> const & B)
{
  return dot(A, B);
}

namespace {

template<typename S>
bool
greater_than(S const & a, S const & b)
{
  return a.first > b.first;
}

} // anonymous namespace

//
// Sort and index in descending order. Useful for ordering singular values
// and eigenvalues and corresponding vectors in the respective decompositions.
//
template<typename T, Index N>
std::pair<Vector<T, N>, Tensor<T, N> >
sort_permutation(Vector<T, N> const & u)
{

  Index const
  dimension = u.get_dimension();

  std::vector<std::pair<T, Index > >
  s(dimension);

  for (Index i = 0; i < dimension; ++i) {
    s[i].first = u(i);
    s[i].second = i;
  }

  std::sort(s.begin(), s.end(), greater_than< std::pair<T, Index > > );

  Vector<T, N> v(dimension);

  Tensor<T, N>
  P = zero<T, N>(dimension);

  for (Index i = 0; i < dimension; ++i) {
    v(i) = s[i].first;
    P(s[i].second, i) = 1.0;
  }

  return std::make_pair(v, P);

}

//
// Extract a row as a vector
//
template<typename T, Index N>
Vector<T, N>
row(Tensor<T, N> const & A, Index const i)
{
  Index const
  dimension = A.get_dimension();

  Vector<T, N>
  v(dimension);

  switch (dimension) {
    default:
      for (Index j = 0; j < dimension; ++j) {
        v(j) = A(i,j);
      }
      break;

    case 2:
      v(0) = A(i,0);
      v(1) = A(i,1);
      break;

    case 3:
      v(0) = A(i,0);
      v(1) = A(i,1);
      v(2) = A(i,2);
      break;
  }

  return v;
}

//
// Extract a column as a vector
//
template<typename T, Index N>
Vector<T, N>
col(Tensor<T, N> const & A, Index const j)
{
  Index const
  dimension = A.get_dimension();

  Vector<T, N>
  v(dimension);

  switch (dimension) {
    default:
      for (Index i = 0; i < dimension; ++i) {
        v(i) = A(i,j);
      }
      break;

    case 2:
      v(0) = A(0,j);
      v(1) = A(1,j);
      break;

    case 3:
      v(0) = A(0,j);
      v(1) = A(1,j);
      v(2) = A(2,j);
      break;
  }

  return v;
}

//
// R^N tensor vector product v = A u
// \param A tensor
// \param u vector
// \return \f$ A u \f$
//
template<typename S, typename T, Index N>
inline
Vector<typename Promote<S, T>::type, N>
dot(Tensor<T, N> const & A, Vector<S, N> const & u)
{
  Index const
  dimension = A.get_dimension();

  assert(u.get_dimension() == dimension);

  Vector<typename Promote<S, T>::type, N>
  v(dimension);

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {

        typename Promote<S, T>::type
        s = 0.0;

        for (Index p = 0; p < dimension; ++p) {
          s += A(i, p) * u(p);
        }
        v(i) = s;
      }
      break;

    case 3:
      v(0) = A(0,0)*u(0) + A(0,1)*u(1) + A(0,2)*u(2);
      v(1) = A(1,0)*u(0) + A(1,1)*u(1) + A(1,2)*u(2);
      v(2) = A(2,0)*u(0) + A(2,1)*u(1) + A(2,2)*u(2);
      break;

    case 2:
      v(0) = A(0,0)*u(0) + A(0,1)*u(1);
      v(1) = A(1,0)*u(0) + A(1,1)*u(1);
      break;

  }

  return v;
}

//
// R^N vector tensor product v = u A
// \param A tensor
// \param u vector
// \return \f$ u A = A^T u \f$
//
template<typename S, typename T, Index N>
inline
Vector<typename Promote<S, T>::type, N>
dot(Vector<S, N> const & u, Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();

  assert(u.get_dimension() == dimension);

  Vector<typename Promote<S, T>::type, N>
  v(dimension);

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {

        typename Promote<S, T>::type
        s = 0.0;

        for (Index p = 0; p < dimension; ++p) {
          s += A(p, i) * u(p);
        }
        v(i) = s;
      }
      break;

    case 3:
      v(0) = A(0,0)*u(0) + A(1,0)*u(1) + A(2,0)*u(2);
      v(1) = A(0,1)*u(0) + A(1,1)*u(1) + A(2,1)*u(2);
      v(2) = A(0,2)*u(0) + A(1,2)*u(1) + A(2,2)*u(2);
      break;

    case 2:
      v(0) = A(0,0)*u(0) + A(1,0)*u(1);
      v(1) = A(0,1)*u(0) + A(1,1)*u(1);
      break;

  }

  return v;
}

//
// R^N tensor tensor product C = A B
// \param A tensor
// \param B tensor
// \return a tensor \f$ A \cdot B \f$
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
dot(Tensor<S, N> const & A, Tensor<T, N> const & B)
{
  Index const
  dimension = A.get_dimension();

  assert(B.get_dimension() == dimension);

  Tensor<typename Promote<S, T>::type, N>
  C(dimension);

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {
        for (Index j = 0; j < dimension; ++j) {

          typename Promote<S, T>::type
          s = 0.0;

          for (Index p = 0; p < dimension; ++p) {
            s += A(i, p) * B(p, j);
          }
          C(i, j) = s;
        }
      }
      break;

    case 3:
      C(0,0) = A(0,0)*B(0,0) + A(0,1)*B(1,0) + A(0,2)*B(2,0);
      C(0,1) = A(0,0)*B(0,1) + A(0,1)*B(1,1) + A(0,2)*B(2,1);
      C(0,2) = A(0,0)*B(0,2) + A(0,1)*B(1,2) + A(0,2)*B(2,2);

      C(1,0) = A(1,0)*B(0,0) + A(1,1)*B(1,0) + A(1,2)*B(2,0);
      C(1,1) = A(1,0)*B(0,1) + A(1,1)*B(1,1) + A(1,2)*B(2,1);
      C(1,2) = A(1,0)*B(0,2) + A(1,1)*B(1,2) + A(1,2)*B(2,2);

      C(2,0) = A(2,0)*B(0,0) + A(2,1)*B(1,0) + A(2,2)*B(2,0);
      C(2,1) = A(2,0)*B(0,1) + A(2,1)*B(1,1) + A(2,2)*B(2,1);
      C(2,2) = A(2,0)*B(0,2) + A(2,1)*B(1,2) + A(2,2)*B(2,2);
      break;

    case 2:
      C(0,0) = A(0,0)*B(0,0) + A(0,1)*B(1,0);
      C(0,1) = A(0,0)*B(0,1) + A(0,1)*B(1,1);

      C(1,0) = A(1,0)*B(0,0) + A(1,1)*B(1,0);
      C(1,1) = A(1,0)*B(0,1) + A(1,1)*B(1,1);
      break;

  }

  return C;
}

//
// R^N tensor tensor product C = A^T B
// \param A tensor
// \param B tensor
// \return a tensor \f$ A^T \cdot B \f$
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
t_dot(Tensor<S, N> const & A, Tensor<T, N> const & B)
{
  Index const
  dimension = A.get_dimension();

  assert(B.get_dimension() == dimension);

  Tensor<typename Promote<S, T>::type, N>
  C(dimension);

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {
        for (Index j = 0; j < dimension; ++j) {

          typename Promote<S, T>::type
          s = 0.0;

          for (Index p = 0; p < dimension; ++p) {
            s += A(p, i) * B(p, j);
          }
          C(i, j) = s;
        }
      }
      break;

    case 3:
      C(0,0) = A(0,0)*B(0,0) + A(1,0)*B(1,0) + A(2,0)*B(2,0);
      C(0,1) = A(0,0)*B(0,1) + A(1,0)*B(1,1) + A(2,0)*B(2,1);
      C(0,2) = A(0,0)*B(0,2) + A(1,0)*B(1,2) + A(2,0)*B(2,2);

      C(1,0) = A(0,1)*B(0,0) + A(1,1)*B(1,0) + A(2,1)*B(2,0);
      C(1,1) = A(0,1)*B(0,1) + A(1,1)*B(1,1) + A(2,1)*B(2,1);
      C(1,2) = A(0,1)*B(0,2) + A(1,1)*B(1,2) + A(2,1)*B(2,2);

      C(2,0) = A(0,2)*B(0,0) + A(1,2)*B(1,0) + A(2,2)*B(2,0);
      C(2,1) = A(0,2)*B(0,1) + A(1,2)*B(1,1) + A(2,2)*B(2,1);
      C(2,2) = A(0,2)*B(0,2) + A(1,2)*B(1,2) + A(2,2)*B(2,2);
      break;

    case 2:
      C(0,0) = A(0,0)*B(0,0) + A(1,0)*B(1,0);
      C(0,1) = A(0,0)*B(0,1) + A(1,0)*B(1,1);

      C(1,0) = A(0,1)*B(0,0) + A(1,1)*B(1,0);
      C(1,1) = A(0,1)*B(0,1) + A(1,1)*B(1,1);
      break;

  }

  return C;
}

//
// R^N tensor tensor product C = A B^T
// \param A tensor
// \param B tensor
// \return a tensor \f$ A \cdot B^T \f$
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
dot_t(Tensor<S, N> const & A, Tensor<T, N> const & B)
{
  Index const
  dimension = A.get_dimension();

  assert(B.get_dimension() == dimension);

  Tensor<typename Promote<S, T>::type, N>
  C(dimension);

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {
        for (Index j = 0; j < dimension; ++j) {

          typename Promote<S, T>::type
          s = 0.0;

          for (Index p = 0; p < dimension; ++p) {
            s += A(i, p) * B(j, p);
          }
          C(i, j) = s;
        }
      }
      break;

    case 3:
      C(0,0) = A(0,0)*B(0,0) + A(0,1)*B(0,1) + A(0,2)*B(0,2);
      C(0,1) = A(0,0)*B(1,0) + A(0,1)*B(1,1) + A(0,2)*B(1,2);
      C(0,2) = A(0,0)*B(2,0) + A(0,1)*B(2,1) + A(0,2)*B(2,2);

      C(1,0) = A(1,0)*B(0,0) + A(1,1)*B(0,1) + A(1,2)*B(0,2);
      C(1,1) = A(1,0)*B(1,0) + A(1,1)*B(1,1) + A(1,2)*B(1,2);
      C(1,2) = A(1,0)*B(2,0) + A(1,1)*B(2,1) + A(1,2)*B(2,2);

      C(2,0) = A(2,0)*B(0,0) + A(2,1)*B(0,1) + A(2,2)*B(0,2);
      C(2,1) = A(2,0)*B(1,0) + A(2,1)*B(1,1) + A(2,2)*B(1,2);
      C(2,2) = A(2,0)*B(2,0) + A(2,1)*B(2,1) + A(2,2)*B(2,2);
      break;

    case 2:
      C(0,0) = A(0,0)*B(0,0) + A(0,1)*B(0,1);
      C(0,1) = A(0,0)*B(1,0) + A(0,1)*B(1,1);

      C(1,0) = A(1,0)*B(0,0) + A(1,1)*B(0,1);
      C(1,1) = A(1,0)*B(1,0) + A(1,1)*B(1,1);
      break;

  }

  return C;
}

//
// R^N tensor tensor product C = A^T B^T
// \param A tensor
// \param B tensor
// \return a tensor \f$ A^T \cdot B^T \f$
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
t_dot_t(Tensor<S, N> const & A, Tensor<T, N> const & B)
{
  Index const
  dimension = A.get_dimension();

  assert(B.get_dimension() == dimension);

  Tensor<typename Promote<S, T>::type, N>
  C(dimension);

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {
        for (Index j = 0; j < dimension; ++j) {

          typename Promote<S, T>::type
          s = 0.0;

          for (Index p = 0; p < dimension; ++p) {
            s += A(p, i) * B(j, p);
          }
          C(i, j) = s;
        }
      }
      break;

    case 3:
      C(0,0) = A(0,0)*B(0,0) + A(1,0)*B(0,1) + A(2,0)*B(0,2);
      C(0,1) = A(0,0)*B(1,0) + A(1,0)*B(1,1) + A(2,0)*B(1,2);
      C(0,2) = A(0,0)*B(2,0) + A(1,0)*B(2,1) + A(2,0)*B(2,2);

      C(1,0) = A(0,1)*B(0,0) + A(1,1)*B(0,1) + A(2,1)*B(0,2);
      C(1,1) = A(0,1)*B(1,0) + A(1,1)*B(1,1) + A(2,1)*B(1,2);
      C(1,2) = A(0,1)*B(2,0) + A(1,1)*B(2,1) + A(2,1)*B(2,2);

      C(2,0) = A(0,2)*B(0,0) + A(1,2)*B(0,1) + A(2,2)*B(0,2);
      C(2,1) = A(0,2)*B(1,0) + A(1,2)*B(1,1) + A(2,2)*B(1,2);
      C(2,2) = A(0,2)*B(2,0) + A(1,2)*B(2,1) + A(2,2)*B(2,2);
      break;

    case 2:
      C(0,0) = A(0,0)*B(0,0) + A(1,0)*B(0,1);
      C(0,1) = A(0,0)*B(1,0) + A(1,0)*B(1,1);

      C(1,0) = A(0,1)*B(0,0) + A(1,1)*B(0,1);
      C(1,1) = A(0,1)*B(1,0) + A(1,1)*B(1,1);
      break;

  }

  return C;
}

//
// R^N tensor tensor double dot product (contraction)
// \param A tensor
// \param B tensor
// \return a scalar \f$ A : B \f$
//
template<typename S, typename T, Index N>
inline
typename Promote<S, T>::type
dotdot(Tensor<S, N> const & A, Tensor<T, N> const & B)
{
  Index const
  dimension = A.get_dimension();

  assert(B.get_dimension() == dimension);

  typename Promote<S, T>::type
  s = 0.0;

  switch (dimension) {

    default:
      for (Index p = 0; p < dimension; ++p) {
        for (Index q = 0; q < dimension; ++q) {
          s += A(p, q) * B(p, q);
        }
      }
      break;

    case 3:
      s+= A(0,0)*B(0,0) + A(0,1)*B(0,1) + A(0,2)*B(0,2);
      s+= A(1,0)*B(1,0) + A(1,1)*B(1,1) + A(1,2)*B(1,2);
      s+= A(2,0)*B(2,0) + A(2,1)*B(2,1) + A(2,2)*B(2,2);
      break;

    case 2:
      s+= A(0,0)*B(0,0) + A(0,1)*B(0,1);
      s+= A(1,0)*B(1,0) + A(1,1)*B(1,1);
      break;

  }

  return s;
}

//
// R^N dyad
// \param u vector
// \param v vector
// \return \f$ u \otimes v \f$
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
dyad(Vector<S, N> const & u, Vector<T, N> const & v)
{
  Index const
  dimension = u.get_dimension();

  assert(v.get_dimension() == dimension);

  Tensor<typename Promote<S, T>::type, N>
  A(dimension);

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {

        typename Promote<S, T>::type const
        s = u(i);

        for (Index j = 0; j < dimension; ++j) {
          A(i, j) = s * v(j);
        }
      }
      break;

    case 3:
      A(0,0) = u(0) * v(0);
      A(0,1) = u(0) * v(1);
      A(0,2) = u(0) * v(2);

      A(1,0) = u(1) * v(0);
      A(1,1) = u(1) * v(1);
      A(1,2) = u(1) * v(2);

      A(2,0) = u(2) * v(0);
      A(2,1) = u(2) * v(1);
      A(2,2) = u(2) * v(2);
      break;

    case 2:
      A(0,0) = u(0) * v(0);
      A(0,1) = u(0) * v(1);

      A(1,0) = u(1) * v(0);
      A(1,1) = u(1) * v(1);
      break;

  }

  return A;
}

//
// R^N bun operator, just for Jay, and now Reese too.
// \param u vector
// \param v vector
// \return \f$ u \otimes v \f$
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
bun(Vector<S, N> const & u, Vector<T, N> const & v)
{
  return dyad(u, v);
}

//
// R^N tensor product
// \param u vector
// \param v vector
// \return \f$ u \otimes v \f$
//
template<typename S, typename T, Index N>
inline
Tensor<typename Promote<S, T>::type, N>
tensor(Vector<S, N> const & u, Vector<T, N> const & v)
{
  return dyad(u, v);
}

//
// R^N diagonal tensor from vector
// \param v vector
// \return A = diag(v)
//
template<typename T, Index N>
Tensor<T, N>
diag(Vector<T, N> const & v)
{
  Index const
  dimension = v.get_dimension();

  Tensor<T, N>
  A = zero<T, N>(dimension);

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {
        A(i, i) = v(i);
      }
      break;

    case 3:
      A(0,0) = v(0);
      A(1,1) = v(1);
      A(2,2) = v(2);
      break;

    case 2:
      A(0,0) = v(0);
      A(1,1) = v(1);
      break;

  }

  return A;
}

//
// R^N diagonal of tensor in a vector
// \param A tensor
// \return v = diag(A)
//
template<typename T, Index N>
Vector<T, N>
diag(Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();

  Vector<T, N>
  v(dimension);

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {
        v(i) = A(i, i);
      }
      break;

    case 3:
      v(0) = A(0,0);
      v(1) = A(1,1);
      v(2) = A(2,2);
      break;

    case 2:
      v(0) = A(0,0);
      v(1) = A(1,1);
      break;

  }

  return v;
}

//
// Zero 2nd-order tensor
// All components are zero
//
template<typename T, Index N>
inline
Tensor<T, N> const
zero()
{
  return Tensor<T, N>(N, ZEROS);
}

template<typename T>
inline
Tensor<T, DYNAMIC> const
zero(Index const dimension)
{
  return Tensor<T, DYNAMIC>(dimension, ZEROS);
}

template<typename T, Index N>
inline
Tensor<T, N> const
zero(Index const dimension)
{
  if (N != DYNAMIC) assert(dimension == N);
  return Tensor<T, N>(dimension, ZEROS);
}

//
// R^N 2nd-order identity tensor
//
namespace {

template< typename T, Index N>
inline
void ones_in_diagonal(Tensor<T, N> & A)
{
  Index const
  dimension = A.get_dimension();

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {
        A(i, i) = 1.0;
      }
      break;

    case 3:
      A(0,0) = 1.0;
      A(1,1) = 1.0;
      A(2,2) = 1.0;
      break;

    case 2:
      A(0,0) = 1.0;
      A(1,1) = 1.0;
      break;

  }

  return;
}

} // anonymous namespace

template<typename T, Index N>
inline
Tensor<T, N> const
identity()
{
  Tensor<T, N> A(N, ZEROS);
  ones_in_diagonal(A);
  return A;
}

template<typename T>
inline
Tensor<T, DYNAMIC> const
identity(Index const dimension)
{
  Tensor<T, DYNAMIC> A(dimension, ZEROS);
  ones_in_diagonal(A);
  return A;
}

template<typename T, Index N>
inline
Tensor<T, N> const
identity(Index const dimension)
{
  if (N != DYNAMIC) assert(dimension == N);

  Tensor<T, N> A(dimension, ZEROS);
  ones_in_diagonal(A);
  return A;
}

//
// R^N 2nd-order identity tensor, à la Matlab
//
template<typename T, Index N>
inline
Tensor<T, N> const
eye()
{
  return identity<T, N>();
}

template<typename T>
inline
Tensor<T, DYNAMIC> const
eye(Index const dimension)
{
  return identity<T>(dimension);
}

template<typename T, Index N>
inline
Tensor<T, N> const
eye(Index const dimension)
{
  return identity<T, N>(dimension);
}

//
// R^N 2nd-order tensor transpose
//
template<typename T, Index N>
inline
Tensor<T, N>
transpose(Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();

  Tensor<T, N>
  B = A;

  switch (dimension) {

    default:
      for (Index i = 0; i < dimension; ++i) {
        for (Index j = i + 1; j < dimension; ++j) {
          std::swap(B(i, j), B(j, i));
        }
      }
      break;

    case 3:
      std::swap(B(0,1), B(1,0));
      std::swap(B(0,2), B(2,0));

      std::swap(B(1,2), B(2,1));

      break;

    case 2:
      std::swap(B(0,1), B(1,0));

      break;

  }

  return B;
}

//
// R^N symmetric part of 2nd-order tensor
// \return \f$ \frac{1}{2}(A + A^T) \f$
//
template<typename T, Index N>
inline
Tensor<T, N>
sym(Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();

  Tensor<T, N>
  B(dimension);

  switch (dimension) {

    default:
      B = 0.5 * (A + transpose(A));
      break;

    case 3:
    {
      T const & s00 = A(0,0);
      T const & s11 = A(1,1);
      T const & s22 = A(2,2);

      T const s01 = 0.5 * (A(0,1) + A(1,0));
      T const s02 = 0.5 * (A(0,2) + A(2,0));
      T const s12 = 0.5 * (A(1,2) + A(2,1));

      B(0,0) = s00;
      B(0,1) = s01;
      B(0,2) = s02;

      B(1,0) = s01;
      B(1,1) = s11;
      B(1,2) = s12;

      B(2,0) = s02;
      B(2,1) = s12;
      B(2,2) = s22;
    }
    break;

    case 2:
    {
      T const & s00 = A(0,0);
      T const & s11 = A(1,1);

      T const s01 = 0.5 * (A(0,1) + A(1,0));

      B(0,0) = s00;
      B(0,1) = s01;

      B(1,0) = s01;
      B(1,1) = s11;
    }
    break;

  }

  return B;
}

//
// R^N skew symmetric part of 2nd-order tensor
// \return \f$ \frac{1}{2}(A - A^T) \f$
//
template<typename T, Index N>
inline
Tensor<T, N>
skew(Tensor<T, N> const & A)
{
  Index const
  dimension = A.get_dimension();

  Tensor<T, N>
  B(dimension);

  switch (dimension) {

    default:
      B = 0.5 * (A - transpose(A));
      break;

    case 3:
    {
      T const s01 = 0.5*(A(0,1)-A(1,0));
      T const s02 = 0.5*(A(0,2)-A(2,0));
      T const s12 = 0.5*(A(1,2)-A(2,1));

      B(0,0) = 0.0;
      B(0,1) = s01;
      B(0,2) = s02;

      B(1,0) = -s01;
      B(1,1) = 0.0;
      B(1,2) = s12;

      B(2,0) = -s02;
      B(2,1) = -s12;
      B(2,2) = 0.0;
    }
    break;

    case 2:
    {
      T const s01 = 0.5*(A(0,1)-A(1,0));

      B(0,0) = 0.0;
      B(0,1) = s01;

      B(1,0) = -s01;
      B(1,1) = 0.0;
    }
    break;

  }

  return B;
}

//
// R^N skew symmetric 2nd-order tensor from vector, undefined
// for N!=3.
// \param u vector
//
template<typename T, Index N>
inline
Tensor<T, N>
skew(Vector<T, N> const & u)
{
  Index const
  dimension = u.get_dimension();

  Tensor<T, N>
  A(dimension);

  switch (dimension) {

    case 3:
      A(0,0) = 0.0;
      A(0,1) = -u(2);
      A(0,2) = u(1);

      A(1,0) = u(2);
      A(1,1) = 0.0;
      A(1,2) = -u(0);

      A(2,0) = -u(1);
      A(2,1) = u(0);
      A(2,2) = 0.0;
      break;

    default:
      std::cerr << "ERROR: " << __PRETTY_FUNCTION__;
      std::cerr << std::endl;
      std::cerr << "Skew from vector undefined for R^" << N;
      std::cerr << std::endl;
      exit(1);
      break;

  }

  return A;
}

} // namespace Intrepid

#endif // Intrepid_MiniTensor_Tensor_i_h