derivative_form.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2013 - 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.
//
// ---------------------------------------------------------------------

#ifndef dealii_derivative_form_h
#define dealii_derivative_form_h

#include <deal.II/base/tensor.h>

DEAL_II_NAMESPACE_OPEN

template <int order, int dim, int spacedim, typename Number = double>
class DerivativeForm
{
public:
  DerivativeForm();

  DerivativeForm(const Tensor<order + 1, dim, Number> &);

  Tensor<order, dim, Number> &operator[](const unsigned int i);

  const Tensor<order, dim, Number> &operator[](const unsigned int i) const;

  DerivativeForm &
  operator=(const Tensor<order + 1, dim, Number> &);

  DerivativeForm &
  operator=(const Tensor<1, dim, Number> &);

  operator Tensor<order + 1, dim, Number>() const;

  operator Tensor<1, dim, Number>() const;

  DerivativeForm<1, spacedim, dim, Number>
  transpose() const;

  typename numbers::NumberTraits<Number>::real_type
  norm() const;

  Number
  determinant() const;

  DerivativeForm<1, dim, spacedim, Number>
  covariant_form() const;

  static std::size_t
  memory_consumption();

  DeclException1(ExcInvalidTensorIndex,
                 int,
                 << "Invalid DerivativeForm index " << arg1);

private:
  DerivativeForm<1, dim, spacedim, Number>
  times_T_t(const Tensor<2, dim, Number> &T) const;


  Tensor<order, dim, Number> tensor[spacedim];
};


/*--------------------------- Inline functions -----------------------------*/

#ifndef DOXYGEN



template <int order, int dim, int spacedim, typename Number>
inline DerivativeForm<order, dim, spacedim, Number>::DerivativeForm()
{
  // default constructor. not specifying an initializer list calls
  // the default constructor of the subobjects, which initialize them
  // selves. therefore, the tensor array  is set to zero this way
}



template <int order, int dim, int spacedim, typename Number>
inline DerivativeForm<order, dim, spacedim, Number>::DerivativeForm(
  const Tensor<order + 1, dim, Number> &T)
{
  Assert((dim == spacedim),
         ExcMessage("Only allowed for forms with dim==spacedim."));
  if (dim == spacedim)
    for (unsigned int j = 0; j < dim; ++j)
      (*this)[j] = T[j];
}



template <int order, int dim, int spacedim, typename Number>
inline DerivativeForm<order, dim, spacedim, Number> &
DerivativeForm<order, dim, spacedim, Number>::
operator=(const Tensor<order + 1, dim, Number> &ta)
{
  Assert((dim == spacedim), ExcMessage("Only allowed when dim==spacedim."));

  if (dim == spacedim)
    for (unsigned int j = 0; j < dim; ++j)
      (*this)[j] = ta[j];
  return *this;
}



template <int order, int dim, int spacedim, typename Number>
inline DerivativeForm<order, dim, spacedim, Number> &
DerivativeForm<order, dim, spacedim, Number>::
operator=(const Tensor<1, dim, Number> &T)
{
  Assert((1 == spacedim) && (order == 1),
         ExcMessage("Only allowed for spacedim==1 and order==1."));

  (*this)[0] = T;

  return *this;
}



template <int order, int dim, int spacedim, typename Number>
inline Tensor<order, dim, Number> &
  DerivativeForm<order, dim, spacedim, Number>::operator[](const unsigned int i)
{
  Assert(i < spacedim, ExcIndexRange(i, 0, spacedim));

  return tensor[i];
}



template <int order, int dim, int spacedim, typename Number>
inline const Tensor<order, dim, Number> &
  DerivativeForm<order, dim, spacedim, Number>::
  operator[](const unsigned int i) const
{
  Assert(i < spacedim, ExcIndexRange(i, 0, spacedim));

  return tensor[i];
}



template <int order, int dim, int spacedim, typename Number>
inline DerivativeForm<order, dim, spacedim, Number>::
operator Tensor<1, dim, Number>() const
{
  Assert((1 == spacedim) && (order == 1),
         ExcMessage("Only allowed for spacedim==1."));

  return (*this)[0];
}



template <int order, int dim, int spacedim, typename Number>
inline DerivativeForm<order, dim, spacedim, Number>::
operator Tensor<order + 1, dim, Number>() const
{
  Assert((dim == spacedim), ExcMessage("Only allowed when dim==spacedim."));

  Tensor<order + 1, dim, Number> t;

  if (dim == spacedim)
    for (unsigned int j = 0; j < dim; ++j)
      t[j] = (*this)[j];

  return t;
}



template <int order, int dim, int spacedim, typename Number>
inline DerivativeForm<1, spacedim, dim, Number>
DerivativeForm<order, dim, spacedim, Number>::transpose() const
{
  Assert(order == 1, ExcMessage("Only for rectangular DerivativeForm."));
  DerivativeForm<1, spacedim, dim, Number> tt;

  for (unsigned int i = 0; i < spacedim; ++i)
    for (unsigned int j = 0; j < dim; ++j)
      tt[j][i] = (*this)[i][j];

  return tt;
}



template <int order, int dim, int spacedim, typename Number>
inline DerivativeForm<1, dim, spacedim, Number>
DerivativeForm<order, dim, spacedim, Number>::times_T_t(
  const Tensor<2, dim, Number> &T) const
{
  Assert(order == 1, ExcMessage("Only for order == 1."));
  DerivativeForm<1, dim, spacedim, Number> dest;
  for (unsigned int i = 0; i < spacedim; ++i)
    for (unsigned int j = 0; j < dim; ++j)
      dest[i][j] = (*this)[i] * T[j];

  return dest;
}



template <int order, int dim, int spacedim, typename Number>
inline typename numbers::NumberTraits<Number>::real_type
DerivativeForm<order, dim, spacedim, Number>::norm() const
{
  typename numbers::NumberTraits<Number>::real_type sum_of_squares = 0;
  for (unsigned int i = 0; i < spacedim; ++i)
    sum_of_squares += tensor[i].norm_square();
  return std::sqrt(sum_of_squares);
}



template <int order, int dim, int spacedim, typename Number>
inline Number
DerivativeForm<order, dim, spacedim, Number>::determinant() const
{
  Assert(order == 1, ExcMessage("Only for order == 1."));
  if (dim == spacedim)
    {
      const Tensor<2, dim, Number> T =
        static_cast<Tensor<2, dim, Number>>(*this);
      return ::determinant(T);
    }
  else
    {
      Assert(spacedim > dim, ExcMessage("Only for spacedim>dim."));
      const DerivativeForm<1, spacedim, dim, Number> DF_t = this->transpose();
      Tensor<2, dim, Number> G; // First fundamental form
      for (unsigned int i = 0; i < dim; ++i)
        for (unsigned int j = 0; j < dim; ++j)
          G[i][j] = DF_t[i] * DF_t[j];

      return (sqrt(::determinant(G)));
    }
}



template <int order, int dim, int spacedim, typename Number>
inline DerivativeForm<1, dim, spacedim, Number>
DerivativeForm<order, dim, spacedim, Number>::covariant_form() const
{
  if (dim == spacedim)
    {
      const Tensor<2, dim, Number> DF_t =
        ::transpose(invert(static_cast<Tensor<2, dim, Number>>(*this)));
      return DerivativeForm<1, dim, spacedim>(DF_t);
    }
  else
    {
      const DerivativeForm<1, spacedim, dim> DF_t = this->transpose();
      Tensor<2, dim, Number>                 G; // First fundamental form
      for (unsigned int i = 0; i < dim; ++i)
        for (unsigned int j = 0; j < dim; ++j)
          G[i][j] = DF_t[i] * DF_t[j];

      return (this->times_T_t(invert(G)));
    }
}


template <int order, int dim, int spacedim, typename Number>
inline std::size_t
DerivativeForm<order, dim, spacedim, Number>::memory_consumption()
{
  return sizeof(DerivativeForm<order, dim, spacedim, Number>);
}

#endif // DOXYGEN



template <int spacedim, int dim, typename Number>
inline Tensor<1, spacedim, Number>
apply_transformation(const DerivativeForm<1, dim, spacedim, Number> &DF,
                     const Tensor<1, dim, Number> &                  T)
{
  Tensor<1, spacedim, Number> dest;
  for (unsigned int i = 0; i < spacedim; ++i)
    dest[i] = DF[i] * T;
  return dest;
}



// rank=2
template <int spacedim, int dim, typename Number>
inline DerivativeForm<1, spacedim, dim>
apply_transformation(const DerivativeForm<1, dim, spacedim, Number> &DF,
                     const Tensor<2, dim, Number> &                  T)
{
  DerivativeForm<1, spacedim, dim> dest;
  for (unsigned int i = 0; i < dim; ++i)
    dest[i] = apply_transformation(DF, T[i]);

  return dest;
}

template <int spacedim, int dim, typename Number>
inline Tensor<2, spacedim, Number>
apply_transformation(const DerivativeForm<1, dim, spacedim, Number> &DF1,
                     const DerivativeForm<1, dim, spacedim, Number> &DF2)
{
  Tensor<2, spacedim, Number> dest;

  for (unsigned int i = 0; i < spacedim; ++i)
    dest[i] = apply_transformation(DF1, DF2[i]);

  return dest;
}


template <int dim, int spacedim, typename Number>
inline DerivativeForm<1, spacedim, dim, Number>
transpose(const DerivativeForm<1, dim, spacedim, Number> &DF)
{
  DerivativeForm<1, spacedim, dim, Number> tt;
  tt = DF.transpose();
  return tt;
}


DEAL_II_NAMESPACE_CLOSE

#endif