derivative_form.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2013 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_derivative_form_h
#define dealii_derivative_form_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/tensor.h>
 
#include <cstddef>
 
 
DEAL_II_NAMESPACE_OPEN
 
template <int order, int dim, int spacedim, typename Number = double>
class DerivativeForm
{
public:
  DerivativeForm() = default;
 
  DerivativeForm(const Tensor<order + 1, dim, Number> &);
 
  DerivativeForm(const Tensor<1, spacedim, Tensor<order, 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<order, spacedim, Tensor<1, dim, Number>> &);
 
  DerivativeForm &
  operator=(const Tensor<1, dim, Number> &);
 
  template <typename OtherNumber>
  DerivativeForm &
  operator=(const DerivativeForm<order, dim, spacedim, OtherNumber> &df);
 
  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;
 
 
  Tensor<2, dim, Number>
  first_fundamental_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(
  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(
  const Tensor<1, spacedim, Tensor<order, dim, Number>> &T)
{
  for (unsigned int j = 0; j < spacedim; ++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<order, spacedim, Tensor<1, dim, Number>> &T)
{
  for (unsigned int j = 0; j < spacedim; ++j)
    (*this)[j] = T[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>
template <typename OtherNumber>
inline DerivativeForm<order, dim, spacedim, Number> &
DerivativeForm<order, dim, spacedim, Number>::operator=(
  const DerivativeForm<order, dim, spacedim, OtherNumber> &df)
{
  for (unsigned int j = 0; j < spacedim; ++j)
    (*this)[j] = df[j];
  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)
{
  AssertIndexRange(i, 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
{
  AssertIndexRange(i, 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."));
      return (std::sqrt(::determinant(first_fundamental_form())));
    }
}
 
 
 
template <int order, int dim, int spacedim, typename Number>
inline Tensor<2, dim, Number>
DerivativeForm<order, dim, spacedim, Number>::first_fundamental_form() const
{
  Assert(order == 1, ExcMessage("Only for order == 1."));
  const DerivativeForm<1, spacedim, dim, Number> DF_t = this->transpose();
 
  Tensor<2, dim, Number> G;
  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 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, Number>(DF_t);
    }
  else
    {
      return (this->times_T_t(invert(first_fundamental_form())));
    }
}
 
 
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 order, int dim, int spacedim, typename Number>
inline std::ostream &
operator<<(std::ostream                                       &out,
           const DerivativeForm<order, dim, spacedim, Number> &df)
{
  for (unsigned int i = 0; i < spacedim; ++i)
    {
      out << df[i];
      if (i != spacedim - 1)
        for (unsigned int j = 0; j < order + 1; ++j)
          out << ' ';
    }
 
  return out;
}
 
 
 
template <int spacedim, int dim, typename Number1, typename Number2>
inline Tensor<1, spacedim, typename ProductType<Number1, Number2>::type>
apply_transformation(const DerivativeForm<1, dim, spacedim, Number1> &grad_F,
                     const Tensor<1, dim, Number2>                   &d_x)
{
  Tensor<1, spacedim, typename ProductType<Number1, Number2>::type> dest;
  for (unsigned int i = 0; i < spacedim; ++i)
    dest[i] = grad_F[i] * d_x;
  return dest;
}
 
 
 
// rank=2
template <int spacedim, int dim, typename Number1, typename Number2>
inline DerivativeForm<1,
                      spacedim,
                      dim,
                      typename ProductType<Number1, Number2>::type>
apply_transformation(const DerivativeForm<1, dim, spacedim, Number1> &grad_F,
                     const Tensor<2, dim, Number2>                   &D_X)
{
  DerivativeForm<1, spacedim, dim, typename ProductType<Number1, Number2>::type>
    dest;
  for (unsigned int i = 0; i < dim; ++i)
    dest[i] = apply_transformation(grad_F, D_X[i]);
 
  return dest;
}
 
 
 
// rank=2
template <int dim, typename Number1, typename Number2>
inline Tensor<2, dim, typename ProductType<Number1, Number2>::type>
apply_transformation(const DerivativeForm<1, dim, dim, Number1> &grad_F,
                     const Tensor<2, dim, Number2>              &D_X)
{
  Tensor<2, dim, typename ProductType<Number1, Number2>::type> dest;
  for (unsigned int i = 0; i < dim; ++i)
    dest[i] = apply_transformation(grad_F, D_X[i]);
 
  return dest;
}
 
 
 
template <int spacedim,
          int dim,
          int n_components,
          typename Number1,
          typename Number2>
inline Tensor<1,
              n_components,
              Tensor<1, spacedim, typename ProductType<Number1, Number2>::type>>
apply_transformation(
  const DerivativeForm<1, dim, spacedim, Number1>        &grad_F,
  const Tensor<1, n_components, Tensor<1, dim, Number2>> &D_X)
{
  Tensor<1,
         n_components,
         Tensor<1, spacedim, typename ProductType<Number1, Number2>::type>>
    dest;
  for (unsigned int i = 0; i < n_components; ++i)
    dest[i] = apply_transformation(grad_F, D_X[i]);
 
  return dest;
}
 
 
 
template <int spacedim, int dim, typename Number1, typename Number2>
inline Tensor<2, spacedim, typename ProductType<Number1, Number2>::type>
apply_transformation(const DerivativeForm<1, dim, spacedim, Number1> &DF1,
                     const DerivativeForm<1, dim, spacedim, Number2> &DF2)
{
  Tensor<2, spacedim, typename ProductType<Number1, Number2>::type> 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)
{
  return DF.transpose();
}
 
 
 
template <int spacedim, int dim, typename Number1, typename Number2>
inline Tensor<1, spacedim, typename ProductType<Number1, Number2>::type>
apply_diagonal_transformation(
  const DerivativeForm<1, dim, spacedim, Number1> &grad_F,
  const Tensor<1, dim, Number2>                   &d_x)
{
  Assert(dim == spacedim,
         ExcMessage("Only dim = spacedim allowed for diagonal transformation"));
  Tensor<1, spacedim, typename ProductType<Number1, Number2>::type> dest;
  for (unsigned int i = 0; i < spacedim; ++i)
    dest[i] = grad_F[i][i] * d_x[i];
  return dest;
}
 
 
template <int dim, typename Number1, typename Number2>
inline Tensor<2, dim, typename ProductType<Number1, Number2>::type>
apply_diagonal_transformation(
  const DerivativeForm<1, dim, dim, Number1> &grad_F,
  const Tensor<2, dim, Number2>              &D_X)
{
  Tensor<2, dim, typename ProductType<Number1, Number2>::type> dest;
  for (unsigned int i = 0; i < dim; ++i)
    dest[i] = apply_diagonal_transformation(grad_F, D_X[i]);
 
  return dest;
}
 
 
 
template <int spacedim,
          int dim,
          int n_components,
          typename Number1,
          typename Number2>
inline Tensor<1,
              n_components,
              Tensor<1, spacedim, typename ProductType<Number1, Number2>::type>>
apply_diagonal_transformation(
  const DerivativeForm<1, dim, spacedim, Number1>        &grad_F,
  const Tensor<1, n_components, Tensor<1, dim, Number2>> &D_X)
{
  Tensor<1,
         n_components,
         Tensor<1, spacedim, typename ProductType<Number1, Number2>::type>>
    dest;
  for (unsigned int i = 0; i < n_components; ++i)
    dest[i] = apply_diagonal_transformation(grad_F, D_X[i]);
 
  return dest;
}
 
 
 
// rank=2
template <int spacedim, int dim, typename Number1, typename Number2>
inline DerivativeForm<1,
                      spacedim,
                      dim,
                      typename ProductType<Number1, Number2>::type>
apply_diagonal_transformation(
  const DerivativeForm<1, dim, spacedim, Number1> &grad_F,
  const Tensor<2, dim, Number2>                   &D_X)
{
  DerivativeForm<1, spacedim, dim, typename ProductType<Number1, Number2>::type>
    dest;
  for (unsigned int i = 0; i < dim; ++i)
    dest[i] = apply_diagonal_transformation(grad_F, D_X[i]);
 
  return dest;
}
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif