packaged_operation.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2014 - 2019 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_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 internal
{
  namespace PackagedOperationImplementation
  {
    // 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

      using type = decltype(test<T>(nullptr, nullptr, nullptr));
    }; // namespace
  }    // namespace PackagedOperationImplementation
} // namespace internal


template <typename Range,
          typename = typename std::enable_if<
            internal::PackagedOperationImplementation::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<
            internal::PackagedOperationImplementation::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<
            internal::PackagedOperationImplementation::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<
            internal::PackagedOperationImplementation::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