packaged_operation.h

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

#ifndef dealii_packaged_operation_h
#define dealii_packaged_operation_h

#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/lac/vector_memory.h>

#include <functional>

DEAL_II_NAMESPACE_OPEN

// Forward declarations:
template <typename Number> class Vector;
template <typename Range, typename Domain, typename Payload> class LinearOperator;
template <typename Range = Vector<double> > class PackagedOperation;


template <typename Range> class PackagedOperation
{
public:

  PackagedOperation ()
  {
    apply = [](Range &)
    {
      Assert(false,
             ExcMessage("Uninitialized PackagedOperation<Range>::apply called"));
    };

    apply_add = [](Range &)
    {
      Assert(false,
             ExcMessage("Uninitialized PackagedOperation<Range>::apply_add called"));
    };

    reinit_vector = [](Range &, bool)
    {
      Assert(false,
             ExcMessage("Uninitialized PackagedOperation<Range>::reinit_vector "
                        "method called"));
    };
  }

  PackagedOperation (const PackagedOperation<Range> &) = default;

  PackagedOperation (const Range &u)
  {
    *this = u;
  }

  PackagedOperation<Range> &operator=(const PackagedOperation<Range> &) = default;

  PackagedOperation<Range> &operator=(const Range &u)
  {
    apply = [&u](Range &v)
    {
      v = u;
    };

    apply_add = [&u](Range &v)
    {
      v += u;
    };

    reinit_vector = [&u](Range &v, bool omit_zeroing_entries)
    {
      v.reinit(u, omit_zeroing_entries);
    };

    return *this;
  }

  operator Range() const
  {
    Range result_vector;

    reinit_vector(result_vector, /*bool omit_zeroing_entries=*/ true);
    apply(result_vector);

    return result_vector;
  }


  PackagedOperation<Range> &operator+=(const PackagedOperation<Range> &second_comp)
  {
    *this = *this + second_comp;
    return *this;
  }

  PackagedOperation<Range> &operator-=(const PackagedOperation<Range> &second_comp)
  {
    *this = *this - second_comp;
    return *this;
  }

  PackagedOperation<Range> &operator+=(const Range &offset)
  {
    *this = *this + PackagedOperation<Range>(offset);
    return *this;
  }

  PackagedOperation<Range> &operator-=(const Range &offset)
  {
    *this = *this - PackagedOperation<Range>(offset);
    return *this;
  }

  PackagedOperation<Range> &operator*=(typename Range::value_type number)
  {
    *this = *this * number;
    return *this;
  }

  std::function<void(Range &v)> apply;

  std::function<void(Range &v)> apply_add;

  std::function<void(Range &v, bool omit_zeroing_entries)> reinit_vector;
};



template <typename Range>
PackagedOperation<Range>
operator+(const PackagedOperation<Range> &first_comp,
          const PackagedOperation<Range> &second_comp)
{
  PackagedOperation<Range> return_comp;

  return_comp.reinit_vector = first_comp.reinit_vector;

  // ensure to have valid PackagedOperation objects by catching first_comp and
  // second_comp by value

  return_comp.apply = [first_comp, second_comp](Range &v)
  {
    first_comp.apply(v);
    second_comp.apply_add(v);
  };

  return_comp.apply_add = [first_comp, second_comp](Range &v)
  {
    first_comp.apply_add(v);
    second_comp.apply_add(v);
  };

  return return_comp;
}

template <typename Range>
PackagedOperation<Range>
operator-(const PackagedOperation<Range> &first_comp,
          const PackagedOperation<Range> &second_comp)
{
  PackagedOperation<Range> return_comp;

  return_comp.reinit_vector = first_comp.reinit_vector;

  // ensure to have valid PackagedOperation objects by catching first_comp and
  // second_comp by value

  return_comp.apply = [first_comp, second_comp](Range &v)
  {
    second_comp.apply(v);
    v *= -1.;
    first_comp.apply_add(v);
  };

  return_comp.apply_add = [first_comp, second_comp](Range &v)
  {
    first_comp.apply_add(v);
    v *= -1.;
    second_comp.apply_add(v);
    v *= -1.;
  };

  return return_comp;
}

template <typename Range>
PackagedOperation<Range>
operator*(const PackagedOperation<Range> &comp,
          typename Range::value_type number)
{
  PackagedOperation<Range> return_comp;

  return_comp.reinit_vector = comp.reinit_vector;

  // the trivial case: number is zero
  if (number == 0.)
    {
      return_comp.apply = [](Range &v)
      {
        v = 0.;
      };

      return_comp.apply_add = [](Range &)
      {
      };
    }
  else
    {
      return_comp.apply = [comp, number](Range &v)
      {
        comp.apply(v);
        v *= number;
      };

      return_comp.apply_add = [comp, number](Range &v)
      {
        v /= number;
        comp.apply_add(v);
        v *= number;
      };
    }

  return return_comp;
}

template <typename Range>
PackagedOperation<Range>
operator*(typename Range::value_type number,
          const PackagedOperation<Range> &comp)
{
  return comp * number;
}

template <typename Range>
PackagedOperation<Range> operator+(const PackagedOperation<Range> &comp,
                                   const Range &offset)
{
  return comp + PackagedOperation<Range>(offset);
}

template <typename Range>
PackagedOperation<Range> operator+(const Range &offset,
                                   const PackagedOperation<Range> &comp)
{
  return PackagedOperation<Range>(offset) + comp;
}

template <typename Range>
PackagedOperation<Range> operator-(const PackagedOperation<Range> &comp,
                                   const Range &offset)
{
  return comp - PackagedOperation<Range>(offset);
}


template <typename Range>
PackagedOperation<Range> operator-(const Range &offset,
                                   const PackagedOperation<Range> &comp)
{
  return PackagedOperation<Range>(offset) - comp;
}




namespace
{
  // Poor man's trait class that determines whether type T is a vector:
  // FIXME: Implement this as a proper type trait - similar to
  // isBlockVector

  template <typename T>
  class has_vector_interface
  {
    template <typename C>
    static std::false_type test(...);

    template <typename C>
    static std::true_type test(decltype(&C::operator+=),
                               decltype(&C::operator-=),
                               decltype(&C::l2_norm));

  public:
    // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
    // otherwise it is std::false_type

    typedef decltype(test<T>(nullptr, nullptr, nullptr)) type;
  };
}


template <typename Range,
          typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
PackagedOperation<Range> operator+(const Range &u, const Range &v)
{
  PackagedOperation<Range> return_comp;

  // ensure to have valid PackagedOperation objects by catching op by value
  // u is caught by reference

  return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries)
  {
    x.reinit(u, omit_zeroing_entries);
  };

  return_comp.apply = [&u, &v](Range &x)
  {
    x = u;
    x += v;
  };

  return_comp.apply_add = [&u, &v](Range &x)
  {
    x += u;
    x += v;
  };

  return return_comp;
}


template <typename Range,
          typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
PackagedOperation<Range> operator-(const Range &u, const Range &v)
{
  PackagedOperation<Range> return_comp;

  // ensure to have valid PackagedOperation objects by catching op by value
  // u is caught by reference

  return_comp.reinit_vector = [&u](Range &x, bool omit_zeroing_entries)
  {
    x.reinit(u, omit_zeroing_entries);
  };

  return_comp.apply = [&u, &v](Range &x)
  {
    x = u;
    x -= v;
  };

  return_comp.apply_add = [&u, &v](Range &x)
  {
    x += u;
    x -= v;
  };

  return return_comp;
}


template <typename Range,
          typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
PackagedOperation<Range> operator*(const Range &u,
                                   typename Range::value_type number)
{
  return PackagedOperation<Range>(u) * number;
}


template <typename Range,
          typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
PackagedOperation<Range> operator*(typename Range::value_type number,
                                   const Range &u)
{
  return number * PackagedOperation<Range>(u);
}


template <typename Range, typename Domain, typename Payload>
PackagedOperation<Range>
operator*(const LinearOperator<Range, Domain, Payload> &op,
          const Domain &u)
{
  PackagedOperation<Range> return_comp;

  return_comp.reinit_vector = op.reinit_range_vector;

  // ensure to have valid PackagedOperation objects by catching op by value
  // u is caught by reference

  return_comp.apply = [op, &u](Range &v)
  {
    op.vmult(v, u);
  };

  return_comp.apply_add = [op, &u](Range &v)
  {
    op.vmult_add(v, u);
  };

  return return_comp;
}


template <typename Range, typename Domain, typename Payload>
PackagedOperation<Domain>
operator*(const Range &u,
          const LinearOperator<Range, Domain, Payload> &op)
{
  PackagedOperation<Range> return_comp;

  return_comp.reinit_vector = op.reinit_domain_vector;

  // ensure to have valid PackagedOperation objects by catching op by value
  // u is caught by reference

  return_comp.apply = [op, &u](Domain &v)
  {
    op.Tvmult(v, u);
  };

  return_comp.apply_add = [op, &u](Domain &v)
  {
    op.Tvmult_add(v, u);
  };

  return return_comp;
}


template <typename Range, typename Domain, typename Payload>
PackagedOperation<Range>
operator*(const LinearOperator<Range, Domain, Payload> &op,
          const PackagedOperation<Domain> &comp)
{
  PackagedOperation<Range> return_comp;

  return_comp.reinit_vector = op.reinit_range_vector;

  // ensure to have valid PackagedOperation objects by catching op by value
  // u is caught by reference

  return_comp.apply = [op, comp](Domain &v)
  {
    GrowingVectorMemory<Range> vector_memory;

    typename VectorMemory<Range>::Pointer i (vector_memory);
    op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/ true);

    comp.apply(*i);
    op.vmult(v, *i);
  };

  return_comp.apply_add = [op, comp](Domain &v)
  {
    GrowingVectorMemory<Range> vector_memory;

    typename VectorMemory<Range>::Pointer i (vector_memory);
    op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/ true);

    comp.apply(*i);
    op.vmult_add(v, *i);
  };

  return return_comp;
}


template <typename Range, typename Domain, typename Payload>
PackagedOperation<Domain>
operator*(const PackagedOperation<Range> &comp,
          const LinearOperator<Range, Domain, Payload> &op)
{
  PackagedOperation<Range> return_comp;

  return_comp.reinit_vector = op.reinit_domain_vector;

  // ensure to have valid PackagedOperation objects by catching op by value
  // u is caught by reference

  return_comp.apply = [op, comp](Domain &v)
  {
    GrowingVectorMemory<Range> vector_memory;

    typename VectorMemory<Range>::Pointer i (vector_memory);
    op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/ true);

    comp.apply(*i);
    op.Tvmult(v, *i);
  };

  return_comp.apply_add = [op, comp](Domain &v)
  {
    GrowingVectorMemory<Range> vector_memory;

    typename VectorMemory<Range>::Pointer i (vector_memory);
    op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/ true);

    comp.apply(*i);
    op.Tvmult_add(v, *i);
  };

  return return_comp;
}


DEAL_II_NAMESPACE_CLOSE

#endif