import sympy


def generate_monomials(order):
    """
    Generate trivariate monomials up to given order.
    """
    monoms = []
    x, y, z = sympy.symbols('x y z')
    for i in range(0, order + 1):
        for j in range(0, order + 1):
            for k in range(0, order + 1):
                if i + j + k <= order:
                    monoms.append(x**i * y**j * z**k)
    return monoms