Explicit calculation of gradient of a scalar field

General Discussion
1

So my whole problem consists of two stages:
First a heat equation is solved:

problem Heat(T2,v) = int3d(Th)(T2 * v) + int3d(Th)(dt * nu * Grad(T2)' * Grad(v)) - int3d(Th)(T1 * v)
                    +int3d(Th, 1, 2, 3, 4, 6)(epsilon * stefConst * (pow(T1, 4.) - pow(Tamb + 273., 4.)) * v)
                    +on(5, T2 = Ttable + 273.)
                    ;

After the calculated temperature T2 should be passed to the following problem

problem Navier([u1, u2, u3], [v1, v2, v3]) = int3d(Th)(lambda * Div(u1, u2, u3) * Div(v1, v2, v3))
                                        +int3d(Th)(2. * mu * (Epsilon(u1, u2, u3)' * Epsilon(v1, v2, v3)))
                                        -int3d(Th)(beta * [v1, v2, v3]' * [0., 0., 0.])
                                        +on(5, u1 = 0, u2 = 0, u3 =0)
                                        ;

The array [0., 0., 0.] should be replaced by gradient of T2. Yet I am not completely sure how to explicitly calculate the gradient of a given scalar field. Here

macro Grad(u) [dx(u), dy(u), dz(u)] //EOM
macro Div(u1,u2,u3) (dx(u1) + dy(u2) + dz(u3)) // EOM
macro Epsilon(u1,u2,u3)  [dx(u1), dy(u2), dz(u3), (dz(u2) + dy(u3)) / sqrt2, (dz(u1) + dx(u3)) / sqrt2, (dy(u1) + dx(u2)) / sqrt2] // EOM

Could anyone help?