diff --git a/pyop2/codegen/loopycompat.py b/pyop2/codegen/loopycompat.py
index e74d2ce62..ae3d5feff 100644
--- a/pyop2/codegen/loopycompat.py
+++ b/pyop2/codegen/loopycompat.py
@@ -172,7 +172,6 @@ def _match_caller_callee_argument_dimension_(program, callee_function_name):
         invocations would demand complex renaming logic which is not
         implemented yet.
     """
-
     assert isinstance(program, TranslationUnit)
     assert isinstance(callee_function_name, str)
     assert callee_function_name not in program.entrypoints
diff --git a/pyop2/codegen/rep2loopy.py b/pyop2/codegen/rep2loopy.py
index 7085c324d..dbdfed4b2 100644
--- a/pyop2/codegen/rep2loopy.py
+++ b/pyop2/codegen/rep2loopy.py
@@ -552,22 +552,16 @@ def renamer(expr):
         headers = headers | set(["#include <petsclog.h>"])
     preamble = "\n".join(sorted(headers))
 
-    from coffee.base import Node
-    from loopy.kernel.function_interface import CallableKernel
-
     if isinstance(kernel.code, loopy.TranslationUnit):
         knl = kernel.code
         wrapper = loopy.merge([wrapper, knl])
-        names = knl.callables_table
-        for name in names:
-            if isinstance(wrapper.callables_table[name], CallableKernel):
-                wrapper = _match_caller_callee_argument_dimension_(wrapper, name)
+        # remove the local kernel from the available entrypoints
+        wrapper = wrapper.copy(entrypoints=wrapper.entrypoints-{kernel.name})
+        wrapper = _match_caller_callee_argument_dimension_(wrapper, kernel.name)
     else:
         # kernel is a string, add it to preamble
-        if isinstance(kernel.code, Node):
-            code = kernel.code.gencode()
-        else:
-            code = kernel.code
+        assert isinstance(kernel.code, str)
+        code = kernel.code
         wrapper = loopy.register_callable(
             wrapper,
             kernel.name,
diff --git a/pyop2/local_kernel.py b/pyop2/local_kernel.py
index 4807463b8..493107ba8 100644
--- a/pyop2/local_kernel.py
+++ b/pyop2/local_kernel.py
@@ -3,7 +3,6 @@
 import hashlib
 from typing import Union
 
-import coffee
 import loopy as lp
 from loopy.tools import LoopyKeyBuilder
 import numpy as np
@@ -36,8 +35,6 @@ def Kernel(code, name, **kwargs):
     """
     if isinstance(code, str):
         return CStringLocalKernel(code, name, **kwargs)
-    elif isinstance(code, coffee.base.Node):
-        return CoffeeLocalKernel(code, name, **kwargs)
     elif isinstance(code, (lp.LoopKernel, lp.TranslationUnit)):
         return LoopyLocalKernel(code, name, **kwargs)
     else:
@@ -120,9 +117,7 @@ def cache_key(self):
     @cached_property
     def _immutable_cache_key(self):
         # We need this function because self.accesses is mutable due to legacy support
-        if isinstance(self.code, coffee.base.Node):
-            code = self.code.gencode()
-        elif isinstance(self.code, lp.TranslationUnit):
+        if isinstance(self.code, lp.TranslationUnit):
             key_hash = hashlib.sha256()
             self.code.update_persistent_hash(key_hash, LoopyKeyBuilder())
             code = key_hash.hexdigest()
@@ -160,10 +155,7 @@ def num_flops(self):
         if not configuration["compute_kernel_flops"]:
             return 0
 
-        if isinstance(self.code, coffee.base.Node):
-            v = coffee.visitors.EstimateFlops()
-            return v.visit(self.code)
-        elif isinstance(self.code, lp.TranslationUnit):
+        if isinstance(self.code, lp.TranslationUnit):
             op_map = lp.get_op_map(
                 self.code.copy(options=lp.Options(ignore_boostable_into=True),
                                silenced_warnings=['insn_count_subgroups_upper_bound',
@@ -195,8 +187,8 @@ class CStringLocalKernel(LocalKernel):
     """:class:`LocalKernel` class where `code` is a string of C code.
 
     :kwarg dtypes: Iterable of datatypes (either `np.dtype` or `str`) for
-        each kernel argument. This is not required for :class:`CoffeeLocalKernel`
-        or :class:`LoopyLocalKernel` because it can be inferred.
+        each kernel argument. This is not required for :class:`LoopyLocalKernel`
+        because it can be inferred.
 
     All other `__init__` parameters are the same.
     """
@@ -215,23 +207,6 @@ def dtypes(self, dtypes):
         self._dtypes = dtypes
 
 
-class CoffeeLocalKernel(LocalKernel):
-    """:class:`LocalKernel` class where `code` has type :class:`coffee.base.Node`."""
-
-    @validate_type(("code", coffee.base.Node, TypeError))
-    def __init__(self, code, name, accesses=None, dtypes=None, **kwargs):
-        super().__init__(code, name, accesses, **kwargs)
-        self._dtypes = dtypes
-
-    @property
-    def dtypes(self):
-        return self._dtypes
-
-    @dtypes.setter
-    def dtypes(self, dtypes):
-        self._dtypes = dtypes
-
-
 class LoopyLocalKernel(LocalKernel):
     """:class:`LocalKernel` class where `code` has type :class:`loopy.LoopKernel`
         or :class:`loopy.TranslationUnit`.
diff --git a/pyop2/op2.py b/pyop2/op2.py
index 726168e79..65affa065 100644
--- a/pyop2/op2.py
+++ b/pyop2/op2.py
@@ -48,7 +48,7 @@
 from pyop2.types import (READ, WRITE, RW, INC, MIN, MAX,
                          ON_BOTTOM, ON_TOP, ON_INTERIOR_FACETS, ALL)
 
-from pyop2.local_kernel import CStringLocalKernel, LoopyLocalKernel, CoffeeLocalKernel, Kernel  # noqa: F401
+from pyop2.local_kernel import CStringLocalKernel, LoopyLocalKernel, Kernel  # noqa: F401
 from pyop2.global_kernel import (GlobalKernelArg, DatKernelArg, MixedDatKernelArg,  # noqa: F401
                                  MatKernelArg, MixedMatKernelArg, MapKernelArg, GlobalKernel)
 from pyop2.parloop import (GlobalParloopArg, DatParloopArg, MixedDatParloopArg,  # noqa: F401
diff --git a/pyop2/parloop.py b/pyop2/parloop.py
index ac78e6bda..ef62a1878 100644
--- a/pyop2/parloop.py
+++ b/pyop2/parloop.py
@@ -14,7 +14,7 @@
 from pyop2.exceptions import KernelTypeError, MapValueError, SetTypeError
 from pyop2.global_kernel import (GlobalKernelArg, DatKernelArg, MixedDatKernelArg,
                                  MatKernelArg, MixedMatKernelArg, GlobalKernel)
-from pyop2.local_kernel import LocalKernel, CStringLocalKernel, CoffeeLocalKernel, LoopyLocalKernel
+from pyop2.local_kernel import LocalKernel, CStringLocalKernel, LoopyLocalKernel
 from pyop2.types import (Access, Global, AbstractDat, Dat, DatView, MixedDat, Mat, Set,
                          MixedSet, ExtrudedSet, Subset, Map, ComposedMap, MixedMap)
 from pyop2.utils import cached_property
@@ -624,7 +624,7 @@ def LegacyParloop(local_knl, iterset, *args, **kwargs):
 
     # finish building the local kernel
     local_knl.accesses = tuple(a.access for a in args)
-    if isinstance(local_knl, (CStringLocalKernel, CoffeeLocalKernel)):
+    if isinstance(local_knl, CStringLocalKernel):
         local_knl.dtypes = tuple(a.data.dtype for a in args)
 
     global_knl_args = tuple(a.global_kernel_arg for a in args)
diff --git a/requirements-git.txt b/requirements-git.txt
index 438205043..d6f3d2182 100644
--- a/requirements-git.txt
+++ b/requirements-git.txt
@@ -1,2 +1 @@
-git+https://github.com/coneoproject/COFFEE.git#egg=coffee
 git+https://github.com/firedrakeproject/loopy.git@main#egg=loopy
diff --git a/setup.py b/setup.py
index a271e2415..ad9f7815b 100644
--- a/setup.py
+++ b/setup.py
@@ -89,11 +89,8 @@ def get_petsc_dir():
     'decorator',
     'mpi4py',
     'numpy>=1.6',
-    'COFFEE',
 ]
 
-dep_links = ['git+https://github.com/coneoproject/COFFEE#egg=COFFEE-dev']
-
 version = sys.version_info[:2]
 
 if version < (3, 6):
@@ -141,7 +138,6 @@ def run(self):
           'Programming Language :: Python :: 3.6',
       ],
       install_requires=install_requires + test_requires,
-      dependency_links=dep_links,
       packages=['pyop2', 'pyop2.codegen', 'pyop2.types'],
       package_data={
           'pyop2': ['assets/*', '*.h', '*.pxd', '*.pyx', 'codegen/c/*.c']},
diff --git a/test/unit/test_caching.py b/test/unit/test_caching.py
index f175bc76f..3a02778b6 100644
--- a/test/unit/test_caching.py
+++ b/test/unit/test_caching.py
@@ -40,8 +40,6 @@
 from pyop2 import op2, mpi
 from pyop2.caching import disk_cached
 
-from coffee.base import *
-
 
 def _seed():
     return 0.02041724
@@ -419,11 +417,7 @@ def test_vector_map(self, iterset, x2, iter2ind2):
     def test_same_iteration_space_works(self, iterset, x2, iter2ind2):
         self.cache.clear()
         assert len(self.cache) == 0
-        kernel_code = FunDecl("void", "k",
-                              [Decl("int*", c_sym("x"), qualifiers=["unsigned"])],
-                              c_for("i", 1, ""),
-                              pred=["static"])
-        k = op2.Kernel(kernel_code.gencode(), 'k')
+        k = op2.Kernel("""static void k(void *x) {}""", 'k')
 
         op2.par_loop(k, iterset,
                      x2(op2.INC, iter2ind2))
diff --git a/test/unit/test_extrusion.py b/test/unit/test_extrusion.py
index dfae39b60..2ae507d7d 100644
--- a/test/unit/test_extrusion.py
+++ b/test/unit/test_extrusion.py
@@ -72,8 +72,6 @@ def compute_ind_extr(nums,
     return ind
 
 
-from coffee.base import *
-
 # Data type
 valuetype = numpy.float64
 
@@ -333,45 +331,6 @@ def extrusion_kernel():
     return op2.Kernel(kernel_code, "extrusion")
 
 
-@pytest.fixture
-def vol_comp():
-    init = FlatBlock("""
-double area = x[0][0]*(x[2][1]-x[4][1]) + x[2][0]*(x[4][1]-x[0][1])
-           + x[4][0]*(x[0][1]-x[2][1]);
-if (area < 0)
-area = area * (-1.0);
-""")
-    assembly = Incr(Symbol("A", ("i0", "i1")),
-                    FlatBlock("0.5 * area * (x[1][2] - x[0][2])"))
-    assembly = c_for("i0", 6, c_for("i1", 6, assembly))
-    kernel_code = FunDecl("void", "vol_comp",
-                          [Decl("double", Symbol("A", (6, 6))),
-                           Decl("double", Symbol("x", (6, 3)))],
-                          Block([init, assembly], open_scope=False),
-                          pred=["static"])
-    return op2.Kernel(kernel_code.gencode(), "vol_comp")
-
-
-@pytest.fixture
-def vol_comp_rhs():
-    init = FlatBlock("""
-double area = x[0][0]*(x[2][1]-x[4][1]) + x[2][0]*(x[4][1]-x[0][1])
-           + x[4][0]*(x[0][1]-x[2][1]);
-if (area < 0)
-area = area * (-1.0);
-""")
-    assembly = Incr(Symbol("A", ("i0",)),
-                    FlatBlock("0.5 * area * (x[1][2] - x[0][2]) * y[0]"))
-    assembly = c_for("i0", 6, assembly)
-    kernel_code = FunDecl("void", "vol_comp_rhs",
-                          [Decl("double", Symbol("A", (6,))),
-                           Decl("double", Symbol("x", (6, 3))),
-                           Decl("int", Symbol("y", (1,)))],
-                          Block([init, assembly], open_scope=False),
-                          pred=["static"])
-    return op2.Kernel(kernel_code.gencode(), "vol_comp_rhs")
-
-
 class TestExtrusion:
 
     """
@@ -423,7 +382,7 @@ def test_extruded_layer_arg(self, elements, field_map, dat_f):
                      pass_layer_arg=True)
         end = layers - 1
         start = 0
-        ref = np.arange(start, end)
+        ref = numpy.arange(start, end)
         assert [dat_f.data[end*n:end*(n+1)] == ref
                 for n in range(int(len(dat_f.data)/end) - 1)]
 
diff --git a/test/unit/test_indirect_loop.py b/test/unit/test_indirect_loop.py
index ab77d182b..005506710 100644
--- a/test/unit/test_indirect_loop.py
+++ b/test/unit/test_indirect_loop.py
@@ -38,8 +38,6 @@
 from pyop2 import op2
 from pyop2.exceptions import MapValueError
 
-from coffee.base import *
-
 
 nelems = 4096
 
@@ -265,14 +263,11 @@ def test_mixed_non_mixed_dat(self, mdat, mmap, iterset):
     def test_mixed_non_mixed_dat_itspace(self, mdat, mmap, iterset):
         """Increment into a MixedDat from a Dat using iteration spaces."""
         d = op2.Dat(iterset, np.ones(iterset.size))
-        assembly = Incr(Symbol("d", ("j",)), Symbol("x", (0,)))
-        assembly = c_for("j", 2, assembly)
-        kernel_code = FunDecl("void", "inc",
-                              [Decl("double", c_sym("*d")),
-                               Decl("double", c_sym("*x"))],
-                              Block([assembly], open_scope=False),
-                              pred=["static"])
-        op2.par_loop(op2.Kernel(kernel_code.gencode(), "inc"), iterset,
+        kernel_inc = """static void inc(double *d, double *x) {
+          for (int i=0; i<2; ++i)
+            d[i] += x[0];
+        }"""
+        op2.par_loop(op2.Kernel(kernel_inc, "inc"), iterset,
                      mdat(op2.INC, mmap),
                      d(op2.READ))
         assert all(mdat[0].data == 1.0) and mdat[1].data == 4096.0
diff --git a/test/unit/test_iteration_space_dats.py b/test/unit/test_iteration_space_dats.py
index af5b817ca..96f530279 100644
--- a/test/unit/test_iteration_space_dats.py
+++ b/test/unit/test_iteration_space_dats.py
@@ -37,8 +37,6 @@
 
 from pyop2 import op2
 
-from coffee.base import *
-
 
 def _seed():
     return 0.02041724
@@ -106,16 +104,14 @@ def test_sum_nodes_to_edges(self):
         e_map = numpy.array([(i, i + 1)
                             for i in range(nedges)], dtype=numpy.uint32)
         edge2node = op2.Map(edges, nodes, 2, e_map, "edge2node")
-
-        kernel_sum = FunDecl("void", "sum",
-                             [Decl(
-                                 "int*", c_sym("edge"), qualifiers=["unsigned"]),
-                              Decl(
-                                  "int*", c_sym("nodes"), qualifiers=["unsigned"])],
-                             c_for("i", 2, Incr(c_sym("*edge"), Symbol("nodes", ("i",)))),
-                             pred=["static"])
-
-        op2.par_loop(op2.Kernel(kernel_sum.gencode(), "sum"), edges,
+        kernel_sum = """
+static void sum(unsigned int *edge, unsigned int *nodes) {
+  for (int i=0; i<2; ++i)
+    edge[0] += nodes[i];
+}
+        """
+
+        op2.par_loop(op2.Kernel(kernel_sum, "sum"), edges,
                      edge_vals(op2.INC),
                      node_vals(op2.READ, edge2node))
 
@@ -124,24 +120,27 @@ def test_sum_nodes_to_edges(self):
 
     def test_read_1d_itspace_map(self, node, d1, vd1, node2ele):
         vd1.data[:] = numpy.arange(nele)
-        k = FunDecl("void", "k",
-                    [Decl("int*", c_sym("d")), Decl("int*", c_sym("vd"))],
-                    c_for("i", 1, Assign(Symbol("d", (0,)), Symbol("vd", ("i",)))),
-                    pred=["static"])
-
-        op2.par_loop(op2.Kernel(k.gencode(), 'k'), node,
+        k = """
+static void k(int *d, int *vd) {
+  for (int i=0; i<1; ++i)
+    d[0] = vd[i];
+}
+        """
+        op2.par_loop(op2.Kernel(k, 'k'), node,
                      d1(op2.WRITE),
                      vd1(op2.READ, node2ele))
         assert all(d1.data[::2] == vd1.data)
         assert all(d1.data[1::2] == vd1.data)
 
     def test_write_1d_itspace_map(self, node, vd1, node2ele):
-        k = FunDecl("void", "k",
-                    [Decl("int*", c_sym("vd"))],
-                    c_for("i", 1, Assign(Symbol("vd", ("i",)), c_sym(2))),
-                    pred=["static"])
-
-        op2.par_loop(op2.Kernel(k.gencode(), 'k'), node,
+        k = """
+static void k(int *vd) {
+  for (int i=0; i<1; ++i)
+    vd[i] = 2;
+}
+        """
+
+        op2.par_loop(op2.Kernel(k, 'k'), node,
                      vd1(op2.WRITE, node2ele))
         assert all(vd1.data == 2)
 
@@ -149,11 +148,13 @@ def test_inc_1d_itspace_map(self, node, d1, vd1, node2ele):
         vd1.data[:] = 3
         d1.data[:] = numpy.arange(nnodes).reshape(d1.data.shape)
 
-        k = FunDecl("void", "k",
-                    [Decl("int*", c_sym("vd")), Decl("int*", c_sym("d"))],
-                    c_for("i", 1, Incr(Symbol("vd", ("i",)), c_sym("*d"))),
-                    pred=["static"])
-        op2.par_loop(op2.Kernel(k.gencode(), 'k'), node,
+        k = """
+static void k(int *vd, int *d) {
+  for (int i=0; i<1; ++i)
+    vd[i] += d[0];
+}
+        """
+        op2.par_loop(op2.Kernel(k, 'k'), node,
                      vd1(op2.INC, node2ele),
                      d1(op2.READ))
         expected = numpy.zeros_like(vd1.data)
@@ -166,17 +167,15 @@ def test_inc_1d_itspace_map(self, node, d1, vd1, node2ele):
 
     def test_read_2d_itspace_map(self, d2, vd2, node2ele, node):
         vd2.data[:] = numpy.arange(nele * 2).reshape(nele, 2)
-        reads = Block(
-            [Assign(Symbol("d", (0,)), Symbol("vd", ("i",), ((1, 0),))),
-             Assign(
-                 Symbol(
-                     "d", (1,)), Symbol("vd", ("i",), ((1, 1),)))],
-            open_scope=True)
-        k = FunDecl("void", "k",
-                    [Decl("int*", c_sym("d")), Decl("int*", c_sym("vd"))],
-                    c_for("i", 1, reads),
-                    pred=["static"])
-        op2.par_loop(op2.Kernel(k.gencode(), 'k'), node,
+        k = """
+static void k(int *d, int *vd) {
+  for (int i=0; i<1; ++i) {
+    d[0] = vd[i];
+    d[1] = vd[i+1];
+  }
+}
+        """
+        op2.par_loop(op2.Kernel(k, 'k'), node,
                      d2(op2.WRITE),
                      vd2(op2.READ, node2ele))
         assert all(d2.data[::2, 0] == vd2.data[:, 0])
@@ -185,14 +184,15 @@ def test_read_2d_itspace_map(self, d2, vd2, node2ele, node):
         assert all(d2.data[1::2, 1] == vd2.data[:, 1])
 
     def test_write_2d_itspace_map(self, vd2, node2ele, node):
-        writes = Block([Assign(Symbol("vd", ("i",), ((1, 0),)), c_sym(2)),
-                        Assign(Symbol("vd", ("i",), ((1, 1),)), c_sym(3))],
-                       open_scope=True)
-        k = FunDecl("void", "k",
-                    [Decl("int*", c_sym("vd"))],
-                    c_for("i", 1, writes),
-                    pred=["static"])
-        op2.par_loop(op2.Kernel(k.gencode(), 'k'), node,
+        k = """
+static void k(int *vd) {
+  for (int i=0; i<1; ++i) {
+    vd[i] = 2;
+    vd[i+1] = 3;
+  }
+}
+        """
+        op2.par_loop(op2.Kernel(k, 'k'), node,
                      vd2(op2.WRITE, node2ele))
         assert all(vd2.data[:, 0] == 2)
         assert all(vd2.data[:, 1] == 3)
@@ -202,16 +202,16 @@ def test_inc_2d_itspace_map(self, d2, vd2, node2ele, node):
         vd2.data[:, 1] = 4
         d2.data[:] = numpy.arange(2 * nnodes).reshape(d2.data.shape)
 
-        incs = Block([Incr(Symbol("vd", ("i",), ((1, 0),)), Symbol("d", (0,))),
-                      Incr(
-                          Symbol("vd", ("i",), ((1, 1),)), Symbol("d", (1,)))],
-                     open_scope=True)
-        k = FunDecl("void", "k",
-                    [Decl("int*", c_sym("vd")), Decl("int*", c_sym("d"))],
-                    c_for("i", 1, incs),
-                    pred=["static"])
+        k = """
+static void k(int *vd, int *d) {
+  for (int i=0; i<1; ++i) {
+    vd[i] += d[0];
+    vd[i+1] += d[1];
+  }
+}
+        """
 
-        op2.par_loop(op2.Kernel(k.gencode(), 'k'), node,
+        op2.par_loop(op2.Kernel(k, 'k'), node,
                      vd2(op2.INC, node2ele),
                      d2(op2.READ))
 
diff --git a/test/unit/test_matrices.py b/test/unit/test_matrices.py
index 34b467e21..d7f27ff8b 100644
--- a/test/unit/test_matrices.py
+++ b/test/unit/test_matrices.py
@@ -39,10 +39,9 @@
 from pyop2 import op2
 from pyop2.exceptions import MapValueError, ModeValueError
 
-from coffee.base import *
-
 from petsc4py.PETSc import ScalarType
 
+
 # Data type
 valuetype = ScalarType
 
@@ -164,7 +163,8 @@ def x_vec(dvnodes):
 
 @pytest.fixture
 def mass():
-    init = FlatBlock("""
+    kernel_code = """
+static void mass(double localTensor[3][3], double c0[3][2]) {
   double CG1[3][6] = { {  0.09157621, 0.09157621, 0.81684757,
                                    0.44594849, 0.44594849, 0.10810302 },
                                 {  0.09157621, 0.81684757, 0.09157621,
@@ -206,26 +206,23 @@ def mass():
       };
     };
   };
-  for(int i_g = 0; i_g < 6; i_g++)
-""")
-    assembly = Incr(Symbol("localTensor", ("i_r_0", "i_r_1")),
-                    FlatBlock("ST0 * w[i_g]"))
-    assembly = Block([FlatBlock("double ST0 = 0.0;\nST0 += CG1[i_r_0][i_g] * CG1[i_r_1][i_g] * (c_q0[i_g][0][0] * \
-                                c_q0[i_g][1][1] + -1 * c_q0[i_g][0][1] * c_q0[i_g][1][0]);\n"), assembly], open_scope=True)
-    assembly = c_for("i_r_0", 3, c_for("i_r_1", 3, assembly))
-
-    kernel_code = FunDecl("void", "mass",
-                          [Decl("double", Symbol("localTensor", (3, 3))),
-                           Decl("double", Symbol("c0", (3, 2)))],
-                          Block([init, assembly], open_scope=False),
-                          pred=["static"])
-
-    return op2.Kernel(kernel_code.gencode(), "mass")
+  for(int i_g = 0; i_g < 6; i_g++) {
+    for (int i_r_0=0; i_r_0<3; ++i_r_0) {
+      for (int i_r_1=0; i_r_1<3; ++i_r_1) {
+        double ST0 = 0.0;
+        ST0 += CG1[i_r_0][i_g] * CG1[i_r_1][i_g] * (c_q0[i_g][0][0] * c_q0[i_g][1][1] + -1 * c_q0[i_g][0][1] * c_q0[i_g][1][0]);
+        localTensor[i_r_0][i_r_1] += ST0 * w[i_g];
+      }
+    }
+  }
+}
+    """
+    return op2.Kernel(kernel_code, "mass")
 
 
 @pytest.fixture
 def rhs():
-    kernel_code = FlatBlock("""
+    kernel_code = """
 static void rhs(double* localTensor, double c0[3][2], double* c1)
 {
   double CG1[3][6] = { {  0.09157621, 0.09157621, 0.81684757,
@@ -284,45 +281,40 @@ def rhs():
       localTensor[i_r_0] += ST1 * w[i_g];
     };
   };
-}""")
-    return op2.Kernel(kernel_code.gencode(), "rhs")
+}"""
+    return op2.Kernel(kernel_code, "rhs")
 
 
 @pytest.fixture
 def mass_ffc():
-    init = FlatBlock("""
-double J_00 = x[1][0] - x[0][0];
-double J_01 = x[2][0] - x[0][0];
-double J_10 = x[1][1] - x[0][1];
-double J_11 = x[2][1] - x[0][1];
-
-double detJ = J_00*J_11 - J_01*J_10;
-double det = fabs(detJ);
-
-double W3[3] = {0.166666666666667, 0.166666666666667, 0.166666666666667};
-double FE0[3][3] = \
-{{0.666666666666667, 0.166666666666667, 0.166666666666667},
-{0.166666666666667, 0.166666666666667, 0.666666666666667},
-{0.166666666666667, 0.666666666666667, 0.166666666666667}};
-
-for (unsigned int ip = 0; ip < 3; ip++)
-""")
-    assembly = Incr(Symbol("A", ("j", "k")),
-                    FlatBlock("FE0[ip][j]*FE0[ip][k]*W3[ip]*det"))
-    assembly = c_for("j", 3, c_for("k", 3, assembly))
-
-    kernel_code = FunDecl("void", "mass_ffc",
-                          [Decl("double", Symbol("A", (3, 3))),
-                           Decl("double", Symbol("x", (3, 2)))],
-                          Block([init, assembly], open_scope=False),
-                          pred=["static"])
-
-    return op2.Kernel(kernel_code.gencode(), "mass_ffc")
+    kernel_code = """
+static void mass_ffc(double A[3][3], double x[3][2]) {
+  double J_00 = x[1][0] - x[0][0];
+  double J_01 = x[2][0] - x[0][0];
+  double J_10 = x[1][1] - x[0][1];
+  double J_11 = x[2][1] - x[0][1];
+
+  double detJ = J_00*J_11 - J_01*J_10;
+  double det = fabs(detJ);
+
+  double W3[3] = {0.166666666666667, 0.166666666666667, 0.166666666666667};
+  double FE0[3][3] = \
+  {{0.666666666666667, 0.166666666666667, 0.166666666666667},
+  {0.166666666666667, 0.166666666666667, 0.666666666666667},
+  {0.166666666666667, 0.666666666666667, 0.166666666666667}};
+
+  for (unsigned int ip = 0; ip < 3; ip++)
+    for (int j=0; j<3; ++j)
+      for (int k=0; k<3; ++k)
+        A[j][k] += FE0[ip][j]*FE0[ip][k]*W3[ip]*det;
+}
+    """
+    return op2.Kernel(kernel_code, "mass_ffc")
 
 
 @pytest.fixture
 def rhs_ffc():
-    kernel_code = FlatBlock("""
+    kernel_code = """
 static void rhs_ffc(double *A, double x[3][2], double *w0)
 {
     double J_00 = x[1][0] - x[0][0];
@@ -355,49 +347,39 @@ def rhs_ffc():
       }
     }
 }
-""")
-    return op2.Kernel(kernel_code.gencode(), "rhs_ffc")
+"""
+    return op2.Kernel(kernel_code, "rhs_ffc")
 
 
 @pytest.fixture
 def rhs_ffc_itspace():
-    init = FlatBlock("""
-double J_00 = x[1][0] - x[0][0];
-double J_01 = x[2][0] - x[0][0];
-double J_10 = x[1][1] - x[0][1];
-double J_11 = x[2][1] - x[0][1];
-
-double detJ = J_00*J_11 - J_01*J_10;
-double det = fabs(detJ);
-
-double W3[3] = {0.166666666666667, 0.166666666666667, 0.166666666666667};
-double FE0[3][3] = \
-{{0.666666666666667, 0.166666666666667, 0.166666666666667},
-{0.166666666666667, 0.166666666666667, 0.666666666666667},
-{0.166666666666667, 0.666666666666667, 0.166666666666667}};
-
-for (unsigned int ip = 0; ip < 3; ip++)
-{
-  double F0 = 0.0;
-
-  for (unsigned int r = 0; r < 3; r++)
-  {
-    F0 += FE0[ip][r]*w0[r];
+    kernel_code = """
+static void rhs_ffc_itspace(double A[3], double x[3][2], double *w0) {
+  double J_00 = x[1][0] - x[0][0];
+  double J_01 = x[2][0] - x[0][0];
+  double J_10 = x[1][1] - x[0][1];
+  double J_11 = x[2][1] - x[0][1];
+
+  double detJ = J_00*J_11 - J_01*J_10;
+  double det = fabs(detJ);
+
+  double W3[3] = {0.166666666666667, 0.166666666666667, 0.166666666666667};
+  double FE0[3][3] = \
+  {{0.666666666666667, 0.166666666666667, 0.166666666666667},
+  {0.166666666666667, 0.166666666666667, 0.666666666666667},
+  {0.166666666666667, 0.666666666666667, 0.166666666666667}};
+
+  for (unsigned int ip = 0; ip < 3; ip++) {
+    double F0 = 0.0;
+
+    for (unsigned int r = 0; r < 3; r++)
+      F0 += FE0[ip][r]*w0[r];
+    for (unsigned int j=0; j<3; ++j)
+      A[j] += FE0[ip][j]*F0*W3[ip]*det;
   }
-
-""")
-    assembly = Incr(Symbol("A", ("j",)), FlatBlock("FE0[ip][j]*F0*W3[ip]*det"))
-    assembly = c_for("j", 3, assembly)
-    end = FlatBlock("}")
-
-    kernel_code = FunDecl("void", "rhs_ffc_itspace",
-                          [Decl("double", Symbol("A", (3,))),
-                           Decl("double", Symbol("x", (3, 2))),
-                              Decl("double*", Symbol("w0"))],
-                          Block([init, assembly, end], open_scope=False),
-                          pred=["static"])
-
-    return op2.Kernel(kernel_code.gencode(), "rhs_ffc_itspace")
+}
+    """
+    return op2.Kernel(kernel_code, "rhs_ffc_itspace")
 
 
 @pytest.fixture
@@ -424,32 +406,26 @@ def zero_vec_dat():
 
 @pytest.fixture
 def kernel_inc():
-    code = c_for("i", 3,
-                 c_for("j", 3,
-                       Incr(Symbol("entry", ("i", "j")), c_sym("*g"))))
-
-    kernel_code = FunDecl("void", "inc",
-                          [Decl("double", Symbol("entry", (3, 3))),
-                           Decl("double*", c_sym("g"))],
-                          Block([code], open_scope=False),
-                          pred=["static"])
-
-    return op2.Kernel(kernel_code.gencode(), "inc")
+    kernel_code = """
+static void inc(double entry[3][3], double *g) {
+  for (int i=0; i<3; ++i)
+    for (int j=0; j<3; ++j)
+      entry[i][j] += g[0];
+}
+    """
+    return op2.Kernel(kernel_code, "inc")
 
 
 @pytest.fixture
 def kernel_set():
-    code = c_for("i", 3,
-                 c_for("j", 3,
-                       Assign(Symbol("entry", ("i", "j")), c_sym("*g"))))
-
-    kernel_code = FunDecl("void", "set",
-                          [Decl("double", Symbol("entry", (3, 3))),
-                           Decl("double*", c_sym("g"))],
-                          Block([code], open_scope=False),
-                          pred=["static"])
-
-    return op2.Kernel(kernel_code.gencode(), "set")
+    kernel_code = """
+static void set(double entry[3][3], double *g) {
+  for (int i=0; i<3; ++i)
+    for (int j=0; j<3; ++j)
+      entry[i][j] = g[0];
+}
+    """
+    return op2.Kernel(kernel_code, "set")
 
 
 @pytest.fixture
@@ -642,20 +618,18 @@ def test_mat_always_has_diagonal_space(self):
 
     def test_minimal_zero_mat(self):
         """Assemble a matrix that is all zeros."""
-
-        code = c_for("i", 1,
-                     c_for("j", 1,
-                           Assign(Symbol("local_mat", ("i", "j")), c_sym("0.0"))))
-        zero_mat_code = FunDecl("void", "zero_mat",
-                                [Decl("double", Symbol("local_mat", (1, 1)))],
-                                Block([code], open_scope=False))
+        zero_mat_code = """
+void zero_mat(double local_mat[1][1]) {
+  local_mat[0][0] = 0.0;
+}
+        """
 
         nelems = 128
         set = op2.Set(nelems)
         map = op2.Map(set, set, 1, np.array(list(range(nelems)), np.uint32))
         sparsity = op2.Sparsity((set, set), (map, map))
         mat = op2.Mat(sparsity, np.float64)
-        kernel = op2.Kernel(zero_mat_code.gencode(), "zero_mat")
+        kernel = op2.Kernel(zero_mat_code, "zero_mat")
         op2.par_loop(kernel, set,
                      mat(op2.WRITE, (map, map)))
 
@@ -919,12 +893,13 @@ def mat(self, msparsity, mmap, mdat):
     @pytest.fixture
     def dat(self, mset, mmap, mdat):
         dat = op2.MixedDat(mset)
-        kernel_code = FunDecl("void", "addone_rhs",
-                              [Decl("double", Symbol("v", (3,))),
-                               Decl("double", Symbol("d", (3,)))],
-                              c_for("i", 3, Incr(Symbol("v", ("i")), FlatBlock("d[i]"))),
-                              pred=["static"])
-        addone = op2.Kernel(kernel_code.gencode(), "addone_rhs")
+        kernel_code = """
+static void addone_rhs(double v[3], double d[3]) {
+  for (int i=0; i<3; ++i)
+    v[i] += d[i];
+}
+        """
+        addone = op2.Kernel(kernel_code, "addone_rhs")
         op2.par_loop(addone, mmap.iterset,
                      dat(op2.INC, mmap),
                      mdat(op2.READ, mmap))
@@ -947,15 +922,15 @@ def test_assemble_mixed_rhs(self, dat):
     def test_assemble_mixed_rhs_vector(self, mset, mmap, mvdat):
         """Assemble a simple right-hand side over a mixed space and check result."""
         dat = op2.MixedDat(mset ** 2)
-        assembly = Block(
-            [Incr(Symbol("v", ("i"), ((2, 0),)), FlatBlock("d[i][0]")),
-             Incr(Symbol("v", ("i"), ((2, 1),)), FlatBlock("d[i][1]"))], open_scope=True)
-        kernel_code = FunDecl("void", "addone_rhs_vec",
-                              [Decl("double", Symbol("v", (6,))),
-                               Decl("double", Symbol("d", (3, 2)))],
-                              c_for("i", 3, assembly),
-                              pred=["static"])
-        addone = op2.Kernel(kernel_code.gencode(), "addone_rhs_vec")
+        kernel_code = """
+static void addone_rhs_vec(double v[6], double d[3][2]) {
+  for (int i=0; i<3; ++i) {
+    v[i*2+0] += d[i][0];
+    v[i*2+1] += d[i][1];
+  }
+}
+        """
+        addone = op2.Kernel(kernel_code, "addone_rhs_vec")
         op2.par_loop(addone, mmap.iterset,
                      dat(op2.INC, mmap),
                      mvdat(op2.READ, mmap))
diff --git a/test/unit/test_subset.py b/test/unit/test_subset.py
index a2c2f9f9b..33936df7d 100644
--- a/test/unit/test_subset.py
+++ b/test/unit/test_subset.py
@@ -37,8 +37,6 @@
 
 from pyop2 import op2
 
-from coffee.base import *
-
 nelems = 32
 
 
@@ -224,15 +222,14 @@ def test_matrix(self):
         mat01 = op2.Mat(sparsity, np.float64)
         mat10 = op2.Mat(sparsity, np.float64)
 
-        assembly = c_for("i", 4,
-                         c_for("j", 4,
-                               Incr(Symbol("mat", ("i", "j")), FlatBlock("(*dat)*16+i*4+j"))))
-        kernel_code = FunDecl("void", "unique_id",
-                              [Decl("double", Symbol("mat", (4, 4))),
-                               Decl("double*", c_sym("dat"))],
-                              Block([assembly], open_scope=False),
-                              pred=["static"])
-        k = op2.Kernel(kernel_code.gencode(), "unique_id")
+        kernel_code = """
+static void unique_id(double mat[4][4], double *dat) {
+  for (int i=0; i<4; ++i)
+    for (int j=0; j<4; ++j)
+      mat[i][j] += (*dat)*16+i*4+j;
+}
+        """
+        k = op2.Kernel(kernel_code, "unique_id")
 
         mat.zero()
         mat01.zero()