Why did you make the integration order configurable? As `diffusion_function` is only evaluated at the entity center (line 535), we know that a quadrature of order 2 will be sufficient ...