derivative_form.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2013 - 2017 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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;

  double 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
double
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> 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