Implement a flexible matrix based operator to be used for reassembling

src/linearsolver.hh

@@ -7,6 +7,42 @@
 
 namespace PPS {
 
+  template<class X>
+  class LinearizedOperatorMixin {
+  public:
+    virtual void linearizeAt(const X &) = 0;
+    virtual ~LinearizedOperatorMixin() = default;
+  };
+
+  template<class Matrix>
+  struct AssembledLinearizedOperatorMatrix
+  {
+    Matrix matrix_;
+  };
+
+  template<class GridOperator, class Matrix, class X, class Y>
+  class AssembledLinearizedOperator :
+    private AssembledLinearizedOperatorMatrix<Matrix>,
+    public Dune::MatrixAdapter<Matrix, X, Y>,
+    public LinearizedOperatorMixin<X>
+  {
+    using MemBase = AssembledLinearizedOperatorMatrix<Matrix>;
+    const GridOperator *go_;
+
+  public:
+    AssembledLinearizedOperator(const GridOperator * go) :
+      MemBase{ go->pattern().template make<Matrix>() },
+      Dune::MatrixAdapter<Matrix, X, Y>(this->matrix_),
+      go_(go)
+    { }
+
+    void linearizeAt(const X &linearizationPoint) override {
+      this->matrix_ = 0;
+      go_->jacobian(linearizationPoint, this->matrix_);
+    }
+
+  };
+
   template<class MatrixType, class VectorType,
            template<class,class,class,int> class Preconditioner,
            template<class> class Solver>