diff --git a/tests/test_matvec.py b/tests/test_matvec.py
index 859cc73..baf4aea 100644
--- a/tests/test_matvec.py
+++ b/tests/test_matvec.py
@@ -1,64 +1,98 @@
 # -*- coding: utf-8 -*-
 #
 # 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/>.
 
 import pytest
 import numpy as np
 import tamaas as tm
 
 from numpy.testing import assert_allclose
 from test_dumper import model_fixture  # noqa
 
 
 linalg = pytest.importorskip('scipy.sparse.linalg')
 
 
 def test_matvec_westergaard(model_fixture):
     m = model_fixture
     m.be_engine.registerNeumann()
     op = m.operators['Westergaard::neumann']
     linop = linalg.aslinearoperator(op)
 
     x = [np.linspace(0, L, N, endpoint=False)
          for (L, N) in zip(m.boundary_system_size, m.boundary_shape)]
 
     coords = np.meshgrid(*x, indexing='ij')
 
     t = m.traction
     res = np.zeros_like(m.displacement)
     volume = len(m.shape) != len(m.boundary_shape)
     basic = m.type in (tm.model_type.basic_1d, tm.model_type.basic_2d)
 
     # Adapting to special cases
     if basic:
         t = t[..., np.newaxis]
         res = res[..., np.newaxis]
 
     t[:] = 0
     t[..., -1] = np.sin(2 * np.pi * coords[0])
 
     # Computing using native API
     op(t, res)
 
     # Only look at surface displacement
     if volume:
         res = res[0]
 
     # Computing using LinearOperator API
     y = np.reshape(linop * np.ravel(t), res.shape)
     assert_allclose(res, y)
+
+
+def test_matvec_hooke():
+    m = tm.Model(tm.model_type.volume_2d, [1, 1, 1], [2, 2, 2])
+    op = m.operators['hooke']
+    linop = linalg.aslinearoperator(op)
+
+    strain = np.zeros(m.shape + [6], dtype=tm.dtype)
+    strain[..., 3:] = 1
+
+    stress = np.zeros_like(strain)
+    op(strain, stress)
+
+    res = linop * np.ravel(strain)
+
+    assert_allclose(res, np.ravel(stress))
+    assert_allclose(res, np.ravel(strain) * m.mu * 2)
+
+
+def test_matvec_mindlin():
+    m = tm.Model(tm.model_type.volume_2d, [1, 1, 1], [2, 2, 2])
+    tm.ModelFactory.registerVolumeOperators(m)
+    op = m.operators['mindlin_gradient']
+    linop = linalg.aslinearoperator(op)
+
+    strain = np.zeros(m.shape + [6], dtype=tm.dtype)
+    strain[..., 3:] = 1
+
+    stress = np.zeros_like(strain)
+    op(strain, stress)
+
+    res = linop * np.ravel(strain)
+
+    assert_allclose(res, np.ravel(stress))