packaged_operation.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) 2015 - 2023 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_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:
#ifndef DOXYGEN
template <typename Number>
class Vector;
template <typename Range, typename Domain, typename Payload>
class LinearOperator;
template <typename Range = Vector<double>>
class PackagedOperation;
#endif
 
 
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 = std::enable_if_t<internal::PackagedOperationImplementation::
                                has_vector_interface<Range>::type::value>>
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 = std::enable_if_t<internal::PackagedOperationImplementation::
                                has_vector_interface<Range>::type::value>>
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 = std::enable_if_t<internal::PackagedOperationImplementation::
                                has_vector_interface<Range>::type::value>>
PackagedOperation<Range>
operator*(const Range &u, typename Range::value_type number)
{
  return PackagedOperation<Range>(u) * number;
}
 
 
template <
  typename Range,
  typename = std::enable_if_t<internal::PackagedOperationImplementation::
                                has_vector_interface<Range>::type::value>>
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