How to properly define matrix in the integral of form

I’m using the latest version of dolfinx and I wonder how to define matrix to compute integral of w^{T}Ay, here, y and w are trial and test functions, I use this method to define matrix

from dolfinx import fem
import ufl
import numpy as np
from ufl import as_ufl


def create_coefficient_matrices(p_f, u_f, T_f, domain):
    # Define constants
    R = fem.Constant(domain, np.float64(287.05))  # Specific gas constant for air
    c_v = fem.Constant(domain, np.float64(718.0))  # Specific heat at constant volume
    c_p = fem.Constant(domain, np.float64(1.0054))
    # Define basic thermodynamic expressions directly
    rho_f_expr = p_f / (as_ufl(R) * T_f)
    v_expr = as_ufl(R) * T_f / p_f
    e_expr = as_ufl(c_v) * T_f
    h_expr = e_expr + p_f * v_expr
    k_expr=0.5 * (u_f[0] * u_f[0] + u_f[1] * u_f[1])
    u_squared = 0.5 * (u_f[0] * u_f[0] + u_f[1] * u_f[1])
    e_tot_expr = e_expr + u_squared

    # Define thermodynamic coefficients directly
    beta_T_expr = 1.0 / p_f
    alpha_p_expr = 1.0 / T_f
    e1_p = rho_f_expr *beta_T_expr*(h_expr + k_expr )-alpha_p_expr*T_f
    e2_p = e1_p + as_ufl(fem.Constant(domain, np.float64(1)))
    e3_p = rho_f_expr * (h_expr+k_expr)
    e4_p=-rho_f_expr *alpha_p_expr*(h_expr+k_expr) + rho_f_expr*c_p
    # Create coefficient matrices using ufl.as_matrix
    get_A0 = ufl.as_matrix([
        [rho_f_expr * beta_T_expr, fem.Constant(domain, np.float64(0)), fem.Constant(domain, np.float64(0)),
         -rho_f_expr * alpha_p_expr],
        [rho_f_expr * beta_T_expr * u_f[0], rho_f_expr, fem.Constant(domain, np.float64(0)),
         -rho_f_expr * alpha_p_expr * u_f[0]],
        [rho_f_expr * beta_T_expr * u_f[1], fem.Constant(domain, np.float64(0)), rho_f_expr,
         -rho_f_expr * alpha_p_expr * u_f[1]],
        [e1_p , rho_f_expr * u_f[0], rho_f_expr * u_f[1],
         rho_f_expr * c_v]
    ])

but after debugging, I find out that the return of this method, which is 4*4, doesn’t contain any attributes like .x or .value, and I am using this in my weak form integration, I think this may cause the process very slow, so how to properly do this in dolfinx, could you please recommend some docs? thank you very much

the function is a mixed element space, which contains p,u,T here u is 2 dimensional

If you want to compute the residual of the PDE system in your code, make sure that you first interpolate your Res at first, otherwise it would be very time consuming or even out of memory