Re: [deal.II] Create a LinearOperator with a function instead of a matrix

[Matthias Maier]

Hi,

the following example (which is a simplified version of test
lac/linear_operator_01) should get you started. In short: Simply create
an empty LinearOperator object (with the appropriate template
parameters) and populate the corresponding std::function objects.

Best,
Matthias



#include <deal.II/lac/linear_operator.h>
#include <deal.II/lac/vector_memory.templates.h>

using namespace dealii;

struct LeftVector {
  typedef double value_type;
  value_type value;

  LeftVector &operator=(value_type new_value)
  {
    value = new_value;
    return *this;
  }
  LeftVector &operator*=(value_type scale)
  {
    value *= scale;
    return *this;
  }
  LeftVector &operator/=(value_type scale)
  {
    value /= scale;
    return *this;
  }
  LeftVector &operator+=(const LeftVector &u)
  {
    value += u.value;
    return *this;
  }
  int size() const
  {
    return 1;
  }
  std::size_t memory_consumption() const
  {
    return 1;
  }
};

struct RightVector {
  typedef double value_type;
  value_type value;

  RightVector &operator=(value_type new_value)
  {
    value = new_value;
    return *this;
  }
  RightVector &operator*=(value_type scale)
  {
    value *= scale;
    return *this;
  }
  RightVector &operator/=(value_type scale)
  {
    value /= scale;
    return *this;
  }
  RightVector &operator+=(const RightVector &u)
  {
    value += u.value;
    return *this;
  }
  int size() const
  {
    return 1;
  }
  std::size_t memory_consumption() const
  {
    return 1;
  }
};

int main()
{
  LinearOperator<LeftVector, RightVector> multiply2;

  multiply2.vmult = [](LeftVector &v, const RightVector &u) {
    v.value = 2 * u.value;
  };
  multiply2.vmult_add = [](LeftVector &v, const RightVector &u) {
    v.value += 2 * u.value;
  };
  multiply2.Tvmult = [](RightVector &v, const LeftVector &u) {
    v.value = 2 * u.value;
  };
  multiply2.Tvmult_add = [](RightVector &v, const LeftVector &u) {
    v.value += 2 * u.value;
  };
  multiply2.reinit_range_vector = [](LeftVector &, bool omit_zeroing_values) {
    // do nothing
  };
  multiply2.reinit_domain_vector = [](RightVector &, bool omit_zeroing_values) {
    // do nothing
  };

  // Small test:

  RightVector u = {4.};
  LeftVector v = {0.};

  multiply2.vmult(v, u);
  std::cout << "2 * " << u.value << " = " << v.value << std::endl;

  multiply2.vmult_add(v, u);
  std::cout << "... + 2 * " << u.value << " = " << v.value << std::endl;

  return 0;
}

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/87h87mar6y.fsf%4043-1.org.
For more options, visit https://groups.google.com/d/optout.