Dear all, I would like to solve equations in which one of the unknowns is a vector that I need to multiply to some matrices. However, these matrices are not computed with the assemble function but with a python function and are thus arrays. To put this into context, I am solving 3D equations with a separation of variable between the (x,y) dependency and the z dependency. Then, the z-discretization is made "by hand" and then substituted into the horizontal weak formulations which are solved with Firedrake. Therefore, the discrete vertical matrices are computed separately, but I would like to include them into my weak formulations. The issue is that it does not seem to work when multiplying a function with an array in the variational problem. How could I do that ? The WF and a description of the objective is in the attached PDF. Thanks in advance, Floriane