I guess, the default quadrature order for bilinear elements should be 2: We have a first order scheme, so its fine to assume that `function` is bilinear. Under this assumption, the products which are integrated contain monomials up to x²y², which are exactly integrated by a (tensor) quadrature of order 2.  Did you have any reasons to increase the order to 4?