diff --git a/tests/python_mpi_projection_tests.py b/tests/python_mpi_projection_tests.py
new file mode 100644
index 0000000..aed2693
--- /dev/null
+++ b/tests/python_mpi_projection_tests.py
@@ -0,0 +1,128 @@
+#! /usr/bin/env python3
+# -*- coding:utf-8 -*-
+"""
+file   python_projection_tests.py
+
+@author Till Junge <till.junge@altermail.ch>
+
+@date   18 Jan 2018
+
+@brief  compare µSpectre's projection operators to GooseFFT
+
+@section LICENSE
+
+Copyright © 2018 Till Junge
+
+µSpectre is free software; you can redistribute it and/or
+modify it under the terms of the GNU General Public License as
+published by the Free Software Foundation, either version 3, or (at
+your option) any later version.
+
+µSpectre 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
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with GNU Emacs; see the file COPYING. If not, write to the
+Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.
+"""
+
+import unittest
+import numpy as np
+import itertools
+
+from mpi4py import MPI
+
+from python_test_imports import µ
+from python_goose_ref import (SmallStrainProjectionGooseFFT,
+                              FiniteStrainProjectionGooseFFT)
+from muSpectre import Formulation
+from muSpectre.fft import Projection
+
+
+def build_test_classes(formulation, RefProjection, fft):
+    class ProjectionCheck(unittest.TestCase):
+        def __init__(self, methodName='runTest'):
+            super().__init__(methodName)
+
+        def setUp(self):
+            self.ref = RefProjection
+            self.resolution = self.ref.resolution
+            self.ndim = self.ref.ndim
+            self.shape = list((self.resolution for _ in range(self.ndim)))
+            self.projection = Projection(self.shape, self.shape, formulation,
+                                         fft, MPI.COMM_WORLD)
+            self.projection.initialise()
+            self.tol = 1e-12 * np.prod(self.shape)
+
+        def test_projection_result(self):
+            # create a bogus strain field in GooseFFT format
+            # dim × dim × N × N (× N)
+            strain_shape = (self.ndim, self.ndim, *self.shape)
+            strain = np.arange(np.prod(strain_shape)).reshape(strain_shape)
+            # if we're testing small strain projections, it needs to be symmetric
+            if self.projection.get_formulation() == µ.Formulation.small_strain:
+                strain += strain.transpose(1, 0, *range(2, len(strain.shape)))
+            strain_g = strain.copy()
+            b_g = self.ref.G(strain_g).reshape(strain_g.shape)
+            strain_µ = np.zeros((*self.shape, self.ndim, self.ndim))
+            for ijk in itertools.product(range(self.resolution), repeat=self.ndim):
+                index_µ = tuple((*ijk, slice(None), slice(None)))
+                index_g = tuple((slice(None), slice(None), *ijk))
+                strain_µ[index_µ] = strain_g[index_g].T
+            res = self.projection.get_subdomain_resolutions()
+            loc = self.projection.get_subdomain_locations()
+            if self.ref.ndim == 2:
+                resx, resy = res
+                locx, locy = loc
+                subdomain_strain_µ = strain_µ[locx:locx+resx, locy:locy+resy]
+            else:
+                resx, resy, resz = res
+                locx, locy, locz = loc
+                subdomain_strain_µ = strain_µ[locx:locx+resx, locy:locy+resy,
+                                              locz:locz+resz]
+            b_µ = self.projection.apply_projection(subdomain_strain_µ.reshape(
+                np.prod(res), self.ndim**2).T).T.reshape(subdomain_strain_µ.shape)
+            for l in range(np.prod(res)):
+                ijk = µ.get_domain_ccoord(res, l)
+                index_µ = tuple((*ijk, slice(None), slice(None)))
+                ijk = loc + np.array(ijk)
+                index_g = tuple((slice(None), slice(None), *ijk))
+                b_µ_sl = b_µ[index_µ].T
+                b_g_sl = b_g[index_g]
+                error = np.linalg.norm(b_µ_sl - b_g_sl)
+                condition = error < self.tol
+                slice_printer = lambda tup: "({})".format(
+                    ", ".join("{}".format(":" if val == slice(None) else val) for val in tup))
+                if not condition:
+                    print("error = {}, tol = {}".format(error, self.tol))
+                    print("b_µ{} =\n{}".format(slice_printer(index_µ), b_µ_sl))
+                    print("b_g{} =\n{}".format(slice_printer(index_g), b_g_sl))
+                self.assertTrue(condition)
+
+    return ProjectionCheck
+
+get_goose = lambda ndim, proj_type: proj_type(
+    ndim, 5, 2, 70e9, .33, 3.)
+get_finite_goose = lambda ndim: get_goose(ndim, FiniteStrainProjectionGooseFFT)
+get_small_goose = lambda ndim: get_goose(ndim, SmallStrainProjectionGooseFFT)
+
+
+small_default_fftwmpi_3 = build_test_classes(Formulation.small_strain,
+                                             get_small_goose(3),
+                                             "fftwmpi")
+small_default_fftwmpi_2 = build_test_classes(Formulation.small_strain,
+                                             get_small_goose(2),
+                                             "fftwmpi")
+
+finite_fast_fftwmpi_3 = build_test_classes(Formulation.finite_strain,
+                                           get_finite_goose(3),
+                                           "fftwmpi")
+finite_fast_fftwmpi_2 = build_test_classes(Formulation.finite_strain,
+                                           get_finite_goose(2),
+                                           "fftwmpi")
+
+if __name__ == "__main__":
+    unittest.main()