diff --git a/python/dune/codegen/pdelab/tensors.py b/python/dune/codegen/pdelab/tensors.py
index a1a52a4c4012f404aaf036b5067facb6067f5b4a..360f9a1425ae1d2111ea15868906618b39b2f914 100644
--- a/python/dune/codegen/pdelab/tensors.py
+++ b/python/dune/codegen/pdelab/tensors.py
@@ -24,8 +24,7 @@ def define_determinant(name, matrix, shape, visitor):
matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
if dim == 2:
- expr_determinant = FMA(matrix_entry[0][0], matrix_entry[1][1],
- -1 * prim.Product((matrix_entry[1][0], matrix_entry[0][1])))
+ expr_determinant = FMA(matrix_entry[0][0], matrix_entry[1][1], -1 * matrix_entry[1][0] * matrix_entry[0][1])
elif dim == 3:
fma_A = FMA(matrix_entry[1][1], matrix_entry[2][2], -1 * matrix_entry[1][2] * matrix_entry[2][1])
@@ -73,35 +72,26 @@ def define_matrix_inverse(name, name_inv, shape, visitor):
sign = 1. if i == j else -1.
exprs[i][j] = prim.Product((sign, prim.Variable(det_inv), matrix_entry[1 - i][1 - j]))
elif dim == 3:
- exprs[0][0] = prim.Product((1., prim.Variable(det_inv),
- prim.Sum((prim.Product((matrix_entry[1][1], matrix_entry[2][2])),
- -1 * prim.Product((matrix_entry[1][2], matrix_entry[2][1]))))))
- exprs[1][0] = prim.Product((-1., prim.Variable(det_inv),
- prim.Sum((prim.Product((matrix_entry[0][1], matrix_entry[2][2])),
- -1 * prim.Product((matrix_entry[0][2], matrix_entry[2][1]))))))
- exprs[2][0] = prim.Product((1., prim.Variable(det_inv),
- prim.Sum((prim.Product((matrix_entry[0][1], matrix_entry[1][2])),
- -1 * prim.Product((matrix_entry[0][2], matrix_entry[1][1]))))))
-
- exprs[0][1] = prim.Product((-1., prim.Variable(det_inv),
- prim.Sum((prim.Product((matrix_entry[1][0], matrix_entry[2][2])),
- -1 * prim.Product((matrix_entry[1][2], matrix_entry[2][0]))))))
- exprs[1][1] = prim.Product((1., prim.Variable(det_inv),
- prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[2][2])),
- -1 * prim.Product((matrix_entry[0][2], matrix_entry[2][0]))))))
- exprs[2][1] = prim.Product((-1., prim.Variable(det_inv),
- prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][2])),
- -1 * prim.Product((matrix_entry[0][2], matrix_entry[1][0]))))))
-
- exprs[0][2] = prim.Product((1., prim.Variable(det_inv),
- prim.Sum((prim.Product((matrix_entry[1][0], matrix_entry[2][1])),
- -1 * prim.Product((matrix_entry[1][1], matrix_entry[2][0]))))))
- exprs[1][2] = prim.Product((-1., prim.Variable(det_inv),
- prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[2][1])),
- -1 * prim.Product((matrix_entry[0][1], matrix_entry[2][0]))))))
- exprs[2][2] = prim.Product((1., prim.Variable(det_inv),
- prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1])),
- -1 * prim.Product((matrix_entry[0][1], matrix_entry[1][0]))))))
+ exprs[0][0] = prim.Variable(det_inv) * FMA(matrix_entry[1][1], matrix_entry[2][2],
+ -1 * matrix_entry[1][2] * matrix_entry[2][1])
+ exprs[1][0] = prim.Variable(det_inv) * FMA(matrix_entry[0][1], matrix_entry[2][2],
+ -1 * matrix_entry[0][2] * matrix_entry[2][1]) * -1
+ exprs[2][0] = prim.Variable(det_inv) * FMA(matrix_entry[0][1], matrix_entry[1][2],
+ -1 * matrix_entry[0][2] * matrix_entry[1][1])
+
+ exprs[0][1] = prim.Variable(det_inv) * FMA(matrix_entry[1][0], matrix_entry[2][2],
+ -1 * matrix_entry[1][2] * matrix_entry[2][0]) * -1
+ exprs[1][1] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[2][2],
+ -1 * matrix_entry[0][2] * matrix_entry[2][0])
+ exprs[2][1] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[1][2],
+ -1 * matrix_entry[0][2] * matrix_entry[1][0]) * -1
+
+ exprs[0][2] = prim.Variable(det_inv) * FMA(matrix_entry[1][0], matrix_entry[2][1],
+ -1 * matrix_entry[1][1] * matrix_entry[2][0])
+ exprs[1][2] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[2][1],
+ -1 * matrix_entry[0][1] * matrix_entry[2][0]) * -1
+ exprs[2][2] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[1][1],
+ -1 * matrix_entry[0][1] * matrix_entry[1][0])
else:
raise NotImplementedError
for j in range(dim):