  This paper, as the sequel to previous work, develops numerical schemes for
fractional diffusion equations on a two-dimensional finite domain with
triangular meshes. We adopt the nodal discontinuous Galerkin methods for the
full spatial discretization by the use of high-order nodal basis, employing
multivariate Lagrange polynomials defined on the triangles. Stability analysis
and error estimates are provided, which shows that if polynomials of degree $N$
are used, the methods are (N+1)-th order accurate for general triangulations.
Finally, the performed numerical experiments confirm the optimal order of
convergence.
