Handle plate models
For plate models, we want to integrate the behaviour at NG = nG x mG Gauss points in total where nG is the number of Gauss points on the plate 2D mesh and mG the number of Gauss points across the thickness.
We need a mechanism to handle the case of non-matching values between NG (used to load the behaviour) and nG the number of Gauss points for the 2D mesh. In addition to just reshaping the appropriate arrays, we can have different expressions for the same gradient depending on the z-coordinate, the integration weights in the form may also be different from one z-Gauss point to another.
Possible syntax:
problem.register_gradient("Strain", [sym(grad(u-zi*theta)) for zi in z_Gauss], weights=w_Gauss)