diff --git a/src/model/hooke.cpp b/src/model/hooke.cpp
index 3fcb61b..3a74d8e 100644
--- a/src/model/hooke.cpp
+++ b/src/model/hooke.cpp
@@ -1,75 +1,95 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "hooke.hh"
 #include "influence.hh"
 #include "loop.hh"
 #include "model.hh"
+#include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 template <model_type type>
 void Hooke<type>::apply(GridBase<Real>& strain, GridBase<Real>& stress) const {
   constexpr auto dim = model_type_traits<type>::dimension;
   const influence::ElasticHelper<dim> elasticity(
       this->model->getShearModulus(), this->model->getPoissonRatio());
 
   // Checking incompressibility
   if (elasticity.nu == 0.5)
     TAMAAS_EXCEPTION("Incompressibility error");
 
   // Checking number of components
   if (strain.getNbComponents() == dim * dim) {
     Loop::loop(
         [elasticity] CUDA_LAMBDA(MatrixProxy<Real, dim, dim> sigma,
                                  MatrixProxy<const Real, dim, dim> epsilon) {
           sigma = elasticity(epsilon);
         },
         range<MatrixProxy<Real, dim, dim>>(stress),
         range<MatrixProxy<const Real, dim, dim>>(strain));
   }
 
   // in case we have voigt notation
   else if (strain.getNbComponents() == voigt_size<dim>::value) {
     Loop::loop(
         [elasticity] CUDA_LAMBDA(SymMatrixProxy<Real, dim> sigma,
                                  SymMatrixProxy<const Real, dim> epsilon) {
           sigma = elasticity(epsilon);
         },
         range<SymMatrixProxy<Real, dim>>(stress),
         range<SymMatrixProxy<const Real, dim>>(strain));
   }
 
   else
     TAMAAS_EXCEPTION("Strain components do not match dimension");
 }
 
+/* -------------------------------------------------------------------------- */
+template <model_type type>
+std::pair<UInt, UInt> Hooke<type>::matvecShape() const {
+  auto shape = this->model->getDiscretization();
+  auto N =
+    std::accumulate(shape.begin(), shape.end(), 6, std::multiplies<UInt>());
+  return std::make_pair(N, N);
+}
+
+template <model_type type>
+GridBase<Real> Hooke<type>::matvec(GridBase<Real>& X) const {
+  constexpr auto dim = model_type_traits<type>::dimension;
+  X.setNbComponents(voigt_size<dim>::value);
+  GridBase<Real> x(X.getGlobalNbPoints(), X.getNbComponents());
+  apply(X, x);
+  return x;
+}
+
+
 /* -------------------------------------------------------------------------- */
 template class Hooke<model_type::basic_1d>;
 template class Hooke<model_type::basic_2d>;
 template class Hooke<model_type::surface_1d>;
 template class Hooke<model_type::surface_2d>;
 template class Hooke<model_type::volume_1d>;
 template class Hooke<model_type::volume_2d>;
 
 }  // namespace tamaas
diff --git a/src/model/hooke.hh b/src/model/hooke.hh
index 34a2516..864b576 100644
--- a/src/model/hooke.hh
+++ b/src/model/hooke.hh
@@ -1,57 +1,63 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef HOOKE_HH
 #define HOOKE_HH
 /* -------------------------------------------------------------------------- */
 #include "integral_operator.hh"
 #include "model_type.hh"
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /**
  * @brief Applies Hooke's law of elasticity
  */
 template <model_type type>
 class Hooke : public IntegralOperator {
 public:
   using IntegralOperator::IntegralOperator;
 
   /// Type of underlying model
   model_type getType() const override { return type; }
 
   /// Hooke's law computes stresses ~ Neumann
   IntegralOperator::kind getKind() const override {
     return IntegralOperator::dirac;
   }
 
   /// Does not update
   void updateFromModel() override {}
 
   /// Apply Hooke's tensor
   void apply(GridBase<Real>& strain, GridBase<Real>& stress) const override;
+
+  /// LinearOperator shape
+  std::pair<UInt, UInt> matvecShape() const override;
+
+  /// LinearOperator interface
+  GridBase<Real> matvec(GridBase<Real>& X) const override;
 };
 
 }  // namespace tamaas
 
 #endif  // HOOKE_HH