diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 0129f4d70888be3db7f4ef1dd418956252754fba..2e0f3e195677d1b6834214764500fe1172712f60 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -98,7 +98,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
                 from dune.perftool.pdelab.geometry import dimension_iname
                 from dune.perftool.pdelab.basis import lfs_child
                 idims = tuple(dimension_iname(context='arg', count=i) for i in range(len(subel.value_shape())))
-                lfs = lfs_child(lfs, idims, shape=subel.value_shape())
+                lfs = lfs_child(lfs, idims, shape=subel.value_shape(), symmetry=subel.symmetry)
                 residual_shape[icount] = name_lfs_bound(name_leaf_lfs(subel.sub_elements()[0], ma.restriction))
                 subel = subel.sub_elements()[0]
             else:
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index edbf22a2f1a7feca3d0d08901ae886a8df559675..4c6c40114c52dd3d5d6cd3d3e63646882c3ff356 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -75,14 +75,25 @@ def define_lfs(name, father, child):
     return "auto {} = child({}, _{});".format(name, father, child)
 
 
-def lfs_child(lfs, children, shape=None):
+def lfs_child(lfs, children, shape=None, symmetry=False):
     from pymbolic.primitives import Call, Product, Sum
     # Old pre-TensorElement implementation kept for comaptibility
     if shape is None:
         indices = (Variable(children[0]),)
     else:
-        children = tuple(Variable(c) for c in children)
-        indices = (Sum(tuple(Product((Product(tuple(shape[j] for j in range(i))), child)) for i, child in enumerate(children))),)
+        if symmetry and len(children) == 2:
+            # I do not want to think about tensors of rank > 2 right now
+            i, j = children
+            if i>j:
+                j, i = i, j
+            i = Variable(i)
+            j = Variable(j)
+            n = len(children)
+
+            indices = (Sum((Product((n-1, i)), Product((.5, i, 1-i)), j)),)
+        else:
+            children = tuple(Variable(c) for c in children)
+            indices = (Sum(tuple(Product((Product(tuple(shape[j] for j in range(i))), child)) for i, child in enumerate(children))),)
 
     from dune.perftool.loopy.functions import LFSChild
     return Call(LFSChild(lfs), indices)
@@ -367,7 +378,7 @@ def evaluate_coefficient(element, name, container, restriction, component):
 
     if isinstance(sub_element, (VectorElement, TensorElement)):
         # TODO we need a concept of symmetry here
-        lfs = lfs_child(lfs, idims, shape=shape_as_pymbolic(shape))
+        lfs = lfs_child(lfs, idims, shape=shape_as_pymbolic(shape), symmetry=element.symmetry)
 
     from dune.perftool.pdelab.argument import pymbolic_coefficient
     coeff = pymbolic_coefficient(container, lfs, index)
@@ -406,7 +417,7 @@ def evaluate_coefficient_gradient(element, name, container, restriction, compone
 
     if isinstance(sub_element, (VectorElement, TensorElement)):
         # TODO we need a concept of symmetry here
-        lfs = lfs_child(lfs, idims[:-1], shape=shape_as_pymbolic(shape[:-1]))
+        lfs = lfs_child(lfs, idims[:-1], shape=shape_as_pymbolic(shape[:-1]), symmetry=element.symmetry)
 
     from dune.perftool.pdelab.argument import pymbolic_coefficient
     coeff = pymbolic_coefficient(container, lfs, index)
diff --git a/test/stokes/CMakeLists.txt b/test/stokes/CMakeLists.txt
index b8bc4b0a2c259746f58099d624104a92ca4460e0..0553ee4c48feeb4da48e5fd912f02c91acc4c01e 100644
--- a/test/stokes/CMakeLists.txt
+++ b/test/stokes/CMakeLists.txt
@@ -21,3 +21,8 @@ dune_add_formcompiler_system_test(UFLFILE stokes_stress.ufl
                                   BASENAME stokes_stress
                                   INIFILE stokes_stress.mini
                                   )
+
+dune_add_formcompiler_system_test(UFLFILE stokes_stress_sym.ufl
+                                  BASENAME stokes_stress_sym
+                                  INIFILE stokes_stress_sym.mini
+                                  )
diff --git a/test/stokes/stokes_stress_sym.mini b/test/stokes/stokes_stress_sym.mini
new file mode 100644
index 0000000000000000000000000000000000000000..9038ebfae579e7c43817ad9e5472c37cd7fbeee7
--- /dev/null
+++ b/test/stokes/stokes_stress_sym.mini
@@ -0,0 +1,18 @@
+__name = stokes_stress_sym_{__exec_suffix}
+__exec_suffix = symdiff, numdiff | expand num
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 16 16
+elementType = simplical
+printmatrix = false
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = hagenpoiseuille_ref
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-6
diff --git a/test/stokes/stokes_stress_sym.ufl b/test/stokes/stokes_stress_sym.ufl
new file mode 100644
index 0000000000000000000000000000000000000000..47e41380efddbaca3fc089dbfa9844ab1e44b695
--- /dev/null
+++ b/test/stokes/stokes_stress_sym.ufl
@@ -0,0 +1,25 @@
+v_bctype = Expression("if (x[0] < 1. - 1e-8) return 1; else return 0;", on_intersection=True)
+g_v = Expression(("4*x[1]*(1.-x[1])", "0.0"))
+g_p = Expression("8*(1.-x[0])")
+g_S = Expression(("0.0", "-8*x[1] + 4.", "0.0"))
+g = g_v * g_p * g_S
+
+cell = triangle
+P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=g_v)
+P1 = FiniteElement("Lagrange", cell, 1)
+P2_stress = TensorElement("DG", cell, 1, symmetry={(0, 1): (1, 0)})
+
+TH = MixedElement(P2, P1, P2_stress)
+
+v, q, T  = TestFunctions(TH)
+u, p, S  = TrialFunctions(TH)
+
+n = FacetNormal(triangle)('+')
+ds = ds(subdomain_data=v_bctype, subdomain_id=0)
+
+r = (inner(grad(v), S) + inner(2*sym(grad(u)) - S, T) - div(v)*p - q*div(u))*dx \
+  - inner(grad(u).T*n, v)*ds
+# \
+#  + inner(S.T*n, v)*ds
+
+forms = [r]